Generalized Gross-Neveu Universality Class with Non-Abelian Symmetry

We use the large $N$ critical point formalism to compute $d$-dimensional critical exponents at several orders in $1/N$ in an Ising Gross-Neveu universality class where the core interaction includes a Lie group generator. Specifying a particular symmetry group or taking the abelian limit of the final exponents recovers known results but also provides expressions for any Lie group or fermion representation.


Introduction
One aspect of quantum field theory that has important applications to Nature is the study of fixed points of the renormalization group functions. These are defined to be the non-trivial zeros of the β-function. Using the location of a fixed point one can compute the values of the renormalization group functions there to produce renormalization scheme independent expressions known as critical exponents [48,49,50,51]. These quantities govern the dynamics of the phase transitions in a material. Indeed accurate measurements of the exponents experimentally as well as the symmetry properties of a material, can equally guide one to the underlying quantum field theory or spin model describing the dynamics. One property is that more than one theory can be a valid description of a phase transition. For instance, both a continuum quantum field theory as well as a discrete spin model with a common symmetry can be valid tools to provide numerical exponent estimates. The equivalence of both theoretical techniques at a fixed point is known as universality [50,51]. Away from a phase transition each theory will have different properties and be inequivalent.
One theory that has risen to the fore in this context in recent years is that of the Ising Gross-Neveu model [1,19]. This is primarily due to the belief that it underpins a particular phase transition in graphene. This material is made up of a sheet of carbon atoms arranged in a hexagonal lattice. When the two dimensional sheet is stretched it can undergo a transition from an electrical conductor to what is termed a Mott-insulating phase. On the theoretical side the Ising Gross-Neveu model can be supplemented with quantum electrodynamics (QED) to describe aspects of other phase transitions. Equally if the basic or Ising Gross-Neveu model is endowed with extra symmetries it contributes to the understanding of transitions in other materials. For example, what is termed the chiral Heisenberg Gross-Neveu model is an extension of the Gross-Neveu model to include an SU(2) symmetry. It is the effective theory for electrons This paper is a contribution to the Special Issue on Algebraic Structures in Perturbative Quantum Field Theory in honor of Dirk Kreimer for his 60th birthday. The full collection is available at https://www.emis.de/journals/SIGMA/Kreimer.html on a half-filled honeycomb lattice where there is a phase transition between an anti-ferromagnetic insulating phase and a semi-metallic one [2,20,36]. Its criticality properties were studied in [22]. More recently a variation of this version of the Gross-Neveu model, called the fractionalized Gross-Neveu model, has been developed [35]. It has a novel spectrum that differs from those of other Gross-Neveu models and has an associated SO(3) symmetry. What is clear from this set of emerging variations of the Gross-Neveu universality class is the common theme of the core interaction being modified to include a non-abelian symmetry. In this respect it is completely parallel to the extension of QED to include a non-abelian symmetry that equates to quantum chromodynamics (QCD) or Yang-Mills theory with fermions in the fundamental representation of the Lie group responsible for colour symmetry. The only difference with the Gross-Neveu class of theories is in the Lorentz structure of the core interaction. As QCD has been studied at length using a general Lie group symmetry rather than specifying SU(3) at the outset, which is the group that governs the strong interactions, it seems sensible to develop a programme to calculate in the Gross-Neveu model with a parallel non-abelian symmetry. Then the properties of the various physical applications can be deduced by specifying the appropriate group parameters. This is the main task here. We will consider a generalized Gross-Neveu universality class with a non-abelian symmetry and calculate the critical exponents of the theory. This will be achieved by using the critical point large N formalism pioneered in [42,43,44] for the nonlinear O(N ) σ model. Here N would correspond to the number of quark flavours in the analogy with QCD. The elegance of the approach is such that we can deduce the critical exponents in d spacetime dimensions as a function of the non-abelian symmetry group Casimirs. The fixed point associated with the formalism is the Wilson-Fisher fixed point in d-dimensions [50]. As the exponents are renormalization group invariants their expansion near 2 and 4 dimensions, where measures the difference in these values from d, will agree with the perturbative evaluation of the same exponents in the respective quantum field theories of the universality class. More usefully since the dimension of interest for the materials application is three, one can determine several terms of the 1/N series for each exponent. These provide reasonably accurate estimates for relatively low N when compared to results from other techniques. However, the benefit of taking the general non-abelian universality class approach is that estimates will be readily available if a phase transition with a new symmetry is discovered.
The article is organized as follows. Relevant background concerning the generalized Gross-Neveu universality class is given in Section 2 together with the basic critical point large N formalism. Subsequently in Section 3 we solve the Schwinger-Dyson equations at criticality at O 1/N 2 to produce the fermion anomalous dimension. Various calculational tools that are necessary for this are reviewed as well. To provide the groundwork for finding the next order of this exponent, the anomalous dmension exponent of the bosonic field is determined at O 1/N 2 in Section 4. One of the other basic exponents in critical systems is that relating to the correlation length behaviour and we determine it at O 1/N 2 in Section 5. Equipped with these results, the large N conformal bootstrap formalism at criticality is applied in Section 6 to deduce the fermion anomalous dimension at O 1/N 3 . We review our results in Section 7 and provide concluding remarks in Section 8.

Background
To begin with we recall the Lagrangian of the chiral Heisenberg Gross-Neveu-Yukawa theory is [20] L cHGNY = iψ iI ∂ /ψ iI + 1 2 ∂ µπ a ∂ µπa + g 1π aψiI T a IJ ψ iJ + 1 24 g 2 2 π aπa 2 which is renormalizable in four dimensions, where the three couplings g 1 , g 2 and g 3 are dimensionless. This is a generalization of the Lagrangian studied in [53] and is in the chiral Heisenberg Gross-Neveu model universality class. The renormalizability dimension is also termed the critical dimension of the theory. The scalar-fermion interaction includes the group generator T a of the Lie algebra and the indices take values in the ranges 1 ≤ i ≤ N , 1 ≤ I ≤ N c and 1 ≤ a ≤ N A , where N is the number of flavours of massless fermions and N c and N A are the respective dimensions of the fundamental and adjoint representations of the symmetry group.
We note that in [20] the specific group considered was SU (2). Within our ultimate critical exponents the generators will manifest themselves through various colour Casimirs such as C F and C A defined by where f abc are the structure constants. The scalar fieldπ a plays a subtle role in the construction of the large N expansion but in four dimensions it corresponds to a fundamental propagating field. The main aspect of the large N critical point formalism of [42,43,44] is that in the approach to criticality at the Wilson-Fisher fixed point the dynamics are driven by the core interaction of the universal quantum field theory. For (2.1) this is the cubic interaction together with the fermion kinetic term. These two terms determine the canonical dimensions of both fields by ensuring the action is dimensionless in d-dimensions. In effect the universal Lagrangian at criticality is where the quadratic term inπ a is necessary for large N renormalizability. We will omit the labels i and I for brevity from now on. We say in effect since at criticality there is no coupling constant in the sense it is conventionally used in perturbation theory. So the critical point universal Lagrangian that will be the foundation for applying the large N critical point formalism of [42,43] is where we have rescaled the scalar fieldπ a to introduce π a . This Lagrangian (2.3) is renormalizable in d-dimensions in the large N formalism [33,40,41,53], where 1/N is the dimensionless ordering parameter since the perturbative coupling constant is absent at criticality. Ensuring that the d-dimensional Lagrangian (2.3) has a dimensionless action means that ψ has canonical dimension 1 2 (d − 1) while that of π a is unity. For (2.3) this implies that g is dimensionless in two dimensions after eliminating the auxiliary field π a producing This is similar to the Ising Gross-Neveu model discussed in [34] and to see the equivalence, one takes the abelian limit of (2.4) by replacing the group generators with the unit matrix. This is completely parallel to taking the abelian limit of QCD to produce QED. The dimensionality of π a at criticality plays a key role in the connection of the universal theory and the Lagrangian of (2.1).
In the latter the three couplings are dimensionless in four dimensions similar to the effective coupling of the 3-point interaction of (2.3). Therefore (2.3) would be strictly non-renormalizable in four dimensions and the quadratic term in π a would have a dimensionful coupling which would be the mass. Instead the standard kinetic term and quartic π a interactions of (2.1) would be relevant. In other words we term the quartic interactions to be spectator interactions that would be active solely in four dimensions. Moreover underlying the first two terms of the universal Lagrangian (2.3) there are an infinite number of Lorentz scalar operators built from combinations of π a and its derivatives. A finite subset of these extra operators would become relevant in even integer dimensions and act as interim spectators in the infinite tower of renormalizable quantum field theories that connect to the universal theory in the neighbourhood of their respective critical dimensions. In this context it is worth noting that the quartic spectator interactions of (2.1) play a role in determining the full fixed point structure of the four dimensional theory. One of these fixed points though will correspond to the Wilson-Fisher fixed point of the generalized Gross-Neveu universality class considered here. However, only a perturbative evaluation of the four dimensional theory's renormalization group functions would determine which one it is. While this is beyond the scope of the present article it is likely to be a solution where the critical values of all three couplings are non-zero. One non-trivial check in establishing such a connection though rests in the agreement of the expansion of the various critical exponents with their large N counterparts that will be determined here. More concretely we now summarize the key aspects of the large N critical point formalism for (2.3). In the approach to the fixed point the propagators have the following asymptotic behaviour in coordinate space [11] ψ(x) ∼ Ax / where the name of the corresponding field is used. The dimensionless quantities A and C are the coordinate independent amplitudes that will always occur in the combination y = A 2 C from the 3-point interaction. The next to leading order terms in (2.5) that involve the exponent λ are called the corrections to scaling. Here λ will be identified with the correlation length exponent ν through 1/ν = 2λ. In addition to the canonical dimension the two fields have anomalous contributions and the respective full dimensions of ψ and π a are where we use d = 2µ for shorthand [43] and η and χ are the fermion field and vertex anomalous dimensions respectively. For applications in condensed matter problems the dimension of π a that is conventionally used is When λ corresponds to the correlation length exponent its canonical dimension will then be taken to be (µ − 1). In this respect the leading and next to leading order terms of (2.5) then have different dimensions which is the reason for the second set of independent dimensionless amplitudes A and C . Each of the exponents that we will compute as well as y will depend on N , µ and the Lie group Casimir invariants. The dependence on the former means that each entity has a Taylor series in powers of 1/N that is formally given by for η and y for example and we will determine the first three terms of η for (2.3). These general considerations cover the basic formalism for the technique introduced in [42,43,44]. To determine all bar η 3 we can apply the original method [42,43] that was used for the Ising Gross-Neveu universality class in [8,13,14,15,39,45]. This required solving the skeleton Schwinger-Dyson equations for the ψ and π a 2-point functions. So the scaling forms of these are needed given that (2.5) represents the critical point behaviour of the propagator. However one definition of the 2-point function is that it is the inverse of the propagator in momentum space and this mapping can be carried out through the Fourier transform. Using the convention given in [42,43] which is we transform (2.5) to momentum space carry out the inversion and then apply the inverse Fourier transform. This results in the coordinate space 2-point function asymptotic scaling forms which are [11] (2.10) The presence of the function a(α) in the Fourier transform produces a complicated dependence on µ and the exponents, leading to the compact functions .
While the large N conformal bootstrap formalism of [44] has its origins in these asymptotic scaling functions and was applied to the Ising Gross-Neveu model in [8,15,39], the extraction of an expression for η 3 derives from the scaling behaviour of the 3-point function. We defer to a later section for the required technicalities of that formalism.

2-point Schwinger-Dyson equation
Equipped with the asymptotic scaling forms of the full propagators, (2.5) and (2.10), which represent the behaviour at criticality, we use them to solve the Schwinger-Dyson equations. In conventional perturbation theory one systematically renormalizes the divergent n-point Green's functions in a renormalizable theory order by order in perturbation theory. This principle is respected in the large N technique of [42,43] except that the ordering of graphs in the n-point functions is achieved by the variable 1/N which is dimensionless across all spacetime dimensions unlike the perturbative coupling constant. For (2.3) the first few terms in the respective 2-point functions of the fields are given in Figure 1, where the dotted line represents the fermion and the wiggly line denotes the π a field. The two loop graphs are the O 1/N 2 corrections to the one loop ones. The counting of powers of N arise from closed fermion loops giving a factor of N and the π a field. The expansion of the amplitude variable begins at O(1/N ) and this translates into each π a line in Figure 1 carrying a power of 1/N . One key point worth noting concerns the lack of dressing of lines with self-energy corrections. Contributions from such graphs are already accounted for in the inclusion of a non-zero anomalous dimension in the power of the asymptotic scaling forms.
At leading order the two equations of Figure 1 equate to where we have included the respective group theory factors which derive from the properties of T a and Tr In this coordinate space representation the one loop graphs require no evaluation. This is because one integrates over the coordinate of the internal vertices. As the one loop graphs have external vertices the corresponding terms of (3.1) are the products of the propagators. In this leading order instance any integration has been effected in the derivation of the scaling forms for the full 2-point functions. The algebraic representation of Figure 1 is (3.1) and to O(1/N ) it contains two unknowns. Using (2.7) and the 1/N expansion of r(α − 1) and p(γ) these are η 1 and y 1 and moreover at this order they occur linearly. Thus eliminating y 1 between the equations of (3.1) produces At next order the situation is not as straightforward due to the additional graphs of Figure 1 being divergent which necessitates the introduction of renormalization constants. The formalism for this was provided in [40,41] and requires the introduction of a regularization which is achieved through the shift [42,43] where ∆ is a small parameter. In effect it equates to an analytic regularization of the propagators and we emphasize that the spacetime dimension d does not play any role in the regularization in contrast to coupling constant perturbation theory. Consequently the extension of (3.1) to the next order has to account for this and so the algebraic representation of Figure 1 becomes in coordinate space where Z V is the vertex renormalization constant. It has the Laurent expansion where the residues are after expanding in powers of 1/N . We follow [40,41] and restrict to the MS scheme. In (3.3) the x 2 dependence arises from the dimensionality of the integrals in the regularized Lagrangian.
As it stands the various factors prevent the limit to criticality from being taken smoothly. Moreover the factors associated with the one loop graphs of Figure 1 will give contributions at O 1/N 2 from the expansion of x 2 χπ . This will produce problematic logarithms but these are connected with the simple poles of the values of the two loop graphs denoted by Σ 1 and Π 1 . In particular they have the formal structure where Σ 1 and Π 1 are finite. They were computed previously in [11], where the explicit d-dependent values are available. We note our trace conventions at this stage are the same as [11] and we use 2 × 2 γ-matrices. To adjust for higher dimensional γ-matrix representations one simply redefines N using where d γ is the dimension of the γ-matrix representation.
x y 0 ≡ a(α)a(β−1)a(γ−1) The key integration tool for evaluating the graphs in [11] is shown in Figure 2 and is termed uniqueness or conformal integration for the scalar-Yukawa interaction. There are several ways to establish the relation provided in Figure 2. One is to use Feynman parameters. In that derivation the final integration is over a hypergeometric function and it cannot proceed unless the uniqueness condition of α + β + γ = 2µ + 1 is fulfilled. Setting this allows the final integration to be completed which produces the factor on the right side of Figure 2. A more elegant alternative is to apply a conformal transformation to the integral which is which implies the mapping for instance. The consequence is that when applied to strings of contracted γ-matrices the transformation does not alter the initial string of γ-matrices. In the application to the left hand side of the equation of Figure 2 the exponent of the scalar becomes (2µ + 1 − α − β − γ). Setting this to zero allows the z-integration to proceed resulting in the expression on the right hand side after undoing the initial conformal transformation.
With the availability of the counterterm from the vertex renormalization constant the divergences are removed minimally. However ln x 2 terms remain through two contributions in each Schwinger-Dyson equation. One is from the power of x 2 2∆ in the O 1/N 2 correction. The other arises from the factor x 2 χπ that is present in the one loop graph. Expanding this in powers of 1/N the O(1/N ) term involves ln x 2 . To be able to take the x 2 → 0 limit safely means that the as yet undefined χ π 1 has to be suitably chosen. Doing so to ensure there are no ln x 2 terms in each Schwinger-Dyson equation at O 1/N 2 requires from the explicit values for the residues which satisfy L 1 = − 2K 1 and implies the same value results for both equations. This finally renders the algebraic representation of Figure 1 finite as well as ensuring that it is scale free. Since the two equations have two unknown variables η 2 and y 2 that appear linearly, eliminating the latter leads to where we use the shorthand which involves the Euler ψ function defined by ψ(z) = d ln Γ(z)/dz. 4 π a critical exponent at O 1/N 2 Having established the fermion critical exponent at O 1/N 2 the next stage in the large N formalism is to determine the same quantity for the boson field. In this instance from the definition (2.6) this requires the vertex anomalous dimension at O 1/N 2 . While χ π 1 followed as a corollary to ensuring the 2-point function was finite in the approach to criticality, in order to proceed to the next order to find χ π 2 by the same method is too intractable. Indeed evaluating the analogous exponent in other models has not been achieved that way. Instead a more direct approach suffices which is to examine the scaling behaviour of the 3-point vertex in the critical limit. In other words the O 1/N 2 graphs illustrated in Figures 3 and 4 are computed.
In practical terms the diagrams are more straightforward to evaluate in momentum space than in coordinate space. This is primarily due to simplifications in the application of the uniqueness rule of [11]. However one can connect the values of graphs in both the coordinate and momentum space evaluations through the Fourier transform (2.8). Underlying this one needs to use the momentum space forms of the asymptotic scaling functions which are It is the inverse Fourier transform of these that produce the leading terms of (2.10). With (4.1) it is straightforward to evaluate the graph of Figure 3 and determine an expression for χ π 1 from the scaling behaviour. It is equivalent to (3.6). The key point is that the same procedure of [40,41] can now be applied to determine χ π 2 . This requires in part the evaluation of the O 1/N 2 corrections illustrated in Figure 4, where the fermions can be routed around the closed loop in both directions. The values of the diagrams in the absence of any group theory considerations were given in [13]. With the presence of the group generator in the vertex of (2.3), obtaining the associated group factor of each graph is straightforward in most cases using (2.2) for example and the definition of the Lie algebra which in our conventions is However for the graphs where there is a closed fermion loop with four π a fields attached a higher group Casimir is present. In particular the group factor associated with each graph will contain the fully symmetric rank 4 tensor among others. To treat these so called light-by-light graphs we have used the color.h routine that accompanies the symbolic manipulation language Form [46,37]. The package allows one to manipulate group theory factors associated with Feynman graphs and write them in terms of Casimirs. It is based on the comprehensive analysis provided in [38]. However, rather than use color.h to determine the group factor solely for these light-by-light graphs we have applied it to all the graphs treated throughout this article for consistency. We note that the rank 3 fully symmetric colour tensor d abc given by also occurs in graphs that contain d abcd F . In the final expressions for the exponents, however, these are absent as a result of a cancellation between graphs where the fermions are routed around the closed loop in different directions. In addition to the O 1/N 2 corrections of Figure 4 there are contributions to χ π 2 from the graph of Figure 3. These arise from the correction to the variableỹ which isỹ 2 as well as the parts from the 1/N expansion of the exponents of the propagators in (4.1). In addition one has to include the vertex renormalization constant Z V in the O(1/N ) graph as this cancels the subgraph divergences in the first three graphs of the first row in Figure 4. This cancellation is necessary in order to ensure the approach to criticality is smooth. Assemblying all these components allows us to deduce the vertex critical exponent at the next order which is where Θ(µ) = ψ (µ) − ψ (1), and the contributions from the light-by-light graphs is evident. We close this section by noting that the d abcd F tensor will appear in the same combination as it does in (4.2) in the perturbative renormalization group functions of (2.1) from Feynman diagrams involving interactions with the couplings g 1 and g 3 .

Correlation length exponent at O 1/N 2
Having established the dimensions of the two fields at O 1/N 2 using the leading term of the asymptotic forms of the propagators in the approach to criticality, it is possible to study the corrections to scaling. These are contained in both (2.5) and (2.10) corresponding to the terms involving the coordinate independent additional dimensionless parameters A and C . The extra exponent λ can be regarded as any exponent but to access the correlation length exponent ν we set λ = 1/(2ν) which has the canonical dimension of (µ − 1). To determine the 1/N corrections in d-dimensions to this exponent requires a consistency equation that extends (3.1) and (3.3). To find λ 1 we substitute the various asymptotic scaling functions into the Schwinger-Dyson equations for the 2-point function which produces the representation where we have omitted the factor of Z 4 V in the respective O 1/N 2 and O(1/N ) correction terms. The counterterms from these only come into effect at the next order.
Unlike the computation of η 1 we have included the two loop graphs of Figure 1, where the correction to scaling is included. These are denoted by Σ 1φ and Π 1φ , where φ indicates that the insertion is on either a ψ or π a line according to whether it is A or C respectively. The reason why these all have to be included in principle resides in the leading order N dependence of the 2-point scaling functions. For the two key combinations that appear in the correction to scaling Schwinger-Dyson equation we note that this dependence is [8,45] In fact the constant of proportionality of the latter is the combination (λ 1 − η 1 − χ π 1 ). As η 1 and χ π 1 are both available this leaves λ 1 as the unknown we seek. These terms will form part of the consistency equation that determines λ to O 1/N 2 and emerges from decoupling the x 2 λ terms in (5.1) and (5.2) which follows on dimensional grounds. The resulting two equations are Alternatively the relevant terms that produce an expression for λ 1 with respect to the large N expansion due to (5.3) can be written as a matrix M, where Examining the (2, 2) element both terms are the same order in 1/N as y is O(1/N ). Setting det M = 0 produces the consistency equation for λ 1 which can be solved to give To proceed to the next stage and find λ 2 the higher order 1/N correction graphs have to be added in to the two Schwinger-Dyson equations that produced the O(1/N ) consistency equations. However the ones we neglected to determine (5.4), since their N dependence in the determinant is an order lower that the contribution from Π 1C , now have to be included. These are Σ 1A , Σ 1C and Π 1A while the graph corresponding to Π 1C has to be expanded to the next order in 1/N since there is N dependence in the propagator exponents. To ease the extraction of the expansion of the consistency equation determinant we formally set to clarify this. The main work however resides in including the final contributions to find λ 2 which are illustrated in Figure 5. While the individual d-dependent values have been recorded in [14] for instance, we have had to append the respective group theory factors. Again we have used the Form color.h routine due to the presence of the light-by-light diagrams. Repeating the exercise of setting the determinant of the consistency equations to zero at the next order in 1/N produces the expression We have introduced the additional shorthand notation 6 Large N conformal bootstrap The final part of our study follows a different tack by applying the large N conformal bootstrap programme developed in [44], based on the early insights of [7,28,29,30,31]. In general terms the focus is initially directed towards the 3-point vertex of (2.3) and its behaviour in the critical region. The fields will still obey the asymptotic scaling forms of (2.10) but in treating the Green's functions in the bootstrap formalism not only are there no self-energy corrections on the propagators but there are no vertex corrections. In effect the primitive graphs are the building blocks and are illustrated in Figure 6. In that figure the dotted vertices do not correspond to the vertex of the Lagrangian (2.3). Instead they denote the presence of a Polyakov conformal triangle [30] which includes all vertex corrections at criticality. It is defined in Figure 7 for the general Yukawa type interaction that includes the one of (2.  The graphs of Figure 6, however, are the lowest order contributions to the full vertex function that we will denote by V (ȳ, α, γ; δ, δ ), where the last two arguments are regularizing parameters. These are required in the derivation of one of the two consistency equations defining the large N bootstrap formalism [29,30,31,44]. The first equation represents Figure 6 and is whereȳ is similar to the early amplitude combination y but includes the normalization of the Polyakov conformal triangle of Figure 7. Indeed (6.1) is responsible for determining the terms in the 1/N expansion ofȳ once the first few orders of η have been reproduced in this formalism. This is because the explicit expressions are necessary to extract η 3 from the third order term of the second bootstrap equation which is We note briefly that the regularizations δ and δ that appear here arise because of singularities in the 2-point Schwinger-Dyson equations when all the vertices are replaced by conformal triangles. In other words it was recognised in the original work of [7] that in the absence of any regularization the 2-point functions with dressed vertices would be finite overall. However each of the contributing diagrams were individually divergent. To accommodate this, and similar to the introduction of ∆ earlier, the vertex anomalous dimension has to be continued in a parallel way to (3.2). This is illustrated in Figure 8 for the leading order contribution to V (ȳ, α, γ; δ, δ ), where we have set χ π = 2∆ π for shorthand. This figure indicates the values of the exponents of the internal lines of the Polyakov triangle. Moreover the appearance of both δ and δ on the external legs of the graph on the left hand side indicate the addition of the regularizations to the exponents of the respective fields. This is reflected internally in the conformal triangles in the right hand graph as each of the vertices have to be unique even when there is a regularization. Given the form of (6.2) we can rederive η 2 from knowledge of the value of the graph of Figure 8. This was determined in [15], using the technique given in [44], for the Ising Gross-Neveu model as a function of the exponents of that theory and that result can be translated to (2.3). If we expand V (ȳ, α, γ; δ, δ ) in large N in the same notation as (2.7) then to the order that will eventually be necessary to evaluate η 3 we recall [15] where the order symbol indicates that terms cubic in any combination of the parameters are neglected. The functions B z and C z are defined by [44] in terms of the Euler ψ function. With (6.3) and setting V = V 1 /N in (6.1) and (6.2) it is straightforward to expand both bootstrap equations to O(1/N ) and verify the earlier expressions for η 1 and η 2 .
Having established the formalism reproduces available results the extension to the next order to find η 3 requires several steps. The first is to compute the next term in the 1/N expansion of (6.3) as that will contribute to the O 1/N 2 part of (6.2). However contributions from all the graphs of Figure 6 bar the one loop one have to be determined and included as they correspond to V 2 . The relevant parts of these graphs were calculated in [15]. By this we mean that a computational shortcut was used akin to that used in [44]. From examining the terms in the formal Taylor expansion of (6.2) in powers of 1/N the part involving V 2 occurs in the combination In [15] the contribution to (6.4) from each of the higher order correction graphs in Figure 6 were determined. Indeed in most cases only the value for the difference in derivatives could be found. Appending the group theory values from the color.h routine allows us to finally extract η 3 . We find Again the pieces arising from the light-by-light graphs can be clearly identified. In addition to the Euler ψ function and its derivatives a new function arises which is related to the function I(µ) defined in [44]. It corresponds to the derivative of a two loop self-energy graph, where a regularizing exponent of one of the propagators is differentiated. In [24] this master graph was expressed in terms of an 4 F 3 hypergeometric function where the regularizing parameter appears in the parameter arguments of the function. A compact expression for I(µ) that is valid in d-dimensions was recorded in [4]. Here we have set so that Ξ(µ) has an expansion involving multiple zeta values [3,24,44]. For instance, In strictly three dimensions [44] I 3 2 = 2 ln 2 + 3ψ 1 2 2π 2 and the expansion around three dimensions is known up to ten terms [5].

Results
We focus in this section on general aspects of the critical exponents of (2.3) that we have determined. One of the main reasons for considering a generalized universality class was the fact that known results for specific models could be extracted as well as be of use where a different Lie group underlies the physics problem. To assist with that we have collected electronic expressions in an attached data file. One aspect of the results that needs to be stated is that we have checked that the d-dimensional exponents are in agreement with several models where the results were determined directly. For instance, the original Ising Gross-Neveu model of [48] corresponds to the abelian limit of the symmetry group by specifying while the Mott insulating phase [2,20] that corresponds to taking the symmetry group to be SU (2) takes the values For the more recent application of (2.3) to the fractionalized Gross-Neveu model discussed in [35] the respective values are For each of these cases the expansion of the exponents near four dimensions with d = 4 − 2 are in full agreement with known three and four loop perturbative results [9,21,23,27,32,35,52,53]. In the case of the Ising Gross-Neveu model exponents these also are in accord with two dimensional perturbation theory [10,12,16,18,25,26,47]. Moreover taking the limits for the three cases, all the large N exponents agree with previous work [8,11,13,14,17,32,34,39,45]. One advantage of the arbitrary group approach in d-dimensions means that the structure of the exponents can be studied in various representations. For example if we restrict the fermions to be in the adjoint representation A whence C F = C A and T F = C A . In that case we find The group valued coefficient of the terms involving ζ 3 and π 2 ln 2, which derive from I 3 2 , has an interesting combination of Casimirs. Indeed there might be instances of this factor being zero for certain Lie groups. However, we have computed the value of C 4 A N A − 24d abcd A d abcd A for all the classical and exceptional Lie groups and found that it is always non-zero.

Discussion
As the Ising Gross-Neveu universality class is central to a number of phase transitions in various materials, we have examined a generalized version of the underlying quantum field theory that incorporates the respective condensed matter systems. The key aspect is that the core interaction is endowed with a non-abelian symmetry that has a parallel in gauge theories. There the gauge interaction of QED is extended from an abelian to a non-abelian one to produce QCD by the inclusion of the generators of a Lie group thereby endowing QED with a colour symmetry. The similar extension of the Ising Gross-Neveu model is simpler in some respects. One obvious one is the absence of gauge symmetry. A benefit, however, is that considering (2.3) at the outset means results for specific phase transitions can be quickly deduced by specifying the Lie group. Indeed if a phase transition were discovered in a material that was in the same universality class as the Ising Gross-Neveu model but possessed a new symmetry other than the specific examples we have noted here, then information on the exponents can readily be deduced from our results. Throughout we have focussed on the application of the critical point large N technique developed in [42,43,44] to determine d-dimensional critical exponents. The advantage of this is that results are available for the renormalization group functions of the four dimensional quantum field theories in the same universality class too. By the same token the large N exponents contain a wealth of information on the structure of the same functions. For instance, coefficients in the anomalous dimension beyond the first few known loop orders can be accessed at successive orders in 1/N . This is particularly useful in that our O 1/N 2 and O 1/N 3 exponents can reveal where the new colour group Casimirs, such as d abcd F d abcd F , appear. In indicating the parallel of the QED to QCD generalization, examining the large N O 1/N 2 exponents in this universality class, albeit with a simpler vertex structure, does provide useful insight into what to expect in the calculation of critical exponents in QCD at O 1/N 2 . We have to qualify this comment by noting that while there are similarities, in the QCD large N critical exponent computation for ν for instance, there will be more graphs to consider than those of Figure 5. This is because in (2.3) Feynman diagrams with subgraphs involving three π a lines connecting to a fermion loop are zero after taking the γ-matrix trace. In QCD this would not be the case due to each vertex adding an extra γ-matrix to the trace. While such graphs remain to be computed the associated group theory factor that would result should not involve any Casimir higher than d abcd F d abcd F .