F eb 2 02 1 Gauge-invariant Renormalization Scheme in QCD : Application to fermion bilinears and the energy-momentum tensor

M. Costa, 2, a I. Karpasitis, 3, b G. Panagopoulos, c H. Panagopoulos, d T. Pafitis, 5, e A. Skouroupathis, 6, f and G. Spanoudes g Department of Physics, University of Cyprus, Nicosia, CY-1678, Cyprus Department of Mechanical Engineering and Materials Science and Engineering, Cyprus University of Technology, Limassol, CY-3036, Cyprus Present address: The Cyprus Institute, Nicosia, CY-2121, Cyprus Department of Physics, Stanford University, California, 94305–2004, USA Present address: Department of Physics, Utrecht University, 3508 TC Utrecht, the Netherlands Cyprus Ministry of Education, Culture, Sport and Youth, Nicosia, CY-1434, Cyprus

We consider a gauge-invariant, mass-independent prescription for renormalizing composite operators, regularized on the lattice, in the spirit of the coordinate space (X-space) renormalization scheme. The prescription involves only Green's functions of products of gauge-invariant operators, situated at distinct space-time points, in a way as to avoid potential contact singularities. Such Green's functions can be computed nonperturbatively in numerical simulations, with no need to fix a gauge: thus, renormalization to this "intermediate" scheme can be carried out in a completely nonperturbative manner.
Expressing renormalized operators in the MS scheme requires the calculation of corresponding conversion factors. The latter can only be computed in perturbation theory, by the very nature of MS; however, the computations are greatly simplified by virtue of the following attributes: i) In the absense of operator mixing, they involve only massless, two-point functions; such quantities are calculable to very high perturbative order.
ii) They are gauge invariant; thus, they may be computed in a convenient gauge (or in a general gauge, to verify that the result is gauge-independent).
iii) Where operator mixing may occur, only gauge-invariant operators will appear in the mixing pattern: Unlike other schemes, involving mixing with gauge-variant operators (which may contain ghost fields), the mixing matrices in the present scheme are greatly reduced. Still, computation of some three-point functions may not be altogether avoidable.
We exemplify the procedure by computing, to lowest order, the conversion factors for fermion bilinear operators of the formψΓψ in QCD.
We also employ the gauge-invariant scheme in the study of mixing between gluon and quark energy-momentum tensor operators: We compute to one loop the conversion factors relating the nonperturbative mixing matrix to the MS scheme.

INTRODUCTION
Renormalization of composite operators is essential when studying matrix elements and correlation functions in Hadronic Physics. It relates bare quantities of the theory to the physical ones. In order to extract nonperturbative physical results from numerical simulations on the lattice, the construction of a proper nonperturbative renormalization scheme is needed. A requirement for such scheme is to be applicable in both continuum and lattice regularizations in order to make contact with the continuum schemes. Nowadays, the most widely used renormalization scheme in lattice simulations is the modified regularization-independent scheme (RI ′ ) [1, 2]; it considers gauge-variant Green's functions (GFs) of composite operators with external elementary quantum fields in momentum space. This is not a unique nonperturbative scheme. In this paper, we consider an alternative approach, which involves gauge-invariant correlation functions of composite operators in coordinate space. This approach, called "X-space" scheme, has been Electronic addresses: a kosta.marios@ucy.ac.cy, b i.karpasitis@cyi.ac.cy , c gpanago@stanford.edu, d haris@ucy.ac.cy, e pafitis.theodosis@ucy.ac.cy, f askour02@ucy.ac.cy, g spanoudes.gregoris@ucy.ac.cy considered in [3] in the context of lattice studies some years ago. Older investigations of coordinate-space methods can be also found in, e.g., [4]. To date, there are only limited lattice applications of the X-space scheme, mainly regarding the multiplicative renormalization of fermion bilinear operators. Other applications regarding more complex operators, such as the four-fermion operators, have been studied in, e.g., [5]. X-space is a promising nonperturbative scheme for the lattice simulations, especially when one considers further applications involving operators which mix under renormalization. However, some extensions are necessary in order to deal properly with the error in nonperturbative calculations and, most importantly, with operator mixing. In this paper we implement a number of extensions to this effect; the resulting scheme will be referred to as "Gauge-Invariant Renormalization Scheme (GIRS)" to emphasize the property of gauge invariance, which is essential when one studies the renormalization of gauge-invariant operators in the presence of mixing.
GIRS involves two-point GFs of the following form: where O 1 (x), O 2 (y) are gauge-invariant operators at two different spacetime points. In many cases the renormalization factors of the operators in GIRS can be extracted by studying only two-point functions; however, as we conclude by this work, in the presence of mixing the study of three-point functions is also needed in numerous cases. This scheme has a number of advantages which make easier its implementation in the lattice simulations:  [6]) do not contribute in these GFs. This property is very useful especially when studying the renormalization of gauge-invariant operators nonperturbatively by lattice simulations; gauge-variant operators, typically, contain ghost fields and/or gauge-fixing terms, which are defined in perturbation theory and their study is not obvious in a nonperturbative context. Secondly, no gauge fixing is needed in GIRS. When fixing a covariant gauge on the lattice, one encounters the problem of Gribov copies (see, e.g., [7,8]). Employing this scheme, we avoid such a problem. Given that GFs in GIRS are independent of the gauge-fixing parameter, one can perform perturbative calculations in the Feynman gauge, where momentum-loop integrals are simpler. 18,19]. A novel aspect of our calculation is that we provide alternative ways of implementing GIRS: i.e. using specific values of x, y, or integrating over timeslices. Such choices might help to reduce statistical noise in the nonperturbative evaluation. Given that the MS renormalization functions are independent of the spacetime points x, y, the nonperturbative estimates can be checked by verifying this property. Another aspect is that we also provide results using the t'Hooft-Veltman d-dimensional definition of γ 5 .
In the second part of our work, we extend the application of GIRS in the presence of mixing: we study the renormalization and mixing of the gluon and quark parts of the QCD energy-momentum tensor (EMT); this is a subject of research with an increased interest in recent years [20][21][22][23][24]. EMT is relevant to the calculation of the renormalized gluon and quark average momentum fractions, which are involved in the study of hadron spin decomposition [25]. In this study, we consider only nondiagonal elements of EMT, which give a simpler mixing pattern. However, the procedure can be similarly extended to the renormalization of the diagonal elements. In order to establish the required number of renormalization conditions we must also consider three-point GFs. A number of candidate GFs can be employed for this purpose, such as: where O µν i is the gluon or quark energy-momentum tensor operator and O 1 1 (x) =ψ(x)ψ(x) is the scalar bilinear operator; we compute the corresponding conversion factors for some of the most prominent candidates.
The outline of this paper is as follows: Sec. 2 regards the renormalization of fermion bilinear operators in GIRS, while Sec. 3 is devoted to the renormalization and mixing of the quark and gluon EMT operators. In both sections, we provide details on the calculational procedure and we present our tree-level and one-loop results for the bare GFs of operators under study, as well as for the conversion factors between GIRS and MS schemes. Finally, we conclude in Sec. 4 with a summary of our calculation and possible future extensions of our work. We also include an appendix containing details on technical aspects of the calculation (Appendix A).

Details of the Calculation
Definitions of renormalization factors are related to specific renormalization schemes. From the perturbative point of view, these factors depend on properties of the scheme, such as the renormalization scale, the regulator, and the imposed renormalization conditions on GFs. These conditions connect bare and renormalized quantities. There exists a variety of schemes for defining renormalized quantities. A gauge-invariant renormalization scheme, which is not strictly perturbative as is MS, is GIRS. In order to determine renormalization factors within GIRS, we will examine GFs which contain the product of two gauge-invariant composite operators, defined at different spacetime points in the massless limit. Such GFs, computed in a lattice simulation, will lead to a nonperturbative determination of the renormalized composite operators in this scheme. Having two different schemes for the same regulator, we can calculate conversion factors between these two schemes. The renormalized quantities in one scheme can be obtained as functions of their values in the other scheme, with the bare quantities being the same in both schemes. Perturbatively, we use the dimensional regularization (DR) and we calculate these conversion factors, which can take us from GIRS to the MS scheme. Although we will be presenting results only up to first order beyond the leading contribution, the fact that these conversion factors can be calculated in DR, reinforces the prospect of evaluating higher-order contributions. To this end, we calculate, in QCD, the following GF for the case of two local fermion bilinear operators O X , O Y : where O Γ (x) ≡ψ(x)Γψ(x), and Γ = X, Y denotes products of Dirac matrices given in Eqs. (3 -7). X can, in principle, differ from Y . Note that in order to obtain a nonzero result the flavor of the fermion (antifermion) field in O X must coincide with the flavor of antifermion (fermion) field in O Y . Depending on the choice of Γ, the operators behave under Lorentz transformations and under parity as: O γ5γµ =ψγ 5 γ µ ψ axial vector, O σµν =ψσ µν ψ tensor, where σ µν ≡ [γ µ , γ ν ]/2. The above composite operators, appear frequently in the study of the eigenstates of the spectrum of a theory, i.e. hadrons (see, e.g., [26]), and therefore it is essential to impose an appropriate renormalization scheme for them. The GF of these operators diverges as the fields are brought near each other. The choice of Eq. (2) ensures that both the GF and the renormalized operators are independent of the gauge. On the contrary, the GFs ψ(q) x O Γ (x)ψ(q ′ ) , which are typically used for defining the RI ′ scheme, are gauge-dependent, as they involve fundamental fields, which are not gauge-invariant. One may consider both flavor singlet ( 1 N f fψ f Γψ f ) and nonsinglet operators (ψ f Γψ f ′ , f = f ′ ). Actually, the one-loop results do not differ between the two cases and thus we have omitted flavor indices on ψ,ψ. Higher-loop contributions are expected to be different as shown in the diagram of Fig both scalar flavor singlet, they develop a finite vacuum expectation value, which gives a mixing coefficient with the unit operator. To avoid such issues, we use normal ordered operators, i.e., O Γ − O Γ 1 1. In our work, we extract the renormalization factors of O Γ , up to one loop, in both GIRS and MS schemes. We define the renormalized operators and parameters of the theory using the following convention: where g is the coupling constant, and µ is a momentum scale. The superscript B denotes bare quantities in the B regularization (e.g., B = DR, LR, where DR (LR) denotes dimensional (lattice) regularization) and the superscript R denotes renormalized quantities in the R renormalization scheme (e.g., R = GIRS, MS). The MS renormalization scaleμ is defined in terms of µ:μ where γ E is Euler's gamma.
There exist several prescriptions [27] for defining γ 5 in d dimensions, such as the naïve dimensional regularization (NDR) [28], the t'Hooft-Veltman (HV) [29], the DRED [30] and the DREZ prescriptions (see, e.g., Ref. [31]). They are related among themselves via finite conversion factors [32]. In our calculation, we apply the NDR and HV prescriptions. The latter does not violate Ward identities involving pseudoscalar and axial-vector operators in d ≡ 4 − 2 ε dimensions. The metric tensor, η µν , and the Dirac matrices, γ µ , satisfy the following relations in d dimensions: In NDR, the definition of γ 5 satisfies: whereas in HV it satisfies: The renormalization factors in GIRS can be obtained by imposing the following condition: wherez is the GIRS 4-vector position scale (z = (0, 0, 0, 0)). As we are interested in applying GIRS in lattice simulations, the scalez may be chosen to satisfy the condition a ≪ |z| ≪ Λ −1 QCD , where a is the lattice spacing and Λ QCD is the QCD physical scale; this condition guarantees that discretization effects will be under control and simultaneously we will be able to make contact with (continuum) perturbation theory.
There are additional, alternative ways for extracting renormalization factors in GIRS, using variants of the GFs of Eq. (2). An option is to take a Fourier transform of Eq.(2); however, this is not an optimal choice as contact terms arise. A more promising option is to integrate Eq.(2) over three of the four components of the position vector (x − y), while setting the fourth component equal to a reference scalet. For the scalar and pseudoscalar operators, the direction of the unintegrated component is immaterial; for the other operators, there are two possible options, depending on whether this direction coincides or not with one of the indices carried by the operators. Due to the anisotropic lattice employed in simulations, the temporal direction is a special one. In this sense, a natural choice for the componentt is to be temporal; we call this variant t-GIRS. Without loss of generality, we set x = (x 1 , x 2 , x 3 , 0) and y = (0, 0, 0,t); then the renormalization condition for t-GIRS takes the following form: Although the choice of prescription is not unique, we benefit from the fact that each choice depends on only one reference scale. In case the renormalization prescription involves GFs which are not integrated over any spacetime direction, as in Eq. (13), the quantity (x− y) can be chosen to take any 4-vector reference value, e.g. the "democratic" choicez = a(n 1 , n 2 , n 3 , n 4 ), where n 1 = n 2 = n 3 = n 4 .
In what follows, we will provide, to one-loop order, the appropriate conversion factor to the MS scheme for all the above choices; it is given by: Any GF computed nonperturbatively and renormalized (also nonperturbatively) according to any of the above choices should, upon conversion to the MS scheme, lead to renormalized GFs which coincide, regardless of the prescription used; this then provides a very strong consistency check of the nonperturbative results. The ultimate selection of a prescription will depend on the statistical errors involved in each case. For completeness, we note that there is another variant of GIRS given in [19]; in this approach, the average of the position vector (x − y) is taken over the 3-dimensional surface of a hypersphere with radius |x − y|, centered at the origin.

Tree-level order
The first step in our perturbative procedure is to calculate the tree-level value of the GF, using dimensional regularization. This simple exercise serves to explain the procedure applied beyond tree level. The Feynman diagrams contributing to this expectation value are shown in Fig. 2. However, since we will consider the operators O X and O Y to be normal ordered, we discard the second diagram of Fig. 2. We are interested in a massindependent scheme and thus all quark masses are set to zero. Then the tree-level contribution takes the following form: where N c is the number of colors. Integrating over the momenta p, k, the resulting expression is: where a summation over repeated indices µ, ν is understood. One should observe that when X, Y transform differently under rotations and parity, Eq. (17) gives zero. The results 1 for the nonzero cases are listed in Table I, in terms of an overall factor: In order to apply t-GIRS, i.e. to integrate over spatial components, we consider all 8 possibilities for X, Y involving both time-like and space-like directions of Dirac matrices: The only nonvanishing contribution under integration over the spatial components, stems from the case X = Y and µ = ν in Eq. (17). The results for G tree ( x, t) in 4 dimensions, after integrating over the spatial components x, are given, in terms of an overall factor of N c /(π 2 |t| 3 ), in Table II. We note that the integrated GFs for γ t and γ 5 γ t vanish; the consequence of this fact for the corresponding operators will be discussed in the following subsection. Table II: Values of the tree-level contributions to the Green's functions OX ( x, 0)OY ( 0, t) in 4 dimensions, in terms of an overall factor of Nc/(π 2 |t| 3 ), after integrating over x.

One-loop order
At one-loop level, the Feynman diagrams which contribute to the GFs G(x − y) are given in Fig. 3. The corresponding contributions are shown below: Diagram 1: Diagram 2: Diagram 3: Note that for one-loop calculations, the bare and renormalized coupling constants differ only by the factor µ −ε (see Eq. (8)), as only the tree-level value of Z g = 1 + O(g 2 ) contributes in this case (i.e., Z g → 1). Our next step is to verify that the one-loop contribution is indeed gauge-invariant, i.e., to verify that terms depending on the gauge parameter 2 , β, cancel out when we sum the three diagrams. To this end, we may use the following Ward-like identity: If we implement this identity on the vertex of a diagram, the initial diagram gets split into two new ones that have a fermion propagator removed. The double application of the above identity on the three one-loop diagrams is given diagrammatically in Fig. 4. The tadpole diagrams lead to a scaleless expression, which therefore vanishes in dimensional regularization. The remaining contributions (β-terms) cancel against each other when summed. We can thus focus on the Feynman gauge.
Recall that the massless tree-level GFs were proportional to a trace of the form Tr(/ zY / zX). The divergent part of the one-loop contribution is expected to contain the same traces. After some Dirac trace algebra, it turns out that the finite parts of GFs are proportional to the tree level as well. In order to determine the conversion factors between GIRS and MS, we must compare one-loop GFs to the corresponding tree-level ones. However, when integrating over spatial components, the GFs containing the operators with γ t and γ 5 γ t vanish at tree-level and one-loop order; thus, Figure 4: Diagrammatic representation of the application of the Ward Identity on the three one-loop diagrams.
this type of GF cannot be used to evaluate the renormalization factors of these particular operators. Nevertheless, we can adopt the natural definition that these renormalization factors be the same as those of γ i and γ 5 γ i , respectively. Higher-loop contributions involve continuum integrands over more than three momenta, including exponentials in the numerator; these make the calculation process more complicated. Note, however, that the presence of exponentials in two-point functions amounts to a simple Fourier transform, once integrals over inner momenta have been performed; such integrals involve only one external momentum and massless propagators, and they have been studied in various contexts in the literature to higher loop order. In this study, we limit ourselves up to one-loop computations. For a four-loop evaluation, see Ref. [17].
The procedure of formulating a gauge-invariant renormalization scheme entails performing perturbative calculations in the continuum, while the necessary lattice calculations can be performed in a completely nonperturbative way. Still, calculations in lattice perturbation theory can be used to check the validity of nonperturbative methods. Also, a perturbative calculation may be employed in order to reduce cutoff effects present in the nonperturbative estimates. While perturbative calculations are easier to implement in the continuum (and also unavoidable by the very nature of the MS scheme), they become exceedingly complicated on the lattice, and consequently, calculations beyond two loops are practically unfeasible 3 . In addition, even some two-loop calculations become prohibitive for some "improved" actions, such as the ones used in many large-scale simulations nowadays. Thus, in practice, lattice results are typically limited to one loop, and this can lead to large systematic errors. In GIRS, even the one-loop calculation is not trivial since, as mentioned in Sec. 1, "n-loop" (i.e., order g 2n ) calculations involve (n + 1)-loop Feynman diagrams. Also, the presence of exponentials in Feynman integrals makes their computation more complex. A method for calculating similar integrals on the lattice can be found in Ref. [37], where the renormalization of nonlocal operators is studied; however, in that work only one-loop Feynman diagrams were considered.

Results
In this subsection, we present our results (up to one loop) on the bare GFs G(x−y), as well as the conversion factors between all the variants of GIRS and MS scheme. As mentioned above, we employ both NDR and HV prescriptions; HV is more useful for comparison with experimental determinations and phenomenological estimates, while NDR is applied for comparison with previous calculations. The one-loop conversion factor between NDR and HV prescriptions can be extracted from our results.
Our resulting expressions for the five nonvanishing bare GFs are (z ≡ x − y): where C F = (N 2 c − 1)/(2 N c ) is the quadratic Casimir operator in the fundamental representation and hv = 0 (1) for the NDR (HV) prescription of γ 5 . Our results agree with Ref. [3] (in the limit ε → 0) 4 , in the case hv = 0 and N c = 3. These results can be used to derive the renormalization factors in MS and in any variant of GIRS. In particular, the MS-renormalized GFs are the same as the above, once the 1/ε poles are removed and the naïve limit ε → 0 is taken in the remaining terms. The vector and axial vector cases are free of poles, as expected.
For another crosscheck of our results, we extract the multiplicative renormalization factors Z DR,MS Γ , using the following relation 5 : The renormalization factors can be read directly from the bare GFs (Eqs. (22 -26)): where we use the notation {S,P,V,A,T} for {scalar, pseudoscalar, vector, axial-vector, tensor} operators. These factors are in agreement with well-known results in the literature ( [38] and references therein).
Applying the condition of Eq. (13) to the resulting expressions for the bare GFs [Eqs. (22 -26)], we extract the conversion factors between GIRS and MS schemes: C GIRS,MS C GIRS,MS Note that the one-loop results of the bare GFs are proportional to the tree-level ones, and therefore only the length (not the orientation) of x − y is relevant as a renormalization scale. Integrating Eqs. (22 -26)) over spatial components and applying the condition of Eq. (14), we also extract the conversion factors between t-GIRS and MS schemes. As previously mentioned, the integration over spatial components separates further these cases into 8 possibilities (S, P, V t , V i , A t , A i , T ti , T ij ) which depend on whether (x µ − y µ ) is temporal or not, and they correspond to give a vanishing contribution, both at tree level, and one loop; however, it is natural to impose that V t (A t ) has the same renormalization factor as V i (A i ). Thus, below we present results for the remaining 6 operators. Note that for extracting the correct one-loop renormalization factors in GIRS, it is essential to include O(ε 1 ) terms of the tree-level GFs G tree (x − y); we recall that such terms are also necessary in the evaluation of the MS-renormalized GFs. The conversion factors are:

Details of the Calculation
The gauge-invariant part of the QCD traceless symmetric energy-momentum tensor (EMT) [39] contains two flavorsinglet operators, a gluonic O µν 1 and a fermionic O µν 2 : where a summation over repeated indices is implied. These two operators are involved in the structure functions in nucleons [25,[40][41][42]: the gluon operator appears in the leading-twist approximation of the gluon parton distribution function, while the fermion operator is related to the unpolarized quark parton distribution function. Furthermore, their matrix elements are directly related to the gluon and quark average momentum fraction of a nucleon state [25,43]. Also, these operators are connected to the anomalous magnetic moment of the muon [44].
A proper renormalization of the above gluon and fermion operators is required, before one can relate their matrix elements, as extracted from numerical simulations, to physical observables. A difficulty in calculating these renormalization factors is that mixing is present; these operators mix among themselves as they have the same transformations under Euclidean rotational (or hypercubic, on the lattice) symmetry. They also mix with other operators, including gauge-variant operators (BRST variations and operators which vanish by the equations of motion; see [45]). The latter include ghost and gauge-fixing terms, which are well-defined in perturbation theory, while their nonperturbative extensions are not obvious; thus, a nonperturbative study of such terms by compact lattice simulations is problematic. However, implementing a gauge-invariant renormalization scheme, such as GIRS, which involves only gauge-invariant GFs, the gauge-variant operators do not contribute in the renormalization process and thus, they are excluded. In our study, we consider only the nondiagonal elements (µ = ν) of the above operators, which give a reduced set of operators under mixing on the lattice, containing only O µν 1 and O µν 2 . Thus, the renormalization of the nondiagonal components of EMT operators entails the construction of a 2 × 2 mixing matrix, which relates the bare to the renormalized operators: The calculation of the 2 × 2 mixing matrix requires a total of four conditions involving GFs of O µν 1 and O µν 2 . Three different two-point GFs can be constructed by taking vacuum expectation values between the two mixing operators: where ν 1 = ν 2 and ν 3 = ν 4 . By rotational (or just hypercubic) invariance: ji (x − y). As it turns out, the two-point GFs G ij at one-loop level are proportional to the tree-level values of the diagonal elements G ii , with a proportionality factor which is independent of the values of the Lorentz indices ν i ; as a consequence, Eq. (47) can lead to only three independent renormalization conditions. We calculate the above GFs up to one loop in dimensional regularization. The tree-level contributions come from Feynman diagrams shown in Fig. 5, while the one-loop contributions stem from Feynman diagrams of Fig. 6. Note that G tree 12 = 0. From the above GFs we can get three renormalization conditions by requiring absence of poles (for MS) or equality to the corresponding tree-level values in a renormalization position scale (for GIRS). These are not enough for fully determining the mixing matrix Z ij in a univocal way, in either MS or GIRS. In order to impose a fourth condition, we need to compute additional GFs; a most natural choice involves products of operators O µν i with lower-dimensional operators. The procedure, which is also alluded to in [16], has the form of a bootstrap: One starts by renormalizing lowest dimensional operators (where no mixing issues are present), and then proceeds to renormalize operators of increasingly higher dimensionality by requiring finiteness (or some other normalization condition) in the GFs involving products of operators up to that dimensionality. In this case, the only available lower-dimensional gauge-invariant local operators are the fermion bilinears, studied in the previous section.
The simplest GF that one might consider is a two-point function constructed from the product of an EMT operator O µν i and one fermion bilinear O Γ at two distinct spacetime points: However, such a GF vanishes for any choice of Γ to all perturbative orders: the corresponding two-point GFs with a scalar, pseudoscalar or tensor operator give traces of an odd number of Dirac matrices (for massless fermions), while where the flavor of the fermion (antifermion) field in O X must coincide with the flavor of the antifermion (fermion) field in O Y . The above GF depends on two position vectors: (x − w) and (y − w). This fact increases the complexity of the perturbative calculation. This also means that the renormalization factors defined in GIRS may depend on two renormalization 4-vector scales. A priori, a possible way of addressing these issues could be to adopt a zero-momentum insertion for one of the three operators. To do this one needs to perform a 4-dimensional integration over the position vector x, y, or w depending on which operator carries zero momentum. Then, the resulting GF will depend on only one vector. However, such an integration over the whole spacetime causes additional complications: contact terms arise when any two position vectors among {x, y, w}, coincide, giving additional UV-divergences. To eliminate such divergences, further additive renormalizations would be needed. One possible alternative way of simplifying the calculation of G µν i;XY (x−w, y −w), which does not create any contact term is to choose the vector (y − w) to be parallel (or antiparallel) to (x − w) (but x = y). In this way, G µν i;XY will depend on a single position vector. A particular example, which we apply in our calculations, is (y − w) = −(x − w) and (without loss of generality) w = 0: The tree-level and one-loop Feynman diagrams contributing to G µν i;XY , (i = 1, 2), are given in Figs. 7 and 8 -9, respectively. Note that (G µν 1;XY ) tree = 0. A method for calculating the d-dimensional integrals stemming from these Feynman diagrams is described in appendix A.   A most natural set of four conditions for calculating the mixing matrix in GIRS, involving the above GFs, is: where ν 1 = ν 2 and ν 3 = ν 4 . Alternatively, we can replace the second condition (Eq. (52)) with: We only need to make one convenient choice for X and Y. All other choices should be related by conversion factors; it is useful to check that these factors are indeed finite. Note that some choices of X and Y may give vanishing contributions, depending on the transformation properties of X and Y under rotations, parity and charge conjugation. The conditions (51) -(54) can be written in the following explicit form 6 : The four elements of the mixing matrix Z ij can be obtained by solving the above system of 4 equations, once all GFs on the left-hand sides have been determined via numerical simulations. Note that renormalization factors Z X and Z Y , appearing in Eq. (59), are eliminated and thus they do not contribute to the calculation of Z ij .
A proper extension of t-GIRS, analogous to what was defined for the renormalization of fermion bilinears (see Eq. (14)), can be applied in this case, leading to the following conditions: No summation over ν 1 , ν 2 is implied. Depending on the choice of X and Y , the fourth condition of t-GIRS may involve odd integrals, which give zero. For example, using X = Y = 1 1 will lead to two structures: (i) δ ν1ν2 , which vanishes since we study the nondiagonal elements (ν 1 = ν 2 ) of O ν1ν2 i , and (ii) x ν1 x ν2 , which will vanish upon integration over x for any choices of ν 1 , ν 2 . In such cases, an appropriate variant of the fourth condition can be applied; e.g., for the case X = Y = 1 1, a possible alternative condition, in place of Eq. (63), is: (no summation over ν 1 , ν 2 is implied), or for a fixed choice of the 3-vector p. These variant schemes lead to relations analogous to Eqs. (56 -59) for the determination of Z ij . In this work, we present one-loop results for three-point GFs with (X, Y ) = (1 1, 1 1), (γ 5 , γ 5 ), (γ ν3 , γ ν4 ), (γ 5 γ ν3 , γ 5 γ ν4 ), before performing integration over x. It is straightforward to integrate all these results over x; we do so for some specific cases (see next section), which are likely the most appropriate for nonperturbative investigations.
After calculating the mixing matrix, we extract the conversion factors between the GIRS (or t-GIRS) and the MS scheme for the EMT operators; they have a 2 × 2 matrix form: The Z factors in Eq. (67) can be computed in any regularization "B"; given that we are making contact with MS, the most natural choice for B is dimensional regularization.

Results
In this subsection, we present our results (up to one loop) for the bare GFs O ν1ν2 for i, j = 1, 2 and (X, Y ) = (1 1, 1 1), (γ 5 , γ 5 ), (γ ν3 , γ ν4 ), (γ 5 γ ν3 , γ 5 γ ν4 ), as well as the conversion factors between all variants of GIRS and MS. The results are expressed in terms of the following Lorentz structures : s [4] ν1ν2ν3ν4 (x) ≡ δ ν1ν3 x ν2 x ν4 x 2 + δ ν1ν4 x ν2 x ν3 x 2 + δ ν2ν3 x ν1 x ν4 The resulting expressions 7 for the bare GFs are (z ≡ y − x): O ν1ν2 O ν1ν2 O ν1ν2 O ν1ν2 s [4] ν1ν2ν3ν4 (x) + 2 s [2] ν1ν2ν3ν4 (x) − 8 s [5] ν1ν2ν3ν4 (x) The above GFs lead to the following results for Z DR,MS ij : These are consistent with the results found using GFs with elementary external fields [46]. The corresponding MSrenormalized GFs can be obtained by removing 1/ε terms from Eqs. (73) -(84) and by taking the naïve limit ε → 0 in the remaining terms. By solving the system of four equations (56 -59), we extract the conversion factors between different variants of GIRS and MS. Below, we present results for five specific variants of GIRS. All these variants are expected to lead to the same MS-renormalized operators, but their respective numerical signals may favor one variant over the others.
Similarly, by choosing X = Y = γ 5 , we will arrive at the same one-loop conversion factors, since the "hv" coefficient, which would have made a difference, appears only in the one-loop GF O ν1ν2 , which does not contribute in the calculation of the conversion factors to one loop. Nevertheless, numerical data will be much different; thus, this provides for an interesting comparison of the corresponding MS-renormalized GFs, as gotten from the lattice.
Since the three-point functions O ν1ν2 have more than one Lorentz structures, it is necessary to isolate a structure by using projectors or by making specific choices for the indices ν 1 -ν 4 and/or the components of x. In this variant, we isolate the structure s [5] ν1ν2ν3ν4 (x) by choosing ν 1 , ν 2 , ν 3 , ν 4 to be all different. In this case all four components of x must be nonzero. Similarly, the choice of X = γ 5 γ ν3 , Y = γ 5 γ ν4 will give the same conversion factors to one loop.
This variant is similar to GIRS 2 with the difference of isolating the structure s [2] ν1ν2ν3ν4 (x).
The conversion factors for the above variants of GIRS are given below: where i = 1, 2, 3, j = 1, 2 and coefficients c kl are given in Table III for each variant of GIRS. Use of Eq. (65) as an alternative renormalization condition requires the integration of various expressions of the form: over spatial components of x = ( x, t). All these integrals can be performed by using the following generating integral function: where K ν (z) is the modified Bessel function of the second kind and a may take noninteger values. Then, differentiating with respect to a or to individual components of p, we can calculate all necessary integrals arising in t-GIRS.

SUMMARY
In this paper, we study a gauge-invariant, mass-independent renormalization scheme (GIRS) for composite operators, which is applicable in both perturbative and nonperturbative studies. This is an extended version of the coordinate space (X-space) renormalization scheme studied in, e.g., Refs. [3,4]. This scheme involves vacuum expectation values of products of gauge-invariant operators located at different spacetime points. The expectation values are gauge-independent and thus, gauge fixing is not needed in this scheme. Also, gauge-variant operators, which may mix with gauge-invariant operators, do not contribute in such Green's functions; as a consequence, they can be safely excluded, leading to a reduced set of mixing operators. In this work, we apply GIRS in the renormalization of fermion bilinear operators, as well as in the renormalization and mixing of the gluon and quark parts of the QCD energy-momentum tensor (EMT). We propose different variants of GIRS, e.g., using specific values for the position vectors of the operators under study, or integrating over timeslices (t-GIRS), which may lead to reduced statistical noise in the nonperturbative calculations via lattice simulations. We provide results, up to one loop, for the conversion factors between the different versions of GIRS and the MS scheme.
As future plans, GIRS and our proposed variants (t-GIRS, etc) could be immediately implemented on operators of similar kind, e.g., four-fermi operators and supersymmetric operators (Gluino-Glue, Noether supercurrent).
2. Integral I 2 : We introduce Schwinger parameters: After integrating over p 1 and p 2 , we make a change of variables: x 1 = λ 3 /(λ 2 + λ 3 ), x 2 = 1 − λ 2 λ 3 /(λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 ) and x 3 = λ 1 λ 2 λ 3 /(λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 ). Integrals over x 2 and x 3 can be calculated algebraically, while the remaining integration over x 1 cannot be obtained in a closed form (for general values of α i ). The resulting expression takes the following form: where (α 1 + α 2 + α 3 − d) < 0, (−α 1 + d/2) > 0, α 1 > 0. The next step is to examine whether the integration over x 1 and the limit of vanishing regulator (ε → 0, ε = 2 − d/2) can be safely interchanged without leading to divergences. For this check, it is useful to express the hypergeometric function appearing in (A7) as a power series in x 1 and (1 − x 1 ), by applying an appropriate transformation formula (see [47]). In case the interchange is indeed permissible, the integration over x 1 can be performed after taking the limit ε → 0; in all other cases, it turns out that the hypergeometric function can be expressed in terms of simpler functions, allowing a direct integration over x 1 . The Laurent expansion of the hypergeometric function 2 F 1 over ε = 0 has been performed with the help of the mathematica package "HypExp" introduced in Ref. [48].