Three-loop formula for quark and gluon contributions to the QCD trace anomaly

In the QCD energy-momentum tensor $T^{\mu\nu}$, the terms that contribute to physical matrix elements are expressed as the sum of the gauge-invariant quark part and gluon part. Each part undergoes the renormalization due to the interactions among quarks and gluons, although the total tensor $T^{\mu\nu}$ is not renormalized thanks to the conservation of energy and momentum. Recently it has been shown that, through the renormalization, each of the quark and gluon parts of $T^{\mu\nu}$ receives a definite amount of anomalous trace contribution, such that their sum reproduces the well-known QCD trace anomaly, $T^\mu_\mu= (\beta/2g)F^{\mu\nu}F_{\mu\nu}+ m (1+\gamma_m)\bar{\psi}\psi$, and the corresponding formulas have been derived up to two-loop order. We extend this result to the three-loop order, working out all the relevant three-loop renormalization structure for the quark and gluon energy-momentum tensors in the (modified) minimal subtraction scheme in the dimensional regularization. We apply our three-loop formula of the quark/gluon decomposition of the trace anomaly to calculate the anomaly-induced mass structure of nucleons as well as pions.


I. INTRODUCTION
The QCD energy-momentum tensor T µν is known to receive the trace anomaly as [1][2][3], representing the broken scale invariance due to the quantum loop effects, with the betafunction β for the QCD coupling constant g and the anomalous dimension γ m for the quark mass m. Here, the suffix R denotes the renormalized quantities, and (F µν F µν ) R and ψ ψ R denote the renormalized composite operators. This QCD trace anomaly signals the generation of a nonperturbative mass scale, say, the nucleon mass M, taking the nucleon matrix element of (1) and using the fact that P |T µν |P = 2P µ P ν , where |P is the nucleon state with 4-momentum P µ .
The energy-momentum tensor consists of the quark part and the gluon part, T µν = T µν q + T µν g . Considering the corresponding renormalized quantities as (T µν ) R = T µν q R + T µν g R , we have representing conservation of energy and momentum, while T µν q R and T µν g R are not conserved separately and thus are different from T µν q and T µν g , respectively, by the ultraviolet (UV) subtraction terms to make the divergent bare quantities finite. Such subtraction terms for the second rank symmetric tensors in various quantum field theories appear to produce a finite contribution, when contracted with the metric tensor η µν [4]. The explicit results for the individual quark and gluon parts in QCD have been derived at the two-loop level in the minimal subtraction (MS) scheme in the dimensional regularization in our recent paper [5], and the corresponding one-loop results read, with such that their total sum, η µν T µν q R + η µν T µν g R , reproduces the QCD trace anomaly (1) with the one-loop terms for the renormalization group (RG) functions, being substituted, where C F = N 2 c −1 2Nc = 4 3 and C A = N c = 3 for N c (= 3) color. Therefore, the quark-loop (∝ 2n f 3 ) and gluon-loop (∝ − 11 3 C A ) contributions to the gluon self-energy (in the background field method, see, e.g., [6]) directly determine the (F 2 ) R terms in (3) and (4), respectively. Similar correspondence may be invoked to the relative weights, 4 and 14, of the (mψψ) R terms in (3) and (4) by interpreting that those weights correspond to the respective contributions caused by the quark and gluon internal propagators in the quark self-energy. Such simple correspondence between the loop contributions of the diagrams and the decomposition of the total anomaly demonstrates that the decomposition as (3) and (4) is of physical origin and well-defined; this point is also supported by the fact [5] that actually (3) and (4) are separately RG-invariant to one-loop order, because (mψψ) R is an exactly RGinvariant quantity and α s (F 2 ) R is RG-invariant up to the corrections of order α 2 s . However, similar intuitive correspondence is not obvious in the two-loop results derived in [5], and each of quark/gluon parts exhibits the RG scale dependence. Therefore, the decomposition at two loops is likely to depend on regularization and renormalization schemes to handle the UV divergences. Within the MS-like (MS, MS) schemes in the dimensional regularization, their mutual relation is straightforward, but the results using the other schemes are not known at present.
Physical relevance as well as phenomenological implications of such decomposition is also demonstrated in [5], such that those RG properties of the quark/gluon parts of the trace anomaly allow us to constrain the twist-four gravitational form factorC q,g , which arises as one of the gravitational form factors [7][8][9][10][11] to parametrize hadron matirx elements of the quark/gluon parts of the QCD energy-momentum tensor, P ′ |T µν q,g |P . In particular, the solution of the corresponding two-loop RG equations provides the model-independent determination of the forward (P ′ → P ) value ofC q,g , at the level of accuracy ∼ 10%. Such quantitative constraint could have an impact on the descriptions of the shape deep inside the hadrons reflecting dynamics of quarks and gluons, such as the pressure distributions inside the hadrons [9,12,13]; indeed, the recent results of the pressure distributions inside the nucleon [14] are based on the determination of the gravitational form factors from the behaviors of the generalized parton distributions (GPDs) [15][16][17], which are obtained by experiments like deeply virtual Compton scattering (DVCS) [16][17][18][19][20][21], deeply virtual meson production [22,23], meson-induced Drell-Yan production [24][25][26], etc. (see also [27,28]).
As another phenomenological implication, the cross section of the near-threshold photoproduction of J/ψ in ep scattering [29] is sensitive to the F 2 part of the trace anomaly (1) [30], which can be conveniently handled [31] through the P ′ → P behavior of the gravitational form factors that represent P ′ |η µν T µν q,g R |P . Also, apparently, the decomposed quark and gluon contributions to the trace anomaly should provide a new insight on understanding the origin of the nucleon mass, one of the main objectives of the future Electron-Ion Collider.
The two-loop results corresponding to (3) and (4) have been derived [5] working out the renormalization mixing among the operators arising in the quark and gluon energymomentum tensors, where the two-loop anomalous dimensions for the second moment of the twist-two quark/gluon distribution functions, as well as for the mixing of the twist-four operators F µν F µν , mψψ, have been used as the necessary input informations. Because all such input informations appear to be available at the three-loop order, we extend, in this paper, the quark/gluon decomposition of the QCD trace anomaly to the three-loop level.
We work out the renormalization-mixing structure relevant for the quark and gluon energymomentum tensors in the MS scheme in Sec II, present their explicit results at three loops, and, as their direct consequence, derive the three-loop results of quark/gluon trace anomaly in Sec. III. As an application of our results, the anomaly-induced mass structure of nucleon as well as pion is discussed in Sec. IV. Sec. V is reserved for conclusions.

II. OPERATOR MIXING IN RENORMALIZATION OF ENERGY-MOMENTUM TENSOR
We follow, here and in the following, the notations and conventions of [5]. Our starting point is the gauge-invariant, symmetric QCD energy momentum tensor which is given by [7] (see also the article by Jackiw in [32]) where the quark and gluon parts read, respectively, in terms of the bare quark/gluon fields, Here, we have neglected the ghost and gauge fixing terms as they do not affect our final results. T µν is conserved and therefore it is a finite, scale-independent operator. However, T µν g and T µν q are not conserved separately (see (94), (95) below) and are subject to regularization and renormalization. Their traceless part and trace part are expressed by the twist-two and twist-four operators, respectively. Therefore, in order to derive the quark and gluon contributions to the quantum anomalies in their trace parts, we have to work out the renormalization structure of the corresponding twist-two as well as twist-four operators associated with (9), clarifying their renormalization mixing. For this purpose, let us write and, using this basis of operators O k (k = 1, 2, 3, 4), T µν of (8) is expressed as We introduce the renormalized composite operators O R k and the corresponding renormalization constants as Here, for simplicity, we do not explicitly show the mixing with the equations-of-motion operators as well as the BRST-exact operators; their matrix elements sandwiched between physical states vanish (see e.g., [4,[33][34][35] respectively, and subtracting those traces from both sides of (17) and (15), we obtain, The differentiation of these relations with respect to the renormalization scale µ yields the RG equations of the twist-two, spin-2 quark and gluon operators, which should coincide with the second moment of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations for the flavor-singlet part of the unpolarized parton distribution functions: with the anomalous dimension matrix, where the ellipses denote the two-and higher-loop terms, given as the second moment of the singlet DGLAP kernel. In mass-independent regularization schemes like the dimensional regularization, the anomalous dimensions of (23) depend on the scale µ only through that of the strong coupling constant α s (µ). As a result, (22) where the beta function is defined as and T (αs) means an ordering in the coupling constant such that the couplings increase from right to left (µ < µ ′ ). In particular, employing the dimensional regularization in d = 4 − 2ǫ space-time dimensions, the µ ′ → ∞ limit of (24) leads to the result,  where (25) implies, This result (26) determines the renormalization constants in (20), (21) as in terms of the anomalous dimensions of (23) in the MS-like schemes to a desired accuracy.
The renormalization constants Z M , Z S in (15), as well as Z K , Z B in (17), are still to be determined. For this purpose, we need to analyze the relevant renormalization structure at the twist-four level by treating the explicit form of the "traces" in (19). For the bare operators, it is straightforward to construct the corresponding contributions as Here, in the second formula, we have used the equations of motion for the quark fields, On the other hand, such straightforward manipulation to construct the trace terms is not useful for the renormalized operators in (19), because the trace operation and the renormalization do not commute, i.e., this is a general phenomenon one encounters when renormalizing a certain operator θ µν , that behaves as a symmetric second rank tensor [4,36].
In the dimensional regularization with d = 4 − 2ǫ dimensions, the counter terms to subtract the UV divergences in the bare operator as, may contain a term proportional to η µν /ǫ. Then, the contribution of this type, η µν /ǫ, produces, for the trace of the renormalized operator, η µν (θ µν ) R , in the MS, while it produces 4/ǫ as a counter term to define θ µ µ R , giving rise to a finite difference between η µν (θ µν ) R and θ µ µ R . Finite contributions of this type represent the trace anomaly in quantum field theories. In particular, such finite contributions for η µν −F µλ F ν λ R , the trace of O R 1 , appear to produce not only the contribution of O R 2 but also of another twist-four operator O R 4 , by contract to (29) [5]; similarly, Thus, we express the corresponding contributions as using certain coefficients x 1 , y 1 , x 3 , and y 3 of order α s and higher, which are to be determined below. Substituting (29)-(34) into (15) and (17), we obtain, where the twist-two operators dropped out using (20), (21), and the further substitution of (16) and (18) yields the relations among the twist-four bare operators as leading to the four conditions, Recalling that Z T , Z L , Z Q , and Z ψ are already determined as (28), these conditions allow us to determine the renormalization constants in the LHS, once Z F , Z C , and the coefficients x 1 , y 1 , x 3 , y 3 are given. We note that the Feynman diagram calculation of Z F and Z C is available to the two-loop order in the literature [37].
Now that explicit forms of all six renormalization constants arising in the RHS of (39)- (42) are available to a certain accuracy, we are able to fix also the coefficients x 1 , y 1 , x 3 , and y 3 , invoking the following property of the renormalization constants in the MS-like schemes: in this scheme, the renormalization constants take the form, with X = T, L, Q, ψ, F , and C; here, a X , b X , c X , . . . , are constants depending only on α s , and δ X,X ′ denotes the Kronecker symbol. Laurent expansion of the RHS of (39)-(42) about ǫ = 0 would produce the O(ǫ 0 ) terms as respectively, but these four formulas equal zero because Z M , Z S , Z B , and Z K should also obey the form of (43). These four conditions allow us to determine x 1 , y 1 , x 3 , and y 3 as power series in α s to the accuracy same as the renormalization constants (43). Then, this result of x 1 , y 1 , x 3 , y 3 allows us to determine Z M , Z S , Z B , and Z K to the same accuracy using (39)-(42); furthermore, using (33) and (34), the result of x 1 , y 1 , x 3 , and y 3 allows us to derive the gluon/quark individual contributions to the trace anomaly as, It is worth mentioning that the obtained results of the renormalization constants satisfy the following constraints, such that T µν = (T µν ) R of (2) is obeyed with (8)- (18). Using (35), (36), and (16), the above results (48) and (49) for individual anomalies may be reexpressed as in terms of the bare operators multiplied by the renormalization constants; adding these formulas and using the constraints (50)-(53), we get Here, in the last equality, we have used (16) and that ǫ (F 2 ) R → 0, as ǫ → 0. Note that this final form coincides with the asymptotic (µ → ∞) limit of the RHS of (1) with the use of (27), so that where we have used the RG invariance of T λ λ . This result shows that the consistency conditions (50)-(53) ensure that the total sum of the individual anomaly contributions (48) and (49) correctly reproduces the QCD trace anomaly (1) at arbitrary order in α s . The formulas derived above hold to arbitrary order in perturbation theory. It is demonstrated in [5] that the above logic leading to the formulas (48) and (49) indeed works at the two-loop level: substituting the anomalous dimension matrix γ(α s ) at two loops [38,39] into (26), we are able to obtain the two-loop formulas of Z T , Z L , Z Q , and Z ψ using (28), and, combined with the two-loop result [37] of Z F and Z C , the coefficients x 1 , y 1 , x 3 , and y 3 are uniquely determined to order α 2 s . As a result, the two-loop result for Z M , Z S , Z B , Z K and that for (48) and (49) are obtained, such that their one-loop terms read (3) and (4).
We emphasize that the input information necessary for this method is the anomalous dimension matrix γ(α s ) of (23) and the renormalization constants Z F and Z C of (16).
The anomalous dimension matrix γ(α s ) for the flavor-singlet part of the DGLAP evolution equation is now known to three-loop order [39,40]. On the other hand, for Z F and Z C beyond two loops, there seems no Feynman-diagram calculation performing the renormalization of the higher-twist scalar operators F 2 andψψ. Fortunately, however, the constraints imposed by the RG invariance of the energy-momentum tensor are strong enough to determine the form of Z F as well as Z C : using the RG invariance of (56) and (1), we have, where we have substituted (16) and (18); this yields the relations, where we have substituted (27) and the anomalous dimension for the quark mass in the form, The power series expansion of (59), (60) in α s reads These formulas completely determine Z F and Z C from the beta function β and the mass anomalous dimension γ m in the MS-like scheme. Indeed, these formulas correctly reproduce the two-loop results of Z F and Z C of [37], which were obtained by the Feynman diagram calculation. As another check, substituting (62) and (63) into the µ derivative of (16), the result can be recast into the following RG equation, which, alternatively, may be derived as a direct consequence of the fact that the total trace anomaly (1) is RG invariant, d dµ T λ λ = 0. The mass anomalous dimension (61) is now known to four-loop order in the literature [41,42]. Therefore, combined with the three-loop DGLAP anomalous dimensions for γ(α s ) in (22)-(28), the formulas derived in this section allow us to work out the new three-loop result for the gluon/quark trace anomalies (48) and (49), which is presented in the next section.

III. RESULTS AT THREE LOOPS
We present the three-loop results in the MS or MS schemes, using the general formulas of Sec. II with the three-loop input informations substituted. For this, the beta-function (27) is given with β 0 of (6) and and the mass anomalous dimension (61) reads [41,42] where ζ(s) is the Riemann zeta-function with ζ(3) = 1.202056903 . . .. The three-loop anomalous dimension matrix (23) for the twist-two flavor-singlet operators reads [39,40] γ(α s ) = α s 4π where and Then, using (28), we obtain the corresponding three-loop renormalization constants as On the other hand, Z F and Z C at three loops are given by (62) and the coefficients in (33) and (34) for the trace after renormalization are obtained as Substitution of these results into (48) and (49) leads to the main result of this paper, corresponding to the three-loop extension of (4) and (3), respectively. In accord with the general result (57), we note that the sum of these two equations reproduces the three-loop expression of (1) and is thus RG-invariant. However, each of them exhibits the dependence on the RG scale µ, i.e., due to the contributions of order α 2 s and higher (see the discussion in Sec. I). The formulas (87)-(89) and the other formulas obtained in this section should be evaluated using the three-loop running coupling constant, which is obtained by solving the RG equation (see (27), (66), (67)), where the constant of integration is represented by the the QCD scale parameter Λ QCD according to the definition in [43,44], and this result may be further solved for α s (µ) iteratively, leading to where L ≡ ln µ 2 /Λ 2 QCD . For example, using the RunDec package [45], we obtain, Before ending this section, we mention an immediate consequence of our three-loop results, combined with the exact operator identities in QCD, which read [10,46,47] ∂ ν T µν q =ψgF µν γ ν ψ , up to the terms which vanish using the equations of motion, (i / D − m) ψ = 0, and Note that (94) and (95) are compatible with the condition (2), using the equations of motion for the gluon fields, D α F αν = gψγ ν ψ, and the fact that the equations of motion are preserved under renormalization. Then, the ∂ µ -derivative of (17) gives Here, the last term plays roles for cases beyond one loop, because Z B − Z Q 4 = O(α 2 s ) from the results (76) and (80). It is remarkable that (96), together with the above results of the three-loop renormalization constants, leads to the three-loop evolution equation for the twist-four quark-gluon operator, The value (93) corresponds to Λ QCD ≃ 0.336 GeV, when (91) is used with n f = 3 fixed.

IV. ANOMALY-INDUCED MASS STRUCTURE OF HADRONS
As mentioned in Sec. I, it is well-known that the QCD trace anomaly (1) is related to the nucleon mass M, as This indicates that almost all of the nucleon mass could be attributed to the quantum loop effects in QCD which induce the trace anomaly. Based on this equation, it is also argued frequently that, in the chiral limit, the entire mass comes from gluons. However, the partition of QCD loop effects for the trace anomaly into the gluon and quark contributions, (4) and (3), and their three-loop extension (87) and (88), shows that the latter statement is not correct: evaluating (87)-(89) with N c = 3, n f = 3, one finds, where where mψψ R is independent of the scale µ, while the µ dependence of (F 2 ) R µ is controlled by the RG equation (65); substituting (93), we obtain Because α 3 s (1 GeV) ≃ 0.1, the neglected four-loop contributions are expected to produce corrections less than ten percent.
For simplicity, in the chiral limit, we find so that the gluon-and quark-loop effects make the nucleon mass heavy and light, respectively, with the magnitude of the former being five times larger than that of the latter, and it appears that the µ-dependence of this result for the relative size of the gluon/quark loop effects in the chiral limit is rather weak. It is also worth noting that the total sum (101) of (104) and (105) allows us to constrain the matrix element of F 2 as When taking into account the quark-mass effects, the matrix element of the quark scalar operator, P |mψψ|P , participates; its value may be constrained from the informations on the sigma terms (see, e.g., [48][49][50][51]), but we do not go into the detail here.
Next, we consider the pion case, making the substitution: |P → |π(p) with p 2 = m 2 π . Then, (100) becomes which implies, in the chiral limit m = 0, Since the PCAC relation (f π is the pion decay constant), indicates m 2 π ∼ m as m → 0, (108) gives the relation among the O(m) terms when the substitution π(p) → π(p) 0 + π(p) 1 + . . ., where π(p) 0 ≡ π(p) m=0 and π(p) 1 is the O(m 1 )-term, is made, such that π(p) mψψ R π(p) → 0 π(p) mψψ R π(p) 0 and . On the other hand, the pion mass can also be calculated as the mass shift due to the ordinary first-order perturbation theory in the quark mass term in the QCD Hamiltonian, as [48] m 2 π = 0 π(p) mψψ π(p) 0 , and, combining this with (108), we obtain to O(m) accuracy. Therefore, up to the corrections of O(m 2 ), the terms associated with the F 2 operator and the mψψ operator in the RHS of (108) contribute to m 2 π according to the relative weights, (1 − γ m (g R )) and (1 + γ m (g R )), respectively, where and 1 2m 2 π π(p) η λν T λν g R µ=1 GeV π(p) = 0.521 , which hold to O(m) accuracy. As in the nucleon case, the µ-dependence of the result is rather weak, but this result shows the structure which is completely different from the nucleon case: both the gluon-and quark-loop effects make the pion mass heavy, each of them giving rise to roughly half of the pion mass. This remarkable feature may be eventually attributed to particular nature of the pion as a Nambu-Goldstone boson, but the clarification needs further studies, including those to make connections with the QCD analysis of hadron masses based on the decomposition of the QCD Hamiltonian [49,50,52] into the quark/gluon components.

V. CONCLUSIONS
We have studied the renormalization of the energy-momentum tensor in QCD at threeloop order and derived the analytic formulas for all the relevant renormalization constants in the MS scheme. To achieve this, we did not need any new loop calculation. We explained our method which allows us to determine the structure of the renormalization mixing involving the twist-four as well as twist-two operators, from the knowledge of the basic RG functions, β(g) and γ m (g), and the anomalous dimensions for twist-two, spin-2 operators, which are available to the required accuracy in the literature. This method utilized, in particular, the constraints from the fact that the trace of the total energy-momentum tensor is determined completely by the basic RG functions in the RG-invariant form, and that the renormalization constants in the MS-like schemes obey the divergent pole structures independent of any mass.
Determining the renormalization structure at three loops, we were immediately able to calculate the trace contributions for the quark and gluon parts of the energy-momentum tensor separately, leading to the decomposition of the three-loop QCD trace anomaly into the quark and gluon parts. Our result extends the previous result in the two-loop order into three loops, and the features found in the former case are preserved at the three-loop level, although their analytic expressions become considerably complicated and cumbersome.
Going from the two-loop formula to the three-loop formula, we have found that the accuracy is generally improved from ten percent level to a few percent level, and, therefore, our three-loop formula should provide good control on the uncertainties in the calculations for most practical purposes. Our formula obtained in the MS scheme will be most convenient for various actual applications, but it would be interesting to derive the similar formula using other regularization schemes, e.g., the gradient flow regularization [36,53], and study the scheme dependence.
As an immediate application, we used our three-loop formula to analyze the anomalyinduced mass structure of nucleon as well as of pion. Our three-loop formula allowed us to carry out such analysis at a few percent-level accuracy, and provided a new insight on the hadron masses such that they are generated from the quantum anomaly effects for quarks as well as for gluons. Our analysis revealed that the nature of the masses, in particular, the dominant roles played by the quark part and the gluon part of the energy-momentum tensor, are completely different between nucleon and pion. Another application, which was not discussed in this paper, is to constrain the gravitational form factorC q,g of a hadron, as treated in [5].C q,g receives much attention in connection with the force distribution inside the nucleon [9,11,31,52,54] and the nucleon's transverse spin sum rule [55][56][57].
Quantitative analysis for the constraints onC q,g will be discussed elsewhere.

ACKNOWLEDGMENTS
The author thanks Yoshitaka Hatta and Abha Rajan for stimulating discussions in the collaboration for [5] which triggered this work.