Two-loop renormalization and mixing of gluon and quark Energy-Momentum Tensor operators

In this paper we present one- and two-loop results for the renormalization of the gluon and quark gauge-invariant operators which appear in the definition of the QCD energy-momentum tensor, in dimensional regularization. To this end, we consider a variety of Green's functions with different incoming momenta. We identify the set of twist-2 symmetric traceless and flavor singlet operators which mix among themselves and we calculate the corresponding mixing coefficients for the nondiagonal components. We also provide results for some appropriate RI'-like schemes, which address this mixing, and we discuss their application to nonperturbative studies via lattice simulations. Finally, we extract the one- and two-loop expressions of the conversion factors between the proposed RI' and the MSbar schemes. From our results regarding the MSbar-renormalized Green's functions one can easily derive conversion factors relating numerous variants of RI'-like schemes to MSbar. To make our results easily accessible, we also provide them in electronic form, as a Mathematica input file, which is included in the distribution package of this paper.


I. INTRODUCTION
An important open question in Hadronic Physics is the hadron spin decomposition, i.e., the distribution of hadron spin amongst its constituent particles. It is well known, by recent experiments, that contributions to the hadron spin arise not only from valence quarks, but also from polarized gluons, as well as sea quarks. Therefore, it is understood that the complete picture of the spin content of a hadron requires taking into account its nonperturbative nature, including gluon and quark disconnected contributions. Useful quantities which give important input to the study of hadron spin structure are the quark and gluon average momentum fractions [1]. Their nonperturbative determination in nucleons is currently under investigation by a number of research groups [2][3][4] and so far the outcomes are very promising for the correct extraction of the nucleon spin decomposition. However, there are still many challenges that need to be faced, including the complete renormalization of these quantities.
the presence of mixing with other operators. Vigorous efforts in addressing the former include optimized algorithms and increased statistics. The latter issue has additional difficulties: the operators which mix among themselves are typically defined in perturbation theory and may involve gauge-variant terms and ghost fields; thus, their nonperturbative calculation, by compact lattice simulations, is not feasible. There remains still a number of conceptual questions to be resolved before a viable nonperturbative evaluation of mixing effects can be implemented. Studying the mixing pattern in higher orders of perturbation theory can give important guidance for the corresponding elimination of operator mixing nonperturbatively.
In this work, we study the renormalization and mixing of gluon and quark singlet gaugeinvariant operators appearing in the definition of the QCD Energy-Momentum Tensor (EMT). These operators are employed in the calculation of the quark and gluon average momentum fractions in hadrons. In terms of the gluon field A a µ and quark field ψ f , they are defined as [5] 1 : F a µν ≡ ∂ µ A a ν − ∂ ν A a µ − gf abc A b µ A c ν is the field strength tensor and f abc are the SU(N c ) structure constants; to extract the mixing coefficients in an unambiguous way, we consider Green's functions with different incoming momenta.
The renormalization factors of gluon and quark EMT operators can be extracted by studying either the diagonal or the nondiagonal components of the operators. As the EMT operators are traceless, it becomes difficult to disentangle the signal of the diagonal part in lattice simulations from the corresponding pure trace. The mixing pattern of the nondiagonal components is simpler comparing to the diagonal ones. For this reason, we choose to consider only nondiagonal components.
We investigated possible ways of defining an appropriate RI ′ -type scheme, which can be applied in the nonperturbative studies on the lattice. Green's functions of GV operators are difficult to obtain nonperturbatively on the lattice due to a number of obstacles: GV operators (BRST variations 2 and operators which vanish by the equations of motion) are defined in a perturbative manner, including gauge-fixing terms, which are not well-defined in the Landau gauge (they contain terms proportional to 1/α, where α is the gauge-fixing parameter and α = 0 in the Landau gauge) and ghost fields. Such terms cannot be studied by compact lattice simulations. In our study, we discuss some possible approaches to overcome this issue.
A novel aspect of this calculation is the extraction of the mixing matrix to two-loop order. A number of previous perturbative and nonperturbative studies of EMT have been carried out in both continuum and lattice regularizations. A one-loop calculation of the mixing matrix in the continuum is presented in Ref. [6]. Corresponding one-loop calculations on the lattice can be found in Refs. [7][8][9]. A conserved EMT for lattice gauge theories is constructed in Refs. [5,10] to one-loop level. Recent nonperturbative studies of the renormalization of gluon and quark EMT operators have been performed in lattice QCD simulations [2-4, 11, 12]. A promising investigation for determining a conserved EMT nonperturbatively on the lattice is given in Ref. [13].
The outline of this paper is as follows: In Sec. II we provide a theoretical analysis of the renormalization of gluon and quark EMT operators based on the Joglekar-Lee theorems and Ward identities of BRST-invariant operators and of conserved quantities. Sec. III contains the calculation setup including details on the calculated Green's functions, description of the proposed renormalization schemes and the conversion to the MS scheme. Our main results are presented in Sec. IV for the MS-renormalized Green's functions, the renormalization functions and the conversion factors between the RI ′ and the MS schemes. In Sec. V we discuss the application of the proposed RI ′ schemes in the nonperturbative studies on the lattice, while in Sec. VI we conclude.

II. THEORETICAL ANALYSIS
According to the Joglekar-Lee theorems [14], a gauge-invariant operator O can mix with three classes of operators which have the same transformations under global symmetries (e.g. Lorentz, or hypercubic on the lattice, global SU(N c ) transformations, etc.) and whose dimension is lower or equal to that of O: Any other operators which respect the same global symmetries, but do not belong to the above classes, can at most have finite mixing with O [14]. In this respect and given that gluon and quark EMT operators are two-index traceless symmetric of dimension 4, the full set of twist-2 operators which mix among themselves, compatibly with Euclidean rotational symmetry, is the following [5]: where a summation over repeated indices is implied; c a andc a are the ghost and anti-ghost fields, respectively, and S is the QCD action: O 1µν and O 2µν are class G operators, O 3µν and O 4µν belong to class A, and O 5µν is a class B operator. Due to the fact that the above operators are even under charge conjugation (C), they cannot mix with the following C-odd traceless symmetric class G operator: By the same token, operator O 2µν with the symmetrized covariant derivatives replaced by only right or left covariant derivatives, is not considered since it is not a pure eigenstate of charge conjugation. We note that there are no operator candidates containing fermions in class A, as any two-index operator with fermion and anti-fermion fields will lead, under BRST transformations, to an operator of dimension at least 5. Moreover, the only potential class B operator stemming from the fermion EOM is pure trace, and thus it is excluded.
On the lattice, where Lorentz symmetry is replaced by hypercubic symmetry, diagonal (µ = ν) and nondiagonal (µ = ν) components of traceless symmetric operators belong to different representations of the hypercubic group, and thus, they renormalize differently. As we are interested in constructing a renormalization scheme applicable to the lattice, we must renormalize diagonal and nondiagonal components separately. However, their corresponding renormalized Green's functions will be components of a common multiplet in the continuum limit, as it happens in continuum regularizations. In this study, we focus on the renormalization of the nondiagonal components of the EMT operators, because they give more accurate results in lattice simulations when inserted in matrix elements between physical states [3]. From now on, when we refer to O iµν , (i = 1, 2, . . . , 5), it is meant that µ = ν.
Operators O 1µν , O 2µν , . . . , O 5µν have some interesting properties which give us an important input in the study of their renormalization. Let us define the mixing matrix Z as follows: where Here, to simplify the notation we omit the Lorentz indices µ, ν. The sum: in dimensional regularization gives the nondiagonal Belinfante symmetrized EMT [15], which is a conserved quantity 3 . As a consequence, this combination of operators has zero anomalous dimension and thus, it is finite. This is also true for the class B operator O DR 5 . This means that (at least) in the MS scheme, we have: Replacing Eq. (10) into Eqs. (11,12), the following relations between the elements of the mixing matrix are extracted [16]: Furthermore, according to the Joglekar-Lee theorems [14], the mixing matrix (at least) in DR and the MS scheme is block triangular, i.e. class A operators cannot mix with class G operators, and class B operators cannot mix with class G and class A operators; thus, Additional relations between the elements of the mixing matrix can be extracted by studying Ward identities (WIs) which contain operators O i . Let us consider the following WI: where δ BRST is the BRST operator. Because of the BRST invariance of both action and class G, A and B operators (modulo equations of motion), Eq. (22) takes the following form: In momentum space, they read: By replacing the bare operators with the renormalized ones, the above relations also hold (at least) in the MS scheme. This is proved by the following arguments. Let us consider, e.g., the Green's function of operator O 1 in the Y renormalization scheme: Using Eqs. (10,26,27,28), the Green's function take the following form: Operators  Further WIs are derived for 1PI Green's functions with conserved quantities. For example, let us consider the following relation: where diagonal components of O DR i are also involved; σ µν ≡ [γ µ , γ ν ]/2. The quantityT µν is conserved in the limit D → 4. Inserting the above equation under the functional integral of the effective action Γ, a master equation is extracted which is suitable for generating WIs: After some operations, two useful WIs are produced for zero momentum transfer 4 : where (D −1 (q)) ab µν and S −1 (q) are the inverse gluon and quark propagators, respectively. Note that in Eqs. (38, 39), indices µ and ν are taken to be different. The above relations can be useful for the construction of the nondiagonal elements of EMT on the lattice.

III. CALCULATION SETUP
In this section we briefly introduce the setup of our calculation. We provide details on the calculated Green's functions, on the renormalization prescriptions that we use in the presence of operator mixing, and on the conversion factors.

A. Green's functions
In order to study the renormalization of the 5 operators defined in Eqs. (3 -7), we must consider a variety of Green's functions (GFs) with different external elementary fields and different incoming momenta. We consider a total of 5 GFs with external gluon fields for two different choices of incoming momenta, and 5 GFs with external fermion fields for one choice of incoming momenta. Based on the different Lorentz and Dirac structures of the pole terms appearing in each GF, this is the minimum number of GFs, which enable us to extract 25 renormalization conditions for the full determination of the mixing matrix. In particular, the GFs that we consider are 5 : 1. Amputated GFs with two external gluon fields and zero-momentum operator insertion: 2. Amputated GFs with two external gluon fields and nonzero-momentum operator insertion. For simplicity, we may set to zero the momentum of one of the two external gluons: These GFs are needed to disentangle operator O 3 from O 4 as they only differ by a total derivative.

Amputated GFs with a pair of external quark and anti-quark fields and zero-momentum operator insertion:
where a f , b f are color indices in the fundamental representation. These GFs are needed to disentangle the fermion operator O 2 from the remaining gluon operators.
Clearly, the above choices of GFs are not unique; for example, one can choose to consider GFs with external ghost fields. However, such a choice is not optimal for studying these operators in compact lattice simulations 6 .
As we are interested in calculating GFs with external gluon and quark fields, we also need to compute the renormalization functions of the external fields. To this end, the gluon and quark propagators must also be calculated up to two-loops: Explicit results for these GFs can be found in the literature up to four loops [19]. Also, fiveloop results for the renormalization functions of the gluon and quark fields are presented in Ref. [20]. For completeness, we calculate these GFs up to two loops and we make the crosscheck. A difference between these studies and our work is that we present the conversion factors of the gluon and quark fields between RI ′ and MS schemes using independent momentum scales.
There are 1 one-loop and 7 two-loop Feynman diagrams contributing to G q (q), shown in Figs. 1, 2 and 4 one-loop and 23 two-loop Feynman diagrams contributing to G g (q), shown in Figs. 3, 4. The diagrams contributing to G qi (q, q) can be produced by inserting the operator O i in the vertices or in the propagators of the diagrams of Figs. 1, 2. Similarly, the diagrams of G gi (q, −q) and G gi (q, 0) can be produced from the diagrams of Figs. 3, 4 using the same procedure. There is a total of 132, 382, 421 diagrams contributing to G qi (q, q), G gi (q, −q), G gi (q, 0), respectively. Note that a number of duplicate diagrams may arise and must not be double-counted. As is standard practice, we apply the IBP (integration by parts) method to reduce two-loop integrals into nested one-loop master integrals, which are evaluated by a well-known one-loop formula (see Ref. [21]). The most difficult part of this calculation regards the nonscalar integrands stemming from the "diamond"-type diagrams (2-3 of Fig. 2 and 5-11 of Fig. 4); we apply an extension of the scalar recursion formula of Ref. [21], including tensor structures: is a function of k, and q is an external momentum D-vector. For n = 0, Eq. (45) reduces to the scalar formula: Another possibility is to express all integrals in terms of scalar functions of the external momentum by multiplying each integral with the appropriate projectors.

B. Renormalization schemes and conversion factors
In our study we adopt two different renormalization schemes: the MS scheme, which is typically used in phenomenology for the analysis of experimental data, and a regularizationindependent (RI ′ ) scheme, which is more immediate for a lattice regularized theory. The latter scheme is appropriate for renormalizing nonperturbative data taken by lattice simulations. Given that MS is defined in a perturbative manner, the best theoretical approach for taking nonperturbative results in MS is to make use of an intermediate scheme, which is applicable in both perturbative and nonperturbative regularizations, and to match the nonperturbative results from this scheme to MS; RI ′ is an example of such an intermediate scheme. RI ′ -renormalized quantities, calculated on the lattice nonperturbatively, can be converted to the MS counterparts through perturbative "conversion" factors between RI ′ and MS schemes; the conversion factors are regularization-independent and thus, calculable in DR.
Below, we provide our conventions for the definition of the renormalization functions, which relate bare to renormalized fields and parameters of the theory: where A µ is the gluon field, ψ f is the quark field of flavor f , g is the coupling constant, α is the gauge-fixing parameter (α = 0 in the Landau gauge), and µ is a momentum scale.
The index X denotes bare quantities in the X regularization and the index Y denotes renormalized quantities in the Y renormalization scheme. The MS renormalization scaleμ is defined in terms of µ:μ where γ E is Euler's gamma. The renormalization functions for the operators under study have been already defined in Eq. (10), in a 5 × 5 matrix form. The renormalized Green's functions of operators and fields, which are defined in the previous subsection, are given by: In the MS scheme, the renormalization condition is defined (in DR) by imposing that the renormalized Green's functions are finite, when the renormalization functions include only negative powers of ε ≡ (4 − D)/2. In a RI ′ -like scheme, there is, a priori, wide flexibility in defining normalization conditions in Green's functions, especially when operator mixing is present. The possible variants differ only by finite terms. Therefore, it is natural to adopt a minimal prescription, which involves the smallest possible set of operators which can mix; this is usually the mixing set found in MS. Of course, the conditions must be regularization-independent and thus, they must also include any possible additional finite or power-divergent mixing, which is present, e.g., in the lattice regularization. Examples of operators with additional mixing on the lattice are the scalar glueball operator, the scalar quark-antiquark operator, as well as the nonlocal quark bilinears [22][23][24] studied in a chiralsymmetry breaking action. In the present case, such admixtures on the lattice are excluded by hypercubic invariance.
A choice of definition for a RI ′ -like scheme, compatibly with MS, is to consider a 5 × 5 mixing matrix. The elements of the mixing matrix are obtained by imposing 5 × 5 = 25 conditions on Green's functions. This can be done by isolating different Lorentz and Dirac structures of each Green's function. Given that the operators under study are two-index (µ, ν) symmetric, the possible structures for the Green's functions under study with external Similarly, the possible structures for the fermionic Green's functions under study are (for µ = ν): We isolate some of these structures, including those with poles, by selecting specific values for the external momentum and/or the Lorentz components of the external fields; for these specific values we impose that: and similarly, The proposed renormalization conditions for this variant of RI ′ , dubbed RI ′ 1 , are: where the trace in Eqs. (64 -67) is taken over color space (in the adjoint representation), and the trace in Eq. (68) is taken over Dirac and color spaces (in the fundamental representation); the 4-vectorq is the RI ′ renormalization scale.
The above prescription is not a minimal one. From our two-loop results, one can observe that the mixing pattern in the MS scheme reduces to a set of three operators: This was expected from the theoretical analysis presented in Sec. II. Thus, a second choice of definition for a RI ′ -like scheme is to consider a 3 × 3 mixing matrix. Now, we only need nine conditions to identify the renormalization factors. The first two and the last condition of the RI ′ 1 scheme (Eqs. 64, 65, 68) taken for the three operators {O 1 , O 2 , O 6 } can be also the conditions for the RI ′ 2 scheme: This scheme has the advantage of not involving GFs with nonzero momentum operator insertions.
A third choice for defining a RI ′ -like scheme is to impose that the EMT, which is constructed by O 1 , O 2 , O 4 and/or O 5 , is still a conserved quantity after its renormalization in RI ′ scheme. In DR, the conservation gives: and Eqs. (13 -17) will also hold to this version of RI ′ . As we insert five new conditions, we must exclude five conditions from the previous definition of RI ′ 1 scheme. For example, we exclude the operator O 4 from each condition (Eqs. (64 -68)). Similarly, we can define the "conserved" version of the RI ′ 2 scheme. On the lattice, the construction of a conserved EMT is more complex due to the presence of discretization effects, which violate translational invariance. A discussion about the possible ways of applying the conservation properties of EMT on the lattice is given in Sec. V.
To complete the renormalization prescription, we also provide the conditions for the RI ′ renormalization factors of gluon and fermion fields: where the trace in Eq. (72) is taken over color space (in the adjoint representation), and the trace in Eq. (73) is taken over Dirac and color spaces (in the fundamental representation).
Finally, the passage to the MS scheme can be achieved by using the conversion factors between the different versions of RI ′ and the MS scheme, defined as: for the set of mixing operators, the gluon and the quark field, respectively. Note that in Eq.

IV. RESULTS
In this section, we present our one-and two-loop results for the MS-renormalized Green's functions of the operators under study, the renormalization factors and the conversion factors between the different RI ′ versions and the MS scheme, which are all described in the previous section. To facilitate the use of all these results, we provide them also in electronic form, in an accompanying Mathematica input file: "Greens Functions and Conversion Factors.m".
In what follows, is the Casimir operator in the fundamental representation and ζ(n) is the Riemann zeta function.
The expressions for the Green's functions with two external gluon fields and zeromomentum operator insertion (G MS gi (q, −q)) are: The expressions for the Green's functions with two external gluon fields and nonzeromomentum operator insertion (G MS gi (q, 0)) are: • G MS g4 (q, 0) = • G MS g5 (q, 0) = The expressions for the Green's functions with a pair of external quark and anti-quark fields and zero-momentum operator insertion (G MS qi (q, q)) are: • G MS q2 (q, q) =

B. Renormalization factors in the MS scheme
Here, we provide our results for the renormalization factors of operators O i , (i = 1, 2, . . . , 5) in (DR, MS), as a Laurent series in ε ≡ (4 − D)/2: where z MS,DR and As was expected, the mixing matrix is block triangular. Also, there is no mixing between

C. Conversion factors
Here, we present our results for the conversion factors C MS,RI ′ ij between the different versions of RI ′ and MS scheme. For the sake of brevity, we provide only our resulting expressions for the RI ′ 2 scheme and its "conserved" version, RI ′ 2 cons , while the conversion factors for RI ′ 1 can be extracted from Eqs. (86 -98). Our results depend on two renormalization scales: the RI ′ scaleq and the MS scaleμ; we have chosen to keep these two scales distinct, for wider applicability. Note that, whereas the matrix Z MS,DR is necessarily block triangular, no such condition applies to the matrix C MS,RI ′ .
The expressions for the conversion factors C MS,RI ′ 2 ij (i, j = 1, 2, 6) between RI ′ 2 and MS are: • C

MS,RI
The construction of a complete nonperturbative renormalization program, which can eliminate operator-mixing effects, is a difficult task; some well-known complications involve power-divergent mixing of lower dimensional operators, as well as additional, finite mixing contributions associated with the reduction of rotational to hypercubic invariance.
Additional complications arise when gauge-variant operators (BRST variations and EOM operators) are included in the set of operators which mix. Such 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.
There are various approaches, used in the literature, for the study of operator mixing on the lattice. The first one is the perturbative approach, where the renormalization factors are extracted by lattice perturbation theory (see e.g., [7,8] for previous application to the EMT operators and [26] for a general setup). In this approach, an intermediate scheme between lattice and MS is not needed; the derivation of the renormalization factors can be obtained directly in the MS scheme by comparing the lattice bare Green's functions with the corresponding MS-renormalized Green's functions calculated in DR. This approach can give reliable results only when higher-loop terms are negligible. The technical complexity of this approach effectively limits the applicability to one-loop order in most cases. A second approach regards the nonperturbative calculation of the mixing matrix by neglecting gauge-variant operators. These operators do not contribute to the calculation of physical quantities. However, they contribute to the correct extraction of the renormalization factors, which relate the bare to the renormalized operators. This approach can give reliable results only when mixing effects by gauge-variant operators are small enough. A third approach is the combination of approaches 1 and 2 (e.g., [3,27]), where some elements of the mixing matrix are calculated nonperturbatively (e.g. the diagonal elements, or those related to lower-dimensional operators) while the remaining elements are calculated in perturbation theory. The mixing with gauge-variant operators is also omitted.
In order to address the effects of gauge-variant operators, we propose an extension of the above approaches, including a semi-nonperturbative determination of the gauge-variant operators' contributions to the renormalization factors: The gluonic and fermionic part of the gauge-variant operators can be calculated by lattice simulations while the ghost part and/or the gauge-fixing terms can be obtained by lattice perturbation theory.
Our proposed method can be applied in the present study of EMT operators and the nonperturbative calculation of their mixing matrix. The RI ′ 1 scheme, defined in Eqs. (64 -68), is not the optimal one, as it contains three operators with ghost and gauge-fixing terms and it entails the nonperturbative calculation of GFs with nonzero momentum transfer. Such calculation requires the use of two distinct momentum scales for the two external fields and the extrapolation of one momentum to zero, before calculating any renormalization factor. On the contrary, the RI (where T a are the generators of the su(N c ) algebra), the first two terms can be investigated nonperturbatively by lattice simulations, while for the last term we content ourselves with its perturbative study.
We note that the conditions of the RI ′ 2 scheme make use of amputated GFs. This may cause worry for the calculation of the gluonic GFs, where the inverse gluon propagator is needed in the process of the amputation; the lattice simulations commonly employ the Landau gauge, in which the gluon propagator is not invertible. However, setting to zero those components of the renormalization scale, which are parallel to the directions of the two external gluons, the amputation can be performed without inverting the whole gluon propagator.
To explain in more detail the previous argument about the amputation of gluonic GFs in the Landau gauge, we consider the following amputated Green's function of the generic operator O ν 3 ν 4 : where is the gluon propagator in a general gauge (Π T (q 2 ), Π L (q 2 ) are scalar functions of q 2 ). In the Landau gauge (α = 0) the propagator is not invertible; however, if q µ = q ν = 0, then the µ th and ν th rows and columns of the propagator matrix take the value: D µν | qµ=qν =0 = (1/q 2 )δ µν Π T (q 2 ). Of course this is not true for the remaining rows and columns. Thus, the propagator takes a block-diagonal form (e.g., for µ = 1 and ν = 2): The propagator is still not invertible. However, the upper block is invertible and can be inverted separately from the lower block. The latter can be inverted only in a general gauge α = 0. Now, going back to Eq. (128) we observe that we do not need to calculate all the matrix elements of the inverse gluon propagator but only the ν 1 th row (for the calculation of (D −1 ) ν 1 ρ , ρ = 1, 2, 3, 4) and the ν 2 th column (for the calculation of (D −1 ) σν 2 , σ = 1, 2, 3, 4). Thus, we do not need to invert the whole propagator matrix, but only the block containing ν 1 , ν 2 components if the propagator matrix is block diagonal. Choosing q ν 1 = q ν 2 = 0, the propagator is indeed block diagonal, and thus, the amputation can be done successfully without inverting the whole gluon propagator. In passing, we note that the Green's function A µ (q)O ν 3 ν 4 A ν (−q) | q=q in a momentum scaleq withq ν 1 =q ν 2 = 0, cannot be generally amputated in the Landau gauge; it can be amputated only in the case of (µ = ν 1 or µ = ν 2 ) and (ν = ν 1 or ν = ν 2 ). Also, a "democratic" momentum renormalization scale cannot be applied in this case as the amputation cannot be implemented in the Landau gauge for this specific choice.
An alternative choice, "RI ′ 3 ", is to consider nonamputated instead of amputated Green's functions: Then the conditions of Eqs. (69 -71) are replaced by: The second condition (Eq. 134) can be alternatively replaced by: where the renormalization 4-vector scale has one (instead of two) zero component. The third condition (Eq. 135) employing fermionic GFs could also involve amputated GFs, as they have no issues in the amputation process. The conversion factors from RI ′ 3 to the MS scheme coincide with those from RI ′ 2 to the MS scheme. For any other variant of RI ′ scheme (e.g., Eq. (136)) the conversion factors can be easily extracted from our expressions of the MS-renormalized amputated GFs given in Eqs. (86 -98).
Another possibility is to modify the RI ′ 2 renormalization scheme in a way that the sum O 1 + O 2 + O 6 is a conserved quantity. In DR, the sum of the bare operators is conserved. However, this is not true on the lattice, where discretization effects violate translational invariance. A proper definition of the RI ′ renormalization scheme can lead to a conserved sum of the renormalized operators even on the lattice. In the continuum this is simple, as we explained in previous section; it requires the sum of RI ′ -renormalized operators to be equal to the sum of the bare operators. The corresponding lattice condition can be obtained by considering the WIs given in Eqs. (38, 39). These WIs are extracted in DR; however, we can impose their validity also to the RI ′ -renormalized operators on the lattice. To avoid any issues regarding Landau-gauge fixing, these relations will give us three conditions by studying the specific choices of Lorentz and Dirac structures, obtained by the conditions of Eqs.
In these conditions, the nonperturbative calculation of the discretized derivatives of gluon and quark propagators w.r.t. external momentum is needed. As we insert three new conditions, we must exclude three conditions from the previous definition of RI ′ 2 scheme. For example, we exclude the operator O 6 from each condition (Eqs. (69 -71)). In this version of RI ′ , an operator with ghost fields is still involved and thus, a combination of perturbative and nonperturbative results is also needed.
The proposed approach does not completely overcome the mixing effects stemming from gauge-variant operators. There are, in the literature, alternative methods for addressing this mixing. One method entails nonperturbative studies of BRST transformations and GFs with ghost fields implemented in the lattice simulations (see Refs. [17,18] and references therein.). Another method investigates the nonperturbative renormalization of EMT on the lattice (Ref. [13]), in which WIs stemming from the conserved properties of the EMT are used in the framework of thermal QCD with a nonzero imaginary chemical potential. Finally, a gauge-invariant renormalization scheme, such as the X-space scheme [28], which considers gauge-invariant GFs in coordinate space, can be applied without the need of involving any gauge-variant operator; in such a scheme gauge-variant operators do not contribute to the extraction of the renormalization factors of the gauge-invariant operators. This scheme has not been applied before in the calculation of the mixing matrix of EMT operators. At the perturbative level, there is a work in progress by our group [29] in this direction. In this case, the gauge-invariant GFs are constructed using only the gluon and quark EMT operators O 1 and O 2 . Complications arise in this method. In order to calculate the 2 × 2 mixing matrix for the renormalization of O 1 and O 2 , we need a total of four conditions. Three conditions can be obtained by studying two-point GFs between the two mixing operators (between themselves and between each other). However, a complete solution needs a fourth condition which cannot be obtained by any other two-point function. More details can be found in our forthcoming paper [29].

VI. SUMMARY
In this paper, we study the two-loop renormalization and mixing of the gluon and quark EMT operators in dimensional regularization. To this end, we compute a set of two-point Green's functions, renormalized in MS; from our results one may directly deduce the conversion factors between MS and a large variety of RI ′ -like schemes which are appropriate for a nonperturbative extraction of renormalization functions through lattice simulations. We provide the conversion factors relating a number of specific versions of the RI ′ scheme to MS.
We discuss in detail the application of our proposed schemes on the lattice and the construction of a nonperturbative renormalization program for the elimination of the operatormixing effects. In particular, we propose a semi-nonperturbative approach, where perturbative and nonperturbative results are combined; the gluonic and fermionic contributions of gauge-variant operators, which mix with the gauge-invariant EMT operators, can be calculated nonperturbatively, while contributions from the ghost parts can be evaluated by lattice perturbation theory. Also, a different version of RI ′ scheme is proposed, which leads to the determination of a conserved RI ′ -renormalized EMT on the lattice. Our approach, along with the results produced in this paper, can be applied in lattice simulations with the expectation of giving more reliable estimates. A complete elimination of mixing effects is currently under investigation by our group using the X-space renormalization scheme.