Critical exponent $\eta$ at $O(1/N^3)$ in the chiral XY model using the large $N$ conformal bootstrap

We compute the $O(1/N^3)$ correction to the critical exponent $\eta$ in the chiral XY or chiral Gross-Neveu model in $d$-dimensions. As the leading order vertex anomalous dimension vanishes, the direct application of the large $N$ conformal bootstrap formalism is not immediately possible. To circumvent this we consider the more general Nambu-Jona-Lasinio model for a general non-abelian Lie group. Taking the abelian limit of the exponents of this model produces those of the chiral XY model. Subsequently we provide improved estimates for $\eta$ in the three dimensional chiral XY model for various values of $N$.


Introduction.
The Gross-Neveu model, [1,2], is one of the core quantum field theories that is on a similar level to scalar φ 4 theory since it too has a quartic self-interaction but is purely fermionic. By contrast although it is renormalizable in two space-time dimensions rather than four, more interestingly it is asymptotically free, [1]. It shares this feature with four dimensional non-abelian gauge theories such as Quantum Chromodynamics (QCD) [3,4] and, moreover, has the property of dynamical symmetry breaking leading to non-perturbative mass generation, [1]. In more recent years the original Gross-Neveu model of [1] and its extensions have been studied in the context of ultraviolet completion, [5,6]. This is where the two dimensional theory is extended to become a renormalizable one in four dimensions which is achieved in several stages. The first is to introduce a Hubbard-Stratonovich transformation which produces an auxiliary scalar field. In two dimensions this field becomes dynamical and corresponds to a bound state of two fermions, [1]. In the ultraviolet completion to four dimensions, however, this auxiliary field develops a fundamental propagator, [5]. So to ensure four dimensional renormalizability it is necessary for the scalar field to have a quartic self-interaction. This in essence describes the structure of a Gross-Neveu-Yukawa theory which is the ultimate point of the ultraviolet completion exercise, [5]. Indeed this particular theory has been of interest for many years due to its potential to model a composite Higgs field in extensions of the Standard Model, [7], for example.
Our description of this connection between two and four dimensions is primarily for the original Gross-Neveu model which has a discrete chiral symmetry. In this context it is sometimes referred to as the Ising Gross-Neveu model. It is not the only version of a Gross-Neveu model in that alternatively one can endow the basic quartic fermion interaction with a continuous chiral symmetry. In two dimensions this is known as the chiral Gross-Neveu model, [1], and that has an ultraviolet completion to what is termed the chiral XY model, [6]. This model as well as the Ising case have been the subject of considerable interest in the last decade or so due to its underlying connection with models of phase transitions in graphene. Stretching a sheet of carbon atoms can change the electrical properties of the material from a semi-conductor to a Mott insulator, [8,9]. Indeed the problem of evaluating the critical exponents describing the phase transition in the dimension of interest, which is three, has generated an immense amount of activity. The main theoretical approaches include the application of the functional renormalization group, [10,11], conformal bootstrap, [12,13], Monte Carlo methods, [14], the expansion, [15,16], and the large N expansion, [17]. We note that this is a representative set of recent articles rather than a definitive one. The connection between the expansion and the large N approach is quite close. For the former the renormalization group functions are calculated to high loop order in the ultraviolet completed four dimensional model. Then the critical exponents at the d = 4 − 2 dimensional Wilson-Fisher fixed point are computed before the expansion of an exponent is summed up prior to setting = 1 2 . In some Gross-Neveu universality classes the same exercise has been additionally performed in the original two dimensional theory. In that case the exponents are determined at the Wilson-Fisher fixed point in an¯ expansion near d = 2 + 2¯ dimensions even though the theory is only perturbatively renormalizable in strictly two dimensions. Despite these respective theories not being in the neighbourhood of three dimensions, accurate estimates for exponents have been extracted in certain cases. For instance, the current situation in the Ising Gross-Neveu model has been summarized in [16].
By contrast in the large N critical point approach pioneered in scalar O(N ) models in [18,19,20], the first three terms of critical exponents have been computed in d-dimensional space-time as a function of d. The large N expansion is renormalizable, [21,22], primarily due to 1/N being a dimensionless parameter. Using the first three terms in the series for relatively low N , estimates of exponents in three dimensions have been shown to be competitive with those of the other techniques mentioned. It is these large N exponents depending on d that drive the close connection with the exponents determined from the expansion of the explicit perturbative renormalization group functions. This is due to the general observation that the exponents depend on both variables N and . So expanding those in powers of 1/N and the expressions derived using either technique separately should be in agreement in the area where terms in the double expansion overlap. Indeed this is the case in many models. Moreover the large N exponents can be expanded to all orders in at each order in 1/N thereby providing non-trivial checks on future perturbative results.
While the main focus so far has been on the Ising and chiral Gross-Neveu models and their connections with graphene phase transitions, the chiral XY model itself has been explored in relation to certain phases in matter. For instance, the correlated symmetry protected phases were investigated in depth in [23]. As this model clearly is of importance to these areas it is worth summarizing the current status of the analytic evaluation of the exponents. For instance, the renormalization group functions of the chiral XY model are available to four loops, [6,15]. Equally the expansion of the large N critical exponents are known to be in agreement with those high order results. More specifically the exponents η, η φ and 1/ν corresponding to the fermion and scalar field dimensions as well as the correlation length exponent respectively are each known to O(1/N 2 ), [17,24]. However, for η this means only the first two terms in 1/N are available since the fermion has zero canonical dimension. However the large N conformal bootstrap formalism of [20] allows one to determine the O(1/N 3 ) term of η. This has been achieved for the Ising Gross-Neveu model in [25,26,27], as well as a variant known as the chiral Heisenberg Gross-Neveu model in [28] but not yet for the chiral XY case. This is due to a technical reason. In the large N conformal bootstrap technique, that is founded on the early work of [29,30,31,32], the expansion is in powers of the vertex anomalous dimensions. In particular in one of the amplitude variables the vertex anomalous dimension appears in the denominator. However in the chiral XY Gross-Neveu case the leading order large N vertex exponent is zero. So one cannot immediately effect the formalism at the outset. In reality one has a mathematically ill-defined quantity akin to 0/0. However since the application of the original formalism of [28,29] to the chiral XY model can indeed produce non-singularly derived expressions for η to O(1/N 2 ), [24], there ought to be a way of applying the general large N conformal bootstrap technology of [20] as well.
Given the need for having as accurate as possible data on the phase transition exponents for the chiral XY model we have found an indirect method to determine η at O(1/N 3 ) in ddimensions. This is the main topic of the article and it is worth outlining our strategy at the outset. Instead of considering the chiral XY model itself we will carry out a computation in a generalized theory which possesses a non-abelian symmetry. In essence this is the Nambu-Jona-Lasinio model, [33]. Critical exponents have been computed in the large N expansion to high order in several cases previously. For instance, the exponents of the fields for the group SU (N c ) have been determined at O(1/N 2 ) as well as that for the correlation length, [24,34]. A calculation of η at O(1/N 3 ) for the SU (2) Nambu-Jona-Lasinio model was given in [35]. However we will consider the case of a general Lie group rather than specific unitary ones. The reason for this is that then the exponents will depend on the general colour Casimir group invariants such as the standard ones of C F and C A as well as higher rank ones. In this way previous results should be accessible in various limits. En route though we will identify several errors in earlier results. Unlike the present situation, when the earlier exponents were determined no high order -expansion perturbative results were available to check against. That is not the case now since four loop chiral XY model results are available, [15]. Large N exponents for that model are subsequently deduced from the general non-abelian Nambu-Jona-Lasinio model by simply taking the abelian limit. It will be evident from the expression for η at O(1/N 3 ) that there are no singularities meaning that a smooth abelian limit can be taken. Indeed we will be able to show the origin of the problem in the direct large N conformal bootstrap application to the chiral XY model. As a corollary we will provide an improved estimate for η in three dimensions to the same number of terms as β φ and ν.
The article is organized as follows. We discuss the Lagrangians of the various theories we will consider for the large N analysis in Section 2 where the large N critical point formalism is briefly reviewed. This is applied in the following section to determine the fermion critical exponent at O(1/N 2 ) in the Nambu-Jona-Lasinio model. Subsequently we complete the evaluation of both scalar field anomalous dimensions at O(1/N 2 ) in Section 4 by computing the vertex anomalous dimensions. These and the other exponents are necessary for the O(1/N 2 ) evaluation of the correlation length exponent 1/ν which is carried out in Section 5. The focus then alters to the large N conformal bootstrap formalism with the derivation of the formal underlying equations provided in Section 6. The details are necessary due to the presence of more than one vertex unlike previous large N applications of the technique. The actual evaluation of η at O(1/N 3 ) is given in Section 7 prior to extracting the corresponding expression for the chiral XY model in the abelian limit in Section 8. Estimates of the exponents in three dimensions are also provided there. Finally, we provide concluding remarks in Section 9.

Background.
To begin with we introduce the theories we will focus on. First the four dimensional chiral XY model has the renormalizable Lagrangian where 1 ≤ i ≤ N introduces our expansion parameter. The renormalization group functions of (2.1) are available in [15] and the expansion of the associated critical exponents at the Wilson-Fisher fixed point can be summed up to obtain estimates of the corresponding theory in three dimensions. Here our convention for is defined by d = 4 − 2 . Underlying the Wilson-Fisher fixed point is a universal theory which is driven by a core interaction that is relevant at criticality in all dimensions. It is this universal theory that one applies the large N formalism of [18,19,20] to as the expansion parameter 1/N is dimensionless in all space-time dimensions.
In the case of (2.1) the universal Lagrangian is The scalar fields σ and π are related toσ andπ of (2.1) respectively by a coupling constant rescaling. This is because at criticality the coupling constant is in effect fixed and plays no role due to the scaling and conformal symmetry. In (2.2) each field has a dimension that comprises a canonical part and an anomalous piece. The former is determined from ensuring the action is dimensionless in d-dimensions and we define the full dimensions of the three fields by where η is the fermion anomalous dimension with χ σ and χ π corresponding to the anomalous dimensions of the respective σ and π a vertex operators. At this stage we assume that the latter two are unequal prior to their determination.
The problem, however, is that at leading order in 1/N both χ σ and χ π are zero, [24]. As the large N conformal bootstrap programme is based on an expansion in the vertex anomalous dimension, for (2.2) the expansion cannot begin. This is because the factor associated with the Polyakov conformal triangle, that is at the core of each vertex in the large N bootstrap formalism, depends on 1/χ σ and 1/χ π which we will demonstrate later. While previous work produced various large N exponents to O(1/N 2 ), [24,27], and showed that the vertex exponents are non-zero at O(1/N 2 ) for (2.2), a strategy needs to be devised to avoid the problem of this singularity in the bootstrap construction at leading order. The direction we have chosen to go in is to consider a more general theory which contains the universality class (2.2) in a specific limit and this is the non-abelian Nambu-Jona-Lasinio model which has the universal Lagrangian, [34], Here the fermion field has a non-abelian symmetry which we take to be any Lie group, classical or exceptional, with generators T a and structure constants f abc while the index ranges are 1 ≤ a ≤ N A and 1 ≤ I ≤ N c where N c and N A are the respective dimensions of the fundamental and adjoint representations and π a is a vector in group space. The four dimensional equivalent to (2.1) is which has a similar coupling structure to (2.1).
Comparing (2.2) and (2.4) the main difference is in the decoration of fields with indices leading to Feynman rules with Lie group dependence. Structurally in terms of Feynman graphs that will arise the only difference is the appearance of group factors. These will depend on the usual group Casimirs which are defined by where I Nc is the N c × N c unit matrix. From the equations defining C F and C A we have the relation However at high loop order it is known that new higher rank Lie group Casimirs can appear. For instance, these first occur in a four dimensional gauge theory in the four loop β-function of QCD, [36]. For (2.4) it transpires that the same feature is present in higher order large N exponents and in particular the square of the fully symmetric tensor d abcd F arises where As will be evident from our final expressions the large N exponents of (2.4) will depend on the group entities T F , C F , C A , N c , N A and d abcd F d abcd F . In light of this and observing that it is possible to extract (2.2) from (2.4) by formally setting T a → 1 then the critical exponents of (2.2) should emerge from those of (2.4) in the abelian limit This then completes our computational strategy. We will focus in the main throughout on determining the large N critical exponents of (2.4) before finding those of (2.2) as a corollary via (2.9).
Having justified why our main focus will be on (2.4) we recall the basic formalism for the large N critical point computations, [18,19]. It is centred on the asymptotic behaviour of the propagators in the limit to the Wilson-Fisher fixed point in d-dimensions where d does not play the role of a regulator. Representing the asymptotic scaling form of each propagator by the name of its field, the coordinate space forms are, [37], where A ψ , B σ and B π are x-independent but d-dependent amplitudes. These will appear in two combinations which we define by and together with η and the other exponents can be expanded in powers of 1/N through for example. In addition to the leading terms of (2.10) we have included corrections to scaling corresponding to the terms with (x 2 ) λ and associated amplitudes A ψ , B σ and B π . The exponent λ is usually used to determine the exponent 1/ν which is related to the correlation length and we will consider it here too. In that case the canonical value of λ is (µ − 1). In order to find the first few orders of each exponent requires the solution of skeleton Schwinger-Dyson 2-point functions which require the asymptotic scaling counterparts to the propagators of (2.10). They were derived in [37] by inverting their momentum space forms using the Fourier transform which is given by, [18], (2.13) in general to set notation where (2.14) Consequently we have where the various functions are defined by for arbitrary α, β and λ.
The asymptotic scaling forms of the propagators are necessary in order to algebraically represent the behaviour of the Schwinger-Dyson equations in the critical region. To determine η at O(1/N 2 ) we consider the 2-point functions of the three fields of (2.2) and to the order we are interested the relevant graphs that contribute are provided in Figures 1, 2 and 3. In each there are no self-energy corrections on any of the propagators. This is because the propagator powers include the non-zero anomalous dimensions η, χ σ and χ π that account for such effects. In terms of counting with respect to the ordering parameter 1/N , each closed fermion loop contributes a power of N whereas a σ or π a field has a factor of 1/N associated with it. This is accounted for through the N dependence of the amplitude combination variables z and y defined in (2.12). So, for example, all the two loop graphs differ from the respective one loop graphs of their Schwinger-Dyson equations by a power of 1/N . For the scalar O(N ) theory of [18,19] various three loop graphs also contribute at O(1/N 2 ) but the corresponding graphs are absent here. This is because the extra diagrams have closed fermion loops with three scalars external to that loop. Such graphs are zero because of the trace over an odd number of γ-matrices. Substituting the propagators into the skeleton Schwinger-Dyson equation for ψ equates algebraically to where we have included the various group theory factors associated with the non-abelian symmetry, as well as the vertex renormalization constants Z Vσ and Z Vπ , and effected the γ-algebra. These will have poles in the analytic regularization ∆ that is introduced through the replacement, [18,19], Once the group and amplitude dependence are factored off the remaining contribution is the Feynman integral itself. Its value is the same irrespective of which of the two interactions are involved. Expressions for both Σ 1 and Π 1 are provided in [37] where the latter arises in Figures  2 and 3. The representation of the respective skeleton Schwinger-Dyson equations analogous to (3.1) are In all three cases we have not included vertex renormalization constants in the higher order terms since those counterterms would contribute to the next order in the respective expansions. Also we have cancelled off a common factor of x 2 whose power is related to the canonical dimension.
What remains are factors of x 2 whose exponent involves a linear combination of the regulator and the vertex anomalous dimensions. These cannot be neglected despite being small. This is because in the approach to criticality there is no overall scale. For instance as to the order of interest. If the ln(x 2 ) terms remain then they would give singularities in the limit as x 2 → 0 or ∞. However the simple pole in each of Σ 1 and Π 1 is removed by the vertex counterterm. Once this is fixed the remaining terms involve χ σ 1 or χ π 1 as well as a term dependent on the counterterm but where each is multiplied by ln(x 2 ). These terms are removed by defining the respective vertex anomalous dimensions. A consistency check on this is that the same solution for each of χ σ 1 and χ π 1 emerges from all three equations. Eliminating z and y at successive orders in 1/N from (3.1), (3.3) and (3.4) and expanding the scaling functions for the 2-point functions of (2.16), we find at leading order. From the O(1/N 2 ) parts of the equation, after renormalization we deduce to ensure no ln(x 2 ) terms remain resulting in from the finite part of the consistency equations. Here we use the shorthand notation where ψ(z) is the Euler ψ function.

Vertex anomalous dimensions at O(1/N 2 ).
One of the key ingredients to proceeding to higher order in the large N expansion is the provision of the vertex anomalous dimensions. The general approach is similar to conventional perturbation theory in that at successive loop orders the wave function renormalization constants are needed first prior to renormalizing the mass and coupling constants. While the fermion anomalous dimension is associated with the critical exponent η those for the scalar fields require χ σ and χ π and we summarize their determination in this section. While χ σ 1 and χ π 1 were deduced as a corollary to finding η 2 by excluding ln(x 2 ) terms thereby ensuring the scaling limit could be taken smoothly, the vertex dimensions can also be computed directly from the two vertex 3-point functions. The leading order 1/N corrections for the σψψ vertex are shown in Figure 4. Those for the π aψ γ 5 T a ψ are virtually the same but with the modification that the external σ line is replaced by one representing π a . It is for this reason that we have not illustrated them.
To extract expressions for χ σ 1 and χ π 1 from Figure 4 we follow the method of [18,19,38,39] which is the large N critical point renormalization formalism. One uses the asymptotic scaling forms of the propagators (2.10) but ∆ regularized. This produces an expansion in powers of 1/N where at each order the terms are a Laurent series in ∆. The poles corresponding to the divergences are absorbed into the vertex renormalization constants Z Vσ and Z Vπ in general although the O(1/N ) terms of each are already available. Once the divergences are removed the remaining finite part does not have a non-singular limit to criticality. By this we mean that in coordinate space ln(x 2 ) parts are present. Equally carrying out the computation in momentum space the partner ln(p 2 ) dependence remains. In each case there are terms which involve χ σ 1 ln(x 2 ) and χ π 1 ln(x 2 ). Thus to ensure a scale free finite Green's function these 0 = π −1 + + +  unknowns are chosen to remove the residual x 2 or p 2 dependence. In the way we have presented the larger formalism in the previous section, the values that we extract for χ σ 1 and χ π 1 from the graphs of Figure 4 fully agree with (3.7). At O(1/N 2 ) there are a considerably larger number of diagrams to evaluate. For instance, each of the three vertices in both graphs of Figure 4 will gain vertex corrections. These are displayed in Figure 5 where again we will only illustrate the graphs for the σψψ 3-point function. Those for the π aψ γ 5 T a ψ are obtained by swapping the external scalar in each graph. In addition to these there will be contributions from expanding the graphs of Figure 4 out to O(1/N 2 ) due to the N dependence in the anomalous part of the propagator exponents in each graph. More significantly the set of graphs in Figure 6 need to be included. These are shown separately, partly for a reason that will become evident later, but mainly because they can be regarded as primitive in the sense that the highest order of pole in ∆ is simple. One aspect of the evaluation of the graphs of Figure 6 requiring care is that of determining the group theory factor associated with the three loop light-by-light graphs. This is particularly more involved for the π aψ γ 5 T a ψ vertex function when all the scalar fields are π a ones. In this instance one has to rationalize traces of the form Tr T a T b T c T d for a general Lie group. To handle this we have used the color.h module that is written in the symbolic manipulation language Form, [40,41]. Indeed we have used Form to handle the algebra associated with solving all the various Schwinger-Dyson self-consistency equations. The color.h routine is an encoding of the Lie group Casimir analysis of [42] that goes beyond the rank 2 ones of (2.6). In particular the earlier trace can be rewritten in terms of d abcd F d abcd F and d abc d abc where is fully symmetric. Though for all the light-by-light graphs of Figure 6 the combination d abc d abc cancels in the final expressions for not only χ σ 2 and χ π 2 but also all the O(1/N 2 ) and higher order exponents where d abcd F d abcd Assembling the values of the various diagrams with their associated group theory factor and extracting the piece associated with the separate vertex exponents we find from the graphs of Figure 5 and 6 where Similarly we find Both (4.2) and (4.4) are more involved than the expression for η 2 primarily due to the lightby-light contributions with the term involving Θ(µ) arising from the non-planar ones. We note though that similar to [24] the vertex anomalous dimensions are equal at O(1/N 2 ) for the SU (2) group. Having (4.2) and (4.4) now means the scalar wave function anomalous dimensions are known at O(1/N 2 ). This was achieved by using the leading terms of the asymptotic scaling forms of the propagators (2.10). To deduce other exponents one considers the correction to scaling terms and to find ν in particular at O(1/N 2 ) we take λ = 1/(2ν) where the canonical value of λ is (µ − 1). For (2.4) we therefore include the correction term on the propagators in the graphs of Figures 1, 2 and 3 retaining only those contributions with one insertion only, [19]. Unlike the scalar O(N ) model of [18,19] for the Gross-Neveu universality class there is a reordering of the diagrams with corrections on the σ and π a lines in both their self-consistency equations, [25,26]. This stems from the correction term of the 2-point functions of (2.16) and in particular q(β σ ) and q(β π ). Specifically the leading order large N term of a(β − µ + λ) of (2.16) is O(N ) for both β σ and β π rather than O(1) for the analogous function in the scalar theory of [19]. Consequently to find ν at O(1/N ) the two graphs of Figures 2 and 3 with corrections on the σ and π a internal propagators have to be included. This is evident when the self-consistency equation for λ is constructed. Clearly the correction to scaling terms in (2.10) are of different length dimension to the leading terms. So the leading and correction terms with regard to x 2 in the algebraic representation of the Schwinger-Dyson equations decouple. The leading term that allowed us to deduce η 2 has been dealt with but three equations are left each involving the correction amplitudes A ψ , B σ and B π . The three equations are combined into one equation involving a 3 × 3 matrix containing λ 1 as the only unknown at leading order. However examining the N dependence of the terms in the matrix it is evident that it is different in different rows. Omitting the higher order graphs for the σ and π a Schwinger-Dyson equations would have left an inconsistent set of equations. The final step to find λ 1 is to set the determinant of this matrix to zero whence we find The consequence of the reordering due to the singular nature of q(β σ ) and q(β π ) is that at next order the higher order graphs of Figure 7 have to be included where the dot on a scalar line indicates the propagator with the correction to scaling term of (2.10). Here we have introduced a shorthand notation to compress the number of graphs that actually contribute. In Figure  7 the zigzag line represents both σ and π a fields. So, for instance, this means that in reality there are four graphs from each possible choice zigzag line in the three loop graphs and eight for the four loop ones. In total that is 48 graphs for each scalar field Schwinger-Dyson equation. The values of each graph independent of the group theory factor are available in [43]. However appending the group values the solution to the determinant of the self-consistency equation to next order in 1/N gives The contributions from the light-by-light graphs of Figure 7 are evident.
Finally we can note that the expressions for the exponents to this stage are in agreement with previously determined exponents in the Gross-Neveu model, [25,26,37,43], in the limit since this excludes graphs involving the π aψ γ 5 T a ψ vertex. For the Gross-Neveu XY model taking the limit of (2.9) we also find agreement with [17,24]. These provide important checks on our exponents for (2.4). 6 Large N conformal bootstrap equations.
Before concentrating on the O(1/N 3 ) evaluation of η using the conformal bootstrap we devote this section to recalling one of the key features of the formalism. This is the Polyakov conformal triangle, [29,44], which represents the full vertex function in the theory at criticality. Ordinarily a vertex operator has a non-zero anomalous dimension. In other words in coordinate space the sum of the exponents of the propagators joining to a 3-point vertex has a non-zero part that means the overall dimension does not equal the canonical dimension. This is in fact the reason why the conformal integration rule referred to as uniqueness, [19,31], cannot directly be applied to the two loop graphs of Figures 1, 2 and 3. The canonical dimension of each vertex in the theory matches the value of (2µ + 1) satisfied by the uniqueness rule of [37] but χ σ and χ π have small but non-zero values. The Polyakov conformal triangle construction allows one to exploit uniqueness as a computational tool by replacing the full vertex with a one loop graph where the exponents of the internal lines are arbitrary but are chosen to make the new internal vertices unique. For a Yukawa vertex the most general conformal triangle is shown in Figure 8, [20,25,26]. Specifically a 1 + a 2 + α 3 = 2µ + 1 a 2 + a 3 + α 1 = 2µ + 1 a 3 + a 1 + α 2 = 2µ + 1 (6.1) where α i are the exponents of the original bare vertex of the underlying Lagrangian and a i are the arbitrary internal exponents. By way of orientation, for example, one could have α 1 = α 2 = α and α 3 = β σ for the σψψ vertex of (2.4). The overall function f (α i , a i ) represents the normalization. While it is straightforward to solve the simultaneous equations (6.1) to find the a i in terms of α i in general we will leave this to its explicit use of the construction for (2.4).
One of the reasons for determining the various exponents for (2.4) at O(1/N 2 ) using the skeleton Schwinger-Dyson equations was in part to establish η, χ σ and χ π at this order for the computation of η at O(1/N 3 ). Since this will use a different formalism knowledge of η 2 will serve as an independent check while χ σ 2 and χ π 2 are necessary to extract η 3 . The main difference in the two formalisms rests in how the class of Feynman diagrams are analysed in the critical region. In the conformal bootstrap the key graphs are the primitive ones in the sense that not only are there no self-energy corrections but there are no vertex corrections. In the former case this was because the propagators are dressed but now in the bootstrap approach the vertices are also dressed. So the focus is not on 2-point functions per se but instead on the vertex or 3-point functions which we denote by V σ and V π for the respective vertices of (2.4). The method to construct the consistency equations for each here has been discussed previously, [25,26,37], based on the early work of [29,30,44,45,46]. As that work involved Lagrangians with a single coupling constant we have repeated those derivations for (2.4) which has two independent vertices and will now summarize.  Figure 11, for example, are used to substitute for the bare vertices in the ψ Schwinger-Dyson equation of Figure 9. In both equations for V σ and V π one has to take account of the two interactions in (2.2) rather than one as in previous Gross-Neveu bootstrap computations. The appearance of the regularizing parameter δ is due to the fact that in the formalism of [45,46] there is at least one divergent graph in the final result for the conformal bootstrap consistency equation for each field. The regularization is introduced in the exponent of the external line where the designation of 2δ and 2δ appears. Unlike [20] we will denote the regularizing parameters by δ and δ instead of and in order to avoid confusion with the parameter that is usually associated with the regulator of dimensional regularization. By contrast we recall that here the regularization is analytic.
The situation after following the procedure of [20,45,46] is illustrated in Figure 13. In that equation the Lorentz index indicates that the original full Schwinger-Dyson equation of Figure  9 has been multiplied by the vector x µ , [45,46], where x is the location of one of the external vertices with the other at the origin. In following the manipulations to arrive at this equation in Figure 13 the contribution from the graphs involving the 4-point boxes have been eliminated. In terms of notation the directed solid line represents a coordinate space vector joining the endpoints of the respective full 3-point vertices. At this point we divide the graphs into those with the directed line joining the external points of one full vertex subgraph and those where there is a directed line between two different full vertices. The former set are divergent while the latter are finite. To see this explicitly one replaces the full vertex by the equivalent Polyakov triangle defined by (6.1). As each of the vertices of the resulting three loop graphs are unique it is a simple exercise to apply the Yukawa uniqueness relation to find explicit expressions for all the graphs. For the graphs where the evaluation produces a finite expression then there is a cancellation due to the relative minus sign. This is because the regularized vertex functions V σ (δ ) and V π (δ ) are unity in the limit δ → 0, [20,45,46]. By contrast for the divergent graphs one has to retain the next term of the Taylor series in δ before taking the δ → 0 limit.
To assist with the evaluation we note that the distribution of the regularizing parameters within a conformal triangle is shown in Figure 14 where these are added or subtracted to the extant internal exponents. In defining the regularization with a shift of 2δ rather than δ, [20], we are merely avoiding the appearance of 1 2 δ internally in the conformal triangle. There is no ambiguity as ultimately the regularization will be lifted. We note that in Figure 11 and other figures we only include the relevant argument of the vertex functions for brevity. The full dependence on the various exponents and parameters in fact is V σ = V σ (z,ȳ; α ψ , β σ , β π ; δ, δ ) , V π = V π (z,ȳ; α ψ , β σ , β π ; δ, δ ) . (6. 2) The parametersz andȳ are similar in origin to those of the skeleton Schwinger-Dyson equations of the earlier approach in that they depend on the amplitude combinations A 2 ψ B σ and A 2 ψ B π . They differ in value, however, due to the normalization factor in the definition of the conformal triangle shown in Figure 8. The result of the exercise is the conformal bootstrap equation for ψ which is where we have included the group factors associated with (2.4) and These are related to the residue with respect to δ of the earlier divergent graphs of Figure 9. In the expression for t σ the notation is for t π which are the solutions to (6.1) for the internal exponents. In arriving at (6.3) we have also made use of another result of [20,45,46] in the derivation of the bootstrap equations for the vertex functions which is 1 = V σ (z,ȳ; α ψ , β σ , β π ; 0, 0) , 1 = V π (z,ȳ; α ψ , β σ , β π ; 0, 0) . (6.7) These reflect the sum of all the graphs with dressed propagators and vertices contributing to the two 3-point vertices. The 1/N leading order graphs for both V σ and V π are shown in Figures  15 and 16 which contain the graphs of Figures 4 and 5. The regularization is not necessary for these bootstrap equations since none of the contributing graphs are divergent.
The final part of the exercise in constructing the bootstrap equations is to repeat the procedure that gave (6.3) for the remaining two equations of Figure 10. The process is similar except that now in these graphs there is a δ regularization on the external σ and π a fields. So the initial part of the derivation of the equations uses the regularized results of Figure 12 where the 4-point boxes will contain graphs with both σ and π a fields. This is the reason for the simpler forms compared with Figure 11. The result of implementing the regularized vertex of Figure 12 is shown in Figure 17 for σ where the Lorentz index indicates that the scaling function has been multiplied by x µ too. The corresponding graphs for π a are similar to Figure 17 with the σ external legs replaced with π a fields in addition to the replacement of σ −1 µ with π −1 µ . While the four graphs of Figure 17 differ from the corresponding ones of Figure 13 it is still the case that those with the directed line connecting two conformal triangles are finite. Evaluating the remaining where the same factors t σ and t π occur. Eliminating these and the amplitudesz andȳ finally produces the consistency equation for η which is where the restriction indicates that both regularizations are lifted at the end of the evaluation of the underlying graphs. At this stage (6.9) is exact and its solution would determine η to all orders. However to achieve this would require the explicit values from all the contributing graphs. The consistency equation (6.9) can be refined though by recalling the observation of in [20]. For the moment we consider a general situation and denote a regularized n-point function by Γ (n) (∆ i , δ, δ ) where each vertex is represented by a conformal triangle, the jth vertex has anomalous dimension χ j ≡ 2∆ j and two of the external legs have δ and δ regularizations. Then the Green's function will have a simple structure. In the unregularized case it will formally be Γ(∆ i , 0, 0)/( n j=1 ∆ j ) but in the regularized case where legs 1 and 2 are chosen to have the respective regularizations δ and δ and Γ (n) (∆ i , δ, δ ) is analytic in all the ∆ i and both regulators. Then for instance. The first term depends only on ∆ 1 and ∆ 2 respectively since the equivalent of (6.7) has been employed and the relation can be written more compactly as, [20], where res indicates the contribution from the derivative of the residue of Γ (n) (∆ i , δ, δ ) only. This was for a general scenario but returning to (2.4) and applying this procedure to (6.9) we have where we have now set χ σ = 2∆ σ , χ π = 2∆ π (6.14) and the arguments of the two vertex functions have been suppressed for brevity but are the same as the corresponding terms in (6.13). For practical purposes it is best to expand the two vertex functions by setting although we will only need V σ 1 , V σ 2 , V π 1 and V π 2 to determine η 3 . If one focuses on the leading order large N piece of (6.13) it is equivalent to that which produced η 1 earlier. So, for instance, a contribution to the O(1/N 3 ) part from the first term on the right hand side of (6.13) can be simplified to for the V σ vertex while the contribution from the second term involving V π is similar. 7 Evaluation of η 3 .
Having derived the key equation (6.13) satisfied by η 3 we now turn to details of its determination. This entails computing the corrections to the vertex functions of the conformal bootstrap equations (6.7) and (6.13). The leading order contributions are shown in Figures 15 and 16 and correspond to V σ 1 and V π 1 respectively where At next order the graphs that comprise V σ 2 and V π 2 are illustrated in Figures 18 and 19 and are clearly primitive. While the leading order 1/N terms of both V σ 1 and V π 1 are needed to verify η 2 the next term in the series is required for η 3 . This N dependence arises through the exponents α ψ , β σ and β π . Therefore what are sometimes termed 3-gamma graphs, which are illustrated in Figures 15 and 16, have to be computed. This is achieved via the method developed in [20] using the properties of the Polyakov conformal triangle. For instance, the explicit designation of the exponents on the lines of both V 1 σσσ and V 1 σππ are shown in Figures 20 and 21 where is set for shorthand due to β σ and β π differing only in the piece involving the respective vertex anomalous dimensions. The factor of 2 associated with ∆ σ and ∆ π is a consequence of (6.5) and (6.6). The graph of Figure 20 was evaluated in [27]. If we set which is consistent with the general form of Γ (n) then we recall that with the shorthand notation, [20], and We note that in [27] not all the terms quadratic in combinations of the parameters δ or δ were recorded. We have included them here for completeness and as an aid to an interested reader although they do not contribute to the consistency equation. We have rederived (7.4) that appeared in [27] for V 1 σσσ in order to extend it for V 1 σππ . The latter is more general than the former since it will depend on χ π as well as χ σ . This leads to a check in that formally taking the limit χ π → χ σ in V 1 σππ should produce the expression for V 1 σσσ . If we define and repeat the procedure that resulted in (7.4) we find Clearly this is consistent with (7.4) since F σππ (δ, δ , ∆ σ , ∆ σ ) = F σσσ (δ, δ , ∆ σ ). These two expressions are sufficient to calculate the contribution to (6.13) from V σ 1 . For the second term of (6.13) involving V π 1 the expressions for the corresponding 3-gamma graphs can be deduced from those of V σ 1 by formally swapping ∆ σ and ∆ π in (7.4) and (7.8). To complete the evaluation of η 3 requires the contribution to (6.13) from V σ 2 and V π 2 . The values of the four topologies underlying the graphs in Figures 18 and 19 have been computed previously, [25,26,27]. Strictly though it is the combination of (6.16) and the π a counterpart that have been calculated. The values of the individual terms cannot be extracted as the underlying three loop integral does not have unique vertices to allow the integration to proceed. While it would appear from the figures that there are three distinct topologies this would be the case if there was no regularization. With non-zero δ and δ through specified external legs there are two independent contributions from the three loop non-planar graphs. Indeed the contribution to (6.7) from these non-planar graphs is the same but not for (6.13). For each of the graphs of Figures 20 and 21 we have computed the associated group factor by again using the color.h package based on [42]. The higher rank colour Casimir will arise in the graphs of Figure 19 where four π a fields connect to the fermion loop.
The next stage is to assemble all the components contributing to (6.7) and (6.13). The purpose of the two equations of (6.7) is to fix the values ofz andȳ at O(1/N 2 ). While these were eliminated to arrive at (6.13) combinations of them are present as factors in each of the graphs of Figures 15, 16, 18 and 19. However we can now illuminate the difficulties with a direct evaluation in the U (1) case. The formal solution forz in (2.4) at leading order is While C A vanishes in the abelian limit producing a singularity, when the leading order values for the vertex dimensions are included these are proportional to C A . So that overall the leading order variablez is finite in the abelian limit. For a direct evaluation one could not determinez 1 since it would equate to the mathematically ill-defined quantity of 0/0. The situation forȳ 1 is similar. It was this that forced us to consider the more general Lagrangian. For (2.4) we find for instance which is finite in the abelian limit. With these values and those for the next order we ultimately arrive at our main result for (2.2) which is which is considerably more involved than the other exponents. We note that electronic versions of (7.11) and all the large N exponents of (2.4) are available in an attached data file.
Like the η 3 result of [20] (7.11) contains the function Ξ(µ) which is related to the function I(µ) of that paper via The expansion of I(µ) around d = D − 2 is known for integer D. In the case of D = 2 and 4 this was derived in [47] and built on the earlier expansion provided in [20]. When D = 3 the expansion was provided in [48] where the leading order term is which was derived in [20] using uniqueness methods. So, for instance, we can record the first three terms of the large N expansion of each exponent in three dimensions for (2.4) which are from which we will be able to deduce the analogous expressions for the chiral XY model. Finally we note that (7.11) agrees with known expressions in the limit of (5.3).  [15] and the functional renormalization group [11].

XY model exponents.
Having established the critical exponents for the non-abelian Nambu-Jona-Lasinio model, (2.4), we can now derive those for the chiral XY model in the abelian limit. It is clear that there are no singular colour group factors for the results of (3.6), (3.7), (3.8), (4.2), (4.4), (5.1), (5.2) and (7.11). Therefore it is safe to take the abelian limit of (2.9) giving at leading order which clearly illustrates the issue surrounding the vanishing leading order vertex anomalous dimension for both Yukawa interactions. At next order we have We have checked that these expressions are in agreement with earlier work, [17,24]. Finally we obtain the main result for this model which is The expression is considerably simpler than its non-abelian counterpart similar to the O(1/N 2 ) exponents. One aspect of this is that the function Ξ(µ) is absent which is not unrelated to the vanishing of χ σ and χ π at leading order.
While the exponents to O(1/N 2 ) agree with previous derivations, it is possible to check them against perturbative computations in four dimensions. In [15] the β-function and anomalous dimensions of the chiral XY model were determined to four loops in the MS scheme. Moreover the expansion of the corresponding critical exponents were computed to O( 4 ). Therefore if we expand our expressions in the same expansion they ought to be in agreement. Setting d = 4 − 2 in (8.1), (8.2) and (8.3) we find Figure 22: Padé approximants of η XY for N t = 2, 4, 6 and 8.
where the combination tallies with the definition used in [15,16]. We have also set N = 4N t to adjust for the different trace and group conventions used here in the underlying master integrals for the skeleton Dyson-Schwinger and conformal bootstrap equations. So N t corresponds to the N used in [15,16]. It is straightforward to verify that each expansion is in full agreement with [15] after allowing for a different convention in the definition of d in terms of . In light of this agreement with high order perturbation theory we can obtain expressions for the exponents in three dimensions. numerically.
The next stage is to extract estimates for the three dimensional exponents for the case of interest in [15] which corresponds to N t = 2 here. To improve the convergence of the expansion we have followed a similar approach to [28] and evaluated several Padé approximants for each of the exponents. These have been plotted in Figures 22,23 and 24 for N t = 2, 4, 6 and 8 for 2 < d < 4 for η XY , η XY φ and 1/ν XY respectively with numerical values for three dimensions recorded in Table 1. The latter three values of N t are selected to gauge the effect the higher order 1/N corrections have. In the case of η it is worth noting that we can plot the behaviour of η 3 unlike [28]. This is because Ξ(µ) is absent from η 3 and it was not possible to write the function in an analytic form that could be passed to the plotting routine in that case. For η XY the [2, 1] Padé had a singularity at around d = 3.6. Equally for N t = 2 the [0, 2] Padé of 1/ν XY had a maximum near a similar value of d. We have omitted plots of both Padé's for each exponent.
What is evident in the plots of η XY for various N t is that there is convergence of the O(1/N 3 ) values down to N t = 4. For N t = 2 there is a difference between the O(1/N 2 ) and O(1/N 3 ) curves but the latter is bounded by the leading two orders. This is similar to higher values of N t . Given this it might not be unreasonable to give 0.083 as a rough estimate of η XY for N t = 2 being the average of the two Padé estimates. This is not inconsistent with the N t = 2 estimates from other approaches which are summarized in Table 2. In the case of η XY φ it is clear from Figure 23 that there is very little difference for the value of the exponent between successive orders in 1/N . Indeed comparing the values for N t = 2 with other estimates there is a degree of consistency. Finally the situation for 1/ν XY does not appear to be as good especially at leading order for low values of N t due to the kink that is evident in Figure 24 in its behaviour across d-dimensions. Although it appears that this feature washes out for larger values of N t , it is not clear whether this would be the case for N t = 2 if an O(1/N 3 ) expression were available. Again one could assume that the two large N Padé approximants for 1/ν XY were bounds of the true result. In this instance one would arrive at an estimate of around 0.845 which is not dissimilar to the perturbative Padé or functional renormalization group estimates of [11,15].

Discussion.
We have provided the large N critical exponents for the Nambu-Jona-Lasinio universality class with Lagrangian (2.4) for a general Lie group using the large N critical point formalism pioneered in [18,19,20]. This included the large N conformal bootstrap formalism that allowed us to extract η 3 . To achieve this we had to extend the formalism to the situation where there are two independent 3-point interactions. By independent we mean the vertex anomalous dimensions of the separate vertices were inequivalent. The derivation of the underlying bootstrap equations was straightforward but was discussed at length. Indeed from that it is evident that they can be readily generalized to the case of more independent 3-point interactions. Taking the more general approach with a general Lie group means that results for exponents in various models can be extracted in certain limits. Previously η 3 was calculated for a limited number of fermionic models including the SU (2) Nambu-Jona-Lasinio model, [35], as the anomalous dimensions of each vertex were the same. To produce the general Lie group results we were able to benefit from the Form encoded color.h package which was based on [42] that allowed us to express our exponents in terms of the general group Casimir invariants. This was particularly important for the light-by-light graphs that contribute to several exponents at O(1/N 2 ) and η 3 . One of the main consequences was that the abelian limit could be taken smoothly to deduce exponents for the chiral XY model. This circumvented the application of the bootstrap formalism directly to (2.2) due to an ill-defined intermediate quantity that formally obstructs the direct derivation in the abelian case. The ultimate general expressions are analytic in the colour group invariants allowing the chiral XY model exponents to emerge smoothly. Finally this more general group approach for the large N formalism offers a more compact way of extracting critical exponents for other universality classes. One benefit for instance is that it gives access to additional information on the terms in the explicit perturbative series of the underlying renormalization group functions.