Large $N$ critical exponents for the chiral Heisenberg Gross-Neveu universality class

We compute the large $N$ critical exponents $\eta$, $\eta_\phi$ and $1/\nu$ in $d$-dimensions in the chiral Heisenberg Gross-Neveu model to several orders in powers of $1/N$. For instance, the large $N$ conformal bootstrap method is used to determine $\eta$ at $O(1/N^3)$ while the other exponents are computed to $O(1/N^2)$. Estimates of the exponents for a phase transition in graphene are given which are shown to be commensurate with other approaches. In particular the behaviour of the exponents in $2$ $<$ $d$ $<$ $4$ is in qualitative agreement with a functional renormalization group analysis. The $\epsilon$-expansion of the exponents near four dimensions are in agreement with recent four loop perturbation theory.

One of the more remarkable fundamental quantum field theories is the Gross-Neveu or Ashkin-Teller model, [1,2]. It is a purely fermionic theory with a sole quartic self-interaction. In some sense it is the parallel to the more widely studied purely scalar quartic field theory which is renormalizable in four spacetime dimensions. The spontaneous symmetry breaking phase in that model is the basis for the Higgs mechanism of the Standard Model. By contrast the O(N ) Gross-Neveu model is asymptotically free, [1], but appears to have little physical application since it is only renormalizable in two dimensions. Put another way prior to the discovery of the W and Z vector bosons of the Standard Model the physics of the weak interactions was described by an effective field theory involving 4-point fermion interactions. The restriction to being effective meant that such interactions could only be reliable as a model of Nature up to a specific momentum scale. However, the general Gross-Neveu class of quantum field theories, several of which were introduced in [1], have enjoyed a renaissance over recent years. This is in the main due to the fact that certain phase transitions in graphene can be described by specific universality classes based on the Gross-Neveu model, [1], where the classes derive from the underlying symmetry of the transition. For instance, stretching a graphene sheet can bring about a transition from a conductor to a Mott-insulating phase, [3,4]. It has been suggested that the physics of this transition can be described by what is termed the chiral Heisenberg Gross-Neveu model, [5,6,7,8,9]. Therefore there has been activity in computing and estimating the fundamental critical exponents of the transition in this model. Prior to 2014 there were virtually no deep computations where the theoretical status can be appreciated in the right panels of Figures 1, 2 and 3 given in [9].
More specifically in [9] the critical exponents η, η φ and 1/ν were estimated in spacetime dimension d where 2 < d < 4 using functional renormalization group techniques. Hence values were given for the three dimensional case of physical interest to graphene. By contrast at that time only two loop expansion results from four dimensions were available to compare with, [10,11]. Since then this four dimensional perturbative work has been extended to three and four loops respectively in [12,13,1]. The situation can be compared to the more widely studied Ising Gross-Neveu model which was given in the left panels of Figures 1, 2 and 3 of [9]. By contrast with the chiral Heisenberg Gross-Neveu model in the Ising case there are now two dimensional expansion estimates to four loops, [15], as well as four loop four dimensional perturbative information, [11]. These higher loop results postdate the lower loop results of [1,16,17,18,19,20,21] which were the state of the art at the time of [9] for the Ising Gross-Neveu universality class. On top of this, Monte Carlo data, [11], as well as estimates from several orders in the large N expansion, [22,23,24,25,26,27,28], were available to give independent analyses of the three dimensional exponents. Overall there was solid agreement for the exponents η φ and 1/ν from all these methods as well as the functional renormalization group analysis itself which was given in [9]. However those from the expansions near two and four dimensions for η were not in close agreement with the other methods, [9]. One simple reason for this is clear from the plot for η in Figure 1 of [9]. It is because η has to vanish in the limits down to two and up to four dimensions separately. So significantly higher orders in as well as knowledge of the asymptotic behaviour of the series would be required even to get close to the three dimensional estimates from the other methods.
A deeper underlying reason for this rests in the various theories which populate a universality class. Viewing the class as driven by a fundamental interaction at the Wilson-Fisher fixed point in d-dimensions, in the neighbourhood of an even spacetime dimension there is a quantum field theory which is renormalizable at that critical (even) dimension. Clearly it would not be renormalizable at another critical dimension. However a different theory would be relevant at that critical dimension but which will have the same fundamental interaction. In the case of the Ising Gross-Neveu universality class the theory which is renormalizable in two dimensions is the Gross-Neveu model of [1]. By contrast in four dimensions the next theory in the tower of theories at the Wilson-Fisher fixed point is what is termed the Gross-Neveu-Yukawa model, [29]. The respective Lagrangians differ only in the terms involving a scalar field. While we noted earlier that the Gross-Neveu model is purely fermionic, the quartic self-interaction can be rewritten via a scalar auxiliary field σ to produce a σψ i ψ i interaction where ψ i are the fermions. This interaction is common to the Gross-Neveu-Yukawa model but the σ field now has a canonical kinetic term and a scalar quartic self-interaction, [29]. It is this field σ which has anomalous dimension η φ . Structurally the Gross-Neveu-Yukawa model has many similarities to the weak sector of the Standard Model itself. As noted in [9], for instance, this has led to the hope that certain phase transitions in graphene could mimic various symmetry breakings in the Standard Model itself such as chiral or spontaneous symmetry breaking. Currently experimental data from such graphene transitions are far from being available but the potential to have a simple laboratory to test Standard Model related phase transitions is a fascinating prospect. Other variations of the Gross-Neveu universality class are therefore becoming important to study so that precise theoretical predictions will be available.
Given the current sparcity of theoretical critical exponent estimates for the chiral Heisenberg Gross-Neveu model it is crucial to bring the analysis up to the same level of precision as that of the Ising Gross-Neveu universality class. By this we mean that the amount of precise data of the two classes as represented in Figures 1, 2 and 3 of [9] should be the same for the different techniques used. This is the purpose of this article where we will compute the critical exponents η, η φ and 1/ν to the same order in the large N expansion as those in the Ising universality class. This will be achieved by applying the critical point large N formalism developed by Vasiliev et al in [30,31,32] for the O(N ) nonlinear σ model or equivalently the O(N ) φ 4 theory. Both these theories lie in the same universality class and the exponents of [30,31,32] are expressed as functions of d. That original approach was later extended to the Ising Gross-Neveu universality class in a series of articles, [22,23,24,25,26,27,28], and it is on a subset of these papers that the large N computations presented here are based. For instance, for the most part we will analyse the basic 2-point functions of the chiral Heisenberg Gross-Neveu theory by solving the skeleton Schwinger-Dyson equations algebraically in the limit as one approaches the ddimensional Wilson-Fisher fixed point. While this will determine η, η φ and 1/ν to O(1/N 2 ), the original method given in [32] pointed the way to determining η at O(1/N 3 ). This was extended to the Ising Gross-Neveu class in [24] and independently in [28]. Therefore, we will apply what we now term the large N conformal bootstrap method to the chiral Heisenberg Gross-Neveu universality class. Given modern usage of the term conformal bootstrap method to mean a largely numerical technique based on solving the crossing symmetric n-point functions with manifest conformal symmetry, we have appended the term large N to the naming of the earlier technique given in [32]. The work of [32] was inspired by the direct three dimensional conformal bootstrap construction of Parisi, [33], upon which the modern bootstrap, [34,35,36,37], is also based.
The advantage of the approach of [30,31,32] is that the expression for η at O(1/N 3 ) as well as the other exponents at O(1/N 2 ) are determined as exact functions of d. It is important to appreciate that this arbitrary spacetime dimension is not the same as the one used for dimensional regularization in perturbation theory. In [30,31,32] the skeleton Schwinger-Dyson equations are analytically and not dimensionally regularized. Therefore the exponents derived in the large N method correspond to the exponents of the true theory at the Wilson-Fisher fixed point in arbitrary (non-integer) dimensions. Having said this the expansion of the large N exponents around two or four dimensions are in exact agreement with the expansion of the perturbative critical renormalization group functions, evaluated at the fixed point, of the underlying theory with a critical dimension of two or four respectively. One benefit of having the exponents as functions of d is that we will be able to plot their behaviour in 2 < d < 4 dimensions which will provide an independent insight into whether the dependence on d is consistent with that given by the functional renormalization group approach of [9] which used a sharp regulator. While three dimensional estimates are clearly and ultimately of physical interest the non-integer dimension structure of a field theory in effect is indicating the properties of the underlying universality class. In the context of the Standard Model these simple Ising Gross-Neveu and chiral Heisenberg Gross-Neveu classes may be pointing us to a novel way of viewing four dimensional physics. For instance, the observation that supersymmetry can emerge at a fixed point in a multicoupling quantum field theory, [38,39,40], may be directing us to a new era of beyond the Standard Model analyses. In other words the couplings of the relevant operators with Standard Model symmetries may flow to a fixed point which is accessed by an effective quantum field theory in the short term. However this could be superseded by a new one in much the same way as the Lagrangian introducing the W and Z vector bosons was eventually constructed. Therefore, exploring this new chiral Heisenberg Gross-Neveu universality class and understanding it in the graphene context offers new exciting avenues to explore.
The article is organized as follows. The following section introduces the chiral Heisenberg Gross-Neveu universality class together with the basic large N formalism which will be used to extract the critical exponents to several orders in d-dimensions. In section 3 the exponent η is determined at O(1/N 2 ) and the evaluation of η φ and 1/ν to the same order are provided in sections 4 and 5 respectively. While these computations in essence follow the same critical point approach using the skeleton Schwinger-Dyson equations, the computation of η at O(1/N 3 ) uses the large N conformal bootstrap method which is outlined in section 6. The analysis of our results is given in section 7 where estimates are shown for the case of graphene and plots of the large N exponents in d-dimensions are given for the Ising, XY and chiral Heisenberg Gross-Neveu models. Finally, we give our conclusions in section 8.

Background.
The chiral Heisenberg Gross-Neveu model (cHGN) is a generalization of the original two dimensional O(N ) or SU (N ) Gross-Neveu model of [1] but where the 4-point interaction is modified to include the three SU (2) Pauli spin matrices σ a where a = 1, 2 and 3. The Lagrangian is where 1 ≤ i ≤ N and g is the dimensionless coupling constant. Throughout our large N computations we will use the spinor trace convention that TrI = 2 rather than 4 which is ordinarily used when extending four dimensional perturbation theory to three dimensions. In other words to circumvent this convention we can rewrite our definition of N using where d γ corresponds to the trace convention of the γ-matrices, [9]. The theory (2.1) is only perturbatively renormalizable in two dimensions. In four dimensions a quartic fermion interaction is non-renormalizable but (2.1) can be reformulated in terms of an auxiliary fieldπ a In this formulation one can connect the original quartic interaction in the two dimensional theory with one in four dimensions which we will term the chiral Heisenberg Gross-Neveu-Yukawa (cHGNY) theory which has the Lagrangian The sharing of the common Yukawa interaction in (2.3) and (2.4) is related to the fact that these theories are in the same universality class at their respective Wilson-Fisher fixed points. The bosonic sectors are different but they play a role in ensuring each Lagrangian is renormalizable in their respective critical dimensions. While both theories lie in the same universality class there is an underlying universal theory within which one can carry out computations in the large N expansion. Specifically d-dimensional critical exponents can be determined to several orders in 1/N , when N is large, using methods based on [30,31,32]. These exponents encode all orders perturbative information on renormalization group functions in d-dimensions and hence contain information on all theories, such as (2.3) and (2.4), in the universality class in their critical dimension. The Lagrangian for the application of the large N method to the underlying universal theory is where the fieldπ a has been rescaled so that there is no coupling constant in the interaction. The method developed in [30,31,32] relies on the fact that at criticality it is the interaction which drives the dynamics at all scales. In comparison to (2.4) there is no quartic interaction. This is not an issue since the structure of the universal theory as such generates the requisite information corresponding to the critical renormalization group functions of (2.4) via π a 4-point subgraphs in Feynman diagrams. This has been elucidated in [40] in the context of the large N f expansion of Quantum Chromodynamics (QCD) where N f is the number of massless quark flavours. Equally for other theories in the tower of this universality class the higher π a n-point subgraphs will contain the vertex structure of the higher n-point interactions relevant to the theory in a particular higher critical dimension.
The large N critical point formalism of [30,31,32] exploits properties of the renormalization group at criticality to provide a method of solving the Schwinger-Dyson equations algebraically. In particular in the approach to a fixed point there is no scale aside from the correlation length. So, for instance, the propagators take an asymptotic scaling form whose behaviour is governed solely by the full dimension of the field. Specifically in coordinate space the scaling forms of the fermion and boson fields of (2.5) are, [22], at leading order in the limit as x → 0 where we use the name of the field for each form. The quantities A and C are x-independent amplitudes and α and γ are the full dimensions of the actual fields in the d-dimensional Lagrangian. These comprise the canonical dimension and the anomalous dimension. The former is deduced from ensuring that the dimension of each term in the Lagrangian is consistent with the action being dimensionless and we define in parallel with [22] where η is the anomalous critical exponent of the fermion. The anomalous dimension of π a includes η and is derived from the Yukawa interaction where χ π is the anomalous dimension of the Yukawa vertex itself. Here we use the convention of [30,31] that the spacetime dimension d is written as d = 2µ for shorthand. In order to be able to compare with the critical exponent of the analogous bosonic field of [12,14] we note that the corresponding exponent, η φ , is related to that of the π a field here by In terms of the connection with the renormalization group functions, η corresponds to the fermion anomalous dimension evaluated at the Wilson-Fisher fixed point. As there is only one interaction the amplitudes A and C always appear in the same combination which we will denote by y = A 2 C. It like all the critical exponents will only depend on µ and N at criticality. Hence they can be expanded in powers of 1/N via and a similar notation will be used for other exponents. In (2.9) and similar expressions one can restore other spinor trace conventions by replacing N with the definition (2.2). As we will be providing various critical exponents we need to recall the background formalism for those. While (2.6) gives the dominant behaviour of the propagators in the approach to criticality corrections to scaling can be included and for the propagators these take the form, [22], where an independent exponent λ governs the scaling of the correction with A and C the respective associated x-independent amplitudes. In principle λ can be any exponent but for the moment we will let it correspond to the critical slope of the β-function of the theory whose critical dimension in this universality class is two and hence corresponds to (2.1) or (2.3). The connection with the exponents of [14] is that 1/ν = 2λ and we note that the leading order canonical dimension of λ is λ 0 = µ − 1.
One of the main methods to determine the explicit values of the exponents to several orders in 1/N is to systematically solve order by order the skeleton Schwinger-Dyson 2-point function for both fields in the approach to criticality. This will also include the corrections to scaling. In order to do this one needs the asymptotic scaling forms of the respective 2-point functions. These are derived from the propagator forms by inverting the momentum space propagators. To do this we use the Fourier transform given in [30,31] which is This produces the coordinate space 2-point function asymptotic scaling forms which include the scaling corrections. The amplitudes differ from the respective ones in the propagator forms and the various functions of the exponents are given by . (2.14) In summary we have introduced the basic structure of the propagators and 2-point functions for the large N critical point analysis. This is sufficient to solve the skeleton Schwinger-Dyson equation at O(1/N 2 ) and determine d-dimensional expressions for η 1 , η 2 , χ π 1 , λ 1 and λ 2 . The values of χ π 2 and η 3 are deduced from solving 3-point functions at criticality and the formalism for each will be discussed later. 3 Evaluation of η 2 .
As aspects of the solution of the 2-and 3-point functions for each of the various exponents are common we illustrate these features in more detail for the derivation of η 2 . The first stage for this is to represent the skeleton Schwinger-Dyson 2-point functions of Figure 1 as algebraic equations. In coordinate space to O(1/N 2 ), which is sufficient to deduce η 1 and η 2 , we have where we have substituted for the asymptotic scaling forms of the inverse propagators. The values of the respective two loop integrals Σ 1 and Π 1 are present and these not only depend on N through the presence of the exponents on the internal lines but also on the regularization ∆. This is introduced into the critical theory by shifting the anomalous dimension of the vertex via, [30,31], corresponding to an analytic regularization. Both graphs are divergent and can be represented by where the O(∆) terms of the integrals would only become relevant for the η 3 derivation using this approach. As discussed in [30,31] and [41] the core vertex is renormalized and this feature is present via the vertex renormalization constant Z V which also depends on N and ∆. It can be written as a double expansion via where the counterterms m ln are then expanded in powers of 1/N Next an overall common factor of the length of the position vector x 2 has cancelled in (3.1). However, from the dimensionality of the Feynman graphs there is a remnant x 2 dependence deriving from the vertex anomalous dimension and the regularization. As it stands (3.1) is a formal representation of Figure 1 but is not applicable in the scaling region, given by x 2 → 0, where it ought to be independent of the length scale x 2 . Equally the renormalization constant has to be determined as well as χ π . It turns out that all these issues are resolved together. To achieve this we formally substitute (3.3) and expand the factors involving x 2 to the appropriate orders in 1/N and ∆ neglecting terms O(1/N 3 ) and O(∆) respectively. To ensure that the algebraic Schwinger-Dyson equation is finite when the regularization is removed leads to and the absence of ln x 2 terms requires We have expressed these conditions in terms of the simple pole of the graph Π 1 rather than Σ 1 . However similar relations emerge from the ψ i Schwinger-Dyson equation and these are not inconsistent since from explicit computation we have which ensures there are no ambiguous expressions for the vertex counterterm and χ π . The resulting Schwinger-Dyson equations are now finite and independent of x 2 to the order in 1/N necessary to deduce η 1 and η 2 .
To proceed one first concentrates on the leading order terms of both equations and eliminates y 1 to produce Setting the leading order value for γ and α = µ + η 1 /N produces With this value established as well as that for the amplitude combination y 1 then η 2 is deduced by including the values of the finite parts of the two loop graphs and eliminating the unknown y 2 . First, with from [22] we find which is needed for the next term in the expansion of p(γ). Then using from [22], where Σ 1 satisfies the same ratio with Π 1 as (3.8), we find where we use the shorthand notation and ψ(z) = d ln Γ(z)/dz is the Euler ψ function.
With the determination of the exponents of the last section we have obtained the critical behaviour of the ψ i -field to O(1/N 2 ) and that of π a to O(1/N ). Since the asymptotic scaling forms of the inverse scaling functions, when there is correction to scaling, involve γ then in order to find λ at O(1/N 2 ) we need to find χ π 2 first. One way is to determine the corrections to the skeleton 2-point functions at next order in 1/N . However this would involve a large amount of computation which would not be necessary since there is an alternative. The extraction of χ π 1 from ensuring there are no terms involving ln x 2 in the limit to criticality is reminiscent of a similar way of extracting the renormalization constant for a mass m by renormalizing that part of a 2-point function in conventional perturbation theory which corresponds to the wave function renormalization. By this we mean that in a renormalizable theory the one loop mass counterterm can be deduced from the two loop 2-point wave function computation by ensuring that there are no non-renormalizable terms involving ln m 2 . While this is not the conventional way to extract such a one loop counterterm from a two loop evaluation it can be regarded in one way as a useful shortcut but more importantly as being consistent with the underlying renormalizability. The same process has in effect been played out in the finding χ π at O(1/N ) through the computation of η at O(1/N 2 ). While our discussion in the previous section followed the original algorithm outlined in [30,31], the critical point large N formalism was subsequently put in a parallel context to conventional perturbation theory in [41,42]. The core aspects of perturbation theory are an ordering of the diagrams constituting a Green's function and a regularization to facilitate their evaluation. With the ordering, such as the power of the coupling constant, the renormalization constants are introduced by multiplicative rescaling and determined with respect to a subtraction criterion. The situation developed in [41,42] is completely parallel. Graphs are ordered with respect to the counting given by the power of 1/N with the regularization introduced by shifting the vertex anomalous dimension, [30,31]. One major difference is that the residues of the poles in the renormalization constants are functions of the spacetime dimension d which does not play a regularizing role. In this process the critical exponents are directly extracted in a renormalization group invariant way. In introducing the vertex renormalization constant Z V previously we were in effect following the 1/N renormalization formalism developed in [41,42] which demonstrated that the ordering of graphs with a suitable regularization in effect was a complete parallel to conventional perturbation theory. Moreover it was shown in [41,42] that the vertex anomalous dimensions, as well as operator renormalization, could be extracted from the critical point evaluation of n-point functions for n ≥ 3. Therefore to deduce χ π 2 we follow that approach. The first stage is to reproduce the value of χ π 1 which requires the evaluation of the graph of Figure 2. In this instance we use the momentum space asymptotic scaling forms for the propagators which are where p is the momentum andÃ andC are the amplitudes in momentum space. However the combinationỹ =Ã 2C will always arise in the graphs. Its O(1/N 2 ) value is known via the Fourier transform (2.11) and we already havẽ where from earlier calculations we can deduce The approach to find the corrections to χ π is similar to that of the 2-point function. Once the graphs have been evaluated in momentum space to the finite part in ∆ the contribution to the respective term of the exponent in 1/N is given by the coefficient of the ln p 2 part. Such terms have to be absent in the limit to the critical point which then fixes the unknown exponent. Therefore following the prescription given in [41] we find the same expression for χ π 1 as before.
The computation to deduce χ π 2 is more involved. The main part is the inclusion of the O(1/N 2 ) graphs of Figure 3. These were evaluated in [23] for the parallel calculation of the vertex critical exponent in the O(N ) Gross-Neveu model. So for (2.5) the values of the master integrals are appended with the corresponding group theory factor. In addition there are next order contributions from the one loop graph of Figure 2. This is because of the presence of η and χ π in the power of each of the propagators. In addition there are vertex counterterms from Z V for each vertex which have to be included in the same graph. They relate to the subgraph divergences arising from the first three graphs of the top row in Figure 3. Piecing all the relevant contributions together and isolating the ln p 2 part we find where Θ(µ) = ψ (µ) − ψ (1) . 5 Evaluation of λ 2 .
The next stage in the evaluation of the large N critical exponents again mimics that of conventional perturbation theory in that with the wave function anomalous dimensions at O(1/N 2 ) one can establish the β-function to the same order. As noted earlier this is achieved by considering corrections to scaling and evaluating the corresponding 2-point Schwinger-Dyson equations of Figure 1 to the next order. There are various aspects of extracting the expansion for λ which is not straightforward. Basic features are best illustrated by considering the equation for ψ i which can be represented by prior to renormalization. The two loop graph denoted by Σ 1 in Figure 1 has been expanded to include the values where there is a correction to scaling in the ψ i and π a fields. These are Σ 1A and Σ 1C respectively and their values have been computed in [27]. The key difference with the inclusion of the corrections is that terms do not have the same dependence on x 2 . This means that the equation decouples into two pieces. One is the piece which contributes to the determination of η while the other involves the unknown λ and is Unlike the equation which determines η the associated correction amplitudes are present linearly. If one ignores the contribution from the two loop graphs then the first two terms are the relevant ones for finding the value of λ 1 . However there is a new feature in comparison with the corresponding equation to determine η. This is to do with the large N dependence in that the first term is O(1) while the second one is O(1/N ) which has an impact on the structure for the π a Schwinger-Dyson equation contributing to λ.
Turning to the second equation of Figure 1 we can represent it in a similar way to that of the first equation and have The values of the two loop graphs with corrections to scaling on the ψ i and π a lines are Π 1A and Π 1C respectively. Again the equation decouples into that which determines η and which is formally similar to (5.2) and also includes the correction amplitudes A and C . Between these two equations at leading order the only unknown is λ 1 . The consistency equation which is formed to solve for it is deduced by first representing (5.2) and (5.4) as a matrix M given by after renormalization where indicates the finite part of Π 1C with respect to ∆. It is worth noting that we have not retained all the terms in the Schwinger-Dyson equations in the matrix. This is because we want to concentrate on finding the value for λ 1 initially. Therefore we have retained the leading terms except for Π 1C which would ordinarily be regarded as next order in the same way that the parent graph only contributed to η 2 and not η 1 . The reason why this graph cannot be omitted resides in the powers of N in the matrix, [27], and the fact that the consistency equation is given by requiring det(M) = 0. Therefore in computing the determinant the two resultant terms have to be the same order in 1/N . The complication arises from the scaling functions since Therefore both terms of the (22) element are the same order in 1/N and the omission of Π 1C would lead to an incorrect value of λ 1 . Therefore setting we have the leading order consistency equation to leading order. With, [27], we find Having outlined in detail the formalism and issues behind the derivation of the leading order consistency equation for λ 1 we extend the analysis to the next order in 1/N . Virtually all aspects of the Schwinger-Dyson equations discussed so far are sufficient to formally deduce the correction to the consistency equation defined by det(M) = 0. For instance, the contributions from the two loop graph of the ψ i equation in Figure 1 need to be included. However, for the π a equation there are two main additions. The first is that the appearance of Π 1C at leading order means that the next term in its large N expansion needs to be included as noted in (5.7) and the correction has the value, [27], Also as another direct consequence of the presence of Π 1C at leading the O(1/N 2 ) corrections to the π a Schwinger-Dyson equation have to be included. These are illustrated in Figure 4 where the dot on a π a line denotes which line has the correction to scaling term. The values of the integrals were given in [27] but in this case the group factors have to be appended. We use a labelling of the graphs which is parallel to that used in [27]. Moreover their contribution has to be included in the extension of (5.4). If we denote the sum of the contributions from the graphs of Figure 4 by Π 2C then the extension of (5.5) is Summing all the explicit contributions to Π C2 we have With this value together with the corrections to the asymptotic scaling functions it is straightforward to solve det(M) = 0 at O(1/N 2 ) and find We note that the coefficients in the Taylor expansion of this function near two dimensions, together with Ψ(µ) and Θ(µ), all involve ζ n for integer n ≥ 3 where ζ z is the Riemann zeta function.
Finally we close this section by briefly mentioning that we have repeated the procedure at leading order to find the critical exponent which relates to one part of the β-functions of (2.4). To do this instead of using the correction to scaling of (2.10) we use the alternative correction where ω 0 = µ − 2. The consistency equation for ω 1 is formally the same as that for λ 1 except that λ is replaced by ω. The same reordering at leading order occurs and the analogous value to Π 1C1 is required. Denoting this by Π 1C1 ω we note that, [43], and find As in [43] this corresponds to one of the eigen-critical exponents of the 2 × 2 matrix of derivatives of the two β-functions of (2.4). While the O(1/N 2 ) corrections are known to the (2.4) βfunctions, [44], the method used in that approach differed from that used here. Specifically the relevant critical exponents were determined by examining 3-point functions. So the values of the master integrals computed in [44] cannot be immediately translated to the extension of our 2-point Schwinger-Dyson equation to next order in 1/N . 6 Large N conformal bootstrap.
The provision of λ 2 completes the evaluation of the three basic exponents η, η φ and 1/ν to O(1/N 2 ) in the large N expansion. The next stage is to proceed to O(1/N 3 ) which is possible in the case of η. However this is not by the evaluation of the two 2-point function Schwinger-Dyson equations. While in principle one should be able to extend the η 2 calculation it transpires that the evaluation of the higher order graphs is not straightforward. Instead we follow the method developed in [32] based on the earlier work of [33] which we will term the large N conformal bootstrap method. In this approach one analyses the skeleton Schwinger-Dyson equation for the 3-point vertex using the dressed propagators (2.6) but with additionally no vertex subgraphs. The focus therefore is on what would be called the primitive graphs in the 3-point function.
These are illustrated in Figure 5 where there is a difference in the representation of the vertices in comparison with the earlier Figures. The dot at each vertex represents what is termed a conformal triangle, [32,33]. In other words the original vertex is replaced by a one loop triangle graph where the exponents of the new internal edges are at this stage arbitrary. Their values are fixed by the criterion that each new vertex in the triangle is unique. A vertex is said to be unique if the sum of the exponents of the lines joining a 3-point vertex is equal to the spacetime dimension d. The concept of uniqueness was introduced in three dimensions in [45] and extended to d-dimensions in [30,31]. One consequence of uniqueness is that in applying a conformal transformation to any of the graphs of Figure 5 immediately reduces it to a 2-point function. While this means that as these graphs stand they can be computed, the question still remains as to how to extract η 3 . + + + + Figure 5: Graphs contributing to the vertex function which determines η 3 using the large N conformal bootstrap method.
The original procedure to carry this out for the nonlinear σ model in d-dimensions was given in [32] and involved extending the work of [33] which was specific to three dimensions. The first stage is to construct the consistency equations whose solution gives η 3 . This has been given in [24,25,46] for the O(N ) Gross-Neveu model and that construction straightforwardly translates to the case of (2.5). The only major difference aside from the different labelling of the exponents is to append the factors deriving from the Pauli matrices in the vertex. Therefore we focus on outlining the formalism. The key ingredient is the vertex function which is denoted by V (y, α, γ; δ, δ ). Here y is the same combination of amplitudes as before and the function depends on two regularizing parameters δ and δ . These are required since in the derivation of the equation for η given in [24,25,32,46] there are divergent 2-point functions. The divergences arise in the same context as those in the earlier 2-point function analyses which required the introduction of the analytic regularization controlled by ∆. As the conformal bootstrap also is a perturbative expansion in the vertex anomalous dimension the vertices of the graphs with conformal triangles also have to regularized. This is achieved by requiring that the external π a and one of the external ψ i legs of the graphs in Figure 5 have their dimension shifted by δ and δ respectively, [24,25,46]. Having outlined features of the vertex function the formalism presented in [24,25,46] leads to the consistency equations. First the vertex function V (y, α, γ; δ, δ ) is defined by the sum of δ-regularized graphs given in Figure 5. The basic equation which in effect determines the hidden amplitide of the vertex order by order in large N is 1 = V (y, α, γ; 0, 0) .
In other words the sum of all the graphs is unity and it turns out that each graph is finite when evaluated after applying conformal techniques. Therefore the regularization for this was unnecessary, [24,25,46]. However, the graphs which contribute to the value of η 3 and equally to the lower order exponents are divergent which means the regularization cannot be neglected in their determination. Consequently the second consistency equation of the set is where the same scaling functions as before are present. At leading order the right hand side is unity which means that the same value for η 1 is recovered as expected. At next order the one loop graph of Figure 5 has to be included. However, it has to be evaluated with the regularized external vertices and conformal triangles. In order to illustrate the earlier points about the graphs in the large N conformal bootstrap construction we have given the explicit allocation of exponents on the lines in Figure 6 for the full graph corresponding to the one loop graph of Figure 5. There for space considerations we have set It is clear that there is a non-zero sum of the exponents at the top and bottom left vertices which are proportional to δ and δ respectively. Equally the sum of all the exponents at internal vertices are unique. The value of the graph in Figure 6 has been determined in the context of the usual Gross-Neveu model as a function of the exponents of that model, [46]. So that formal expression can be adapted to this computation as the group factor is separate for the graph in Figure 6. Exporting that function and inserting it into (6.1) and (6.2) we recover the earlier value of η 2 which provides a useful check on our conformal bootstrap construction.
The final stage is to include the remaining graphs of Figure 5 in the two consistency equations. Again the values of the contributing integrals have been computed in [28] in the context of the usual Gross-Neveu model. So those master integral values need only be decorated with the group factors for the specific case we are interested in. Unlike the one loop graph of Figure 5 some of the higher order δ-regularized graphs cannot be evaluated completely. However for these cases it transpires that the difference can be determined. This is all that is necessary since this difference for the higher order graphs is sufficient to deduce η 3 in the 1/N expansion of the right hand side of (6.2). If we denote the contributions to the right hand side of (6.2) from all the graphs in Figure 5 except the one loop triangle by V 2 then This involves a new function Ξ(µ) which is related to a particular two loop graph introduced in [32] and denoted there by I(µ). In terms of I(µ) we have so that the expansion of Ξ(µ) near two dimensions only involves multiple zeta values, [47]. For instance, the first few terms are In three dimensions the integral I(µ) is known exactly, [32], since which will be needed for our three dimensional estimates later. Finally including V 2 in (6.2) and expanding the one loop contribution from Figure 6 to O(1/N 2 ) we find that which completes the evaluation of all the critical exponents. To assist with analyses we have provided an attached data file where electronic versions of all the exponents computed here are given.

Results.
We devote this section to discussing our results and give estimates of critical exponents in three dimensions. As a first stage we must indicate that all the expressions derived here are consistent with known perturbation theory near four dimensions. Recently the three and four loop renormalization group functions have been provided for the chiral Heisenberg Gross-Neveu-Yukawa theory, [12,14], which built on the early loop work of [10,11]. Setting d = 4 − 2 in the d-dimensional exponents we find Comparing with the results of [10,12,13,14] and allowing for the different expansion convention we find exact agreement with the exponents of [10,12,13,14]. For future five loop computations we have also given the O( 5 ) terms. We have also checked that the expansion of ω 1 near four dimensions is in agreement with one of the eigen-critical exponents of the Hessian defined by the derivatives of the two β-functions of [14]. This check is similar to that carried out in the Ising Gross-Neveu model [43].
Having established the consistency of our d-dimensional critical exponents with known four dimensional perturbation theory it is a simple exercise to determine the values in three dimensions. We have where the O(1/N ) correction to η φ is zero or numerically. Equipped with these three dimensional expressions we can now provide estimates for the exponents in the case which is of interest to graphene problems. Therefore recalling our convention for the spinor trace (2.2) the value of N we use in order to compare with the numerical estimates of critical exponents of other methods for (2.5) is N = 4. Given this we have determined numerical estimates for the three critical exponents using Padé approximants. These are given in the last line of Table 1 together with results quoted in [14] by other approaches for comparison. In the table bracketed values for 1/ν represent the value derived from ν which was computed directly. For the large N estimates we have quoted the [1,2], [0, 2] and [1,1] Padé approximants for η, η φ and 1/ν respectively. In general the values for large N are similar but larger than those of the functional renormalization group values of [48] for η φ and 1/ν. For η the situation is much different with virtually no overlap for any of the different analyses. As [14] also considered the XY Gross-Neveu model we have also provided the parallel results for N = 4 in Table 2 in order to compare results for the same methods in a different context. Again it is the case that the estimates for η φ and 1/ν are slightly larger than those of the functional renormalization group values in [14] with again no clear consensus for η. Although it is worth noting that the large N estimate for η is a [1,1] Padé since the value of η 3 for that model has not been computed. Perhaps a more instructive way of viewing the situation with the exponent estimates is to use a graphical representation motivated by the functional renormalization group approach [9]. One aspect of that method is that it is not limited to a discrete spacetime dimension. In other words the critical exponents can be determined as functions of d and plots of the three exponents in 2 < d < 4 were given in [9]. As our expansion parameter, 1/N , is dimensionless in d-dimensions it is possible to provide similar plots for the same exponents. By this we mean that we can plot the various Padé approximants used to obtain the estimates in Table 1 for the chiral Heisenberg Gross-Neveu model as functions of d. These are given in Figure 7 for N = 4 except that the [1,1] Padé is used for η. This is because we do not have a closed analytic form for the Ξ(µ) as a function of µ. Instead we have used the first two terms of η. However in the plot of η we have included the [2, 1] Padé estimate, indicated by an open circle, and the [1, 2] estimate indicated by a solid circle. These are meant to guide roughly where the d-dimensional line would intersect if η 3 was known as an analytic function in d-dimensions. One interesting aspect of the three curves in Figure 7 is that they are in good qualitative agreement with those of the right hand panels of Figures 1, 2 and 3 of [9]. By this we mean that the shape in terms of concavity and convexity for η φ and 1/ν are very similar as well as the offset of the peak for η which is not in the neighbourhood of d = 3. Finally, in order to assist this comparison we have included similar plots for the XY Gross-Neveu model when N = 4 and the ordinary or Ising Gross-Neveu model for N = 8 in Figures 8 and 9 respectively. The latter value of N in that case is the parallel one to compare with [9]. In Figure 8 only the first two terms of η were available which why the Padé estimate lies on the line. While the shapes of the plots for the Ising Gross-Neveu model are also qualitatively similar to those of [9] it is worth noting that the one for η φ differs in concavity to that for the chiral Heisenberg Gross-Neveu model which is consistent with the functional renormalization group approach.    We have completed the evaluation of the three core critical exponents η, η φ and 1/ν to several orders in the large N expansion as a function of the spacetime dimension d for the chiral Heisenberg Gross-Neveu universality class. The expansion of the expressions near four dimensions agrees exactly with the recent explicit four loop renormalization group functions of (2.4) at the Wilson-Fisher fixed point, given in [14]. Indeed our large N results are an important independent check on that work. Moreover we have provided the next terms in the series as a future check for any five loop renormalization of (2.4). As the exponents depend on d, plots of the Padé approximant of each exponent in dimensions 2 < d < 4 have been given in order to compare with the functional renormalization group approach of [9] where a sharp regulator was used. All the large N plots for not only the chiral Heisenberg Gross-Neveu universality class but the Ising and XY Gross-Neveu classes in d-dimensions are in good qualitative agreement with [9]. Again this consistency provides independent evidence that these methods are capturing the proper and general behaviour of the exponents of the universality class across the dimensions. At the outset we drew attention to Figures 1, 2 and 3 of [9] in relation to the status of results available for the Ising and chiral Heisenberg Gross-Neveu universality classes. In addition to the three and four loop results of [12,13,14] the results here are now edging towards the latter class having commensurate data with the former. What is lacking is Monte Carlo results and higher order two dimensional perturbative renormalization group functions. The latter should be possible to obtain to four loops with the recent derivation of the β-function for the Ising Gross-Neveu model in two dimensions, [15]. This is not as straightforward a computation as that for the four dimensional case, [12,13,14]. In two dimensions quartic fermion self-interactions are not multiplicatively renormalizable when the Lagrangian is dimensionally regularized. Instead additional evanescent quartic interactions are generated and their presence means that one has to be careful in extracting the true renormalization group functions after the regularization is lifted. Aside from this technical issue it should be possible in the future to add this extra information to the analysis of the chiral Heisenberg universality class.