Symmetric point flavour singlet axial vector current renormalization at two loops

We calculate the two loop correction to the quark 2-point function with the non-zero momentum insertion of the flavour singlet axial vector current at the fully symmetric subtraction point for massless quarks in the modified minimal subtraction (MSbar) scheme. The Larin method is used to handle $\gamma^5$ within dimensional regularization at this loop order ensuring that the effect of the chiral anomaly is properly included within the construction.

One of the more curious experimental results over a generation ago was that of the EMC collaboration, [1]. They measured the origin of the proton spin and discovered that against expectations it was not due in a major part to the valence quarks. As the proton is a bound state of three quarks it was widely assumed that the combination of their quark spins would be the source of the overall spin- 1 2 of the proton. Instead the experiment observed that the gluons binding the quarks together give a sizeable contribution. This was surprising due to the fact that in some sense the gluons are sea partons. While the original experiment was subsequently refined and improved to confirm the original observation, [2,3,4,5], a clear theoretical understanding was sought to explain the phenomenon. As such a venture requires the use of the strong sector of the Standard Model described by Quantum Chromodynamics (QCD), tools had to be developed and refined to tackle the problem. Moreover, to do so one has to study an energy regime which is in the infrared and hence outside the region where perturbation theory is valid. Therefore the only viable approach was the application of lattice gauge theory which can access the non-perturbative structure of the proton through heavy use of supercomputers. Clearly such an exercise required new methods such as the inclusion of dynamical fermions and the field is better placed now to answer the theoretical question of the source of the proton spin. This is not an isolated exercise for lattice studies. As an aside it is worth mentioning that recently the breakdown of the proton mass in terms of its constituent entities such as quark, gluon, weak sector and anomaly contributions has been accurately estimated on the lattice, [6]. This entailed measuring the diagonal components of the energy-momentum tensor. Indeed the study given in [6] has indicated that more accurate knowledge of the internal proton structure can be adduced theoretically in the near future. Parenthetically it is also worth noting the related problem of the pressure inside the proton. Experimentally this can be deduced accurately now as demonstrated in [7]. In terms of theoretical studies and in particular those for the lattice such a pressure problem translates into requiring precise measurements of the off-diagonal components of the energy-momentum tensor. Such an exercise has been carried out recently in [8,9]. However progress on studying the source of the proton spin on the lattice over the last few years can be seen in a non-exhaustive set of articles, [10,11,12,13,14], while a more detailed overview of this and the status of future directions of hadron physics computed on the lattice can be found in [15,16,17,18].
Concerning proton spin measurements on the lattice the quantum field theory formalism behind such potential calculations originate in the work of Ji, [19]. There the relevant, in the sense of important, operators were identified and it was shown how their expectation values relate to the overall proton spin. Indeed central to the three properties of the proton mentioned already is the need to study matrix elements of various key operators for each observable. While the energy-momentum tensor provides the main operator in relation to mass and pressure, in the spin case it is the flavour singlet axial vector current orψγ 5 γ µ ψ where ψ andψ are the respective quark and anti-quark of the same flavour. Indeed one of the motivations for studying the singlet axial vector operator is that the difference between the singlet and non-singlet cases provides a way of quantifying the strange quark contribution to the proton spin, [20,21,22]. Over the years similar quark bilinear operators have been widely studied on the lattice, [20,23,24,25,26,27,28]. Aside from the associated matrix element which determines the non-perturbative structure, the operator renormalization has to be understood in the lattice regularized theory in order to implement the renormalization group running over momentum scales. In indicating the similar operators we mean the flavour non-singlet quark bilinear operators which are termed scalar, vector, tensor, axial vector and pseudoscalar depending on their Lorentz properties. Moreover, the matrix elements have been computed for a variety of configurations. These break into two classes known as forward and non-forward where the former has the operator inserted at zero momentum in a quark 2-point function. In the latter case it has a momentum flowing through it and the square of that momentum and those of the two external quarks can take different non-zero values. This non-forward set-up allows more freedom to probe detailed structure within nucleons. Most lattice studies of the flavour axial singlet current have been for the forward case, [20,23,24,25,26,27,28], and in the associated lattice renormalization scheme termed the modified regularization invariant (RI ) scheme introduced in [29,30].
One aspect of making lattice measurements in general is that of ensuring the continuum limit is taken accurately. In recent years to assist this the process has been adopted where the relevant matrix element is evolved to the ultraviolet region and matched to the continuum perturbative expansion of the same matrix element or Green's function. Clearly the more loop orders that are available in the perturbative expansion means the matching will be more accurate and hence the lattice error estimates can be improved. For quark bilinear operators the early work in this direction was provided in [29,30]. In the context of lattice spin measurements there have been studies, [23,24,25,26,27,28,31] where the nucleon isovector scalar, axial vector and tensor charges were measured. For example [31] built on an earlier parallel study in [27] where the nucleon axial form factors were computed and the issues centered on the chiral anomaly taken into account through a nonperturbative treatment. However in the matching to the continuum in the work of [31] only one loop information was available for the flavour singlet axial vector current. This is because at that loop order the matrix element for the flavour singlet and non-singlet axial vector currents are the same. The difference in these only occur at two loop order. The main reason for this is that the flavour singlet axial vector current is not conserved due to the chiral anomaly and its effect in the matrix element becomes present at two loops. In Feynman graph language there are graphs that are zero in the flavour non-singlet case but non-zero for the flavour singlet operator. By contrast in lattice language there are disconnected contributions to the proton correlation function when probed with the axial vector current in the flavour singlet case. Therefore while the latter have been incorporated in the lattice simulations, [27,31], their omission in the matching to the continuum in [31], albeit due to taking only one loop data, means that the effect of the chiral anomaly has not been taken into account. In other words the central values measured on the lattice for the hadronic matrix element will not have accommodated the discrepancy in the continuum difference between the flavour non-singlet and singlet operators. This is due to the lack of a two loop computation of the relevant Green's function containing the flavour singlet axial vector current.
Therefore it is the purpose of this article to close this particular gap. By doing so we will bring all the quark bilinear operator Green's functions to the same loop level for the non-forward case. Specifically we will compute the two loop quark 2-point function with the singlet axial vector current at non-zero momentum insertion for the fully symmetric momentum configuration. Such a configuration is a non-exceptional one and hence should avoid infrared problems. The extension of the early lattice work on operator renormalization at an exceptional point of [29,30] was extended to the non-exceptional case in [32] at one loop and later to two loops in a variety of articles, [33,34,35,36,37]. More recently in the lattice study of [38] the operator renormalization constants for all the flavour non-singlet quark bilinear operators were measured and it was demonstrated that for massless quarks those for the non-exceptional configuration were much more reliable in the infrared limit. Therefore in focussing on the fully symmetric point we are aiming at minimizing other potential sources of avoidable error for matching to lattice data. While straightforward to state, the two loop computation we will undertake in dimensional regularization is also fraught with technical complications. One obvious one concerns the treatment of γ 5 but we will use the Larin approach, [39], which is valid at least to the loop orders we are interested in. Moreover, the chiral anomaly has been correctly treated in that approach. As part of our study we will extend the Larin construction in the sense that the nonforward matrix elements are computed. So we will show how the same finite renormalization constant emerges as that of [39] to ensure chiral symmetry is correctly present after renormalization and the subsequent lifting of the dimensional regularization. This is a non-trivial task and in some sense our study reinforces the elegance of Larin's approach. In saying this we will also check the same finite renormalization constant arises for the pseudoscalar current in the non-forward configuration to substantiate that a consistent picture emerges. En route we will discuss a minor modification of the Larin approach for flavour non-singlet operators. In [39] the criterion to define the finite renormalization constant needed to ensure four dimensional properties correctly emerge for the two spin-1 quark bilinear operators was to implement equality of the two relevant Green's functions. By contrast we show that this can also be achieved from ensuring the currents are conserved in four dimensions after renormalization.
The article is organized as follows. In the following section we introduce the formalism used to carry out the two loop computation as well as discuss the treatment of γ 5 in dimensional regularization using the Larin method, [39]. Included in this section is the algorithm we implement to evaluate the singlet axial vector current Green's function at the symmetric point. Section 3 is devoted to the discussion of our results where we quantify the difference between the flavour non-singlet and singlet axial vector operator Green's functions prior to providing concluding remarks in Section 4.

Background.
We devote this section to describing the details of the computation and en route review previous renormalizations of the operators in question at the symmetric point. In order to appreciate the subtleties between the flavour non-singlet and singlet quark bilinear operators we consider we define the non-singlet operators as S ns ≡ψ i ψ j , V ns ≡ψ i γ µ ψ j , T ns ≡ψ i σ µν ψ j , A ns ≡ψ i γ 5 γ µ ψ j , P ns ≡ψ i γ 5 ψ j (2.1) and for the singlet case where i and j are flavour indices and there is no sum over i in the latter set. Given that our main interest is in the perturbative structure of a specific Green's function for a particular external momentum configuration we define the general case, which includes the above operators as well, by where the operator O corresponds to any one of (2.1) and (2.2). We use a similar notation to [35,36,37] and follow the general approach provided there. The two independent external momenta, p and q, satisfy the condition for the symmetric subtraction point which is, [40,41], where µ is a mass scale. This setup is a nonexceptional configuration and therefore there are no infrared issues. The Green's function (2.3) is illustrated in Figure 1 where an operator O I of (2.1) or (2.2) is indicated by the crossed circle.
Having defined the Green's function the first step in its perturbative evaluation is the generation of the two loop Feynman graphs. To do this we have used the Qgraf package, [42], which produces 1 one loop and 13 two loop graphs. At the former order the expressions for the respective flavour non-singlet and singlet Green's functions are the same irrespective of the subtraction point. The first place where any discrepancy will appear as a result of flavour symmetry in the chiral limit is at two loops and is due to the two graphs of Figure 2. This is because these graphs are zero for the flavour non-singlet case as the insertion of such an operator into the closed fermion loop gives a trace over a traceless flavour group generator. So such graphs will not contribute to Γ I ns . By contrast for the flavour singlet case the graphs of Figure 2 will not be zero by this particular flavour trace argument. However for the operators S s , T s and P s the graphs of Figure 2 are zero since there are an odd number of γ-matrices in the closed loop in the chiral limit. These graphs would correspond to the disconnected graphs in the proton correlation function on the lattice. For V s and A s there will be an even number of γ-matrices in the loop in the massless case considered here. So these are the two possible instances of the flavour symmetry producing different two loop Green's functions (2.3). As we will be using dimensional regularization in our calculations one concern that could be raised is whether this argument applies for operators containing γ 5 . In the case of A ns and P ns we are permitted to use the naive anticommuting γ 5 since it always appears in Feynman integrals in an open string of γ-matrices. It is its presence, or an odd number of them, within a closed fermion loop that requires special treatment in dimensional regularization. We will discuss this in detail later aside from noting that when γ 5 , or an odd number of them, is present in a closed loop with an odd number of ordinary γ-matrices then the spinor trace is still zero as in four dimensions. Having taken flavour and Lorentz symmetries into account there is one final constraint to be considered however which is that corresponding to the colour vector space. All operators of (2.1) and (2.2) are colour singlets. In other words they do not include a colour group generator. So the colour trace is not the same as the parallel graphs contributing to the quark-gluon vertex function. In that instance if the operator insertion was replaced by a gluon then the sum of the fermion one loop subgraphs would be proportional to the structure constants f abc . In the case of V s the colour trace with one fewer group generator produces δ ab which together with the relative minus sign arising from the respective γ-matrix traces means that the two graphs of Figure 2 sum to zero for this operator. For A s by contrast while the colour argument applies equally, the presence of γ 5 in the spinor trace prevents the same procedure giving zero since the traces sum and are not equal and opposite. Therefore the only singlet operator that needs to be considered at any subtraction point, and not just the symmetric one, is A s as it will produce contributions additional to its non-singlet partner. Having isolated the only case we have to consider for singlet operators then in order to evaluate Γ A s we note that we have taken a different path to the partner computation of Γ A ns of [35,36,37]. To appreciate the contrast it is first best to summarize the earlier approach. In [35,36,37] a projection method was used whereby the Green's function was decomposed into a basis of tensors consistent with the symmetries of the operator. For instance, for V this was

⊗ ⊗
and this Lorentz basis for each bilinear operator Green's function is the same whether the operator is flavour non-singlet or singlet. Here we have introduced the generalized γ-matrices of [43,44,45,46,47] which are denoted by Γ µ 1 ...µn and defined by for integers n with 0 ≤ n < ∞. Each matrix is fully antisymmetric in the Lorentz indices for n ≥ 2 and the full set span the infinite dimensional spinor space when dimensional regularization is implemented. We note that Γ µ 1 µ 2 µ 3 µ 4 µ 5 is not related to γ 5 and in four dimensions Γ µ 1 ...µn (n) = 0 for n ≥ 5. One advantage of these matrices is that they provide a natural partition since with I µ 1 ...µmν 1 ...νn denoting the unit matrix in Γ-space. Here we use the convention that when a momentum is contracted with a Lorentz index in Γ µ 1 ...µn then the momentum itself appears in place of the index.
Once the basis for the Green's function has been chosen then the coefficients of each tensor is deduced by multiplying it by a d-dimensional linear combination of the tensors to produce a sum of scalar Feynman integrals. These are then evaluated by applying the Laporta algorithm, [48], which systematically integrates by parts all the contributing graphs. This produces a linear combination of a small set of master integrals with d-dependent rational polynomial coefficients. For the symmetric point configuration explicit expressions for the one and two loop master integrals have been available for many years, [49,50,51,52]. Though in more recent years they have been understood in the language of cyclotomic polynomials, [53]. To be explicit to two loops the Green's functions with the kinematic configuration of (2.4) involve different linear combinations of numbers from the set Here ζ n is the Riemann zeta-function, s 2 (z) and s 3 (z) are defined in terms of the polylogarithm function Li n (z) and ψ(z) is the derivative of the logarithm of the Euler Γ-function. A final important aspect of the computation was the extensive use of symbolic manipulation which was facilitated through the use of the language Form and its threaded version Tform, [54,55]. Using the Reduze package, [56,57], written in C++ which implements the Laporta reduction we inserted the relations generated by the package for the required Feynman integrals via a Form module so that all our computations were carried out automatically. In addition we followed the prescription of [58] to implement the operator renormalization within the same automatic framework.
As we will be considering Green's functions involving operators with γ 5 present we need to be careful in its treatment within dimensional regularization. Moreover the strategy we will choose to follow has to be robust in the sense that it should reproduce the known two loop symmetric point results of earlier work, [32,33,34,35]. For the non-singlet example of A ns and P ns there are several ways one can proceed. First though we need to discuss the details of the Larin approach as general features of that computation will be used. In doing so it is important to note that in [39] the calculations were carried out at an exceptional momentum configuration where the operator is inserted at zero momentum. In this case and also for the symmetric point, which is non-exceptional, for open strings of γ-matrices which include γ 5 matrices their naive anti-commutation with d-dimensional γ-matrices is valid. However, in [59] an alternative point of view was introduced in d-dimensions where a generalized set of quark bilinear operators were renormalized in the dimensionally regularized theory. These were (2.10) In four dimensions they reduce to S ns , V ns and T ns respectively for n = 0, 1 and 2. For n = 3 and 4 we have where µνσρ is the totally antisymmetric strictly four dimensional pseudo-tensor. Like γ 5 it has no existence outside strictly four dimensions. For n ≥ 5 the operators of (2.10) are evanescent in the sense that they are not present in strictly four dimensions due to their being more free Lorentz indices than the spacetime dimension which cannot be possible for a fully antisymmetric object.
In the dimensionally regularized theory all these operators exist and can be renormalized multiplicatively. Therefore given the relation (2.11) between the generalized operators there ought to be a connection with not only the renormalization of all the operators of (2.11) but the respective Green's functions should be in full agreement in strictly four dimensions after renormalization. It was one of the aims of [39] to substantiate this given that one in principle already knows the answer from the naive anticommutation of γ 5 . Therefore in [39] a two stage process to renormalize A ns and P ns from the operators (2.10) was proposed. The first part was to renormalize the operators O µνσ ns (3) and O µνσρ ns (4) separately in the usual way in the modified minimal subtraction (MS) scheme to produce what is termed the naive renormalization constant for the d-dimensional operator. Aside from the one loop scheme independent term of the operator anomalous dimension the renormalization constants for the respective pairs O ns (0) and O µνσρ ns (4) , and O µ ns (1) and O µνσ ns disagreed. We note that while Γ µ (1) and γ µ are equivalent we will retain the former notation when we discuss the renormalization of the set of operators of the form (2.10). This disagreement in the naive renormalization constants is clearly not consistent with expectations from the naive anticommutation of γ 5 . However the second stage of the Larin process was to recognise that it is not possible to retain the symmetry properties of a purely four dimensional entity, γ 5 , in the dimensionally regularized theory. So to circumvent the absence of the four dimensional properties in the dimensionally regularized theory, one has to augment the naive renormalization of O µνσ ns (3) and O µνσρ ns (4) by including a finite renormalization constant pertinent to each operator. The condition to define this for the flavour non-singlet case is to impose the constraint that after the naive renormalization of the Green's function with the nindex operator insertion the result is equivalent to that when n is replaced by (4 − n), [39]. Here 4 is chosen as it is the critical dimension of QCD. There is a practical caveat to this in that the Lorentz tensor basis has to be written in terms of purely four dimensional objects first before the finite renormalization can be made. Also the finite renormalization constant in MS derives from the basis tensor corresponding to the tree term of the Green's function.
For the symmetric point case we consider here we have adapted this approach but first checked that the previous non-singlet results are reproduced for n = 3 and 4. However for the Green's functions for both operators the basis of tensors is larger than at the exceptional point. As a first step for n = 3 and 4 we have taken a slightly different tack to [35,36,37] and avoided using a projection method on the Green's function. Instead to evaluate the contributing Feynman graphs we first removed all the γ-algebra from the Feynman integrals and evaluated the underlying tensor integrals. One reason for this is that it bypasses the complication of constructing a decomposition into a Lorentz basis with 3 and 4 free indices in the case of A ns and P ns respectively. In the latter case for instance that basis would include Γ µνσρpq (6) as one example. While this is evanescent it would have to be included to ensure the tensor basis was complete. So in avoiding a direct projection and consequently reproducing the same results as previous computations, [32,33,34,35], will ensure that we have established a valid algorithm for handling γ 5 in the non-singlet case. This therefore will eventually be our strategy for the vector and axial vector operators in the singlet case when the complications due to the graphs of Figure 2 have to be taken into account. While the correct tensor basis will emerge from this integral projection the four dimensional basis (2.5) will not be the relevant one for the naive renormalization of A ns . Using relations such as (2.12) in four dimensions, for instance, the analogous basis will be in four dimensions. Next we recall that to properly renormalize the operator V ns requires some care. As the non-singlet vector current is conserved in the chiral limit ∂ µ ψ i γ µ ψ j = 0 (2.14) its renormalization constant is unity to all orders in perturbation theory and in all renormalization schemes. In the MS context this means that the Green's function for V ns is completely finite. However in a kinematic renormalization scheme where renormalization constants can have non-zero finite parts the renormalization constant for V ns could mistakenly be chosen by demanding that the coefficient of P V (1)µ (p, q) is unity for instance. This would clearly contradict the general result that physical operators do not get renormalized. In practical terms the conservation of the currents is encoded in the quantum theory by a Ward-Takahashi identity. In the case of V ns this corresponds to Green's function of ∂ µ ψ i γ µ ψ j being related to the quark 2-point function. Therefore we have checked that for our integral projection construction that this is indeed the case in MS for the fully symmetric and hence non-exceptional configuration. We note that we found full agreement with the one and two loop results of [32,33,34,35]. This was also checked in [59] for the direct projection on the operator Green's function.
Since the non-singlet axial vector current is also conserved similar general reasoning also applies. Calculating the Green's function of ∂ µ ψ i γ 5 γ µ ψ j , however, it does not agree with γ 5 times the quark 2-point function. This is consistent with Larin's observation that treating the naive renormalization of the generalized operator O µνσ ns (3) as being equivalent to A ns is incorrect. Instead an additional finite renormalization constant is required to ensure the consistency with symmetry properties of the strictly four dimensional theory. Ensuring the consistency with the quark 2-point function in this case we find that to two loops at the symmetric point the finite renormalization is in full agreement with [39] where a = g 2 /(16π 2 ) and g is the gauge coupling constant. It should be stressed that we have verified that the finite renormalization is independent of the subtraction point which is a non-trivial observation. As can be seen in earlier work [32,34,35,36,37] the two loop explicit expression for (2.3) for each operator involves linear combinations of (2.8).
These in principle could have been present in the two loop finite renormalization constant with our choice of (2.4). In the momentum configuration used in [39] only rational numbers or ζ 3 were present in the finite parts of the Green's function at three loops. So it is reassuring that the polylogarithms, for instance, of (2.8) are absent. We note that although the three loop finite renormalization is also known in the MS scheme, [39], in order to verify the next term of (2.16) at the symmetric point would require a new three loop computation. At present the three loop master integrals are not known in order to be able to perform this calculation. This comparison of the n = 1 and 3 Green's function is in keeping with the spirit of [39,59]. However in the singlet case the presence of the extra graphs of Figure 2 will not make this a viable way to proceed as was indicated in [39]. In the non-singlet case the corresponding Ward-Takahashi identity for the axial vector operator provides a more field theoretic alternative method to determine the finite renormalization. In other words the Green's function of ∂ µ ψ i γ 5 γ µ ψ j has to be equivalent to the quark 2-point function multiplied by γ 5 . To ensure this we have checked that the same gauge parameter independent finite renormalization constant, (2.16), is required. It should also be noted that we have checked that once this finite renormalization is determined from ensuring that the current conservation is preserved then the expressions for both Green's functions of V ns and A ns are also in full agreement. By this we mean that the coefficients of P I (i)µ (p, q) for O ns = V ns and A ns are equal for each i = 1 to 6 and thus establishes the equivalence with Larin's strategy in strictly four dimensions after renormalization.
One of the reasons for reviewing the non-singlet operator renormalization and indicating an alternative way of defining the finite renormalization constant is that it avoids the need to connect Green's functions for different operators. For the singlet axial vector operator there is clearly no analogous partner. Instead as highlighted in [39] the associated finite renormalization is derived from ensuring that the chiral anomaly is correctly restored in four dimensional Green's function after renormalization. The non-singlet axial vector current is non-anomalous. Specifically the anomaly for the singlet case is given by where the right hand side can be written as a derivative of a gluonic current and (2.14) and (2.15) become partners in this alternative view of [39]. In [39] the two loop finite renormalization constant was computed and is given by which differs from (2.16) in the final two loop term. We have not reproduced this expression beyond one loop here due to a subtlety in (2.17). This is to do with the fact that both operators of (2.17) have first to be renormalized and the triangular mixing matrix of the operators determined, [39]. Then to find Z fin A s the anomaly equation itself, (2.17), has to be inserted in a Green's function and evaluated in strictly four dimensions. In [39] this was carried out by inserting in a gluon 2-point function where the momentum flowing into one of the external gluon legs vanishes. With a total derivative present due to having to take the divergence of the current, the momentum flow into the operator cannot be nullified. This momentum configuration allows one to use the Mincer algorithm, [60], which evaluates three loop massless 2-point functions in dimensional regularization. A three loop calculation is in fact required to determine (2.18) due to the presence of the coupling constant a on the right hand side of (2.17). This means that the loop order of the anomaly Green's function and its coupling constant expansion are mismatched by one. Since this will also be the case for the symmetric point setup, the absence of the three loop master integrals for such Green's function evaluations means that that calculation cannot be carried out at present. This also means that for the same reasons one cannot compute the two loop corrections to the singlet axial vector current conversion or matching function calculated at one loop in [32]. This function allows one to translate results between two different renormalization schemes and in the case of [32] the respective schemes were MS and the symmetric momentum subtraction (SMOM) scheme. As the conversion function is computed from the singlet axial vector current renormalization constants in the two schemes then the finite renormalization constant for restoring four dimensional symmetries through the Larin method will also be needed for this. In turn this requires the finite parts of the three loop Green's function in both schemes. The absence of the three loop symmetric masters means that this cannot be carried out here. Having demonstrated, however, that the non-singlet axial vector finite MS renormalization constant consistently emerges in the symmetric point configuration, instead we merely accept the result (2.18) and use it for our MS computations. In other words we include the extra graphs of Figure 2 where O A s is inserted. The γ-algebra is removed and the Lorentz tensor integrals evaluated by projection and the naive axial vector operator renormalization constant determined which is This is in full agreement with the expression given in [39]. Reproducing this is a check on our different projection strategy since the O(a 2 ) term arises purely from the graphs of  3 Results.
Having discussed the computational strategy in detail we now present our results for the singlet axial vector Green's function. In order to give an indication of the structure of the final expression in four dimensions after the implementation of the Larin method we note that in the Landau gauge for N f = 3 we have for SU (3) where the restriction also means that the symmetric point conditions of (2.4) have been implemented and α is the covariant gauge parameter. As a further check on the computation we note that the symmetry due to the interchange of the external legs of Figure 1 is present. This is the reason why pairs of basis tensors appear in the second and third terms and is consistent with the structure that emerged in the projection method used to evaluate Γ A ns . The full expressions for the Green Clearly the coefficient of the γ 5 γ µ term of the two loop Landau gauge Yang-Mills expression has the smallest magnitude in the MS scheme. In order to gauge the effect the inclusion of the graphs of Figure 2 have in comparison with the flavour non-singlet axial vector Green's function it is instructive to compute the difference (Γ A ns µ − Γ A s µ ). For example when N f = 3 then for SU (3) we have which is independent of α or numerically.
Given this in order to appreciate the effect of the finite renormalization associated with the chiral anomaly in a clear way, it is instructive to plot (3.3) as a function of a dimensionless momentum variable. To achieve this we recall that solving the two loop β-function as a function of the squared momentum Q 2 and the QCD Λ parameter gives where the β-function coefficients are, [61,62,63,64], in the MS scheme and we use the shorthand for the logarithm which is present. In Figure 3 we have plotted the two loop coefficient of γ 5 γ µ for Γ A ns and Γ A s for several values of N f where the dimensionless variable x is defined by x = Q 2 /Λ 2 . At x = 3 for example, the difference between the coefficients of this particular tensor range from around 2% to 7% at x = 3 from N f = 3 to 6 respectively. As an alternative to the explicit values Figure 4 shows the difference for the same values of N f against x. Clearly at very high energy the discrepancy tends to zero.

Discussion.
We have evaluated the two loop quark 2-point function with the flavour singlet axial vector current inserted at non-zero momentum at the non-exceptional symmetric point. The main feature of this result is that unlike the one loop Green's function we have had to ensure that our construction is consistent with the non-conservation of the singlet axial current due to the chiral anomaly. This was a non-trivial task since the use of dimensional regularization means that the purely four dimensional γ 5 -matrix has to reconciled. To achieve this we implemented the Larin method, [39], which was originally developed for the quark bilinear operators. However, the finite renormalization constant associated with each γ 5 dependent operator was determined in [39] for an exceptional momentum configuration. Here we have confirmed that these expressions for the flavour non-singlet operators are independent of the subtraction point for the renormalization scheme that we have used which is MS. If one were to use another scheme such as the momentum subtraction scheme of [40,41] then a different finite renormalization constant would emerge. For the flavour singlet axial vector case one would have to extend our two loop computation to three loops at the symmetric point to verify Larin's two loop finite renormalization constant for that operator. This is due to the non-conserving part of the anomaly equation being proportional to the coupling constant.
Concerning the comparison of the axial vector flavour non-singlet and singlet corrections at two loops, the difference at a representative momentum scale is of the order of a few percent depending on the value of N f . It is not clear whether such a value would make a significant difference to the analysis already carried out in [27,32] for instance. It would depend on whether there is a clear signal in the lattice measurements to differentiate between the various Green's functions. Of course this situation could be improved by extending to three loops which is the next natural step in the process. Aside from the absence of the three loop master integrals at the symmetric point, one would still have to be careful with the treatment of γ 5 in dimensional regularization. For instance, the finite renormalization constant for the axial vector current is known to three loops for the non-singlet case but only two loops for the singlet case. To put the latter on the same level as the former would require a four loop evaluation. While the computational tools are available through the development of the Forcer package, [65,66], the treatment of γ 5 at four loops possibly requires detailed care. Recently there has been progress in this direction through the careful determination of the four loop Standard Model β-functions, [67]. This was achieved by ensuring general quantum field theory consistency conditions were satisfied for the γ 5 sector. By contrast for any treatment of γ 5 that uses the Larin approach where spinor traces involve a large number of γ-matrices it has been noted in [68] that one has always to be fully aware of potential hidden evanescence issues. In our case since we use (2.6) in dimensional regularization the latter point would need to be accommodated at three and higher loops.