Renormalization of QCD in the interpolating momentum subtraction scheme at three loops

We introduce a more general set of kinematic renormalization schemes than the original momentum (MOM) subtraction schemes of Celmaster and Gonsalves. These new schemes will depend on a parameter $\omega$ which tags the external momentum of one of the legs of the $3$-point vertex functions in Quantum Chromodynamics (QCD). In each of the three new schemes we renormalize QCD in the Landau and maximal abelian gauges and establish the three loop renormalization group functions in each gauge. As an application we evaluate two critical exponents at the Banks-Zaks fixed point and demonstrate that their values appear to be numerically scheme independent in a subrange of the conformal window.


I. INTRODUCTION
The renormalization of a renormalizable quantum field theory is a technical exercise which first requires the evaluation of the Feynman diagrams of the relevant divergent n-point Green's functions of the theory to a specific order in perturbation theory. The machinery to subsequently render the theory finite is well established and is completed by encoding the computed renormalization constants in the fundamental renormalization group functions. These functions, such as the β function, which relates to the renormalization of the coupling constant, allows one to determine the behavior and value of all the Green's functions with the renormalization scale and energy. Of course, this has to be tempered by noting that any values can only be reliable within the confines of the perturbative approximation made or equivalently the range of validity of the loop expansion. However, in principle, with sufficient orders in perturbation theory, any estimates should be reasonably reliable. Indeed, the last known term of the perturbative series can sometimes be used as a way of estimating errors. In outlining the general process of renormalization, several more technical issues lurk within the procedure. One of these major areas is that of how the renormalization constants are determined. There are two main aspects to achieving this. First, one has to specify the point where the renormalization constants are to be defined. By this, we mean the momentum configuration of the external legs of the divergent n-point functions. In other words, the values of the square of each external momentum have to be specified. Clearly, there are infinitely many possibilities for such momenta values, but there is a subset which has to be avoided. These are where the sum of a subset of the external momenta is zero. Termed an exceptional momentum configuration, such momenta values can lead to infrared problems in the evaluation or running of the final value of the Green's function. The second general feature of renormalization is that once the renormalization point is specified one has to specify the prescription to defining the renormalization constant associated with each Green's function. This is known as the renormalization scheme. Again, there are infinitely many of ways of achieving this. The most commonly used scheme is the modified minimal subtraction (MS) scheme [1]. It is a variation on the original minimal subtraction scheme denoted by MS, [2,3]. In the MS scheme, the renormalization constants are determined by removing only the divergences with respect to the regulator. The MS scheme is a variant on this where not only are the poles removed but also a specific finite part which is lnð4πe −γ Þ, where γ is the Euler-Mascheroni constant, [1]. The removal of this extra piece appears to improve the convergence of the series for the Green's function, [1].
One major benefit of the MS scheme is that the evaluation of (massless) Feynman diagrams can be completed to very high-loop order. Particularly impressive has been the progress in determining the β function and other renormalization group functions of QCD to five-loop order [4][5][6][7] in the MS scheme. This, together with other five-loop results [8][9][10][11], built on the earlier one-to four-loop MS β-function results of Refs. [12][13][14][15][16][17][18] over a period of around 40 years. At this stage, we make a side remark in relation to the general renormalization process. We noted that the renormalization scheme involves the divergent part of a Green's function, which is quantified by a regulator. The specific regularization was not stated at that point as its nature is irrelevant to the scheme definition. For instance, there are several main regularizations such as Paul-Villars cutoff, spacetime lattice regularization, and dimensional regularization. Each has particular attributes best suited to the problem of interest. For instance, spacetime lattice regularization is particularly useful for studying infrared physics numerically where perturbation theory is not applicable. Although that regularization breaks Lorentz invariance, this can be addressed in order to obtain reliable physics results. Also, it is not an easy regularization to use for analytic perturbative results, and only a few loop orders have ever been determined analytically. Equally, a standard cutoff approach is only useful for a few loop orders and has the drawback of breaking gauge invariance. To circumvent these technical constraints, high-order loop computations are efficiently carried out using dimensional regularization in which the critical (integer) dimension of the quantum field theory is replaced by a continuum spacetime dimension d. This is an analytic continuation with the regulator introduced as a small perturbation from the critical dimension. For QCD, we then have d ¼ 4 − 2ϵ, where ϵ is a complex variable which has magnitude very much less than unity and is the regulator. Unlike the other regularizations, gauge and Lorentz symmetry are not broken.
One feature of the MS scheme that is maybe not immediately obvious but that is exploited in the higherloop computations is that the correct MS renormalization of a Green's function emerges even at an exceptional momentum configuration. For instance, the determination of the coupling constant renormalization requires the evaluation of one of the three three-point functions in QCD. The MS divergences can be extracted by nullifying one of the external momenta. This relegates the evaluation of threepoint functions to the level of a two-point function, which is significantly easier to determine. This observation has been beneficial to finding the QCD β function at various loop orders. In this sense, MS is regarded as a nonkinematical renormalization scheme. It carries no information within the renormalization constants with respect to the location of the subtraction point. By contrast, the renormalization constants of a kinematic scheme contain data corresponding to that point. Several kinematic renormalization schemes have been used to study QCD. For example, there are an on-shell scheme of Ref. [19] and the momentum subtraction (MOM) schemes of Refs. [20,21]. In the latter scheme, a nonexceptional momentum configuration is chosen to evaluate the three-point vertex functions of QCD. The second aspect of the MOM scheme definition is that at that subtraction point the divergences as well as the finite parts of the two-and three-point functions are absorbed into the renormalization constants. Specifically, the subtraction point for the MOM scheme is defined as the point where the squares of the external momenta are all equal [20,21]. This is known as the symmetric (subtraction) point. For QCD, it leads to three separate MOM schemes known as MOMg, MOMh, and MOMq, corresponding to schemes based on the triple gluon, ghost-gluon, and quarkgluon vertices, respectively. In Refs. [20,21], the two-loop renormalization group functions were determined to two loops. More recently, this was extended to three loops in Ref. [22] due to the advance in the determination of the two-loop three-point integrals for nonexceptional momenta configurations [23][24][25][26].
At this point, one natural question arises, which has to do with convergence. It could be the case that the value of a Green's function, for instance, appears to converge more quickly in one scheme than another at the same loop order. If one knew the full series, then there would be no difference in the value of the Green's function at the same evaluation point. However, with a truncated series, the numerical values of the coefficients of the coupling constant differ in different schemes. Also, a Green's function itself is not a physical quantity, and a more proficient way of seeing the scheme dependence is by computing a renormalization group invariant. One simple and accessible example of such a quantity is the set of critical exponents at a phase transition. Critical exponents are related to the underlying renormalization group functions themselves since they are the values of the latter at the fixed point or zero of the β function. In general, the main fixed point is the Wilson-Fisher fixed point in the d-dimensional version of the quantum field theory [27][28][29][30]. For kinematic schemes, these have been studied in QCD in detail at three loops in Ref. [31], for example. A second fixed point, however, is available in QCD, which is the Banks-Zaks fixed point [32]. It has been of intense interest recently [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50], from the point of view of studying various nonkinematical schemes and responding to the new higher-loop QCD data. The Banks-Zaks fixed point is a strictly four-dimensional phenomenon and exists in the conformal window defined by the two-loop β function. For the SUð3Þ color group, the window is 9 ≤ N f ≤ 16, [32], where N f is the number of massless quarks. At the Banks-Zaks fixed point, one can determine critical exponents, and there has been theoretical analysis into extracting accurate estimates which have been shown to be scheme independent [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. Detailed analyses in the main used the MS data of Refs. [4][5][6][7][8][9]11]. However, in Ref. [31], the critical exponents for the three MOM schemes of QCD were examined with the aim of seeing to what extent the renormalization group invariance held as a function of N f in the conformal window. This is not a trivial exercise because of the different structure of the renormalization group functions. For instance, the classes of numbers appearing in the MS scheme renormalization group functions are the rationals and the Riemann zeta function evaluated for integers n ≥ 3. By contrast, the MOM scheme functions in addition to rationals and Riemann zetas involve polylogarithms reflecting the kinematical information of the subtraction point. Therefore, analytically, it is difficult to ascertain the true scheme independence of the results for the Banks-Zaks critical exponents. While this was resolved numerically in Ref. [31] for the region 12 < N f ≤ 16, it was not clear if the consistency for the MOM schemes was a coincidence or not. For instance, there may have been something implicitly related to the MS scheme in the choice of renormalizing at the fully symmetric point for the vertex functions.
Therefore, we have chosen to reexamine the problem of scheme independence of the Banks-Zaks critical exponents in a new set of MOM-related schemes, which we will term the interpolating MOM (iMOM) schemes. There will be an iMOM scheme for each of the three three-point vertices of QCD, which will be termed iMOMg, iMOMh, and iMOMq in direct parallel to the earlier MOM ones with the set denoted by iMOMi. The renormalization group functions in the iMOM schemes will depend on a parameter ω, which is restricted to 0 < ω < 4. It tags one of the external momenta of the three-point vertices, and the concept was introduced in Refs. [51,52] for the specific case of the quark mass operator renormalization only. The QCD Lagrangian itself was not treated in our iMOM scheme defined here since the application in Refs. [51,52] was to assist with matching to a lattice gauge theory computation in which the coupling constant was renormalized in the MS scheme. By contrast, here, we will actually renormalize the QCD Lagrangian itself by defining the scheme originally and determining all the renormalization group functions. The earlier MOM construction of Refs. [20,21] will correspond to ω ¼ 1. By allowing for a parameter, we will be able to quantify where and when the divergence from renormalization group invariance of the Banks-Zaks critical exponents begins in the conformal window. In practical terms, our approach for the exponents will still be numerical and focus on the specific values of ω ¼ 1 2 and 2, although the full analytic renormalization will be for arbitrary ω. These two values will be sufficient to band the MOM ω ¼ 1 value and gauge the tolerance on the exponents. Alternatively, one could regard these two values as a method of error estimation for truncated renormalization group invariant or physical quantities. One popular way of assessing the effect of higher-loop behavior in physical quantities is to evaluate the quantity as a function of the scale μ. Then, the values at 1 2 μ and 2μ are regarded as the error bounds. One disadvantage of this is that the lower value may be beyond the region of perturbative validity and hence unreliable for perturbative error estimates. We will introduce a new approach here, which will have a scheme motivation. With the parameter ω acting as a variation on the subtraction point, its variation between 1 2 and 2 would be a better measure of the errors. In other words, it tracks the effect of the vertex subtraction within the graphs constituting the truncated series of the quantity of interest. The exponents we will compute will demonstrate the effectiveness of this different interpolating approach, and we suggest it would be useful to investigate other QCD quantities with it.
One property of the β function in kinematic schemes is that it is gauge dependent [20,21]. However, in general, the gauge parameter of a linear covariant gauge can be regarded as a second coupling constant. So, at criticality, the renormalization group function of the gauge parameter has to be zero, which corresponds to the Landau gauge. Therefore, all the MOM and iMOM scheme data will be in that gauge. However, there is a second covariant gauge which is of interest, called the maximal Abelian gauge (MAG). It is based on gauge fixing the gluon in the Abelian subgroup of the color group differently from the other gluons. The MAG was introduced in Refs. [53][54][55] to study Abelian monopole condensation since this was believed to be a potential mechanism for color confinement [56][57][58]. Subsequently, it has been studied in that context [56][57][58] and shown to be renormalizable [59][60][61][62][63][64]. The renormalization group functions are available at three loops in MS [65] and the MOM schemes [66,67]. Therefore, it seems natural to examine the Banks-Zaks fixed point in the iMOM schemes in the MAG context as well, and we will carry out that analysis in parallel with the Landau gauge. The aim is to quantify how far the gauge independence of the Banks-Zaks exponents extends into the conformal window. This will also be carried out to the same threeloop order as the linear covariant gauge in order to have as comprehensive overview of the scheme and gauge dependence of the strictly four-dimensional exponents.
The article is organized as follows. We review the formalism required to evaluate and renormalize three-point functions at the nonexceptional momentum configuration in Sec. II for the two gauges we will consider. The definition of the iMOM scheme is given there as well before we record all the renormalization group functions in Sec. III. Subsequently, Sec. IV is devoted to the determination of the two critical exponents of interest at the Banks-Zaks fixed point. We present our conclusions in Sec. V. An Appendix records the tensor bases and projection matrices for each of the three vertices at the interpolating substraction point.

II. FORMALISM
To renormalize QCD in the interpolating MOM schemes in both gauges, we have to introduce a wide body of formalism, such as notation and conventions, as well as the computational tools required for the whole process. There are common aspects of the renormalization for both the Landau and MAG gauges, which can be outlined together. These will center on the two-and three-point functions or self-energy and vertex Green's functions, respectively.
The two-point functions are relatively straightforward to treat in the sense that with the massless fields we use here there is only one scale which is the external momentum. Therefore, in each of the iMOMi schemes, we choose to define the wave function renormalization constants at a point p 2 ¼ −μ 2 , where μ is the mass scale introduced when we dimensionally regularize in d ¼ 4 − 2ϵ dimensions. This scale is necessary to ensure the coupling constant is dimensionless in d dimensions. As the evaluation of the two-point functions is a straightforward exercise for which we use the MINCER algorithm [68,69], we devote the remainder of the present discussion to the vertex function computation. For these, the situation is more involved in that there are two independent external momenta, and so one has to be careful in specifying the point where the three Green's functions are renormalized. For the present discussion, we will focus on the canonical linear covariant gauge and then indicate the modification to the formalism to accommodate the MAG. First, to be more concrete, the three vertex functions we consider are where A a μ , c a , and ψ i are the gluon, Faddeev-Popov ghost, and quark fields, respectively. The color group generators are T a , and the group has structure constants f abc . The indices have the ranges 1 ≤ a ≤ N A and 1 ≤ I ≤ N F for the linear covariant gauge with N F and N A corresponding to the dimension of the fundamental and adjoint representations, respectively. Throughout, we will use notation similar to Refs. [22,67]. In (2.1), we have indicated the momenta of the external legs, which are p, q, and r, but we take the first two as the independent ones and set ð2:2Þ For simplicity, we have factored off the color group structure in each case in (2.1) since to two loops there are no other color tensors. The definition of the iMOMi schemes does not depend on these tensors. The amplitude which remains is a Lorentz tensor, which, although it depends in general on the two independent external momenta, is restricted as indicated to the particular external momentum configuration of the iMOMi subtraction point.
In particular, the squared external momenta are constrained to satisfy where ω is our interpolating parameter. These lead to The latter relations introduce restrictions on the range of validity of ω such that 0 < ω < 4. The lower bound would correspond to an infrared divergence, and the supremum leads to collinear singularities. The original MOM configuration of Celmaster and Gonsalves [20,21] corresponds to ω ¼ 1, which will be used as an internal check throughout. To carry out the renormalization of a Green's function, two aspects have to be considered. The first is the specification of the values of the external momenta. For the iMOMi schemes, we have already indicated this with (2.3) and (2.4). The second is the prescription to define the renormalization constants. We have not phrased this as the way to remove the divergences with respect to the regularization as any scheme has to do this at the very least. The crucial part is the treatment of the nondivergent pieces. As the subtraction point is specified, the complicated function of the external momenta which ordinarily is present within a Green's function is reduced to a particular value of this function. While in our case there are no internal masses upon which this function could also depend, there will be dependence on the parameter ω. Therefore, the finite part of any renormalization constant would correspond to a particular number or for iMOM involve a parameter which sweeps over a range of numbers. For the iMOMi schemes, the subtraction prescription is that the renormalization constants for the two-and three-point functions are chosen so that at the subtraction point there are no OðaÞ corrections where a ¼ g 2 =ð16π 2 Þ. This is within the spirit of the original Celmaster and Gonsalves scheme for the symmetric point case in which ω ¼ 1 [20,21]. Therefore, as noted, there are three different iMOMi schemes for both the linear covariant gauge and MAG which derive from the respective 3-point vertices of the two gauges. The next stage of the vertex function evaluation is the evaluation of the Lorentz tensor amplitudes. For an iMOMi renormalization, this means that the amplitudes have to be computed to the finite parts. The first stage is to decompose these into a set of scalar amplitudes for each vertex by the projection method discussed in Refs. [22,67] where P ggg ðkÞμνσ ðp; qÞ, P qqg ðkÞσ ðp; qÞ, and P ccg ðkÞσ ðp; qÞ are the independent Lorentz tensors which can be built out of the independent external momenta and tensors such as η μν for each vertex. The evaluation at (2.3) and (2.4) is understood on the right-hand side of (2.5). In addition, for the quarkgluon vertex, the spinor indices of the quark fields have to be attached to γ matrices. This means that we include combinations of γ matrices in the set of objects from which the basis tensors of the vertex decomposition are selected. As we dimensionally regularize in order to evaluate the Feynman integrals, we use generalized γ matrices [70][71][72][73][74], which are denoted by Γ ðnÞμ 1 …μ n and defined by Γ μ 1 …μ n ðnÞ ¼ γ ½μ 1 …γ μ n : ð2:6Þ An overall factor of 1=n! is understood in the definition of this totally antisymmetric object. The benefits of using this is that these matrices span the spinor space of the dimensionally regularized theory as well as give a partition since and is a symmetric matrix with d-and ω-dependent entries. Finally, the symmetric matrix M ggg kl is the inverse of N ggg kl [22,67]. The process for the remaining vertices is similar, and we have provided the projection matrices in the Appendix. The only caveat is that for the construction of the quark-gluon projection matrix a trace over the spinor indices is also taken. As the decomposition (2.5) is based on the Lorentz sector, we will use the same projection matrices for the MAG. The final Green's function we will consider is that from which we can extract the renormalization of the quark mass operatorψψ. In particular, we will evaluate hψðpÞ½ψψðrÞψðqÞij ω also to the finite part. The procedure to achieve this is completely parallel to that outlined for the vertex functions and in particular the quark-gluon vertex. The reason for the connection with that specific vertex is the open spinor indices. For the projection, there are two independent tensors in the basis, which actually partitions into two separate sectors due to (2.7). Aside from the unit matrix Γ ð0Þ , the second basis element involves Γ ð2Þμν and is given in the Appendix together with the ω-dependent projection matrix.
The outcome of the projection process is to relegate the Green's functions to a sum over Lorentz scalar amplitudes for each gauge. So far, this process has been general and did not appeal to the the specific Feynman diagrams which comprise the vertex functions. In practice, we generate all the graphs by the FORTRAN-based QGRAF package [75]. The electronic representations of the graphs are then individually passed through the projection algorithm once the color, spinor, flavor, and Lorentz indices have been appended. The consequence is that the amplitude for each Feynman graph is a sum of Feynman integrals, which has scalar products of the external and internal momenta. To two loops, the majority of these scalar products can be rewritten in terms of the propagators of a topology, which allows one to naively reduce the powers of various propagators in an integral. However, at two loops for a three-point function, there is no more than one irreducible scalar product. While this would ordinarily be a difficulty in carrying out the integration, there is an elegant algorithm designed to accommodate this problem, which was developed by Laporta [76]. In essence, the set of propagators of a topology is completed by including an additional propagator in the two-loop three-point case, but more can be added for other Green's functions and loops, which is a process termed completion. Although this extension puts the irreducible numerator propagator on the same footing as the other propagators of the original integral, the new integral cannot correspond to a bona fide Feynman graph. This is not a problem as the Laporta algorithm [76] at large involves integrating by parts completed Feynman integrals. Then, the huge set of resultant linear equations, which as it turns out are redundant, is solved. The end product is that all the integrals contributing to a Feynman graph of the original Green's function can be written as a sum over a relatively small set of what is termed master integrals. Their ϵ expansion has to be determined by explicit evaluation. For our two-loop three-point function, the master integrals have been known for some time [23][24][25][26]. In terms of practically implementing the procedure we have outlined, we have used both versions of the REDUZE package [77,78], which encodes the Laporta algorithm. One advantage is that the output of the reduction to masters can be converted into the symbolic manipulation language FORM [79,80]. This was used intensively with the whole process of evaluating each of our Green's functions to two loops prior to the renormalization in a fully automatic way. The latter is achieved by the method of Ref. [81] by computing all the Green's functions in terms of the bare coupling constant and gauge parameter. Their renormalized counterparts are introduced by the canonical rescaling. The renormalization constants are fixed in a particular scheme by the procedure given earlier in which all the one-loop Green's functions are converted to renormalized variables first before repeating the exercise at two loops. This is carried out first for all the wave function, coupling constant, and gauge parameter renormalizations. Then, to extract the renormalization constant for the quark mass operatorψψ, the Green's function hψðpÞ½ψψðrÞψðqÞij ω is also evaluated in terms of the bare coupling and gauge parameter. These then are rescaled in a particular iMOMi scheme, and the operator renormalization constant for that scheme is determined by ensuring that there are no OðaÞ corrections at the subtraction point. This process is repeated separately for the other two schemes.
To two loops, the iMOM master integrals have been deduced from Refs. [23][24][25][26] and discussed in Ref. [67] for the renormalization of the quark mass operator as a function of ω used for lattice matching. Therefore, we summarize key aspects required for the present work. It transpires that several complicated functions arise in the finite parts of the one-and two-loop master integrals. These are [24,25] which are symmetric in x and y. Here, Li n ðzÞ is the polylogarithm function, and the other functions are defined by The variables x and y are specific to the three-point function with nonzero external p, q, and r since [51,52] x ¼ p 2 r 2 ; y¼ q 2 r 2 ; r 2 ¼ −ωμ 2 : ð2:12Þ Within the interpolating setup, two specific argument combinations emerge for Φ i ðx; yÞ and Ω 2 ðx; yÞ with ðx; yÞ being either ð1; ωÞ or ðω; ωÞ.
As our focus will be on the renormalization of QCD in various schemes, we need to recall several aspects of the formalism. If the subscript o denotes a bare quantity, then the relation between such an object and its renormalized counterpart for a linear covariant gauge is given by where the constant of proportionality is the renormalization constant, α is the gauge parameter of the linear covariant gauge fixing, and the mass scale μ ensures the dimensionlessness of the coupling constant. The definition of the renormalization constants for the MAG is somewhat different and can be found in Ref. [64]. If we introduce the shorthand notation O ¼ψψ ð2:14Þ for the quark mass operator, then we have for both gauges. However, in defining these relations, it is worth stressing that there is an infinite number of different renormalized fields and variables. This is because the renormalized objects are in a particular scheme. With the interpolating momentum subtraction scheme, the parameter ω will play the role of potentially running over a range of different possible schemes. Although we will focus primarily on two values, the key point is that any set of renormalization constants will depend on variables, such as the coupling constant, defined with respect to a scheme. Another set will depend on the variables in a different scheme. The two sets, however, can be related through properties of the renormalization group. We can illustrate this formally for the linear covariant gauges in the iMOMi schemes. From the defining relations of renormalized variables, we can deduce for the coupling constant and gauge parameter. We have been careful to label the quantities in the iMOMi schemes, and this label means that the variables in the particular object are in that scheme. We take the convention that unlabeled quantities are in the MS scheme, which is regarded as the reference scheme. There is no specific reason why we take MS as the reference aside from the fact that it is the most widely used as well as being the scheme for which the renormalization group functions are known to the highest loop order. For completeness, we note that the anomalous dimensions of the fields and gauge parameter in any scheme are related to their renormalization constants by where the same conventions as in Refs. [22,66,67] are used. We have included an α dependence in the β function as it is gauge dependent in general. In the MS scheme, it can be shown that the gauge dependence is absent [2,3]. Our apparently nonstandard relation for γ α ða; αÞ is because we have not made any assumption on the form of Z α . In a linear covariant gauge in our convention, Z α ¼ 1, but this is not true in general. In particular, in the MAG, the corresponding parameter of the off-diagonal gauge fixing is not unity [59][60][61][62][63][64].
Once the renormalization constants have been determined in a scheme, for example, then the renormalization group functions can be determined in that scheme from (2.17). However, there is a second way in which expressions can be deduced without evaluating (2.17). This requires the renormalization group functions in another, or what we will term a base, scheme. Here, this will be the MS scheme, and the iMOMi schemes are the ones for which we wish to deduce the renormalization group functions. To achieve this, we define the respective conversion functions by where ϕ will be used to represent one of the three fields or quark mass operator. The arguments of the conversion functions are in the base or MS scheme. This is an important observation since the iMOMi renormalization constants on the right-hand side have their variables in that scheme. Naively evaluating these functions by taking the explicit renormalization constants will lead to a divergent function with respect to ϵ. To obtain finite expressions, the variables of the iMOMi scheme have to be converted to their MS counterparts. To find the relations between the coupling constant and gauge parameter in the different schemes, the relations (2.16) are solved recursively order by order in perturbation theory. Once the maps are known to the necessary loop order, then one can explicitly find the conversion functions to the same order. Equipped with these, the renormalization group functions in the iMOMi scheme for the linear covariant gauge are given by from the renormalization group equation. The restriction on each equation indicates that once the right-hand side has been computed, which necessarily will be in the base scheme variables, they have to be mapped to the iMOMi variables. This mapping is the (perturbative) inverse of the one used in deriving the conversion functions themselves. While this process is equivalent to the direct evaluation of the renormalization group functions using (2.17), the advantage of the conversion function approach is that if the conversion functions are known at L loops and the renormalization group functions of the base theory are available at (L þ 1) loops then the latter formalism produces the renormalization group functions of the new scheme at (L þ 1) loops without having to perform an explicit (L þ 1)-loop renormalization in the new scheme. We will benefit from this observation by using it to construct the three-loop iMOMi renormalization group functions. We close this section by briefly discussing the definition of the MAG and the relevant aspects of the renormalization of QCD in this gauge with respect to the iMOM schemes. More background details to our conventions can be found in Refs. [65][66][67]. In essence, the gauge is defined by treating the diagonal color fields separately from the offdiagonal ones [59][60][61][62][63][64]. First, we write the group valued gluon field A as where the group generators T a are treated differently. Our notation varies slightly from that of [65], given that we are treating QCD in two different gauges.  [59][60][61][62][63][64]. By contrast, the off-diagonal gluons have a gauge fixing similar to the linear covariant gauge. However, with the different treatment of the gluons, the resultant MAG gauge-fixed Lagrangian involves a significantly larger number of interactions, which includes gluon-ghost and purely quartic ghost interactions. In addition to the two sets of gluons being gauge fixed differently, there are two sets of Faddeev-Popov ghosts, cā and c¯i, which are responsible for the additional type of interactions. Consequently, one has a much larger number of Feynman graphs to evaluate when examining the three-point vertices. One of the complications in dealing with the enlarged basis of terms in the Lagrangian is that the group theory used in evaluating the Feynman graphs becomes more involved. This has been detailed in Ref. [65], and we refer the interested reader to that paper. We have used the same group theory routines in carrying out our MAG computations here as were used in Refs. [65][66][67].
While the MAG has additional interactions and more structure to handle for a computation, its renormalization has several interesting features deriving from the Slavnov-Taylor identities constructed in Ref. [64]. The first is that the diagonal gluons effectively act as a background field similar to the background field gauge [82][83][84][85][86]. For instance, the renormalization of the diagonal gluon is in one-to-one correspondence with the coupling constant renormalization, which is parallel to the situation for the background gluon of the background field gauge. However, for the MAG in the context of the iMOM renormalization, we derive each coupling constant renormalization from the three-point functions directly and extend the MOM approach of Refs. [66,67]. This therefore retains the spirit of the MOM approach in that the iMOM vertex renormalization incorporates kinematical information which will involve ω and hence allows us to explore the dependence on ω. Another consequence of the Slavnov-Taylor identities is that we need only to consider the three vertices defined by the purely off-diagonal gluons as the basis for our MAG iMOM schemes as in the MOM study of Refs. [66,67]. Vertex functions with a diagonal gluon will not lead to any iMOM schemes. Therefore, our focus in the MAG will be on the Green's function parallel to (2.1), except that the adjoint indices of the linear covariant gauge are replaced by their barred counterparts. The color and Lorentz decomposition are completely similar in the MAG. Therefore, we do not need to add anything to this part of the discussion except to note that we have applied the same projection algorithms to the respective MAG vertices. Therefore, after this point, we will not need to distinguish between barred and unbarred adjoint color indices in the discussion. Once all the renormalization constants have been extracted for the MAG, we use the parallel definitions to (2.17) for the off-diagonal gluon and ghosts as well as the other anomalous dimensions. One complication is that, like the MOM results of Refs. [66,67], the explicit expressions will depend on N o A and N d A , for instance. Equally, there are parallel definitions for the corresponding conversion functions (2.18) and (2.19). The only caveat is that we need to distinguish our linear covariant gauge formalism from the corresponding MAG expressions. We do this by denoting the iMOM schemes in the MAG in general by iMOMmi, where the three MAG schemes will be denoted by iMOMmg, iMOMmh, and iMOMmq for the triple off-diagonal gluon, off-diagonal gluon-ghost, and off-diagonal gluon-quark schemes, respectively.

III. RENORMALIZATION GROUP FUNCTIONS
We now discuss the results of our renormalization in both gauges. In this context, it is worthwhile displaying one of the renormalization group functions in analytic form as it illustrates several important features. For instance, the β function for the iMOMh scheme for the SUð3Þ group in the Landau gauge is ð3:1Þ where we have introduced the shorthand notation Analytic expressions for this β function for an arbitrary color group together with all the other renormalization group functions and conversion functions in both gauges are included in the Supplemental Material [87], which also contains the six MS vertex functions. In addition, the decomposition of the vertex functions into the tensor basis are also provided in the data file. While we have presented the SUð3Þ expression, the arbitrary color result in the Landau gauge shares the same property that the one-and two-loop terms are in agreement with the schemeindependent parts which were first computed in Refs. [12][13][14][15]. For a nonzero α, the two-loop term is in fact α dependent, as one would expect. This is because α can be regarded as a second coupling constant and the theorem which shows that the β function is scheme independent at two loops only applies to a theory with a single coupling constant. The main difference of the β function, however, is in the scheme-dependent part, which is clearly ω dependent. Moreover, it contains poles at ω ¼ 0 and 4 as expected from the kinematics of the external momenta of the three-point functions. The same feature occurs in all the other renormalization group functions in this and the other schemes in the scheme-dependent terms for an arbitrary color group. Therefore, in terms of validity, our results are restricted to the range 0 < ω < 4. The remainder of our results will be given in numerical form for both gauges as this is the most straightforward way to compare the same expressions in different schemes as well as to see the effect of varying ω within a renormalization group function. As noted earlier, we will focus specifically on the cases of ω ¼ 1 2 and 2 so that one can gauge how scheme independent the renormalization group invariant critical exponents are. While the difference in these values of ω from the MOM case of unity are not the same, we note that the difference in j ln ωj from zero for these ω values is equivalent. To obtain numerical information, we evaluate each of the special functions which appear at both values of ω, and these are for ω ¼ 2. The first stage in the construction of the three-loop renormalization group functions is the determination of the various conversion functions in different schemes. As our main focus will be on two specific values of ω, we concentrate on these cases. Moreover, all the numerical results presented in the remainder of this section are for the SUð3Þ color group. The respective coupling constant mapping functions for arbitrary α are for ω ¼ 2. We recall that in all the conversion functions the coupling constant and gauge parameter variables are the MS ones. The corresponding expressions for the quark mass operator are more compact in that The same feature emerges here as in the ω ¼ 1 situation in that in each scheme the conversion function is the same for a particular value of ω. This is also the case for the various wave function conversion functions, which are þ Oða 3 Þ: ð3:8Þ They are totally equivalent to the MOM conversion functions of Refs. [20][21][22]. This is not unexpected as the renormalization of the two-point functions here was carried out using the MOM renormalization criterion. As the corresponding conversion functions in the MAG are more involved but are also equivalent to those of the MOM scheme [66,67], we record their explicit iMOMmi numerical values, which are Equipped with the coupling constant and mass operator conversion functions, we can now determine the respective threeloop iMOMi and iMOMmi renormalization group functions. First, in order to compare, we recall both the MS and MOM scheme β functions in numerical form, which are [12][13][14][15][16][20][21][22]66,67] where we have restricted to the Landau gauge for the linear covariant gauge and α ¼ 0 for the MAG for the MOM schemes. Our notation for the renormalization group functions is that the second argument in a MOM or iMOM scheme is the gauge parameter. Such an argument is not needed for the MS β function. In the MS, MOM, and iMOMi scheme renormalization group functions, we recall that the coupling constant and gauge parameter correspond to the variable in the specific case given by the label on the left-hand side. This is to avoid encumbrancing the variables themselves with labels. The two sets of β functions are In the Landau gauge, the correct scheme-independent part emerges in each case with the scheme dependence only present at three loops. However, we can now quantify the effect of varying the parameter ω between 1 2 and 2 in comparison with the symmetric point MOM scheme of Refs. [20,21]. If one examines the N f -independent part of each three-loop term, for example, it is evident that the corresponding coefficients in the ω case lie roughly halfway between the coefficients for ω ¼ 1 2 and 2 in each of the respective iMOMi schemes. This is consistent with our expectations. However, the comparison with, say, the iMOMg ω ¼ 2 β function with the MS case is not appropriate as the coupling constants are not as similar as within an ω comparison. Moreover, β functions are not physically meaningful quantities, and the running of the coupling constant with scale is not the same at three loops in each of the schemes. This is one of the reasons why it is more beneficial to examine critical exponents since they are renormalization group invariants.
We repeat the exercise for the quark mass operator. First, the known MS and MOM scheme results for comparison are [31,[88][89][90][91][92] where we have computed the ω ¼ 1 expressions for the MAG as a corollary of the iMOMmi calculation. Again, the second argument corresponds to the gauge parameter in the non-MS schemes. We recall that the scheme dependence begins at two loops for the quark mass anomalous dimension. The corresponding iMOMi results are for ω ¼ 2. Now, if we compare the N f ¼ 0 two-and three-loop terms within the same vertex scheme, we find the same general trend that was apparent in the β function. First, at two loops, the ω ¼ 1 coefficient lies roughly halfway between the ω ¼ 1 2 and 2 values. In this case, for the schemes based on the ghost and quark vertices, there is a sign change in the coefficient across the range which does not affect this overall observation. At three loops, there is a slackening of the feature of the ω ¼ 1 coefficients lying roughly halfway between the other two values. This is perhaps not surprising as the higherloop expressions are teasing out the effective asymmetry in the range. It will be interesting to see if this is evident in the critical exponent case and how pronounced it is.
Our final results-oriented remarks concern the internal checks on our computation. We have determined all the twoloop renormalization group functions in two ways for each class of gauges. The first is the direct evaluation and renormalization of the two-and three-point functions in the respective schemes by using the various vertex functions as well as the iMOMi renormalization conditions to determine the renormalization constants. From these, we have produced the respective renormalization group functions. There is also an internal check at that stage in that the nonsimple poles in ϵ in the renormalization constants are predetermined by the simple poles at the previous loop orders. An error in these nonsimple poles would have resulted in ϵ divergent renormalization group functions. With the various parameters such as the two separate gauge parameters and the color group Casimirs present, this is a useful check. The second way we have constructed the renormalization group functions is via the conversion function route using the formalism of (2.19) once the various renormalization constants were available. Therefore, we are able to verify that the renormalization group functions computed directly are consistent. This method has the advantage that it automatically produces the three-loop terms of the renormalization group functions. The threeloop iMOMi coupling constant renormalization constants cannot be adduced from these as there is no point of contact with the finite part of the corresponding three-loop iMOMi coupling constant renormalization constants. However, for the wave function renormalizations, we have directly renormalized the respective two-point functions at three loops and checked that they agree precisely with the three-loop expressions constructed using the conversion functions. This represents a useful verification of the conversion function formalism. The final check on our results is that the ω → 1 limit correctly emerges.

IV. CRITICAL EXPONENTS
As an application of our results, we now turn to the evaluation of various critical exponents, which are renormalization group invariants. Therefore, the values of the exponents in any scheme will be the same, provided the renormalization group functions are known to all orders. However, when one has a truncated perturbative expansion, the values for exponents in different schemes will be different with the hope that the discrepancy reduces at high loop order. This is part of our motivation for determining the renormalization group functions in these new schemes here. Our choice of the two specific values of the parameter ω can be used to quantify the variation in some way. First, we summarize the formalism we will apply [31], concentrating on the Landau gauge for illustration. We define the β function in a scheme S by where the coupling constant is understood to be in the scheme S. The associated partial sums or truncated β functions are where there is no OðaÞ term since we are only considering the Banks-Zaks fixed points and not d-dimensional Wilson-Fisher critical points. We denote the critical coupling constant at the Lth -loop order by a L and define it as the solution of the Lth partial sum β S L ða L ; 0Þ ¼ 0 ð4:3Þ in scheme S. As we will be considering the critical exponent associated with the quark mass renormalization, we formally define the anomalous dimension in the scheme S by [31] γ S ψψ ða; 0Þ ¼ in the Landau gauge and the corresponding partial sums by γ S ψψn ða; 0Þ ¼ X n r¼1 γ S r a r : ð4:5Þ The same formalism will also apply to the case of the MAG. Then, from each partial sum, the truncated critical exponents we will evaluate in each of the iMOMi schemes arẽ ð4:6Þ in the notation of Ref. [31]. Here, we use the notation ofω as the exponent corresponding to correction to scaling in order not to confuse it with our choice of the interpolating parameter, which was introduced in Refs. [51,52] in the study of the quark mass renormalization. In defining the exponentsω and ρ with the factors specified in (4.6), we have the same definition of Refs. [31,43]. However, since we used the β-function conventions of Ref. [65], comparing the location of the critical couplings with Ref. [43], there will be a difference of a factor of 4π. This has been absorbed into our coupling constant.
Having introduced the partial sum formalism, we have solved (4.3) for the Banks-Zaks fixed point at two and three loops in each of the iMOMi schemes. We concentrate on the values of ω ¼ 1 2 and 2 and provide numerical results in a series of tables. The critical couplings at two and three loops for both values of ω are given in Table I. Since the β function is scheme independent to two loops inclusive, we have chosen to present the values forω 3 in Table II again for each of the three schemes together. The subsequent three tables (Tables III-V) contain the results for the quark mass exponent in the iMOMq, iMOMh, and   IV. Exponent ρ for the iMOMh scheme for ω ¼ 1 2 , 1, and 2 at two loops (left bank) and three loops (right bank). .   Critical couplings for the MAG in the iMOMmh scheme for ω ¼ 1 2 , 1, and 2 at two loops (left bank) and three loops (right bank).        [31]. We have included them here for comparison with the new values and in order to gauge, for instance, what the range of the exponent is when ω is varied. For all these tables, we have concentrated on the SUð2Þ and SUð3Þ color groups for their two-loop conformal windows, which are 6 ≤ N f ≤ 10 and 9 ≤ N f ≤ 16, respectively. The subsequent tables contain the same data but for the case of the MAG. For instance, the critical couplings for each of iMOMmq, iMOMmh, and iMOMmg schemes are given in Tables VI, VII, and VIII, respectively, while the corresponding values ofω are given in Tables IX, X, and XI for the same three schemes. Finally, Tables XII, XIII, and XIV record the parallel two-and three-loop estimates for the quark mass exponent ρ also at two and three loops with ω values of 1 2 , 1, and 2 for each iMOMm scheme. In several of the iMOMmi schemes, the lower end of the conformal window is at N f ¼ 8 rather than 9 for SUð3Þ, and we have included some data for these schemes. In these cases, we omitted recording two-loop exponents purely for the reason that the values were several orders of magnitude greater than either the subsequent N f estimates or the three-loop value. This no doubt reflects the fact that perturbation theory is probably not truly reliable at that point. However, as the three-loop data for N f ¼ 8 are not unreasonable compared to the N f ¼ 9 value, we have included those for guidance.
To ascertain how the exponentω depends on the scheme, we have plotted the data of Table I for the Landau gauge in Fig. 1 for both SUð2Þ and SUð3Þ. While we determined the exponents for discrete values of N f , we have chosen to TABLE XIII. Critical exponent ρ for the MAG in the iMOMmh scheme for ω ¼ 1 2 , 1, and 2 at two loops (left bank) and three loops (right bank).  present piecewise linear connections between the spot values for this and the other tables, similar to Ref. [50], in order to see any trends more clearly. In Fig. 1, the three schemes are shown in the order of iMOMq, iMOMh, and iMOMg, and on each plot, we have included the corresponding values for the MS and mini-MOM (mMOM) schemes, which were computed in Ref. [31]. As these are the same in each plot, they provide a guide point for comparing the scheme results for each color group. The definition of the mMOM scheme [93] centers on the fact that the ghost-gluon vertex in the Landau gauge does not get renormalized [94]. Therefore, the coupling constant renormalization constant in the mMOM scheme is determined by evaluating the ghost-gluon vertex function with one external ghost leg nullified and then requiring that the nonrenormalization condition is preserved away for all covariant gauges [93]. While the scheme is motivated by a specific property of the gauge theory, the fact that the subtraction point is at an exceptional momentum configuration means that there may be infrared problems unlike the iMOMi schemes. Several general themes appear to emerge for the Landau gauge. First, for each of the three schemes comparing the SUð2Þ to the respective SUð3Þ plot, the MS values diverge from the iMOMi values at about the midpoint of the conformal window as N f decreases. This is not unexpected, as this is around the place where one would expect perturbation theory to become less reliable. The effect is most pronounced for the iMOMg scheme, whereas for the iMOMh scheme, there is a smaller spread across the schemes for relatively low N f . For the iMOMi schemes, the spread over the range of ω is relatively small, which is perhaps surprising for low values of N f but in keeping with our expectations for the higher values where one is in the perturbative region. For SUð3Þ regarding N f ¼ 12 as a rough boundary of this region then estimates forω would appear to be in agreement from all the schemes except the MS one.
One of the features of the three-loop plots forω is the relatively small spread for the range of ω we took. However, the two-loop value for this exponent is scheme independent, and therefore we cannot say whether the momentum subtraction based schemes have any marked difference with the nonkinematic schemes. By contrast, the quark mass anomalous dimension is scheme dependent at two loops, and so we can examine scheme features over several loop orders. We have given plots of ρ 2 and ρ 3 in Figs. 2 and 3 for both SUð2Þ and SUð3Þ. In both figures, the left set is the two-loop values, and the right set is the threeloop ones with Fig. 2 giving the SUð2Þ values. The order of the schemes is the same as Fig. 1. Some general comments are in order at the outset. First, in both cases, the two-loop results at the lower boundary of the conformal window for both groups are clearly unreliable. While this is more marked than for the exponentω, it is clear that there is a huge difference at this end of the window when one compares with the three-loop plots. Even for values of N f above the lower boundary, there is still a large discrepancy between the two-and three-loop cases as the large vertical scale at two loops camouflages the disparity. The other general feature is one shared with the exponentω in that the MS scheme, and to a lesser extent the mMOM scheme, has different behavior compared to the iMOMi schemes as N f decreases. This should also be tempered by the fact that the discrepancy becomes apparent at around N f ¼ 12 for SUð3Þ, which is where perturbation theory is perhaps on the limit of credibility. In terms of the different schemes, there is clearly a parallel structure when comparing each scheme for both groups, which is reassuring. Equally, for the three-loop plots aside from the iMOMg scheme, there is a slight discrepancy between the MS and mMOM scheme estimates and the iMOMi ones, which has a more dramatic low N f limit for the iMOMh case. Clearly, for the iMOMi schemes, there is a significant difference in the limit to the lower window boundary, and therefore no significance should be placed on any estimate of ρ 3 in this case. However, compared with the dramatic change from two loops, it would suggest that a four-loop evaluation could improve the picture for lower N f in the kinematic schemes. One interesting feature emerges if one examines the three-loop plots for both exponents. For the most part, the mMOM plots appear to faithfully track the MS ones. Both schemes are defined in closely similar ways. For instance, the MS scheme is a mass-independent scheme and can be defined at an exceptional momentum configuration as a consequence. Equally, the mMOM scheme has its origin in preserving a property of a vertex function at specific exceptional momentum configuration. However, the exception to parallel behavior for mMOM exponents compared to the MS scheme is the exponentω in the iMOMh scheme. For both color groups, the mMOM exponents are virtually on top of each of the iMOMi schemes for the whole range of the conformal window. It is premature to say that this is a general feature ahead of a four-loop analysis. In other words, it may be the fact that as the mMOM scheme preserves by definition a property of the ghost-gluon vertex then this is reflected in the agreement with the kinematical scheme behavior. Indeed, of the three schemes, the iMOMhω exponents have minimal spread for all N f . Again, this observation needs to be balanced by noting that the iMOMg behavior of ρ 3 is parallel to the MS and mMOM schemes for low N f .
For the MAG, we have provided similar plots for both exponents at three loops in Fig. 4 for the iMOMmq, iMOMmh, and iMOMmg schemes, respectively, for SUð3Þ. Fewer plots for this gauge have been included as there is a strong general similarity with the Landau gauge plots at two loops. Also, we have omitted points for N f ¼ 8 for the two schemes in which there is a window as the relatively large values of the respective exponents would skew the analysis. For instance, if the data from that value of N f were included, then the plots would make it appear that for N f ≥ 9 all the scheme estimates were equivalent. Taking the range as 9 ≤ N f ≤ 16 allows the finer detail to be seen. The MS and mMOM three-loop estimates are again provided for comparison. First, forω, the three-loop plots are virtually the same as for the Landau gauge. Down to around N f ¼ 13, there is little difference between the two reference schemes and the kinematic ones. For the border point of N f ¼ 12, the iMOMmq and iMOMmh schemes are practically the same but back up the earlier observation that this is probably the place where higherorder corrections could remove scheme ambiguity. For the iMOMmg, there is clearly a discrepancy at N f ¼ 13 which is slightly larger than the Landau gauge. This could be due in part to the nature of the MAG in which a subset of gluon fields is isolated in the definition of the gauge itself. However, the closer agreement for the other two schemes would suggest that with higher orders this discrepancy could wash out. The situation for ρ is somewhat different with a different functional behavior for each scheme below N f ¼ 16. However, the general behavior of the three iMOMmi schemes is not dissimilar to that of the Landau gauge plots. It may be that the different behavior lies in the nature of the quantity plotted, which is the quark mass operator. As an operator, it does not have any gluon content in which the split color group property would be significant. However, the plots may be misleading in that the difference between ρ exponent estimates between N f ¼ 13 and 15 range from 5% to 8% Finally, what is noticeable in both gauges is that the behavior of the schemes based on the triple gluon vertex is different from the other two schemes at the lower end of the conformal window. That this is the case in the MAG as well as the Landau gauge suggests that it is a feature of the particular vertex, which has significantly more graphs at two loops, and these are predominantly gluonic. It will not be until three loops that there would be a commensurate number of gluonic contributions to the quark-and ghost-gluon vertex functions with which to compare. It may be then that the behavior at the lower end of the window becomes similar across all three iMOMi schemes.

V. DISCUSSION
We conclude with various remarks. First, we have completed the full renormalization of QCD in a new set of kinematic schemes called iMOM, which extend the MOM schemes of Celmaster and Gonsalves [20,21]. In particular, we have derived all the renormalization group functions at three loops in the three iMOM schemes for both the Landau and maximal Abelian gauges. The schemes depend on a parameter ω, which is restricted to the range 0 < ω < 4. One motivation for introducing the iMOM schemes was to provide a testing ground for evaluating quantities of physical interest in truncated perturbation theory and seeing how far the scheme independence was apparent. The major application of this idea here was to the Banks-Zaks fixed point in the conformal window of QCD. As critical exponents are the evaluation of the renormalization group functions at a fixed point and therefore physical quantities, it was important to study the exponents in the iMOM scheme. The conformal window is such that for values of N f near the upper limit perturbation theory should be a good tool for reliable information. By contrast, as N f reduces inside the window, perturbation will cease to be a reliable guide. However, where the breakdown occurs is not immediately obvious without a numerical analysis. Overall, for both gauges, it seems that from the three-loop results one cannot fully rely on the estimates at N f ¼ 12. This should be qualified by noting that this is from the raw results without resummation to improve convergence.
One aspect of our results which is worth remarking on is in regard to agreement between MS and iMOMi scheme results for the upper end of the window. Numerically, the data in the plots for both sets of schemes lie on top of each other. This strongly suggests bona fide scheme independence. However, this needs to be balanced by the fact that the numerology of the MS and iMOMi schemes is different with the differences appearing first in the schemedependent terms. Therefore, this ought to motivate a future analytic study in order to see if this can be established beyond numerical evidence. Of course, one could extend the iMOMi schemes beyond the appearance of one parameter. For instance, a more general set of schemes could involve two parameters related to the dimensionless variables x and y appearing in the underlying polylogarithms of the master one-and two-loop integrals. While we have not studied this, we would expect the outcome to be the same. In other words, there would be scheme independence. Such a more general set of schemes might be useful in the extension of these ideas to other quantities of physical interest such as the R ratio. There, one has experimental data for which the truncated perturbative expansion is also available but has been computed primarily in the MS scheme. By recomputing in the iMOMi schemes, one could systematically provide bounds on the measured value at a particular energy scale by using the tolerance from the values at ω ¼ 1 2 and 2. This would appear to be a more quantum field theory motivated approach as ω tracks the effect the scheme has through the Feynman diagrams underlying the quantity of interest. This is in contrast to what is currently used in terms of varying the scale itself of where the measurement is made. The toy example of the critical exponents suggests that the scheme variation would be a more robust procedure. Finally, in completing the derivation of the three-loop renormalization group functions in these classes of kinematic schemes and gauges, the natural extension is to proceed to the next-loop order in future work.

ACKNOWLEDGMENTS
This work was carried out in part with the support of the STFC through the Consolidated Grant No. ST/L000431/1 and a studentship (R. M. S.). The authors thank J. Hackett, Professor T. A. Ryttov, and Professor R. Shrock for useful conversations.