The two-loop energy-momentum tensor within the gradient-flow formalism

The gradient-flow formulation of the energy-momentum tensor of QCD is extended to NNLO perturbation theory. This means that the Wilson coefficients which multiply the flowed operators in the corresponding expression for the regular energy-momentum tensor are calculated to this order. The result has been obtained by applying modern tools of regular perturbation theory, reducing the occurring two-loop integrals, which also include flow-time integrations, to a small set of master integrals which can be calculated analytically.


Introduction
The gradient-flow formalism as introduced by Lüscher [1] and further formalized by Lüscher and Weisz [2] has proven useful in lattice QCD in many respects. One of its main virtues is that composite operators at finite flow time t do not require ultra-violet (UV) renormalization beyond the one of the involved parameters and fields. This means that the operators do not mix under renormalization-group running, which makes it particularly simple to combine results from different regularization schemes. This feature opens promising prospects for a cross-fertilization of lattice and perturbative calculations, such as a possible lattice determination of α s (M Z ), for example [3].
A particularly powerful way to exhibit this possible interplay is obtained by considering the expansion of composite operators in the limit of small flow time, which expresses flowed operators in terms of QCD operators at t = 0, with t-dependent Wilson coefficients [2]. This method has been used by Makino and Suzuki [4,5] to derive a regularization-independent formula for the energya e-mail: harlander@physik.rwth-aachen.de b e-mail: kluth@physik.rwth-aachen.de c e-mail: flange@physik.rwth-aachen.de momentum tensor (EMT) T μν which has already led to promising results (see, e.g., Refs. [6][7][8][9][10]).
The universal Wilson coefficients that occur in the formula of Ref. [5] for the EMT have been calculated through next-toleading order (NLO) in perturbation theory [4,5]. This corresponds to a one-loop calculation in the sense that it involves integrals over a single D-dimensional momentum. In this paper, we will carry this calculation to the next perturbative order. 1 It is important to note at this point that the integrals which occur in the gradient-flow formalism are of a more general type than in regular QCD. They involve additional exponential factors which depend on loop and external momenta, as well as on flow-time variables, some of which are also integrated over. Nevertheless, the first two-loop result was already obtained in Ref. [1], even in analytic form. The extension to the three-loop level required significant aid from computer algebra and numerical tools [3]. From the quantumfield theoretical point of view, it closely followed the steps of Ref. [1] by directly expressing the Green's functions in terms of integrals with the help of Wick's theorem. The integrals themselves were evaluated using sector-decomposition [12,13] in order to isolate the poles in D − 4, whose coefficients were determined using high-precision numerical methods [14,15].
In the current calculation, we apply a completely independent setup. On the one hand, it applies the gradient-flow formalism described in terms of a five-dimensional quantum field theory [2], which leads to well-defined, albeit nonstandard Feynman rules. On the other hand, rather than evaluating the resulting integrals numerically, we express them in terms of master integrals using the integration-by-parts method of Chetyrkin and Tkachov [16]. This reduces the NLO calculation of the Wilson coefficients of the EMT to a single one-loop integral without flow-time integration. The next-to- next-to-leading order (NNLO) calculation leads to four twoloop master integrals without flow-time integration, and two two-loop master integrals with a single flow-time integration. All master integrals can be calculated analytically by standard means for general values of D, the number of space-time dimensions. By suitable renormalization, the Wilson coefficients of the EMT can be defined in such a way that they are formally renormalization-scale independent. For a fixed-order perturbative result, this means that the renormalization-scale dependence is formally of higher order. This allows one to estimate the perturbative uncertainty on the Wilson coefficients through their residual dependence on the renormalization scale μ around a particular "central" value. Based on the form of the analytical result, we argue for a specific choice of this central value. Our numerical study shows that the higher-order terms indeed lead to an appreciable reduction of the μ-variation. However, by comparison of the successive higher-order terms, it appears that the uncertainty estimate from a variation within μ ∈ [μ 0 /2, 2μ 0 ], as it is common practice in regular perturbative QCD calculation, might be too optimistic.
While we consider the NNLO expressions for the Wilson coefficients of the EMT as our main result, our calculation allows us to obtain a number of additional results that might be useful in a broader context. Among these are the flowed quark-field renormalization constant Z χ , and the matrix of anomalous dimensions for the set of operators which form the energy-momentum tensor in regular (non-flowed) QCD through NNLO. The remainder of the paper is structured as follows. After briefly introducing the perturbative gradient-flow formalism in order to define our notation in Sect. 2, we outline the approach of Refs. [4,5] for using this formalism to define the EMT in Sect. 3. Technical details of our calculation are described in Sect. 4. Section 5 contains our main result, the Wilson coefficients through NNLO QCD in the MS scheme. As pointed out in Ref. [5], the trace anomaly of the EMT allows for a welcome check of the calculation; we briefly describe the derivation of the resulting relations among the coefficient functions in Sect. 6. Finally, in Sect. 7, we use the finiteness condition of the flowed operators in order to derive the anomalous-dimension matrix for the set of operators which occur in the EMT in regular QCD. Section 8 presents our conclusions.

QCD gradient flow in perturbation theory
In the following, we will work in D-dimensional Euclidean space-time with D = 4 − 2ε. The gradient-flow formalism continues the gluon and quark fields A a μ (x) and ψ(x) of reg-ular 2 QCD to (D +1)-dimensional fields B a μ (t, x) and χ(t, x) through the boundary conditions and the flow equations where the "flow time" t is a parameter of mass dimension minus two, and κ is an additional gauge parameter which drops out of physical observables. The (D + 1)-dimensional field-strength tensor is defined as the covariant derivative in the adjoint representation is given by and As usual, the color indices of the adjoint representation are denoted by a, b, c, . . ., while μ, ν, ρ, . . . are D-dimensional Lorentz indices. Color indices of the fundamental representation are suppressed throughout this paper, unless required by clarity. The symmetry generators T a obey the commutation relation with the structure constants f abc . The flow-field equation leads to a smearing of gauge-field configurations at finite flow time t > 0. As a consequence, composite operators at t > 0 do not require renormalization beyond the renormalization of the parameters and fields of the Lagrangian. For the strong coupling and the quark mass, the renormalization constants are identical to those at t = 0; the flowed gluon fields do not require renormalization at finite flow time as was pointed out in Ref. [2]. The renormalization constant for the flowed quark field through NNLO is a byproduct of this paper and will be given below.

Energy-momentum tensor
In a continuous D-dimensional space-time, the gauge invariant part of the EMT reads where g 0 is the bare coupling constant of QCD. The operators are defined as , , , where f labels the n F different quark flavors, m f,0 is the bare quark mass, and The notation i f ∈ {i 1 , . . . , i n F } for the indices which label different flavors will be useful later on in this paper. In general, T μν may contain gauge-dependent operators which vanish when evaluating physical matrix elements [17].
Here and in what follows, we implicitly assume that the vacuum expectations values of all composite operators have been subtracted 3 so that O i,μν (x) ≡ 0 ∀i. In this paper, we will focus on the case where physical matrix elements of the EMT itself are considered, i.e. no other operator multiplies the EMT at the same space-time point. In this case, the equations-of-motion (EOM) render the set of operators in Eq. (8) redundant. In particular, the EOM for the quark fields in regular QCD implies which allows us to eliminate O 5,μν from the set of operators in Eq. (8). Note that due to this relation, the last two terms in Eq. (7) cancel.
We will further assume all quark masses to be equal to each other 3 In other words, the precise definition of O 1,μν , for example, would be given by Therefore, the different quarks are indistinguishable and the mixing between two different quark flavors cannot depend on the flavors. Defining the analogous operators of Eq. (8) for flowed fields, we writẽ , , where Z χ f ≡ Z χ is the renormalization constant for the flowed quark fields, and Since we have eliminated O 5,μν from the set of operators by using Eq. (10), we do not need to include a flowed version of this operator in Eq. (12). Similar to the composite operators of regular QCD, we assume that the vacuum expectation values of the flowed composite operators have been subtracted, i.e.
We can now use the expansion in small flow time [2] to get a relation between flowed and regular QCD operators. In Eq. (14), and similarly in what follows, a sum 4 j=1 is understood. The ellipsis denotes terms that vanish as t → 0 which will be neglected throughout this paper. As discussed above, matrix elements of the l.h.s. of this equation are finite after renormalization of the QCD parameters, while those of the regular QCD operators on the r.h.s. are in general divergent. The mixing matrix ζ i j (t) will therefore be divergent as well.
Inverting Eq. (14) and using it to re-express the regular QCD operators in the energy-momentum tensor in terms of flowed fields, one arrives at where Since matrix elements of theÕ i as well as the energymomentum tensor itself are finite (after mass and charge renormalization), the universal coefficients c i (t) of Eq. (16) are finite as well. In Ref. [5], they have been calculated in perturbation theory through NLO QCD. The goal of the current paper is to evaluate them through NNLO QCD.

Method of projectors
To compute the coefficients ζ i j (t) we use the so-called "method of projectors" [18,19], which consists of constructing external states |k and differential operators D k for which where we have dropped the Lorentz indices for convenience, and we define the matrix element to include only diagrams which are one-particle irreducible (1PI) with respect to (w.r.t.) QCD particles. Applying P k on both sides of Eq. (14), one obtains Since the ζ i j (t) only depend on the flow time t and the renormalization scale μ, we can choose arbitrary values for all other dimensional parameters in this equation. Setting them to zero turns all higher-order corrections on the r.h.s. into massless tadpoles, so that Eq. (17) is only required to hold at tree-level. One thus obtains where m and p collectively denote all masses and external momenta. The right-hand side thus results in vacuum diagrams whose only dimensional scale is t.
In order to find suitable projectors, we first derive the Feynman rules for the relevant terms of the operators. For example, 4 (20) where the momenta are defined to be outgoing. This suggests to use where N A is the dimension of the adjoint representation of the gauge group; for SU(N c ), it is N A = N 2 c − 1. The projector 4 Feynman diagrams in this paper were drawn using TikZ-Feynman [20].
onto the Lorentz structure is defined by for any other linearly independent Lorentz tensor.
In the appendix, one can find the relevant parts of the Feynman rules for the other operators, which in a similar way lead to the projectors , , where N c is the dimension of the fundamental representation of the gauge group, i.e. the number of colors, and i and j are the corresponding indices. The trace appearing in the projectors P 3 f and P 4 f is taken w.r.t. the spinor indices of the Green's function, and f denotes the associated quark flavor. Note that P 4 f is constructed such that in order to ensure that only those fermionic operators are taken into account which do not vanish according to the EOM, see Eq. (10).
With this procedure, we get a mixing matrix which distinguishes between different quark flavors. To avoid confusion with ζ i j (t), which is the mixing matrix between operators summed over all flavors, we will call it Ω i j (t). This matrix is then defined bỹ where double indices in this expression are summed over {1, 2, 3 1 , . . . , 3 n F , 4 1 , . . . , 4 n F } (see Eq. (8)), in contrast to Eq. (15), where double indices are summed over {1, 2, 3, 4}. Its general structure is given by where an element Ω i j represents the mixing between O i andÕ j , taking into account individual flavors. Therefore, an underlined and a double underlined element denotes an n F -dimensional vector and an n F × n F dimensional matrix, respectively. Their elements describe the mixing between different flavors. As the quarks are indistinguishable, the n F and n 2 F dimensional objects appearing in Eq. (26) can each be described by two independent parameters, named ω i j and ω i j : Summing over the different flavors occurring in Eq. (26), the relation between Ω(t) and ζ(t) can be easily established by

Computational methods
The gradient-flow formalism in perturbation theory can be formulated in terms of a Lagrangian field theory, where the flow equations (2) are implemented with the help of Lagrange-multiplier fields [2]. The crucial difference between the regular QCD Feynman rules and those in the gradient-flow formalism is the occurrence of exponential factors exp(−sp 2 ), where s is a "flow-time variable", and p the linear combination of D-dimensional external and/or loop momenta. Vertices involving flowed fields induce an integration over all positive values of the corresponding flowtime variable, which is, however, bounded from above by "propagators" of the Lagrange-multiplier fields, since they introduce step functions of the flow-time variables.
We have implemented the Feynman rules into the program qgraf [21,22], which generates the Feynman diagrams for the desired matrix elements. Its output is then transformed to FORM [23,24] notation by q2e/exp [25,26]. An in-house set of FORM routines inserts the Feynman rules, performs the projections onto the relevant color and Lorentz structures according to the P j of Eq. (21) and (23), and evaluates the Dirac traces. The result is then expressed in terms of a linear combination of integrals whose general form is as follows: ((d 1 , . . . , d f ), (b 1 , . . . , b n ), (a 1 , . . . , a n )) where the a i and d i are integers (d i ≥ 0), f and l is the number of flow-time and loop integrations, respectively, the b j are polynomials in (rescaled) flow-time parameters u i and the q i are linear combinations of the loop momenta k j . For the problem and the perturbative order under consideration, it is 0 ≤ f ≤ 4, 1 ≤ l ≤ 2, and 0 ≤ n ≤ 3, respectively. Note that the projectors defined in Eqs. (21) and (23) eliminate all dependence on external momenta and masses, so that, after making a suitable ansatz for the index structure of the integrals, we only have to evaluate scalar vacuum integrals. Using the identities [16] and similar ones for the flow-time integrations, one can derive relations among these integrals by explicitly performing the derivatives in the integrand on the l.h.s. These so-called "integration-by-parts (IBP) relations" were fed to Kira [27] which allowed us to reduce all occurring integrals to a single master integral at one-loop level, and six master integrals at two-loop level using the Laporta algorithm [28]. 5 Their analytical evaluation is possible along the lines of Ref.

Coefficient functions through NNLO QCD
The strong coupling and the quark mass require the regular QCD renormalization according to where we write the renormalization constants Z g and Z m as C F and C A are the quadratic Casimir eigenvalues of the fundamental and the adjoint representation of the gauge group, respectively. Furthermore, T F = T n F , with n F the number of quark flavors, and T the trace normalization in the fundamental representation. For SU(N c ), it is C F = (N 2 c − 1)/(2N c ), C A = N c , and T = 1/2. In addition, the flowed quark fields also require renormalization according to leading to the factor Z χ in the definition of the operatorsÕ 3,4 in Eq. (12). Z χ differs from the quark-field renormalization of regular QCD. Through NLO, the MS result has been evaluated in Ref. [33]. To determine the renormalization constant at NNLO as required by our calculation, we may use the fact that the coefficient c 3 (t) must be finite after renormalization. Writing the MS expression as we find 7 This allows us to evaluate the coefficients of the energymomentum tensor in the MS scheme through NNLO QCD: 8   (47) where we introduced the parameter 9 with the Euler-Mascheroni constant γ E = 0.57721 . . .. Through NLO, these results are in full agreement with those of Refs. [4,5]. We have carried out the calculation in the general R ξ gauge of regular QCD; the fact that the gaugeparameter dependence cancels in the final result serves as another welcome check. The gauge parameter κ of Eq. (2) has been set to 1. While the energy-momentum tensor T μν (x) is renormalization-scheme independent, this is not necessarily the case for the operatorsÕ i,μν (t, x) and the coefficient functions c i (t). SinceÕ 1,μν andÕ 2,μν do not require operator renormalization, their matrix elements as well as the coefficient function are indeed renormalization-scheme independent. On the other hand, using the quark-field renormalization Z χ of Eq. (43) in the MS scheme, matrix elements ofÕ i,μν and coefficient functions c i (t) become explicitly dependent on the renormalization scale μ for i ∈ {3, 4}.
However, this renormalization-scheme dependence can be avoided by introducing so-called "ringed" quark fields as suggested in Ref. [5]. This corresponds to replacing Z χ in Eq. (12) bẙ Eur. Phys. J. C (2018) 78:944 Currently,Z χ (t) is available only through NLO QCD. Its explicitly μ-dependent terms can be reconstructed from the requirement that c i (t) must be finite and μ-independent for i ∈ {3, 4} though. In this way we find for the ratio to the MS quark-field renormalization constant of Eq. (43): The constant C 2 cannot be determined in this way, but requires a dedicated three-loop calculation of the two-point function occurring in the denominator of Eq. (50). A detailed outline of this calculation is beyond the scope of this paper; it will be presented together with a more complete description of our setup in a forthcoming publication [32]. At this point, we simply quote the numerical value of this result up to three significant digits, which is more than sufficient in the light of the theoretical uncertainties to be discussed below: 10 Multiplication of c 3 (t) and c 4 (t) in Eqs. (47) and (48) by this ratio makes also these coefficients formally μ-independent, i.e., As in any perturbative calculation, the μ-independence only holds up to higher orders in g. The decrease of the residual μdependence is thus commonly used as a qualitative check of the perturbation expansion for the specific observable under consideration. We thus study the μ-dependence of the four coefficients after dividing c 3 and c 4 by the ratio ζ χ defined in Eq. (51). We fix a characteristic value for the flow time t and vary the renormalization scale μ around the central value μ 0 , which we define such that L(μ 0 , t) = 0, cf. Eq. (49), i.e. and g (n F =5) (130 GeV) = 1.19. The μ-variation of the strong coupling constant g(μ) is determined by numerically solving the corresponding renormalization-group equation with the help of RunDec [34,35] at one-, two-, and three-loop level for the LO, NLO, and the NNLO curve, respectively. In Fig. 1, the value of t is chosen such that the central scale of Eq. (54) is μ 0 = 3 GeV. At this central scale, the NNLO corrections increase the modulus of the coefficients c 1 and c 2 by 10% and 13% relative to NLO, respectively. This is within twice the NLO uncertainty due to missing higher-order effects as estimated by varying μ/μ 0 between 1/2 and 2, where one finds 7.3% for c 1 , and 8.0% for c 2 . We are therefore confident that the NNLO uncertainty estimated in the same way is rather reliable: it is given by 5.7% for c 1 and 7.2% for c 2 . Note that the dominant contribution to these numbers comes from the downward variation of μ, where g(μ) starts to become sensitive to the non-perturbative region. The behavior of c 1 and c 2 towards larger values of μ seems to suggest that this uncertainty estimate may actually be too conservative.
As opposed to the gluonic coefficients c 1 and c 2 , the coefficients of the fermionic operators c 3 and c 4 exhibit a residual scale dependence only from the NLO term onwards. One therefore expects a stronger μ-dependence at NNLO for these terms. Nevertheless, forc 3 , the estimate of the theory uncertainty due to scale variation still decreases from 9.8% to 8.1%. The increase of the result due to the NNLO effects is 8.3% relative to the NLO result at μ = μ 0 .
The behavior ofc 4 , on the other hand, is less satisfactory at μ 0 = 3 GeV. The NNLO effects more than double the NLO result in this case, and the uncertainty estimate due to scale variation actually increases from 47% to 71% when going from NLO to NNLO. Note, however, that c 4 = 0 at LO, which means that this coefficient is numerically sub-dominant.
As one would expect, for μ 0 = 130 GeV, the perturbative behavior of all coefficients is significantly improved, cf. Fig. 2. For c 1 , c 2 , andc 3 , the scale uncertainty is at the subpercent level already at NLO; at NNLO, it amounts to less than 0.2% in all three cases. The effect of the NNLO corrections relative to the NLO result is about 2% for c 1 and c 2 , and 0.8% forc 3 . Also forc 4 , the situation improves significantly: the NNLO terms add 38% to the NLO result, and the uncertainty goes down from 10% at NLO to 5.8% at NNLO. It is also worth pointing out that the choice of the central scale μ 0 as defined in Eq. (54) seems justified by the behavior of the successively higher orders. In almost all cases, the NLO and the NNLO corrections are both relatively small at μ = μ 0 . At the same time, the NNLO corrections relative to the NLO result are always smaller than the NLO corrections compared to the LO result. The only exception to this isc 4 at μ 0 = 3 GeV, where, however, no choice of μ seems to stand out over any other. In summary, we conclude that the NNLO terms lead to a significant improvement of the perturbative accuracy of the Wilson coefficients.

Trace anomaly
As a test of our result, we use the trace anomaly of the EMT. As suggested in Ref. [17], a simple derivation consists of taking the trace of the EMT in D = 4 − 2ε dimensions. By use of the equations of motion, this gives for the gauge invariant part where we have used Eq. (10) in the last step. Using the mixing matrix ζ i j (t), we can rewrite this in terms of flowed operators: Note thatc 1 (t) =c 3 (t) = 0, asÕ 1,μν (t, x) andÕ 3,μν (t, x) have a non-trivial index structure and therefore O 2,μν (x) and O 4,μν (x) cannot mix with them. SinceÕ 2,μμ = DÕ 1,μμ and 2Õ 4,μμ = DÕ 3,μμ , we cannot equate coefficients with Eq. (15) for all i individually. Instead, only the weaker conditions can be derived. We checked that these equations are indeed fulfilled by our result.

Operator renormalization
Using the fact that flowed operators are finite after mass and field renormalization, we can also compute the renor- This multiplication of O 1,μν and O 2,μν by 1/g 2 0 ensures that the mass dimension of all operatorsÔ i,μν is equal to D. The renormalization matrix is then defined as Expressing theÔ i,μν (x) in terms of flowed operators, one can determine its entries in the MS scheme by demanding that be finite. In analogy to Eqs. (40) and (43), we write Our result for the coefficients of the anomalous dimension at NLO is in agreement with Ref. [5]: The renormalization matrix Z i j and the mixing matrix ζ i j are provided in ancillary files with this paper. In this way, one obtains the following expression for the energy-density operator in terms of flowed operators, for example: × g 2 (4π) 2 C F (5 + 6L(μ, t)) + g 4 (4π) 4 × γ χ,0 β 0 + γ χ,0 2 L 2 (μ, t) Through NLO, this result agrees with Refs. [5,36]. Similar relations can be derived for all other operators of Eq. (8), of course.

Conclusions
We have presented the universal Wilson coefficients for the gradient-flow definition of the energy-momentum tensor through NNLO QCD. The NNLO corrections modify the three numerically dominant coefficients c 1 , c 2 ,c 3 at the level of 10% (1-2%) for a central scale of μ 0 = 3 GeV (μ 0 = 130 GeV), where μ 0 is related to the flow time t according to Eq. (54). We observe a reduction of the theoretical uncertainty relative to the NLO result as derived from varying the renormalization scale by a factor of two around its central value. The behavior of the fourth coefficientc 4 is less satisfactory, but its impact is expected to be numerically suppressed.
Aside from this main outcome, new results presented in this paper include the flowed quark-field renormalization constant to NNLO in the MS scheme, and the anomalous dimension matrix for the regular QCD operators which make up the EMT. In conclusion, we hope that our results will help to improve the studies of the EMT on the lattice. They are the first outcome of a systematic setup for higher-order perturbative calculations within the gradient-flow formalism [32], which should prove useful also in other applications of this theoretical framework.