Correction exponents in the Gross - Neveu - Yukawa model at $1/N^2$

We calculate the critical exponents $\omega_\pm$ in the $d$-dimensional Gross-Neveu model in $1/N$ expansion with $1/N^2$ accuracy. These exponents are related to the slopes of the $\beta$-functions at the critical point in the Gross - Neveu - Yukawa model. They have been computed recently to four loops accuracy. We checked that our results are in complete agreement with the results of the perturbative calculations.


Introduction
The Gross -Neveu (GN) model, originally introduced in [1] in 1974 as a toy model for a study of dynamical symmetry breaking, found its application in many areas of physics. The model describes a system of fermion fields with a quartic interaction. It is renormalizable in two dimensions, asymptotically free and, moreover, admits an exact solution [2]. Above two dimensions the N-component GN -model is renormalizable within the 1/N expansion technique. It possesses a nontrivial Wilson -Fisher fixed point and gives an example of a conformal field theory (CFT) in d-dimensions. Basic critical indices of the GN model are available at 1/N 2 order [3][4][5][6][7] and the index η -the anomalous dimension of the fermion field -at 1/N 3 [8,9]. These results were obtained by methods of the self-consistency equation and conformal bootstrap developed in [10][11][12]. A recent revival of interest to this model is due to a relevance of fermionic systems with quartic interactions to the description of the phase transition in graphene [13].
The UV completion of the GN model contains an additional scalar field with a quartic self-interaction and is known as the Gross-Neveu-Yukawa (GNY) model [14]. Moreover, the chiral extension of the GNY model -the Nambu-Jona-Lasinio (NJL) model -describes, for a low number of fermion a e-mail: alexander.manashov@desy.de b e-mail: matthias.strohmaier@ur.de flavors, a system with an emergent supersymmetry [15]. The recent calculation of the corresponding renormalization group (RG) functions with four-loop accuracy [16][17][18] confirms the emergence of supersymmetry at the critical point.
The four-loop calculations are quite involved and carried out with the help of computer algebra. The results available from 1/N expansion provide an additional check for the perturbative calculations. The four loop RG functions obtained in [18] are in a perfect agreement with the results of 1/N calculations [3][4][5][6][7][8][9]. In this paper we present 1/N 2 expressions for the other two indices -the so-called correction exponents that are related to the slopes of β functions at the critical point and were known only with 1/N accuracy [19].
The 1/N calculations are usually done with the help of the methods of the self-consistency equations or conformal bootstrap [10][11][12]. These methods are very effective for the calculation of the critical indices of the basic and auxiliary fields, but are not very suitable in the case of local operators. Therefore we use a different method developed in [20][21][22]. A detailed description of the method can be found in [22,23] and an application to the GN model in [24].
The paper is organized as follows: In sect. 2 we recall the formulation of the GNY model and show that the slopes of β functions at the critical point coincide with the critical dimensions of certain operators of dimension four. In sect. 3 we review briefly the 1/N expansion technique for the GN model and present the rules for calculating the anomalous dimensions of local operators. Section 4 contains the details of the calculation of the correction exponents at order 1/N 2 . In Appendix A we collect expressions for relevant basic indices. Appendix B contains details of the calculation of the self-energy and vertex correction diagrams and in Appendix C we collect results for the individual diagrams. The Lagrangian of the GNY model in d = 4 − 2ε Euclidean space takes the form where q is the N-component fermion field and σ is a scalar field. The model is multiplicatively renormalized [14] Here M is the renormalization scale and Z i are the renormalization constants. In MS-like schemes the renormalization factors Z i are given by a series in 1/ε, Z i = 1 + ∑ k z k i /ε k . The anomalous dimensions of the fields and the β functions are defined as usual where Z σ = Z The GNY model possesses non-trivial fixed points One of this points is infrared stable, i.e. the eigenvalues of the matrix are both positive. These eigenvalues, ω ± , are called correction indices. The indices ω ± can be identified with the anomalous dimensions of certain composite operators that will be discussed in the rest of this section. The partition function Z(J) defined by the path integral where J = { j, η,η} and JΦ = jσ +ηq+qη, is finite in perturbation theory. Moreover, derivatives of (7) with respect to the renormalized couplings g i are finite, too. Thus one concludes that where the sum goes over some set of renormalized operators of the canonical dimension four with finite coefficients c a . By ∂ B we denote terms (not necessarily finite) which vanish after integration. Since ∂ g 1 L R ≡ M ε σqq + singular 1/ε terms (9) and similar for another derivative one concludes that where we omitted terms with total derivatives. Applying the operator (4) to both sides of Eqs. (10) and taking into account that D M S R = D M S 0 = 0 one derives the RG equation for the operators where the anomalous dimension matrix γ has the following form Using the known lowest order expression for the beta-functions (see e.g. ref. [18]) one finds that at the Wilson-Fisher critical point the anomalous dimensions take the form resulting in the following scaling dimensions for the operators O ± , Diagonalizing the matrix ∆ one constructs two operators whose scaling dimensions are given by the eigenvalues of ∆ . Scaling dimensions of operators are physical observables and do not depend on a regularization, expansion, subtraction scheme, etc. We use this property to calculate the anomalous dimension of the operators [O ± ] in 1/N expansion.

Large N expansion for the GN model
The GNY model is critically equivalent to the d-dimensional N-component GN model, i.e. critical indices in both models coincide. In the GN model the indices can be calculated using the 1/N expansion. Below we briefly review the method that we use for the calculation. A more extensive review of the 1/N techniques including the self-consistency equations and conformal bootstrap can be found in the book [25]. Let us write the action of the GN model in the form suitable for generating the 1/N expansion Here q = {q i , i = 1, . . . , N} is the N-component fermion field and σ is an auxiliary scalar field which can be excluded by the equation of motion (EOM). At a certain (critical) value of the coupling, g = g * the system undergoes a second order phase transition. At this point correlators of the basic and auxiliary fields q,q, σ exhibit power law behaviour and, as it can be shown, the model enjoys scale and conformal invariance [5].
It can be shown that in the infrared region (IR) the dominant contribution to the propagator of the σ -field comes from the fermion loop [14] Here n = N tr1l, where the trace is taken in the space of ddimensional spinors and the normalization factor is For practical calculations it is convenient to use a simplified (massless) version of the GN model which is critically equivalent to (15). The action of the model is given by the following expression [20, 25] The kernel L has the form where D q is the fermion propagator The regularized kernel L ∆ is The first two terms in (18) can be viewed as the free part of the action, S 0 , and the remaining ones -as an interaction. The last term in (18) cancels diagrams with insertions of simple fermion loops in the σ -lines. As a consequence, in the leading order the propagator of the σ field is given by the inverse kernel L ∆ . It this work we fix the constant C(∆ ), which is arbitrary save the condition C(0) = 1, by the requirement for the propagator D σ to have the following form The parameter ∆ should be considered as a regularization parameter. The divergences in the correlators appear as poles in ∆ and are removed by adding the counterterms to the action (18). The renormalized action takes the form Though the model is renormalizable, the renormalization is not multiplicative, i.e. S R (q, σ ) = S (q 0 , σ 0 ). Due to the nonmultiplicative character of the renormalization the anomalous dimensions of the fields and operators cannot be, in general, determined via the corresponding renormalization factors. However, as it was shown in [22] this problem appears starting from 1/n 3 order only. Up to the order 1/n 2 the anomalous dimension of the operators can be obtained via the renormalization factors as follows. First, one has to modify the σ -propagator by multiplying it by some parameter u, D σ → uD σ and determine the renormalization factors calculating the corresponding diagrams with the modified propagator. Let O i be a system of operators mixing under renormalization, where Φ = {q,q, σ }, and O B i are bare operators. Up to the order 1/n 2 the anomalous dimensions for the operators can be obtained as follows [22] We use this formula for the calculation of the anomalous dimensions of the operators in question.

Correction exponents at 1/n 2
In the perturbative expansion the correction exponents to the scaling are related to the scaling dimensions of the operators σ 4 and σqq. The scaling dimensions of these operators are ∆ a = 4 + O(1/n). Let us first identify the corresponding operators in the 1/n expansion. There exist three scalar operators of dimension four, The last one is a total derivative and can be neglected. Note here that the operator σqq in the 1/n expansion has the dimension 4 − 2ε + O(1/n). Therefore it does not mix with the above operators and has no relation to the correction exponents 1 . Thus we need to find the anomalous dimension matrix for the operators O + = σ 4 and O − = σ ∂ 2 σ . It is easy to see that the diagonal entries γ 11 , γ 22 and γ 12 of the anomalous dimension matrix are of order 1/n. At the same time the anomalous dimension γ 21 is of order 1/n 3 . It happens because the only diagram which contributes to this matrix element at order 1/n 2 , σ 3 → σ ∂ 2 σ , see Fig. 1, does not diverge. It can be checked by an explicit calculation. Also, the finiteness of this diagram is guaranteed by the conformal symmetry. The arguments are exactly the same as in the case of σ 2 → ∂ 2 σ transition in the nonlinear σ model. An interested reader can find detailed discussions in Refs. [23,26]. Thus at the order 1/n 2 one can neglect mixing and calculate only the diagonal matrix elements.
We will use either upper or lower index for the coefficients of the 1/n expansion trying to keep notations close to generally accepted. For the reader's convenience we collect 1/n 2 expressions for the ba-

Operator σ
Here we present our result for the anomalous dimension of the operator σ . The mixing affects the anomalous dimension of this operator only starting from 1/n 3 . Analyzing the contributing diagrams it is easy to notice that the anomalous dimension for the operator σ can be written in the form where C k are the binomial coefficients. The anomalous dimension γ σ and γ σ 2 are known with 1/n 2 accuracy [3][4][5][6], see Appendix A. The last contribution starts from 1/n 2 . The diagrams contributing to ∆ γ σ 3 are shown in Fig. 2. Their calculation is rather straightforward so we will not dwell on it. We obtained the following expression for ∆ γ σ 3 ∆ γ Here ψ(x) is the Euler ψ function and an explicit expression for the index η 1 = −B(µ)/2µ can be found in Appendix A. Thus for the correction index ω + we find The explicit expression for the 1/n 2 correction, ω + , is rather lengthy and will not be given here. Instead we present an expansion of the index ω + around µ = 2 − ε up to ε 4 terms ω (0) in complete agreement with the results of the four-loop calculation [18].

Operator σ ∂ 2 σ
In this section we discuss diagrams contributing to the anomalous dimension of the operator O − = σ ∂ 2 σ which, as usual, can be split in two parts In the leading order in 1/n only the two diagrams shown in Fig. 3 contribute to γ. The first diagram has the symmetry coefficient C sym = 2. The calculation of both diagrams is rather straightforward and leads to the following expression where we introduced the notation α ≡ µ − 1. Taking into account that γ that agrees with the result of [19] derived with the help of the self-consistency equations approach. At the order 1/n 2 there are in total 23 different diagrams to be calculated. Namely, Before presenting our final answer for γ we briefly discuss the calculation of the diagrams. -In order to calculate the operator correction diagrams it is sufficient to calculate the 1/n order diagrams up to finite, O(∆ 0 ), terms. For the sum of the two diagrams we obtained 2 D(∆ ) = u 2 D(∆ )/∆ , where where we accepted the notations of [25] B The answer for the operator correction diagrams (D OC ) reads -Instead of calculating separately each diagram of the second group, that is not quite simple, we use the result of [22] where it was shown that the contribution of all such diagrams to the anomalous dimension can be extracted from the 1/n order diagrams with dressed propagators and vertices. We discuss details of the calculation in Appendix B . -The single box diagrams shown in Fig. 5 as well as double box diagrams, Fig. 7, have only a superficial divergence. After integration over the operator insertion coordinate one gets a propagator type diagram, The ∆ pole arises after the Fourier transform. For the anomalous dimension it is sufficient to know A k (∆ ) at ∆ = 0. Since in coordinate space the diagram is finite at ∆ → 0, one can put ∆ = 0 from the beginning. It greatly simplifies calculations since one can use the star -triangle relation for the basic σqq vertex, shown in Fig. 6. Fig. 6 The star-triangle relation. The indices α, β , γ satisfy the uniqueness condition α Next, multiplication by x 2 makes the diagrams logarithmically divergent, i.e. the residue of the ∆ pole does not depend on the external momenta. Thus one can use the freedom to change the momentum flow through a diagram to facilitate calculations [27]. For example, to calculate the second diagram in Fig. 5 it is convenient to choose the momentum flow (after a multiplication by x 2 ) through the vertices connected to the crossed circle (operator insertion) and use the star-triangle relation for the upper σqq vertex. After this transformation the diagram is calculated by consecutive application of chain integration rules.
-The calculation of the first three double box diagrams which follows the same pattern is more involved and one has to use many tricks like the integration by part [28], the tetrahedron transformation [29], etc. A review of the relevant multiloop calculation techniques can be found in [25]. In principle the last double box diagram could be calculated in the same way. However it is more convenient to apply an inversion transformation first, simplify the numerators and use the star-triangle relation for the bosonic triangle. As the result one gets a sum of three-and two-loop diagrams whose calculation is more or less straightforward.
Finally, collecting contributions from all diagrams (see Appendix C) we obtain for γ 2 Thus for the second correction exponent we get Substituting the expressions for γ σ and γ and expanding around µ = 2 − ε up to ε 4 terms one obtains that agrees with the results of [18].

Summary
We have calculated the correction exponents ω ± in the GN model with 1/n 2 accuracy. These exponents are related to the slopes of the beta functions in the critical GNY model. Our results are in complete agreement with the expressions for the perturbative four loop beta functions in the GNY model obtained recently in [18]. For the calculation we used the method developed in [20][21][22]. This method allows one to use the standard RG technique and, at the same time, to resum effectively a certain subset of diagrams reducing the total number of diagrams to be calculated. Values for each individual diagram are presented in Appendix C and will be useful for calculations of the corresponding exponents in the chiral versions of the GN model.

Appendix A: 1/n 2 indices
In this appendix we collected the results for the basic indices η, κ = 2χ and λ related to the anomalous dimensions of the fields q, σ and the composite operator σ 2 that can be found either in the original papers [3][4][5][6][7] or in the book [25]. The leading order coefficients read and the 1/n 2 coefficients take the form We remind here that α = µ − 1, η = ∑ i η i /n i , etc.

Appendix B: SE and Vertex correction diagrams
It was shown in [22] that the contribution to an anomalous dimension coming from diagrams with self-energy insertions and vertex corrections to a LO diagram can be obtained by evaluation the LO diagram with dressed propagators and vertices. The dressed propagators have the form 3 The expressions for dressed propagators given in the Appendix of ref. [24] correspond to a different normalization condition for the propagator of the σ field, and the expression for the fermion propagator contains a typo. and up to O(1/n 2 ) terms The dressed vertex is shown in Fig. 8.
Its form is fixed by the conformal invariance χ = κ/2 is defined in (A.1) and for the factor Z we derived In order to obtain the contribution to an anomalous dimension from 1/n diagrams, 4 see Fig. 3, with all possible SE insertions and vertex corrections (SEV) one has to replace all propagators and vertices in the 1/n diagrams by the dressed ones and introduce a regulator ∆ in one of the lines, e.g. to shift the index of the fermion propagator, ∆ q → ∆ q − ∆ . The corresponding diagrams have only a superficial divergence. Then the contribution to the anomalous dimension reads [22] where r i is the residue at the ∆ pole of the dressed diagram, D i (∆ ) = r i /∆ + O(∆ 0 ) and the sum goes over all such diagrams (in the case under consideration i = 1, 2). Proceeding in this way one reduces the number diagrams to be calculated. Moreover, individual vertex correction diagrams appear to be more complicated: all of them have a divergent subgraph. The presence of the regulator ∆ breaks the uniqueness condition for the basic scalar-fermion vertices that rather complicates a calculation. On contrary, all Fig. 9 An example of the calculation of the dressed diagram.
vertices in a dressed diagram obey the uniqueness condition and a diagram can be transformed with the help of the startriangle relation. In order to further simplify the calculation one can dress only part of the propagators and vertices. Indeed, the residues r i depend on two parameters, η and χ, which we will, for a while, consider as independent parameters. Since r i is needed up to 1/n 2 terms one can write r(η, χ) = r(0) + ηr η (0) + χr χ (0) + . . . , (B.10) where r η = ∂ η r, etc. It means that up to 1/n 2 terms the calculation of the diagram D(η, χ) is equivalent to the calculation of two simpler diagrams: D(η, χ = 0) and D(η = 0, χ). Note, that the parameter χ plays the role of a regulator in the vertex V . It can be checked that in the limit χ → 0 the vertex reduces to the point-like basic vertex with unit coefficient. Thus the diagram D(χ = 0, η) coincides with the 1/n diagram with modified propagators and point-like vertices. Moreover, taking into account that contributions from dressed vertices and propagators to an anomalous dimension are additive (at the order 1/n 2 ) one can, in order to simplify a calculation, modify only part of the vertices and propagators. A more general discussion of this technique can be found in [30] while here we illustrate it on the example of the calculation of the right diagram in Fig. 3.
The corresponding dressed diagram is shown on the l.h.s in Fig. 9 -the fat lines stand for the dressed propagators and the grey triangles for the dressed vertices. One needs to find the corresponding residue, r(χ, η). Since we restrict ourselves to the 1/n 2 accuracy the corresponding result can be extracted from another three diagrams, see Fig. 9. Namely, The last diagram which is D(0, η) depends only on η. The sum of the previous two diagrams gives D(χ, 0)+D(0, 0) up to O(1/n 3 ) terms. The white triangles stand for the dressed vertex V | η=0 and lines without indices -for the bare propagators. Let us start the analysis from the last diagram, which we denote by D 3 (a, b). Since the indices a = µ − 1 + η/2, b = 1 − η obey the uniqueness condition 2a + b = 2µ − 1 the diagram can be easily calculated. For the corresponding residue we get and With the help of the star-triangle relation the second (fourloop) diagram can be transformed (up to a prefactor) into the diagram D 3 (a χ , b χ ), where a χ = α + χ and b χ = 1−2χ, that gives for the residue r 2 where H(χ) = a 2 (1 − χ)a(2µ − 3 + 2χ) and The calculation of the first diagram follows the same lines.
Using the star-triangle relation one can show that It results in the following expression for r 1 where the factor E(χ) = D 3 (α, 1 + χ)/D 3 (α, 1) + O(1/n 2 ) has the form The calculation of the first diagram in Fig. 3 is much simpler. The final answer for both diagrams with self-energy insertions and vertex corrections are given in Appendix C.
The diagrams in Fig. 4 give the following contributions For the single box diagrams we obtain: Finally, our results for double box diagrams (Fig. 7) are ∆ µα(2α − 1)η 1 − 2α 2 − 2α + 1 (α − 1)(α + 1) Here B = B 2 3 −C 3−µ . The contribution to the anomalous dimension from a diagram D a reads a is the simple pole residue, D a = D a /∆ + . . .. Finally we give the contributions to the anomalous dimensions from the self-energy and vertex correction diegrams where γ SEV 1 , γ SEV 2 correspond to the left and right diagrams in Fig. 3, respectively.