Gauge-invariant Renormalization Scheme in QCD: Application to fermion bilinears and the energy-momentum tensor
M. Costa, I. Karpasitis, G. Panagopoulos, H. Panagopoulos, T. Pafitis, A. Skouroupathis, G. Spanoudes
aa r X i v : . [ h e p - l a t ] F e b Gauge-invariant Renormalization Scheme in QCD: Application to fermion bilinearsand the energy-momentum tensor
M. Costa,
1, 2, a I. Karpasitis,
1, 3, b G. Panagopoulos, c H. Panagopoulos, d T. Pafitis,
1, 5, e A. Skouroupathis,
1, 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 op-erators, regularized on the lattice, in the spirit of the coordinate space (X-space) renormalizationscheme. 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. SuchGreen’s functions can be computed nonperturbatively in numerical simulations, with no need to fixa gauge: thus, renormalization to this “intermediate” scheme can be carried out in a completelynonperturbative manner.Expressing renormalized operators in the MS scheme requires the calculation of correspondingconversion factors. The latter can only be computed in perturbation theory, by the very nature ofMS; 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 quan-tities are calculable to very high perturbative order.ii) They are gauge invariant; thus, they may be computed in a convenient gauge (or in a generalgauge, to verify that the result is gauge-independent).iii) Where operator mixing may occur, only gauge-invariant operators will appear in the mixingpattern: Unlike other schemes, involving mixing with gauge-variant operators (which may containghost fields), the mixing matrices in the present scheme are greatly reduced. Still, computation ofsome three-point functions may not be altogether avoidable.We exemplify the procedure by computing, to lowest order, the conversion factors for fermionbilinear operators of the form ¯ ψ Γ ψ in QCD.We also employ the gauge-invariant scheme in the study of mixing between gluon and quarkenergy-momentum tensor operators: We compute to one loop the conversion factors relating thenonperturbative mixing matrix to the MS scheme.
1. INTRODUCTION
Renormalization of composite operators is essential when studying matrix elements and correlation functions inHadronic Physics. It relates bare quantities of the theory to the physical ones. In order to extract nonperturbativephysical results from numerical simulations on the lattice, the construction of a proper nonperturbative renormalizationscheme is needed. A requirement for such scheme is to be applicable in both continuum and lattice regularizationsin order to make contact with the continuum schemes. Nowadays, the most widely used renormalization scheme inlattice simulations is the modified regularization-independent scheme (RI ′ ) [1, 2]; it considers gauge-variant Green’sfunctions (GFs) of composite operators with external elementary quantum fields in momentum space. This is not aunique nonperturbative scheme. In this paper, we consider an alternative approach, which involves gauge-invariantcorrelation functions of composite operators in coordinate space. This approach, called “X-space” scheme, has been Electronic addresses: a [email protected], b [email protected] , c [email protected], d [email protected], e pafi[email protected], f [email protected], g [email protected] considered in [3] in the context of lattice studies some years ago. Older investigations of coordinate-space methods canbe also found in, e.g., [4]. To date, there are only limited lattice applications of the X-space scheme, mainly regardingthe 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 schemefor the lattice simulations, especially when one considers further applications involving operators which mix underrenormalization. However, some extensions are necessary in order to deal properly with the error in nonperturbativecalculations and, most importantly, with operator mixing. In this paper we implement a number of extensions to thiseffect; the resulting scheme will be referred to as “Gauge-Invariant Renormalization Scheme (GIRS)” to emphasizethe property of gauge invariance, which is essential when one studies the renormalization of gauge-invariant operatorsin the presence of mixing.GIRS involves two-point GFs of the following form: hO ( x ) O ( y ) i , ( x = y ) , (1)where O ( x ) , O ( y ) are gauge-invariant operators at two different spacetime points. In many cases the renormalizationfactors of the operators in GIRS can be extracted by studying only two-point functions; however, as we conclude bythis work, in the presence of mixing the study of three-point functions is also needed in numerous cases. This schemehas a number of advantages which make easier its implementation in the lattice simulations:1. The GFs under consideration are gauge-invariant. The benefits from this property are: Firstly, when mixingoccurs, the set of operators under mixing is reduced. Gauge-variant operators [BRST (Becchi-Rouet-Stora-Tyutin) variations and operators which vanish by the equations of motion] which mix with gauge-invariantoperators (according to the Joglekar-Lee theorems [6]) do not contribute in these GFs. This property is veryuseful especially when studying the renormalization of gauge-invariant operators nonperturbatively by latticesimulations; gauge-variant operators, typically, contain ghost fields and/or gauge-fixing terms, which are definedin perturbation theory and their study is not obvious in a nonperturbative context. Secondly, no gauge fixing isneeded 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 thegauge-fixing parameter, one can perform perturbative calculations in the Feynman gauge, where momentum-loopintegrals are simpler.2. Contact terms are automatically excluded ( x = y ) in contrast to the standard renormalization schemes inmomentum space (e.g., RI ′ scheme).3. In the absence of operator mixing, perturbative calculations in GIRS involve diagrams with only one incon-ing/outgoing momentum. Given that one may also adopt a massless renormalization scheme, one can makeuse of techniques for evaluating such diagrams which have been developed to very high perturbative order (see,e.g., [9–15]). Analogous techniques can be used even in the presence of mixing.4. The fact that GIRS renormalization functions can be fully obtained non-perturbatively, without recurrence tolattice perturbation theory, has the consequence that conversion to the modified minimal subtraction (MS)scheme entails only continuum perturbative calculations.There are also some disadvantages of GIRS:1. Computations in GIRS, at a given order in perturbation theory, involve diagrams with one more loop.2. The cost for not generating contact terms is the presence of exponentials in Feynman integrals, which makestheir computation somewhat more complex.3. When mixing occurs, one must often study also ( n > ′ and GIRS are used as intermediate schemes in which renormalization functions can be directly obtainedvia lattice simulations. The ultimate goal is to obtain renormalized GFs in the MS scheme, which is the standardscheme used in the analysis of experimental data. For this purpose, one must compute appropriate conversion factors,which are finite and regularization independent.Our work is divided into two parts: The first part focuses on the employment of GIRS in the multiplicativerenormalization of fermion bilinear operators. This serves both as an example for describing the renormalizationprocedure in GIRS and as a necessary ingredient in possible variants of GIRS appearing in the second part. Thefermion bilinear operators have already been studied in GIRS in both continuum [16, 17] and lattice regularizations [3,18, 19]. A novel aspect of our calculation is that we provide alternative ways of implementing GIRS: i.e. usingspecific values of x, y , or integrating over timeslices. Such choices might help to reduce statistical noise in thenonperturbative 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 provideresults using the t’Hooft-Veltman d-dimensional definition of γ .In the second part of our work, we extend the application of GIRS in the presence of mixing: we study therenormalization and mixing of the gluon and quark parts of the QCD energy-momentum tensor (EMT); this is a subjectof research with an increased interest in recent years [20–24]. EMT is relevant to the calculation of the renormalizedgluon 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, theprocedure can be similarly extended to the renormalization of the diagonal elements. In order to establish the requirednumber of renormalization conditions we must also consider three-point GFs. A number of candidate GFs can beemployed for this purpose, such as: hO µνi ( x ) O y ) O z ) i , where O µνi is the gluon or quark energy-momentum tensoroperator and O x ) = ¯ ψ ( x ) ψ ( x ) is the scalar bilinear operator; we compute the corresponding conversion factors forsome 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 GFsof operators under study, as well as for the conversion factors between GIRS and MS schemes. Finally, we concludein Sec. 4 with a summary of our calculation and possible future extensions of our work. We also include an appendixcontaining details on technical aspects of the calculation (Appendix A).
2. RENORMALIZATION OF FERMION BILINEAR OPERATORS IN GIRS2.1. Details of the Calculation
Definitions of renormalization factors are related to specific renormalization schemes. From the perturbative pointof view, these factors depend on properties of the scheme, such as the renormalization scale, the regulator, and theimposed renormalization conditions on GFs. These conditions connect bare and renormalized quantities. There existsa variety of schemes for defining renormalized quantities. A gauge-invariant renormalization scheme, which is notstrictly perturbative as is MS, is GIRS. In order to determine renormalization factors within GIRS, we will examineGFs which contain the product of two gauge-invariant composite operators, defined at different spacetime pointsin the massless limit. Such GFs, computed in a lattice simulation, will lead to a nonperturbative determination ofthe renormalized composite operators in this scheme. Having two different schemes for the same regulator, we cancalculate conversion factors between these two schemes. The renormalized quantities in one scheme can be obtained asfunctions 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 tothe MS scheme. Although we will be presenting results only up to first order beyond the leading contribution, the factthat 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 : G ( x − y ) ≡ hO X ( x ) O Y ( y ) i , x = y, (2)where O Γ ( x ) ≡ ¯ ψ ( x )Γ ψ ( x ), and Γ = X, Y denotes products of Dirac matrices given in Eqs. (3 - 7). X can, inprinciple, 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 operatorsbehave under Lorentz transformations and under parity as: O ψψ scalar , (3) O γ = ¯ ψγ ψ pseudoscalar , (4) O γ µ = ¯ ψγ µ ψ vector , (5) O γ γ µ = ¯ ψγ γ µ ψ axial vector , (6) O σ µν = ¯ ψσ µν ψ tensor , (7)where σ µν ≡ [ γ µ , γ ν ] /
2. The above composite operators, appear frequently in the study of the eigenstates of thespectrum of a theory, i.e. hadrons (see, e.g., [26]), and therefore it is essential to impose an appropriate renormalizationscheme 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 h ψ ( q ) P x O Γ ( x ) ¯ ψ ( q ′ ) i , which are typically used for defining the RI ′ scheme, are gauge-dependent, as they involvefundamental fields, which are not gauge-invariant.One may consider both flavor singlet ( N f P f ¯ ψ f Γ ψ f ) and nonsinglet operators ( ¯ ψ f Γ ψ f ′ , f = f ′ ). Actually, theone-loop results do not differ between the two cases and thus we have omitted flavor indices on ψ , ¯ ψ . Higher-loopcontributions are expected to be different as shown in the diagram of Fig. 1. Note that if the operators in Eq. (2) are X Y
Figure 1: A two-loop order Feynman diagram contributing to the expectation value hO X ( x ) O Y ( y ) i for flavor singlet operators O X and O Y . A wavy (solid) line represents gluons (quarks). both scalar flavor singlet, they develop a finite vacuum expectation value, which gives a mixing coefficient with theunit operator. To avoid such issues, we use normal ordered operators, i.e., O Γ − hO Γ i O Γ , up to one loop, in both GIRS and MS schemes. Wedefine the renormalized operators and parameters of the theory using the following convention: O R Γ ≡ Z B,R Γ O B Γ , g R ≡ µ ( d − / (cid:0) Z B,Rg (cid:1) − g B , (8)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 renormalizationscale ¯ µ is defined in terms of µ : ¯ µ ≡ µ (cid:18) πe γ E (cid:19) / , (9)where γ E is Euler’s gamma.There exist several prescriptions [27] for defining γ 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 andHV prescriptions. The latter does not violate Ward identities involving pseudoscalar and axial-vector operators in d ≡ − ε dimensions. The metric tensor, η µν , and the Dirac matrices, γ µ , satisfy the following relations in d dimensions: η µν η µν = d, { γ µ , γ ν } = 2 η µν . (10)In NDR, the definition of γ satisfies: { γ , γ µ } = 0 , ∀ µ, (11)whereas in HV it satisfies: { γ , γ µ } = 0 , µ = 1 , , , , [ γ , γ µ ] = 0 , µ > . (12)The renormalization factors in GIRS can be obtained by imposing the following condition: hO GIRS X ( x ) O GIRS Y ( y ) i| x − y =¯ z ≡ Z B, GIRS X Z B, GIRS Y hO BX ( x ) O BY ( y ) i| x − y =¯ z = hO GIRS X ( x ) O GIRS Y ( y ) i tree | x − y =¯ z , (13)where ¯ z is the GIRS 4-vector position scale (¯ z = (0 , , , z may be chosen to satisfy the condition a ≪ | ¯ z | ≪ Λ − , where a is the lattice spacingand Λ QCD is the QCD physical scale; this condition guarantees that discretization effects will be under control andsimultaneously 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 GFsof Eq. (2). An option is to take a Fourier transform of Eq.(2); however, this is not an optimal choice as contactterms 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 scale ¯ t . 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 theanisotropic lattice employed in simulations, the temporal direction is a special one. In this sense, a natural choice forthe component ¯ t is to be temporal; we call this variant t-GIRS. Without loss of generality, we set x = ( x , x , x , y = (0 , , , ¯ t ); then the renormalization condition for t-GIRS takes the following form: Z d ~x h O t − GIRS X ( ~x, O t − GIRS Y ( ~ , ¯ t ) i = Z d ~x h O t − GIRS X ( ~x, O t − GIRS Y ( ~ , ¯ t ) i tree . (14)Although the choice of prescription is not unique, we benefit from the fact that each choice depends on only onereference scale. In case the renormalization prescription involves GFs which are not integrated over any spacetimedirection, as in Eq. (13), the quantity ( x − y ) can be chosen to take any 4-vector reference value, e.g. the “democratic”choice ¯ z = a ( n , n , n , n ), where n = n = n = n .In what follows, we will provide, to one-loop order, the appropriate conversion factor to the MS scheme for all theabove choices; it is given by: C GIRS , MSΓ ≡ Z DR , MSΓ Z DR , GIRSΓ = Z LR , MSΓ Z LR , GIRSΓ . (15)Any GF computed nonperturbatively and renormalized (also nonperturbatively) according to any of the above choicesshould, upon conversion to the MS scheme, lead to renormalized GFs which coincide, regardless of the prescriptionused; this then provides a very strong consistency check of the nonperturbative results. The ultimate selection of aprescription 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 theposition vector ( x − y ) is taken over the 3-dimensional surface of a hypersphere with radius | x − y | , centered at theorigin. The first step in our perturbative procedure is to calculate the tree-level value of the GF, G ( x − y ) ≡ hO X ( x ) O Y ( y ) i ,using dimensional regularization. This simple exercise serves to explain the procedure applied beyond tree level. TheFeynman diagrams contributing to this expectation value are shown in Fig. 2. However, since we will consider the Figure 2: Tree-level Feynman diagrams contributing to the expectation value hO X ( x ) O Y ( y ) i . A solid line represents quarks. operators O X and O Y to be normal ordered, we discard the second diagram of Fig. 2. We are interested in a mass-independent scheme and thus all quark masses are set to zero. Then the tree-level contribution takes the followingform: G tree ( x − y ) = N c Z d d p (2 π ) d Z d d k (2 π ) d e i ( k − p )( x − y ) p k Tr( /pY /kX ) , (16)where N c is the number of colors. Integrating over the momenta p , k , the resulting expression is: G tree ( x − y ) = N c (Γ(2 − ε )) π − ε ( x µ − y µ ) ( x ν − y ν )(( x − y ) ) − ε Tr( γ µ Y γ ν X ) . (17)where a summation over repeated indices µ , ν is understood. One should observe that when X , Y transform differentlyunder rotations and parity, Eq. (17) gives zero. The results for the nonzero cases are listed in Table I, in terms ofan overall factor: N ( x − y ) ≡ N c (Γ(2 − ε )) / (4 π − ε (( x − y ) ) − ε ): The different definitions of γ give identical results at tree level, since only the first four components of the vectors x and y can benonzero. X Y G tree ( z ) / N ( z )1 1 4 z γ γ − z γ ν γ ν z ν z ν − δ ν ν z γ γ ν γ γ ν z ν z ν − δ ν ν z σ ν ν σ ν ν δ ν ν z ν z ν − δ ν ν z ν z ν − δ ν ν z ν z ν + δ ν ν z ν z ν ) − δ ν ν δ ν ν − δ ν ν δ ν ν ) z Table I: Values of the tree-level contributions to the Green’s functions (cid:10) O X ( x ) O Y ( y ) (cid:11) in 4 − ε dimensions in terms of anoverall factor of N ( z ) ≡ N c (Γ(2 − ε )) / (4 π − ε ( z ) − ε ), where z ≡ x − y . In order to apply t-GIRS, i.e. to integrate over spatial components, we consider all 8 possibilities for X , Y involvingboth time-like and space-like directions of Dirac matrices: , γ , γ t , γ i , γ γ t , γ γ i , γ t γ i , γ i γ j , where i = j = t = i .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 , aregiven, in terms of an overall factor of N c / ( π | t | ), in Table II. We note that the integrated GFs for γ t and γ γ t vanish;the consequence of this fact for the corresponding operators will be discussed in the following subsection. X = Y R d ~x G tree ( ~x , t ) / [ N c / ( π | t | )]1 1 / γ − / γ t γ i − / γ γ t γ γ i − / σ ti / σ ij − / (cid:10) O X ( ~x, O Y ( ~ , t ) (cid:11) in 4 dimensions, in terms of anoverall factor of N c / ( π | t | ), after integrating over ~x . At one-loop level, the Feynman diagrams which contribute to the GFs G ( x − y ) are given in Fig. 3. YX d1 Y d2 X X Y d3 Figure 3: One-loop order Feynman diagrams contributing to the expectation value hO X ( x ) O Y ( y ) i . A wavy (solid) line representsgluons (quarks). The corresponding contributions are shown below:Diagram 1: d = Z d d k (2 π ) d Z d d p (2 π ) d Z d d q (2 π ) d g B N c − e ik ( x − y ) (cid:20) p − q ) p ( p − k ) q ( q − k ) T r ( X/pγ µ /qY ( /q − /k ) γ µ ( /p − /k )) (cid:21) , (18)Diagram 2: d = Z d d k (2 π ) d Z d d p (2 π ) d Z d d q (2 π ) d g B N c − e ik ( x − y ) (cid:20) p − q ) ( p ) ( p − k ) q T r ( X/pγ µ /qγ µ /pY ( /p − /k )) (cid:21) , (19)Diagram 3: d = Z d d k (2 π ) d Z d d p (2 π ) d Z d d q (2 π ) d g B N c − e ik ( x − y ) (cid:20) p − q ) ( p ) ( p − k ) q T r ( X ( /p − /k ) Y /pγ µ /qγ µ /p ) (cid:21) . (20)Note that for one-loop calculations, the bare and renormalized coupling constants differ only by the factor µ − ε (seeEq. (8)), as only the tree-level value of Z g = 1 + O ( g ) contributes in this case (i.e., Z g → , β , cancel out when we sum the three diagrams. To this end, we may use the followingWard-like identity: /k/p ( /p + /k ) = 1 /p − /p + /k (21)If we implement this identity on the vertex of a diagram, the initial diagram gets split into two new ones thathave a fermion propagator removed. The double application of the above identity on the three one-loop diagrams isgiven diagrammatically in Fig. 4. The tadpole diagrams lead to a scaleless expression, which therefore vanishes indimensional regularization. The remaining contributions ( β -terms) cancel against each other when summed. We canthus 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 ofthe one-loop contribution is expected to contain the same traces. After some Dirac trace algebra, it turns out thatthe finite parts of GFs are proportional to the tree level as well. In order to determine the conversion factors betweenGIRS and MS, we must compare one-loop GFs to the corresponding tree-level ones. However, when integrating overspatial components, the GFs containing the operators with γ t and γ γ t vanish at tree-level and one-loop order; thus, β = 0 (1) corresponds to the Feynman (Landau) gauge. X Y
X Y X YX YX Y d2 X YX Y
X Y
X YX Y d3 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, wecan adopt the natural definition that these renormalization factors be the same as those of γ i and γ γ i , respectively.Higher-loop contributions involve continuum integrands over more than three momenta, including exponentials inthe numerator; these make the calculation process more complicated. Note, however, that the presence of exponentialsin 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 variouscontexts in the literature to higher loop order. In this study, we limit ourselves up to one-loop computations. For afour-loop evaluation, see Ref. [17].The procedure of formulating a gauge-invariant renormalization scheme entails performing perturbative calculationsin 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 natureof the MS scheme), they become exceedingly complicated on the lattice, and consequently, calculations beyond twoloops are practically unfeasible . 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 typicallylimited to one loop, and this can lead to large systematic errors. In GIRS, even the one-loop calculation is not trivialsince, as mentioned in Sec. 1, “ n -loop” (i.e., order g n ) calculations involve ( n + 1)-loop Feynman diagrams. Also,the presence of exponentials in Feynman integrals makes their computation more complex. A method for calculatingsimilar 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. In this subsection, we present our results (up to one loop) on the bare GFs G ( x − y ), as well as the conversion factorsbetween 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 isapplied for comparison with previous calculations. The one-loop conversion factor between NDR and HV prescriptionscan be extracted from our results. The only three-loop calculations on the lattice existing thus far regard “vacuum” diagrams, that is diagrams without external lines andmomenta [33–35]. In addition, stochastic perturbation theory has been carried out to higher loops [36].
Our resulting expressions for the five nonvanishing bare GFs are ( z ≡ x − y ): hO DR x ) O DR y ) i = N c (Γ(2 − ε )) π − ε ( z ) − ε " g C F π (cid:18) ε + 2 + 6 ln(¯ µ z ) + 12 γ E −
12 ln(2) (cid:19) + O ( ε , g ) + O ( g ) , (22) hO DR γ ( x ) O DR γ ( y ) i = −hO DR x ) O DR y ) i −
16 hv N c (Γ(2 − ε )) π − ε ( z ) − ε g C F π + O ( g ) , (23) hO DR γ ν ( x ) O DR γ ν ( y ) i = N c (Γ(2 − ε )) π − ε ( z ) − ε (cid:16) z ν z ν − δ ν ν z (cid:17) g C F π + O ( ε , g ) + O ( g ) ! , (24) hO DR γ γ ν ( x ) O DR γ γ ν ( y ) i = hO DR γ ν ( x ) O DR γ ν ( y ) i + 8 hv N c (Γ(2 − ε )) π − ε ( z ) − ε (cid:16) z ν z ν − δ ν ν z (cid:17) g C F π + O ( g ) , (25) hO DR σ ν ν ( x ) O DR σ ν ν ( y ) i = N c (Γ(2 − ε )) π − ε ( z ) − ε h δ ν ν z ν z ν − δ ν ν z ν z ν − δ ν ν z ν z ν + δ ν ν z ν z ν ) − ( δ ν ν δ ν ν − δ ν ν δ ν ν ) z i × " g C F π (cid:18) − ε + 6 − µ z ) − γ E + 4 ln(2) (cid:19) + O ( ε , g ) + O ( g ) , (26)where C F = ( N c − / (2 N c ) is the quadratic Casimir operator in the fundamental representation and hv = 0 (1) forthe NDR (HV) prescription of γ . Our results agree with Ref. [3] (in the limit ε → , 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, theMS-renormalized GFs are the same as the above, once the 1 /ε poles are removed and the na¨ıve limit ε → Z DR , MSΓ , using thefollowing relation : (cid:10) O MS X ( x ) O MS Y ( y ) (cid:11) = Z DR , MS X Z DR , MS Y (cid:10) O DR X ( x ) O DR Y ( y ) (cid:11) | ε → . (27)The renormalization factors can be read directly from the bare GFs (Eqs. (22 - 26)): Z DR, MS S = 1 − g C F π ε + O ( g ) , (28) Z DR, MS P = 1 − g C F π ε + O ( g ) , (29) Z DR, MS V = 1 + O ( g ) , (30) Z DR, MS A = 1 + O ( g ) , (31) Z DR, MS T = 1 + g C F π ε + O ( g ) , (32)where we use the notation { S,P,V,A,T } for { scalar, pseudoscalar, vector, axial-vector, tensor } operators. These factorsare in agreement with well-known results in the literature ([38] and references therein). Up to possible typos: An overall minus sign is missing in the pseudoscalar case. Furthermore, a sign must be altered in the definitionof the parameter ˆ ε , as follows: 1 / ˆ ε ≡ /ε + ln(4 π ) − γ E . As usual, perturbative corrections in Z DR , MSΓ may only be proportional to inverse powers of ε , with coefficients chosen in a way as togive a well-defined limit to the right-hand side of Eq. (27). C GIRS , MS S = 1 + g C F π (cid:2) µ ¯ z ) + 6 γ E − (cid:3) + O ( g ) , (33) C GIRS , MS P = 1 + g C F π (cid:2) µ ¯ z ) + 6 γ E − (cid:3) + O ( g ) , (34) C GIRS , MS V = 1 + g C F π (cid:18) (cid:19) + O ( g ) , (35) C GIRS , MS A = 1 + g C F π (cid:18)
32 + 4 hv (cid:19) + O ( g ) , (36) C GIRS , MS T = 1 + g C F π (cid:2) − ln(¯ µ ¯ z ) − γ E + 2 ln(2) (cid:3) + O ( g ) . (37)Note that the one-loop results of the bare GFs are proportional to the tree-level ones, and therefore only the length (notthe orientation) of x − y is relevant as a renormalization scale. Integrating Eqs. (22 - 26)) over spatial componentsand 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 X = Y = { , γ , γ t , γ i , γ γ t , γ γ i , γ t γ i , γ i γ j } . Two of them ( V t , A t ) give a vanishing contribution, both at tree level, andone loop; however, it is natural to impose that V t ( A t ) has the same renormalization factor as V i ( A i ). Thus, below wepresent results for the remaining 6 operators. Note that for extracting the correct one-loop renormalization factorsin GIRS, it is essential to include O ( ε ) terms of the tree-level GFs G tree ( x − y ); we recall that such terms are alsonecessary in the evaluation of the MS-renormalized GFs. The conversion factors are: C t − GIRS , MS S = 1 + g C F π (cid:18) −
12 + 6 ln(¯ µ ¯ t ) + 6 γ E (cid:19) + O ( g ) , (38) C t − GIRS , MS P = 1 + g C F π (cid:18) −
12 + 6 ln(¯ µ ¯ t ) + 6 γ E + 8 hv (cid:19) + O ( g ) , (39) C t − GIRS , MS V i = 1 + g C F π (cid:18) (cid:19) + O ( g ) , (40) C t − GIRS , MS A i = 1 + g C F π (cid:18)
32 + 4 hv (cid:19) + O ( g ) , (41) C t − GIRS , MS T ti = 1 + g C F π (cid:18) − t ¯ µ ) − γ E (cid:19) + O ( g ) , (42) C t − GIRS , MS T ij = 1 + g C F π (cid:18) − t ¯ µ ) − γ E (cid:19) + O ( g ) . (43)
3. RENORMALIZATION AND MIXING OF THE QUARK AND GLUON ENERGY-MOMENTUMTENSOR OPERATORS IN GIRS3.1. Details of the Calculation
The gauge-invariant part of the QCD traceless symmetric energy-momentum tensor (EMT) [39] contains two flavor-singlet operators, a gluonic O µν and a fermionic O µν : O µν = 2Tr[ F µρ F νρ ] − d δ µν F ρσ F ρσ ] , (44) O µν = N f X f =1 (cid:20) (cid:16) ¯ ψ f γ µ ←→ D ν ψ f + ¯ ψ f γ ν ←→ D µ ψ f (cid:17) − d δ µν (cid:16) ¯ ψ f γ ρ ←→ D ρ ψ f (cid:17)(cid:21) , (45)where F µν ≡ ∂ µ A ν − ∂ ν A µ + ig [ A µ , A ν ], ←→ D µ ≡ ( −→ D µ − ←− D µ ) / −→ D µ ≡ −→ ∂ µ + igA µ , and ←− D µ ≡ ←− ∂ µ − igA µ ; a summationover repeated indices is implied. These two operators are involved in the structure functions in nucleons [25, 40–42]:1the gluon operator appears in the leading-twist approximation of the gluon parton distribution function, while thefermion operator is related to the unpolarized quark parton distribution function. Furthermore, their matrix elementsare directly related to the gluon and quark average momentum fraction of a nucleon state [25, 43]. Also, theseoperators 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 matrixelements, as extracted from numerical simulations, to physical observables. A difficulty in calculating these renormal-ization factors is that mixing is present; these operators mix among themselves as they have the same transformationsunder Euclidean rotational (or hypercubic, on the lattice) symmetry. They also mix with other operators, includinggauge-variant operators (BRST variations and operators which vanish by the equations of motion; see [45]). Thelatter include ghost and gauge-fixing terms, which are well-defined in perturbation theory, while their nonpertur-bative extensions are not obvious; thus, a nonperturbative study of such terms by compact lattice simulations isproblematic. However, implementing a gauge-invariant renormalization scheme, such as GIRS, which involves onlygauge-invariant GFs, the gauge-variant operators do not contribute in the renormalization process and thus, theyare excluded. In our study, we consider only the nondiagonal elements ( µ = ν ) of the above operators, which give areduced set of operators under mixing on the lattice, containing only O µν and O µν . Thus, the renormalization of thenondiagonal components of EMT operators entails the construction of a 2 × O µν R O µν R = Z B,R Z B,R Z B,R Z B,R O µν B O µν B . (46)The calculation of the 2 × O µν and O µν . Threedifferent two-point GFs can be constructed by taking vacuum expectation values between the two mixing operators: G ν ν ; ν ν ij ( x − y ) ≡ (cid:10) O ν ν i ( x ) O ν ν j ( y ) (cid:11) , ( i, j ) = (1 , , (1 , , (2 , , (47)where ν = ν and ν = ν . By rotational (or just hypercubic) invariance: G ν ν ; ν ν ij ( x − y ) = G ν ν ; ν ν ji ( x − y ) = G ν ν ; ν ν ji ( x − y ). As it turns out, the two-point GFs G ij at one-loop level are proportional to the tree-level values ofthe 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 aboveGFs up to one loop in dimensional regularization. The tree-level contributions come from Feynman diagrams shownin Fig. 5, while the one-loop contributions stem from Feynman diagrams of Fig. 6. Note that G tree12 = 0. Figure 5: Diagrams which contribute to the tree-level Green’s functions G ν ν ; ν ν and G ν ν ; ν ν . A diamond (square) denotesinsertion of O µν ( O µν ). From the above GFs we can get three renormalization conditions by requiring absence of poles (for MS) or equalityto the corresponding tree-level values in a renormalization position scale (for GIRS). These are not enough for fullydetermining 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-dimensionaloperators. The procedure, which is also alluded to in [16], has the form of a bootstrap: One starts by renormalizinglowest dimensional operators (where no mixing issues are present), and then proceeds to renormalize operators ofincreasingly higher dimensionality by requiring finiteness (or some other normalization condition) in the GFs involvingproducts of operators up to that dimensionality. In this case, the only available lower-dimensional gauge-invariantlocal 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: hO µνi ( x ) O Γ ( y ) i . (48)However, such a GF vanishes for any choice of Γ to all perturbative orders: the corresponding two-point GFs with ascalar, pseudoscalar or tensor operator give traces of an odd number of Dirac matrices (for massless fermions), while2 G G G Figure 6: Diagrams which contribute to the one-loop Green’s functions G ν ν ; ν ν , G ν ν ; ν ν and G ν ν ; ν ν = G ν ν ; ν ν .Wavy (solid, dashed) lines represent gluons (quarks, ghosts). A diamond (square) denotes insertion of O µν ( O µν ). the GF with an axial vector gives traces containing one γ and Dirac matrices with symmetrized indices; even the GFwith a vector operator hO ν ν i ( x ) O γ ν ( y ) i vanishes, since O ν ν i is C-even, while O γ ν is C-odd.The next most “economic” possibility is to consider a three-point function constructed from the product of O µνi with two fermion bilinear operators ( O X , O Y ) at three distinct spacetime points: G µνi ; XY ( x − w, y − w ) ≡ hO µνi ( w ) O X ( x ) O Y ( y ) i , (49)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 complexityof the perturbative calculation. This also means that the renormalization factors defined in GIRS may depend on tworenormalization 4-vector scales. A priori , a possible way of addressing these issues could be to adopt a zero-momentuminsertion for one of the three operators. To do this one needs to perform a 4-dimensional integration over the positionvector x , y , or w depending on which operator carries zero momentum. Then, the resulting GF will depend on onlyone vector. However, such an integration over the whole spacetime causes additional complications: contact termsarise when any two position vectors among { x, y, w } , coincide, giving additional UV-divergences. To eliminate suchdivergences, 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 contactterm is to choose the vector ( y − w ) to be parallel (or antiparallel) to ( x − w ) (but x = y ). In this way, G µνi ; XY willdepend 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: G µνi ; XY ( x ) ≡ hO µνi (0) O X ( x ) O Y ( − x ) i . (50)The tree-level and one-loop Feynman diagrams contributing to G µνi ; XY , ( i = 1 , G µν XY ) tree = 0. A method for calculating the d-dimensional integrals stemming from theseFeynman diagrams is described in appendix A.3 Figure 7: Diagram which contributes to the tree-level Green’s functions G µν XY . A square denotes insertion of O µν . A crossdenotes insertion of a fermion bilinear operator. A diagram having the arrows of the fermion lines in counterclockwise directionmust also be considered.Figure 8: Diagrams which contribute to the one-loop level Green’s functions G µν XY . A diamond denotes insertion of O µν . Across denotes insertion of a fermion bilinear operator.Figure 9: Diagrams which contribute to the one-loop level Green’s functions G µν XY . A square denotes insertion of O µν . Across denotes insertion of a fermion bilinear operator. Diagrams having the arrows of the fermion lines in counterclockwisedirection must also be considered.
4A most natural set of four conditions for calculating the mixing matrix in GIRS, involving the above GFs, is:1 . hO ν ν ( x ) O ν ν ( y ) i| x − y =¯ z = hO ν ν ( x ) O ν ν ( y ) i tree | x − y =¯ z , (51)2 . hO ν ν ( x ) O ν ν ( y ) i| x − y =¯ z = hO ν ν ( x ) O ν ν ( y ) i tree | x − y =¯ z , (52)3 . hO ν ν ( x ) O ν ν ( y ) i| x − y =¯ z = hO ν ν ( x ) O ν ν ( y ) i tree | x − y =¯ z = 0 , (53)4 . hO ν ν (0) O GIRS X ( x ) O GIRS Y ( − x ) i| x =¯ z = hO ν ν (0) O GIRS X ( x ) O GIRS Y ( − x ) i tree | x =¯ z = 0 , (54)where ν = ν and ν = ν . Alternatively, we can replace the second condition (Eq. (52)) with: hO ν ν (0) O GIRS X ( x ) O GIRS Y ( − x ) i| x =¯ z = hO ν ν (0) O GIRS X ( x ) O GIRS Y ( − x ) i tree | x =¯ z . (55)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 vanishingcontributions, 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 : Z G + 2 Z Z G + Z G = G tree11 , (56) Z G + 2 Z Z G + Z G = G tree22 , (57) Z Z G + ( Z Z + Z Z ) G + Z Z G = G tree12 = 0 , (58) Z X Z Y ( Z G XY + Z G XY ) = G tree1; XY = 0 . (59)The four elements of the mixing matrix Z ij can be obtained by solving the above system of 4 equations, once all GFson 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:1 . Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i = Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i tree , (60)2 . Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i = Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i tree , (61)3 . Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i = Z d ~x hO ν ν − GIRS ( ~x, O ν ν − GIRS ( ~ , ¯ t ) i tree = 0 , (62)4 . Z d ~x hO ν ν − GIRS ( ~ , O t − GIRS X ( ~x, ¯ t ) O t − GIRS Y ( − ~x, − ¯ t ) i = Z d ~x hO ν ν − GIRS ( ~ , O t − GIRS X ( ~x, ¯ t ) O t − GIRS Y ( − ~x, − ¯ t ) i tree = 0 . (63)No summation over ν , ν is implied. Depending on the choice of X and Y , the fourth condition of t-GIRS may involveodd integrals, which give zero. For example, using X = Y = δ ν ν , which vanishessince we study the nondiagonal elements ( ν = ν ) of O ν ν i , and (ii) x ν x ν , which will vanish upon integration over ~x for any choices of ν , ν . In such cases, an appropriate variant of the fourth condition can be applied; e.g., for thecase X = Y =
1, a possible alternative condition, in place of Eq. (63), is: Z d ~x x ν x ν x hO ν ν − GIRS ( ~ , O t − GIRS ~x, ¯ t ) O t − GIRS − ~x, − ¯ t ) i = Z d ~x x ν x ν x hO ν ν − GIRS ( ~ , O t − GIRS ~x, ¯ t ) O t − GIRS − ~x, − ¯ t ) i tree = 0 , (64)(no summation over ν , ν is implied), or Z d ~x e i~p · ~x hO ν ν − GIRS ( ~ , O t − GIRS ~x, ¯ t ) O t − GIRS − ~x, − ¯ t ) i = Z d ~x e i~p · ~x hO ν ν − GIRS ( ~ , O t − GIRS ~x, ¯ t ) O t − GIRS − ~x, − ¯ t ) i tree = 0 , (65) For simplicity, we omit superscripts referring to the regularization and renormalization scheme, as well as Lorentz indices. We also omitthe dependence on spacetime coordinates. ~p . These variant schemes lead to relations analogous to Eqs. (56 - 59) for thedetermination of Z ij . In this work, we present one-loop results for three-point GFs with ( X , Y ) = ( γ , γ ),( γ ν , γ ν ), ( γ γ ν , γ γ ν ), before performing integration over ~x . It is straightforward to integrate all these resultsover ~x ; we do so for some specific cases (see next section), which are likely the most appropriate for nonperturbativeinvestigations.After calculating the mixing matrix, we extract the conversion factors between the GIRS (or t-GIRS) and the MSscheme for the EMT operators; they have a 2 × O µν O µν = C GIRS , MS11 C GIRS , MS12 C GIRS , MS21 C GIRS , MS22 · O µν O µν ⇒ (66) C GIRS , MS11 C GIRS , MS12 C GIRS , MS21 C GIRS , MS22 = Z B, MS11 Z B, MS12 Z B, MS21 Z B, MS22 · Z B, GIRS11 Z B, GIRS12 Z B, GIRS21 Z B, GIRS22 − . (67)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. In this subsection, we present our results (up to one loop) for the bare GFs hO ν ν i ( x ) O ν ν j ( y ) i , and hO ν ν i (0) O X ( x ) O Y ( − x ) i for i, j = 1 , X, Y ) = ( γ , γ ), ( γ ν , γ ν ), ( γ γ ν , γ γ ν ), as well as theconversion factors between all variants of GIRS and MS. The results are expressed in terms of the following Lorentzstructures : s [1] ν ν ( x ) ≡ x ν x ν x , (68) s [2] ν ν ν ν ( x ) ≡ x ν x ν x δ ν ν , (69) s [3] ν ν ν ν ≡ δ ν ν δ ν ν + δ ν ν δ ν ν , (70) s [4] ν ν ν ν ( x ) ≡ δ ν ν x ν x ν x + δ ν ν x ν x ν x + δ ν ν x ν x ν x + δ ν ν x ν x ν x , (71) s [5] ν ν ν ν ( x ) ≡ x ν x ν x ν x ν ( x ) . (72)The resulting expressions for the bare GFs are ( z ≡ y − x ): hO ν ν ( x ) O ν ν ( y ) i = (cid:16) s [3] ν ν ν ν − s [4] ν ν ν ν ( z ) + 8 s [5] ν ν ν ν ( z ) (cid:17) N c C F (1 − ε ) (Γ(2 − ε )) π − ε ( z ) − ε × " − g π (cid:16) N f (cid:16) ε + 2 γ E − µ z ) + 16 (cid:17) + 20 N c (cid:17) + O ( ε , g ) + O ( ε , g ) + O ( g ) , (73) hO ν ν ( x ) O ν ν ( y ) i = (cid:16) s [3] ν ν ν ν − s [4] ν ν ν ν ( z ) + 8 s [5] ν ν ν ν ( z ) (cid:17) N c N f (2 − ε ) (Γ(2 − ε )) (1 − ε )2 π − ε ( z ) − ε × g π C F (cid:18) ε + 2 γ E − µ z ) − (cid:19) + O ( ε , g ) + O ( g ) , (74) For brevity, decimal numbers in our results are presented only with six digits after the decimal point; they are known to higher accuracy. hO ν ν ( x ) O ν ν ( y ) i = (cid:16) s [3] ν ν ν ν − s [4] ν ν ν ν ( z ) + 8 s [5] ν ν ν ν ( z ) (cid:17) N c N f (2 − ε ) (Γ(2 − ε )) π − ε ( z ) − ε × (75) " − g π (cid:18) C F (cid:19) (cid:16) ε + 2 γ E − µ z ) − (cid:17) + O ( ε , g ) + O ( ε , g ) + O ( g ) , (76) hO ν ν (0) O x ) O − x ) i = − s [1] ν ν ( x ) N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε g π C F (cid:20) ε + ln(¯ µ x ) − . (cid:21) + O ( ε , g ) + O ( g ) , (77) hO ν ν (0) O γ ( x ) O γ ( − x ) i = −hO ν ν (0) O x ) O − x ) i , (78) hO ν ν (0) O γ ν ( x ) O γ ν ( − x ) i = N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε g π C F × "(cid:16) s [4] ν ν ν ν ( x ) + 2 s [2] ν ν ν ν ( x ) − s [5] ν ν ν ν ( x ) (cid:17) (cid:16) ε + ln(¯ µ x ) − . (cid:17) + 12 s [2] ν ν ν ν ( x ) + 34 s [3] ν ν ν ν − s [4] ν ν ν ν ( x ) + O ( ε , g ) + O ( g ) , (79) hO ν ν (0) O γ γ ν ( x ) O γ γ ν ( − x ) i = hO ν ν (0) O γ ν ( x ) O γ ν ( − x ) i , (80) hO ν ν (0) O x ) O − x ) i = − s [1] ν ν ( x ) N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε × " g π C F (cid:18) ε + ln(¯ µ x ) + 2 . (cid:19) + O ( ε , g ) + O ( g ) , (81) hO ν ν (0) O γ ( x ) O γ ( − x ) i = −hO ν ν (0) O x ) O − x ) i + hv s [1] ν ν ( x ) N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε g π C F + O ( g ) , (82) hO ν ν (0) O γ ν ( x ) O γ ν ( − x ) i = N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε × "(cid:16) s [4] ν ν ν ν ( x ) + 2 s [2] ν ν ν ν ( x ) − s [5] ν ν ν ν ( x ) (cid:17) × − g π C F (cid:16) ε + ln(¯ µ x ) − . (cid:17)! − g π C F (cid:16) s [2] ν ν ν ν ( x ) + 98 s [3] ν ν ν ν − s [4] ν ν ν ν ( x ) ! + O ( ε , g ) + O ( g ) , (83) hO ν ν (0) O γ γ ν ( x ) O γ γ ν ( − x ) i = hO ν ν (0) O γ ν ( x ) O γ ν ( − x ) i + hv N c N f (Γ(2 − ε )) Γ(3 − ε )2 − ε π − ε ( x ) − ε × (cid:16) s [4] ν ν ν ν ( x ) + 2 s [2] ν ν ν ν ( x ) − s [5] ν ν ν ν ( x ) (cid:17) g π C F + O ( g ) . (84)7The above GFs lead to the following results for Z DR , MS ij : Z DR , MS11 = 1 + g π ε (cid:18) N f (cid:19) + O ( g ) , (85) Z DR , MS12 = g π ε (cid:18) − C F (cid:19) + O ( g ) , (86) Z DR , MS21 = g π ε (cid:18) − N f (cid:19) + O ( g ) , (87) Z DR , MS22 = 1 + g π ε (cid:18) C F (cid:19) + O ( g ) . (88)These are consistent with the results found using GFs with elementary external fields [46]. The corresponding MS-renormalized GFs can be obtained by removing 1 /ε terms from Eqs. (73) - (84) and by taking the na¨ıve limit ε → : X = Y =
1; no integration over ~x .Similarly, by choosing X = Y = γ , 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 hO ν ν (0) O γ ( x ) O γ ( − x ) i ,which does not contribute in the calculation of the conversion factors to one loop. Nevertheless, numerical datawill be much different; thus, this provides for an interesting comparison of the corresponding MS-renormalizedGFs, as gotten from the lattice.2. GIRS : X = γ ν , Y = γ ν ; no integration over ~x ; ν , ν , ν , ν are all different.Since the three-point functions hO ν ν i (0) O γ ν ( x ) O γ ν ( − x ) i have more than one Lorentz structures, it is neces-sary to isolate a structure by using projectors or by making specific choices for the indices ν – ν and/or thecomponents of x . In this variant, we isolate the structure s [5] ν ν ν ν ( x ) by choosing ν , ν , ν , ν to be all different.In this case all four components of x must be nonzero. Similarly, the choice of X = γ γ ν , Y = γ γ ν will givethe same conversion factors to one loop.3. GIRS : X = γ ν , Y = γ ν ; no integration over ~x ; ν = ν , ν = ( ν , ν ); and x ν = 0.This variant is similar to GIRS with the difference of isolating the structure s [2] ν ν ν ν ( x ).4. t − GIRS : X = γ ν , Y = γ ν ; integration over ~x ; ν = ν , ν = ν ; ν , ν are both spatial.The integration over the spatial components of x will give zero unless the four indices are paired (both in thetwo-point and three-point GFs), e.g., ν = ν , ν = ν . In principle, there are two distinct possibilities: ν , ν are both spatial or ν is spatial and ν is temporal. However, the latter case will be impossible to satisfy, sincethe combination (cid:16) s [3] ν ν ν ν − s [4] ν ν ν ν ( z ) + 8 s [5] ν ν ν ν ( z ) (cid:17) / ( z ) appearing in the two-point functions [see Eqs.( 73 - 76)] vanishes upon integration over spatial components.5. t − GIRS : X = Y =
1; integration over ~x ; ν = ν , ν = ν ; ν , ν are both spatial; projector: x ν x ν /x (seeEq. (64)).8The conversion factors for the above variants of GIRS are given below: C GIRS i , MS11 = 1 − g π h N c + c N f + 23 N f ln(¯ µ ¯ z ) i + O ( g ) , (89) C GIRS i , MS12 = − g π C F h c −
83 ln(¯ µ ¯ z ) i + O ( g ) , (90) C GIRS i , MS21 = − g π N f h c −
23 ln(¯ µ ¯ z ) i + O ( g ) , (91) C GIRS i , MS22 = 1 − g π C F h c + 83 ln(¯ µ ¯ z ) i + O ( g ) , (92) C t − GIRS j , MS11 = 1 − g π h N c + c N f + 23 N f ln(¯ µ ¯ t ) i + O ( g ) , (93) C t − GIRS j , MS12 = − g π C F h c −
83 ln(¯ µ ¯ t ) i + O ( g ) , (94) C t − GIRS j , MS21 = − g π N f h c −
23 ln(¯ µ ¯ t ) i + O ( g ) , (95) C t − GIRS j , MS22 = 1 − g π C F h c + 83 ln(¯ µ ¯ t ) i + O ( g ) , (96)where i = 1 , , j = 1 , c kl are given in Table III for each variant of GIRS. GIRS GIRS GIRS t − GIRS t − GIRS c -0.043464 -0.043464 -0.043464 0.236288 0.236288 c c c -3.896079 -3.896079 -3.896079 -2.777072 -2.777072Table III: Values of the one-loop coefficient c ij in the definition of the conversion factors of EMT operators, given in Eqs. ( 89- 96). Use of Eq. (65) as an alternative renormalization condition requires the integration of various expressions of theform: e i~p · ~x x µ x µ . . . x µ j (ln( x )) i ( x ) k , ( i, j, k : nonnegative integers) (97)over spatial components of x = ( ~x, t ). All these integrals can be performed by using the following generating integralfunction: Z d x e i~p · ~x ( x ) a = 2 / − a π / Γ( a ) ( | p | / | t | ) − / a K / − a ( | p || t | ) , (98)where K ν ( z ) is the modified Bessel function of the second kind and a may take noninteger values. Then, differentiatingwith respect to a or to individual components of p , we can calculate all necessary integrals arising in t-GIRS.9
4. SUMMARY
In this paper, we study a gauge-invariant, mass-independent renormalization scheme (GIRS) for composite op-erators, which is applicable in both perturbative and nonperturbative studies. This is an extended version of thecoordinate space (X-space) renormalization scheme studied in, e.g., Refs. [3, 4]. This scheme involves vacuum expec-tation values of products of gauge-invariant operators located at different spacetime points. The expectation valuesare gauge-independent and thus, gauge fixing is not needed in this scheme. Also, gauge-variant operators, whichmay mix with gauge-invariant operators, do not contribute in such Green’s functions; as a consequence, they can besafely excluded, leading to a reduced set of mixing operators. In this work, we apply GIRS in the renormalizationof fermion bilinear operators, as well as in the renormalization and mixing of the gluon and quark parts of the QCDenergy-momentum tensor (EMT). We propose different variants of GIRS, e.g., using specific values for the positionvectors of the operators under study, or integrating over timeslices (t-GIRS), which may lead to reduced statisticalnoise in the nonperturbative calculations via lattice simulations. We provide results, up to one loop, for the conversionfactors 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 ofsimilar kind, e.g., four-fermi operators and supersymmetric operators (Gluino-Glue, Noether supercurrent).
Acknowledgments
M.C., H.P. and A.S. acknowledge financial support from the project “Quantum Fields on the Lattice” , fundedby the Cyprus Research and Innovation Foundation (RIF) under contract number EXCELLENCE/0918/0066. G.S.acknowledges financial support by the University of Cyprus, under the research programs entitled “
Quantum Fieldson the Lattice ” and “
Nucleon parton distribution functions using Lattice Quantum Chromodynamics ”. We thank C.Alexandrou, M. Dalla Brida, and K. Hadjiyiannakou for useful comments.
Appendix A: Technical aspects of the calculation
There are three types of scalar Feynman integrals appearing in our calculation: I ( ξ ; α ) ≡ Z d d p (2 π ) d e ip · ξ ( p ) α , (A1) I ( ξ , ξ ; α , α , α ) ≡ Z d d p d d p (2 π ) d e ip · ξ e ip · ξ ( p ) α ( p ) α (( − p + p ) ) α , (A2) I ( ξ , ξ ; α , α , α , α , α ) ≡ Z d d p d d p d d p (2 π ) d e ip · ξ e ip · ξ ( p ) α ( p ) α (( − p + p ) ) α ( p ) α (( − p + p ) ) α . (A3)For simplicity, we write down only scalar integrals for each type; integrands containing additional factors of p µ , p ν can be handled in a similar way, or by taking derivatives of the results with respect to ξ µ , ξ ν , respectively. Below,we briefly describe the procedure for calculating each type of integral:1. Integral I : We introduce Schwinger parameters: 1( p ) α = 1Γ( α ) Z ∞ dλ λ α − e − λp . (A4)After integrating over p and λ , we get: I ( ξ ; α ) = Γ( − α + d/
2) ( ξ ) α − d/ α π d/ Γ( α ) . (A5)02. Integral I : We introduce Schwinger parameters:1( p ) α ( p ) α (( − p + p ) ) α = 1Γ( α ) Γ( α ) Γ( α ) × Z ∞ dλ Z ∞ dλ Z ∞ dλ λ α − λ α − λ α − e − λ p − λ p − λ ( − p + p ) . (A6)After integrating over p and p , we make a change of variables: x = λ / ( λ + λ ), x = 1 − λ λ / ( λ λ + λ λ + λ λ ) and x = λ λ λ / ( λ λ + λ λ + λ λ ). Integrals over x and x can be calculated algebraically,while the remaining integration over x cannot be obtained in a closed form (for general values of α i ). Theresulting expression takes the following form: I ( ξ , ξ ; α , α , α ) = Γ( − α + d/
2) Γ( − α − α − α + d )4 α + α + α π d Γ( α ) Γ( α ) Γ( d/
2) ( ξ ) α + α + α − d × Z dx " (1 − x ) − α + α − d/ x − α + α − d/ × F − α + d/ , − α − α − α + d, d/ − ( ξ + x ξ ) (1 − x ) x ξ ) ! , (A7)where ( α + α + α − d ) <
0, ( − α + d/ > α >
0. The next step is to examine whether the integrationover x and the limit of vanishing regulator ( ε → ε = 2 − d/
2) can be safely interchanged without leadingto divergences. For this check, it is useful to express the hypergeometric function appearing in (A7) as a powerseries in x and (1 − x ), by applying an appropriate transformation formula (see [47]). In case the interchangeis indeed permissible, the integration over x 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 directintegration over x . The Laurent expansion of the hypergeometric function F over ε = 0 has been performedwith the help of the mathematica package “HypExp” introduced in Ref. [48].3. Integral I : We introduce: 1 = 1 d X ρ ∂p ρ ∂p ρ . (A8)After integrating by parts, we get the following recursive relation, which can eliminate inverse powers of p , or p , or p : I ( ξ , ξ ; α , α , α , α , α ) = 1 − α − α − α + d · h α (cid:16) I ( ξ , ξ ; α − , α , α + 1 , α , α ) − I ( ξ , ξ ; α , α − , α + 1 , α , α ) (cid:17) + α (cid:16) I ( ξ , ξ ; α − , α , α , α , α + 1) − I ( ξ , ξ ; α , α , α , α − , α + 1) (cid:17)i . (A9)In the case where α , α , α are positive integers, which is true in the computation at hand, an iterativeimplementation of Eq. (A9) leads to terms with one propagator less. One momentum can then be integratedusing a well-known one-loop formula (see Eqs. (A.1 – A.2) in Ref. [49]); the remaining integrals are of type 1or 2. [1] M. Bochicchio, L. Maiani, G. Martinelli, G. C. Rossi, M. Testa, Chiral Symmetry on the Lattice with Wilson Fermions,Nucl. Phys. B262 (1985) 331. doi:10.1016/0550-3213(85)90290-1 . [2] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa, A. Vladikas, A General method for nonperturbative renormalizationof lattice operators, Nucl. Phys. B 445 (1995) 81–108. arXiv:hep-lat/9411010 , doi:10.1016/0550-3213(95)00126-D .[3] V. Gimenez, L. Giusti, S. Guerriero, V. Lubicz, G. Martinelli, S. Petrarca, J. Reyes, B. Taglienti, E. Trevi-gne, Non-perturbative renormalization of lattice operators in coordinate space, Phys. Lett. B598 (2004) 227–236. arXiv:hep-lat/0406019 , doi:10.1016/j.physletb.2004.07.053 .[4] K. Jansen, C. Liu, M. Luscher, H. Simma, S. Sint, R. Sommer, P. Weisz, U. Wolff, Nonperturbative renormalization of lat-tice QCD at all scales, Phys. Lett. B 372 (1996) 275–282. arXiv:hep-lat/9512009 , doi:10.1016/0370-2693(96)00075-5 .[5] P. Dimopoulos, G. Herdo´ıza, M. Papinutto, C. Pena, D. Preti, A. Vladikas, Non-Perturbative Renormalisation andRunning of BSM Four-Quark Operators in N f = 2 QCD, Eur. Phys. J. C 78 (7) (2018) 579. arXiv:1801.09455 , doi:10.1140/epjc/s10052-018-6002-y .[6] S. D. Joglekar, B. W. Lee, General Theory of Renormalization of Gauge Invariant Operators, Annals Phys. 97 (1976) 160. doi:10.1016/0003-4916(76)90225-6 .[7] A. Maas, Describing gauge bosons at zero and finite temperature, Other thesis (2013). arXiv:1106.3942 , doi:10.1016/j.physrep.2012.11.002 .[8] A. Cucchieri, D. Dudal, T. Mendes, O. Oliveira, M. Roelfs, P. J. Silva, Faddeev-Popov Matrix in Linear Covariant Gauge:First Results, Phys. Rev. D 98 (9) (2018) 091504. arXiv:1809.08224 , doi:10.1103/PhysRevD.98.091504 .[9] K. Chetyrkin, A. Retey, Renormalization and running of quark mass and field in the regularization invariantand MS-bar schemes at three loops and four loops, Nucl. Phys. B 583 (2000) 3–34. arXiv:hep-ph/9910332 , doi:10.1016/S0550-3213(00)00331-X .[10] B. Ruijl, T. Ueda, J. Vermaseren, A. Vogt, Four-loop QCD propagators and vertices with one vanishing external momentum,JHEP 06 (2017) 040. arXiv:1703.08532 , doi:10.1007/JHEP06(2017)040 .[11] T. Luthe, A. Maier, P. Marquard, Y. Schr¨oder, The five-loop Beta function for a general gauge group and anomalousdimensions beyond Feynman gauge, JHEP 10 (2017) 166. arXiv:1709.07718 , doi:10.1007/JHEP10(2017)166 .[12] K. Chetyrkin, G. Falcioni, F. Herzog, J. Vermaseren, Five-loop renormalisation of QCD in covariant gauges, JHEP 10(2017) 179, [Addendum: JHEP 12, 006 (2017)]. arXiv:1709.08541 , doi:10.1007/JHEP10(2017)179 .[13] J. A. Gracey, R. M. Simms, Renormalization of QCD in the interpolating momentum subtraction scheme at three loops,Phys. Rev. D 97 (8) (2018) 085016. arXiv:1801.10415 , doi:10.1103/PhysRevD.97.085016 .[14] F. Herzog, S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren, A. Vogt, Five-loop contributions to low-N non-singletanomalous dimensions in QCD, arXiv e-prints (2018) arXiv:1812.11818 arXiv:1812.11818 .[15] P. Baikov, K. Chetyrkin, Transcendental structure of multiloop massless correlators and anomalous dimensions, JHEP 10(2019) 190. arXiv:1908.03012 , doi:10.1007/JHEP10(2019)190 .[16] J. A. Gracey, Three loop anti-MS operator correlation functions for deep inelastic scattering in the chiral limit, JHEP 04(2009) 127. arXiv:0903.4623 , doi:10.1088/1126-6708/2009/04/127 .[17] K. G. Chetyrkin, A. Maier, Massless correlators of vector, scalar and tensor currents in position space at orders α s and α s :Explicit analytical results, Nucl. Phys. B844 (2011) 266–288. arXiv:1010.1145 , doi:10.1016/j.nuclphysb.2010.11.007 .[18] K. Cichy, K. Jansen, P. Korcyl, Non-perturbative renormalization in coordinate space for N f = 2 maximally twistedmass fermions with tree-level Symanzik improved gauge action, Nucl. Phys. B865 (2012) 268–290. arXiv:1207.0628 , doi:10.1016/j.nuclphysb.2012.08.006 .[19] M. Tomii, N. H. Christ, O (4)-symmetric position-space renormalization of lattice operators, Phys. Rev. D99 (1) (2019)014515. arXiv:1811.11238 , doi:10.1103/PhysRevD.99.014515 .[20] Y.-B. Yang, M. Gong, J. Liang, H.-W. Lin, K.-F. Liu, D. Pefkou, P. Shanahan, Nonperturbatively renormalized gluemomentum fraction at the physical pion mass from lattice QCD, Phys. Rev. D98 (7) (2018) 074506. arXiv:1805.00531 , doi:10.1103/PhysRevD.98.074506 .[21] P. E. Shanahan, W. Detmold, Gluon gravitational form factors of the nucleon and the pion from lattice QCD, Phys. Rev.D99 (1) (2019) 014511. arXiv:1810.04626 , doi:10.1103/PhysRevD.99.014511 .[22] Y.-B. Yang, J. Liang, Y.-J. Bi, Y. Chen, T. Draper, K.-F. Liu, Z. Liu, Proton Mass Decomposition from the QCD EnergyMomentum Tensor, Phys. Rev. Lett. 121 (21) (2018) 212001. arXiv:1808.08677 , doi:10.1103/PhysRevLett.121.212001 .[23] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, H. Panagopoulos,G. Spanoudes, Complete flavor decomposition of the spin and momentum fraction of the proton using lattice QCD simula-tions at physical pion mass, Phys. Rev. D101 (9) (2020) 094513. arXiv:2003.08486 , doi:10.1103/PhysRevD.101.094513 .[24] M. Dalla Brida, L. Giusti, M. Pepe, Non-perturbative definition of the QCD energy-momentum tensor on the lattice,JHEP 04 (2020) 043. arXiv:2002.06897 , doi:10.1007/JHEP04(2020)043 .[25] X.-D. Ji, Off forward parton distributions, J. Phys. G 24 (1998) 1181–1205. arXiv:hep-ph/9807358 , doi:10.1088/0954-3899/24/7/002 .[26] M. Constantinou, Hadron Structure, PoS LATTICE2014 (2015) 001. arXiv:1411.0078 , doi:10.22323/1.214.0001 .[27] S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B303 (1993) 113–118. arXiv:hep-ph/9302240 , doi:10.1016/0370-2693(93)90053-K .[28] M. S. Chanowitz, M. Furman, I. Hinchliffe, The Axial Current in Dimensional Regularization, Nucl. Phys. B159 (1979)225–243. doi:10.1016/0550-3213(79)90333-X .[29] G. ’t Hooft, M. J. G. Veltman, Regularization and Renormalization of Gauge Fields, Nucl. Phys. B44 (1972) 189–213. doi:10.1016/0550-3213(72)90279-9 .[30] W. Siegel, Supersymmetric Dimensional Regularization via Dimensional Reduction, Phys. Lett. 84B (1979) 193–196. doi:10.1016/0370-2693(79)90282-X .[31] A. Patel, S. R. Sharpe, Perturbative corrections for staggered fermion bilinears, Nucl. Phys. B395 (1993) 701–732. arXiv:hep-lat/9210039 , doi:10.1016/0550-3213(93)90054-S .[32] A. J. Buras, P. H. Weisz, QCD Nonleading Corrections to Weak Decays in Dimensional Regularization and ’t Hooft-VeltmanSchemes, Nucl. Phys. B333 (1990) 66–99. doi:10.1016/0550-3213(90)90223-Z .[33] B. All´es, A. Feo, H. Panagopoulos, Asymptotic scaling corrections in QCD with Wilson fermions from the three loop averageplaquette, Phys. Lett. B 426 (1998) 361–366, [Erratum: Phys.Lett.B 553, 337–338 (2003)]. arXiv:hep-lat/9801003 , doi:10.1016/S0370-2693(98)00295-0 .[34] H. Panagopoulos, A. Skouroupathis, A. Tsapalis, Free energy and plaquette expectation value for gluons on the lattice, inthree dimensions, Phys. Rev. D 73 (2006) 054511. arXiv:hep-lat/0601009 , doi:10.1103/PhysRevD.73.054511 .[35] A. Athenodorou, H. Panagopoulos, A. Tsapalis, The Lattice Free Energy of QCD with Clover Fermions, up to Three-Loops,Phys. Lett. B 659 (2008) 252–259. arXiv:0710.3856 , doi:10.1016/j.physletb.2007.11.064 .[36] M. Brambilla, F. Di Renzo, High-loop perturbative renormalization constants for Lattice QCD (II): three-loop quarkcurrents for tree-level Symanzik improved gauge action and n f =2 Wilson fermions, Eur. Phys. J. C 73 (12) (2013) 2666. arXiv:1310.4981 , doi:10.1140/epjc/s10052-013-2666-5 .[37] M. Constantinou, H. Panagopoulos, Perturbative renormalization of quasi-parton distribution functions, Phys. Rev. D96 (5)(2017) 054506. arXiv:1705.11193 , doi:10.1103/PhysRevD.96.054506 .[38] J. A. Gracey, Three loop anomalous dimension of nonsinglet quark currents in the RI-prime scheme, Nucl. Phys. B662(2003) 247–278. arXiv:hep-ph/0304113 , doi:10.1016/S0550-3213(03)00335-3 .[39] D. Z. Freedman, I. J. Muzinich, E. J. Weinberg, On the Energy-Momentum Tensor in Gauge Field Theories, Annals Phys.87 (1974) 95. doi:10.1016/0003-4916(74)90448-5 .[40] X.-D. Ji, Breakup of hadron masses and energy - momentum tensor of QCD, Phys. Rev. D 52 (1995) 271–281. arXiv:hep-ph/9502213 , doi:10.1103/PhysRevD.52.271 .[41] X.-D. Ji, Gauge-Invariant Decomposition of Nucleon Spin, Phys. Rev. Lett. 78 (1997) 610–613. arXiv:hep-ph/9603249 , doi:10.1103/PhysRevLett.78.610 .[42] A. Radyushkin, Nonforward parton distributions, Phys. Rev. D 56 (1997) 5524–5557. arXiv:hep-ph/9704207 , doi:10.1103/PhysRevD.56.5524 .[43] R. Horsley, R. Millo, Y. Nakamura, H. Perlt, D. Pleiter, P. Rakow, G. Schierholz, A. Schiller, F. Winter,J. Zanotti, A Lattice Study of the Glue in the Nucleon, Phys. Lett. B 714 (2012) 312–316. arXiv:1205.6410 , doi:10.1016/j.physletb.2012.07.004 .[44] T. Aoyama, et al., The anomalous magnetic moment of the muon in the Standard Model (6 2020). arXiv:2006.04822 .[45] S. Caracciolo, P. Menotti, A. Pelissetto, One loop analytic computation of the energy momentum tensor for lattice gaugetheories, Nucl. Phys. B 375 (1992) 195–239. doi:10.1016/0550-3213(92)90339-D .[46] G. Panagopoulos, H. Panagopoulos, G. Spanoudes, Two-loop renormalization and mixing of gluon and quark energy-momentum tensor operators, Phys. Rev. D 103 (2021) 014515. arXiv:2010.02062 , doi:10.1103/PhysRevD.103.014515 .[47] I. Gradshteyn, I. M. Ryzhik, Tables of integrals, series, and products, Elsevier, 2007, Equation 9.132.1.[48] T. Huber, D. Maitre, HypExp: A Mathematica package for expanding hypergeometric functions around integer-valuedparameters, Comput. Phys. Commun. 175 (2006) 122–144. arXiv:hep-ph/0507094 , doi:10.1016/j.cpc.2006.01.007 .[49] K. Chetyrkin, F. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B192 (1981) 159–204. doi:10.1016/0550-3213(81)90199-1doi:10.1016/0550-3213(81)90199-1