Lattice energy-momentum tensor from the Yang-Mills gradient flow -- inclusion of fermion fields

Local products of fields deformed by the so-called Yang--Mills gradient flow become renormalized composite operators. This fact has been utilized to construct a correctly normalized conserved energy--momentum tensor in the lattice formulation of the pure Yang--Mills theory. In the present paper, this construction is further generalized for vector-like gauge theories containing fermions.


Introduction and the main result
Energy and momentum are fundamental notions in physics. In lattice field theory (the best-developed non-perturbative formulation of quantum field theory), however, the construction of the corresponding Noether current, the energy-momentum tensor [1][2][3][4], is not straightforward [5,6], because the translational invariance is explicitly broken by the lattice structure. In a recent paper [7], a possible method to avoid this complication being inherent in the energy-momentum tensor on the lattice has been proposed on the basis of the Yang-Mills gradient flow (or the Wilson flow in the context of lattice gauge theory) [8][9][10]. 1 See Ref. [18] for a recent review on the gradient flow and Refs. [19][20][21][22][23][24][25] for its applications in lattice gauge theory.
The basic idea of Ref. [7] is the following: 2 Consider the pure Yang-Mills theory. The gradient flow deforms the bare gauge field A μ (x) according to a flow equation with a flow time t [Eq. (3.1) below] and this makes gauge field configurations "smooth" for positive flow times. It can then be shown to all orders in perturbation theory that any local products of the flowed gauge field B μ (t, x) for any strictly positive flow time t is ultraviolet (UV) finite when expressed in terms of renormalized parameters [9]. In particular, no multiplicative renormalization factor is required to make those local products finite. In other words, they are renormalized composite operators. Such UV-finite quantities should be "universal" in the sense that they are independent of the UV regularization chosen, in the limit in which the regulator is removed. This suggests a possibility that by using the gradient flow PTEP 2014, 063B02 H. Makino and H. Suzuki as an intermediate tool one may bridge composite operators defined with the dimensional regularization, with which the translational invariance is manifest, 3 and those in the lattice regularization with which one may carry out non-perturbative calculations. 4 Following this idea, a formula that expresses the correctly normalized conserved energy-momentum tensor as a t → 0 limit of a certain combination of the flowed gauge field was derived [7]. This formula provides a possible method to compute correlation functions of the energy-momentum tensor by using the lattice regularization because the universal combination in the formula should be independent of the regularization. An interesting point is that the small flow-time behavior of the (universal) coefficients in the formula can be determined by perturbation theory thanks to the asymptotic freedom. This implies that if the lattice spacing is fine enough a further non-perturbative determination of the coefficients is not necessary. (Practically, non-perturbative determination of those coefficients may be quite useful and how this determination can be carried out has been investigated in Ref. [26].) Although the validity of the formula in Ref. [7], especially the restoration of the conservation law in the continuum limit, still remains to be carefully investigated, the measurement of the interaction measure (the trace anomaly) and the entropy density of the Yang-Mills theory at finite temperature on the basis of the formula [27] shows encouraging results; the method appears to be promising even practically.
In Ref. [7], the method was developed only for the pure Yang-Mills theory. It is then natural to ask for wider application if the method can be generalized to gauge theories containing matter fields, especially fermion fields. In the present paper, we work out this generalization. We thus suppose a vector-like gauge theory 5 with a gauge group G that contains N f Dirac fermions in the gauge representation R. (1.1) For simplicity, we assume that all N f fermions possess a common mass; this restriction might be appropriately relaxed. Our main result is Eq. (4.70), and a step-by-step derivation of this master formula is given in subsequent sections. For those who are interested mainly in the final result, here we give a brief explanation of how to read the master formula (4.70): The left-hand side is the correctly normalized conserved energy-momentum tensor (with the vacuum expectation value subtracted); our formula (4.70) holds only when the energy-momentum tensor is separated from other operators in correlation functions in the position space. The combinations in the right-hand side are defined by Eqs. (4.1)-(4.5). There, G a μν (t, x) and D μ are the field strength and the covariant derivative of the flowed gauge field, respectively (the definition of the flowed gauge field is identical to that of Refs. [8][9][10]); our ringed flowed fermion fields • χ(t, x) and • χ(t, x) are, on the other hand, somewhat different from those of Ref. [10], χ(t, x) andχ(t, x), and they are related by Eqs. (3.20) and (3.21). The coefficient functions c i (t) in Eq. (4.70) are given by Eqs. (4.72)-(4.76). There,ḡ(q) andm(q) are the running coupling and the running mass parameter defined by Eqs. (4.20) and (4.21), respectively; throughout this paper, m denotes the (common) renormalized mass of the fermions. Equations (4.72)-(4.76) are for the minimal subtraction (MS) scheme, and the result for the modified minimal subtraction (MS) scheme can be obtained by the replacement (4.77). b 0 and d 0 are the first coefficients of renormalization group functions, Eq. (2.16) and Eq. (2.18), respectively. The energy-momentum tensor is given by the 3 The drawback of the dimensional regularization is, of course, that it is defined only in perturbation theory. 4 We thus implicitly assume that the finiteness of the flowed fields that can rigorously be proven only in perturbation theory persists even in the non-perturbative level. 5 We assume that the theory is asymptotically free.

2/27
Downloaded from https://academic.oup.com/ptep/article-abstract/2014/6/063B02/1561362 by guest on 30 July 2018 t → 0 limit in the right-hand side of Eq. (4.70). As in the pure Yang-Mills case mentioned above, the combination in the right-hand side is UV finite and one may use the lattice regularization to compute the correlation functions of the combination in the right-hand side. In this way, correlation functions of the correctly normalized conserved energy-momentum tensor are obtained. To ensure the "universality," however, the continuum limit has to be taken before the t → 0 limit. Practically, with a finite lattice spacing a, the flow time t cannot be taken arbitrarily small to keep the contact with the continuum physics. Instead, we have a natural constraint, where R denotes a typical physical scale, such as the hadronic scale or the box size. The extrapolation for t → 0 thus generally requires a sufficiently fine lattice.
Here is our definition of the quadratic Casimir operators: We set the normalization of anti-Hermitian generators T a of the representation R as tr R (T a T b ) = −T (R)δ ab and T a T a = −C 2 (R)1. We also denote tr R (1) = dim(R). From the structure constants in [T a , T b ] = f abc T c , we define f acd f bcd = C 2 (G)δ ab . For example, for the fundamental N representation of SU (N ) for which dim(N ) = N , the conventional normalization is

Energy-momentum tensor with dimensional regularization
The description of the energy-momentum tensor in gauge theory [3,4] is particularly simple with the dimensional regularization. 6 This is because this regularization manifestly preserves the (vectorial) gauge symmetry and the translational invariance. Thus, in this section, we briefly recapitulate basic facts concerning the energy-momentum tensor on the basis of the dimensional regularization. The action of the system under consideration in a D dimensional Euclidean space is given by where g 0 and m 0 are bare gauge coupling and the mass parameter, respectively. The field strength is defined by for A μ (x) = A a μ (x)T a and F μν (x) = F a μν (x)T a , and the covariant derivative on the fermion is Here, and in what follows, the summation over N f fermion flavors is always suppressed. Our gamma matrices are Hermitian and for the trace over the spinor index we set tr1 = 4 for any D. We also set Assuming the dimensional regularization, one can derive a Ward-Takahashi relation associated with the translational invariance straightforwardly. We consider the following infinitesimal variation (2.5) Then, since the action changes as 7 where the energy-momentum tensor T μν (x) is defined by where O in (O out ) is a collection of gauge-invariant operators localized inside (outside) the finite integration region D. This relation shows that the energy-momentum tensor generates the infinitesimal translation and, at the same time, the bare quantity (2.7) does not receive the multiplicative renormalization. Thus, we define a renormalized finite energy-momentum tensor by subtracting its (possibly divergent) vacuum expectation value: (2.10) A fundamental property of the energy-momentum tensor is the trace anomaly [29][30][31][32]. One simple way to derive this [34,35] is to set ξ μ (x) ∝ x μ in Eq. (2.5) and compare the resulting relation with the renormalization group equation. After some consideration, this yields where we assume that the renormalized operators in the right-hand side are defined in the MS scheme [28]. 8 Throughout the present paper, we always assume that the vacuum expectation value is subtracted in renormalized operators. In Eq. (2.11), the renormalization group functions β and γ m are defined by where g and m are the renormalized gauge coupling and the renormalized mass, respectively, and the derivative with respect to the renormalization scale μ is taken with all bare quantities kept fixed. 7 Here, to make the perturbation theory well defined, we implicitly assume the existence of the gauge-fixing term and the Faddeev-Popov ghost fields. However, since they do not explicitly appear in correlation functions of gauge-invariant operators, we neglect these elements in what follows. 8 The renormalization in the MS scheme to the one-loop order is summarized in Appendix A. The renormalization constants are defined by

Flow equations and the perturbative expansion
The Yang-Mills gradient flow is a deformation of a gauge field configuration generated by a gradient flow in which the Yang-Mills action integral is regarded as a potential height. To be explicit, for the gauge potential A μ (x), the flow is defined by [8,9] where t is the flow time and G μν (t, x) is the field strength of the flowed field, and the covariant derivative on the gauge field is 3) The first term in the right-hand side of Eq. (3.1) is the "gradient" in the functional space, where S is the Yang-Mills action integral for the flowed field. Note that since is a sort of diffusion equation and the flow for t > 0 effectively suppresses high-frequency modes in the configuration. The second term in the right-hand side of Eq. (3.1) with the parameter α 0 is a "gauge-fixing term" that makes the perturbative expansion well defined. It can be shown, however, that any gauge-invariant quantities are independent of α 0 . In actual perturbative calculation in the next section, we adopt the "Feynman gauge" α 0 = 1 with which the expressions become simplest. The formal solution of Eq. (3.1) is given by [8,9] is the heat kernel and denotes non-linear interaction terms. By iteratively solving Eq. (3.4), we have a perturbative expansion for the flowed field B μ (t, x) in terms of the initial value A ν (y).
A similar flow may also be considered for fermion fields [10]. For our purpose, it is not necessary that the flow of fermion fields is a gradient flow of the original fermion action. A possible choice introduced in Ref. [10] is The formal solutions of the above flow equations are and The initial values for the above flow, A μ (x), ψ(x), andψ(x), are quantum fields being subject to the functional integral. The quantum correlation functions of the flowed fields are thus obtained by expressing the flowed fields in terms of original un-flowed fields (the initial values) and taking 9 Throughout the present paper, we use the abbreviation, the functional average of the latter. Equations (3.4), (3.12), and (3.13) provide an explicit method to carry this out. For example, in the lowest (tree-level) approximation, we have in the "Feynman gauge" in which λ 0 = α 0 = 1, where λ 0 is the conventional gauge-fixing parameter.
Similarly, for the fermion field, in the tree-level approximation, Besides these "quantum propagators," we also have heat kernels, Eqs. We now explain a diagrammatic representation of the perturbative expansion of flowed fields (the flow Feynman diagram). For quantum propagators (3.17) and (3.18), we use the standard convention that the free propagator of the gauge boson is denoted by a wavy line and the free propagator of the fermion is denoted by an arrowed straight line. We stick to these conventions because these are quite natural in a system containing fermions. In Refs. [7] and [9], on the other hand, the arrowed straight line was adopted to represent the gauge boson heat kernel (3.6). Since we already used this in this paper for the fermion propagator, we instead use "doubled lines" to represent heat kernels (3.6) and (3.14). For instance, Figs. 1-5 are one-loop flow Feynman diagrams which contribute to the twopoint function of the flowed fermion field. In these figures, the doubled straight line represents the fermion heat kernel (3.14); the arrow denotes the flow of the fermion number, not the direction of the flow time. Similarly, in Fig. 15, the doubled wavy line is the gauge boson heat kernel (3.6). In the present representation, we thus lose the information of the direction of the flow time, which is represented by an arrow in Refs. [7] and [9]. This information, however, can readily be traced back.
Another element of the flow Feynman diagram is the vertex. The vertices that come from the original action (2.

Ringed fermion fields
A salient feature of the flowed fields is the UV finiteness: Any correlation functions of the flowed gauge field B μ (t, x) with strictly positive t are, when expressed in terms of renormalized parameters, UV finite without the multiplicative (wave function) renormalization [9]. Moreover, this finiteness persists even for the equal-point limit. Thus, any correlation functions of any local products of B μ (t, x) (with t > 0) are UV finite without further renormalization other than the parameter renormalization. In other words, although those local products are given by a certain combination of the bare gauge field A μ (x) through the flow equation, they are renormalized finite quantities. A basic reason for this UV finiteness is that the propagator of the flowed gauge field (3.17) contains the Gaussian dumping factor ∼ e −t p 2 which effectively provides a UV cutoff for t > 0. To prove the above finiteness, however, one has to also utilize a Becchi-Rouet-Stora (BRS) symmetry underlying the present system that is inhomogeneous with respect to the gauge potential [9]. Regrettably, the above finiteness in the first sense does not hold for the flowed fermion field. It requires the wave function renormalization. Although its propagator (3.18) also possesses the dumping factor ∼ e −t p 2 , the BRS transformation is homogeneous on the fermion field (and on general matter fields) and the finiteness proof in Ref. [9] does not apply. In fact, computation of the one-loop diagrams in Figs. 1-5 (and diagrams with opposite arrows) shows that the wave function renormalization makes correlation functions UV finite [10]. Finiteness in the above second sense still holds: Any correlation functions of any local products of χ R (t, x) andχ R (t, x) remain UV finite [10]. Although the finiteness in the second sense is quite useful for our purpose, we still need to incorporate the renormalization factor Z χ in Eq. (3.19). This introduces a complication to our problem, because we have to find a matching factor between Z χ in the dimensional regularization and that in the lattice regularization.
One possible way to avoid this complication is to normalize the fermion fields by the vacuum expectation value of the fermion kinetic operator: 10 where so that the multiplicative renormalization factor Z χ is cancelled out in the new ringed variables. Note that the mass dimension of The vacuum expectation value of the kinetic operator in the lowest (one-loop) order approximation is given by diagram D01 in Fig. 6. For D = 4 − 2 dimensions, we have  Table 1. In total, we have Using Eq. (A1) for the normalization factor in Eqs. (3.20) and (3.21), we have

Energy-momentum tensor constructed from the flowed fields 4.1. Small flow-time expansion and the renormalization group argument
To express the energy-momentum tensor in Eqs. (2.7) and (2.10) in terms of the flowed fields, we consider following second-rank symmetric tensors (which are even under the CP transformation) constructed from the flowed fields: Note that all the above operatorsÕ iμν (t, x) are of dimension 4 for any D. 10 We also introduce corresponding bare operators in the D-dimensional x-space: We now consider the situation in which the flow time t in Eqs. (4.1)-(4.5) is very small. Since a flowed field at position x with a flow time t is a combination of un-flowed fields in the vicinity of x of radius ∼ √ 8t, 11 the operators (4.1)-(4.5) can be regarded as local operators in x-space in the t → 0 limit. Since the bare operators in Eqs. (4.6)-(4.10) span a complete set of symmetric second-rank gauge-invariant local operators of dimension 4 (for D → 4) which are even under CP, we have the following (asymptotic) expansion for small t: where the abbreviated terms are contributions of operators of mass dimension 6 (for D → 4) or higher.
In writing down the small flow-time expansion (4.11), we have assumed that there is no other Ddimensional composite operator at the point x. That the expansion (4.11) cannot necessarily hold when there is another operator at the point x [say, P(x)] can be seen by noting that the product of the left-hand side of Eq. (4.11) with P(x) does not possess any divergence for t > 0, while each term in the right-hand side can have an equal-point singularity with P(x) because the operators in the right-hand side of Eq. (4.11) are D-dimensional (i.e., non-flowed) composite operators. This is a contradiction if Eq. (4.11) holds. This implies that the formula we will derive for the energy-momentum tensor below holds only when the energy-momentum tensor has no overlap with other operators. 12 In particular, we cannot say anything about whether the Ward-Takahashi relation (2.9) is reproduced with our construction. Still, our construction is expected to have a correct normalization because it is determined from matching the energy-momentum tensor in the dimensional regularization which fulfills Eq. (2.9). Our lattice energy-momentum tensor is thus useful to compute correlation functions in which the energy-momentum tensor is separated from other operators. This is the case for correlation functions relevant to the viscosities [40][41][42], for example.
The expansion (4.11) may be inverted as we have the expression (for D = 4) where (4.16) In Eq. (4.14), we have used the fact that the finite operatorÕ 1μν (t, x) − (1/4)Õ 2μν (t, x) is traceless in D = 4 and thus has no vacuum expectation value. Equation (4.14) shows that if one knows the t → 0 behavior of the coefficients c i (t), the energy-momentum tensor can be obtained as the t → 0 limit of the combination in the right-hand side. As already noted, since the composite operators (4.1)-(4.5) constructed from (ringed) flowed fields should be independent of the regularization adopted, one may use the lattice regularization to compute correlation functions of the quantity in the righthand side of Eq. (4.14). This provides a possible method to compute correlation functions of the correctly normalized conserved energy-momentum tensor with the lattice regularization. Thus, we are interested in the t → 0 behavior of the coefficients c i (t) in Eq. (4.14). Quite interestingly, one can argue that the coefficients c i (t) can be evaluated by the perturbation theory for t → 0 thanks to the asymptotic freedom. To see this, we apply on both sides of Eq. (4.14), where μ is the renormalization scale and the subscript 0 implies that the derivative is taken while all bare quantities are kept fixed. Since the energy-momentum tensor is not multiplicatively renormalized, (μ∂/∂μ) 0 (left-hand side of Eq. (4.14)) = 0. On the right-hand 12 These observations imply Then the standard renormalization group argument says that c 1,2,3,4 (t) and mc 5 (t) are independent of the renormalization scale, if the renormalized parameters in these quantities are replaced by running parameters defined by where μ is the original renormalization scale. Thus, since c 1,2,3,4 (t) and mc 5 (t) are independent of the renormalization scale, two possible choices, q = μ and q = 1/ √ 8t, should give an identical result. In this way, we infer that where we have explicitly written the dependence of c i (t) on renormalized parameters and on the renormalization scale. Finally, since the running gauge couplingḡ(1/ √ 8t) → 0 for t → 0 thanks to the asymptotic freedom, we expect that we can compute c i (t) for t → 0 by using the perturbation theory. Although we are interested in low-energy physics for which the perturbation theory is ineffective, the coefficients c i (t) in Eq. (4.14) for t → 0 can be evaluated by perturbation theory; this might be regarded as a sort of factorization.

c i (t) to the one-loop order
We thus evaluate the above coefficients c i (t) in Eq. (4.14) to the one-loop order approximation. For this, we compute the mixing coefficients ζ i j (t) in Eq.   (1) i j (t) + ζ (2) i j (t) + · · · . (4.25) As Ref. [7], it is straightforward to compute ζ (1) i j (t). For example, to obtain ζ (1) i j (t) with j = 1 and 2, we consider the correlation function In the Feynman gauge which we use throughout the present paper, this has the structure βγ (k, l). (4.27) After expanding M μν,βγ (k, l) to O(k, ), we can make use of the following correspondence to read off the operator mixing: In this way, we obtain ζ (1) i j (t) with j = 1 and 2. 13 Similarly, to obtain ζ (1) i j (t) with j = 3, 4, and 5, we consider whose general structure reads We then expand M μν (k, l) to O(k) and O( ) and use the correspondence to read off the operator mixing. A remark on the correspondences in Eqs. (4.28) and (4.31): In our small flow-time expansion (4.11), the operators O jμν (x) in the right-hand side are bare, not renormalized. Such bare operators in the momentum space are represented by tree-level vertices in the left-hand sides of Eqs. (4.28) and (4.31). This is in accord with the definition of the bare operator in the conventional operator renormalization, in which multiplicative renormalization factors to bare operators are chosen so that the insertion of corresponding tree-level vertices with the renormalization factors produces finite correlation functions.
For ζ (1) 1 j (t), diagrams A03, A04, A05, A06, A07, A08, A09, A11, A13, A14, A15, A16, and A19 in Figs. 14-30 contribute. 14 In these and following diagrams, the cross generically represents one of composite operators in Eqs. (4.1)- (4.5) [in the present case,Õ 1μν (t, x)]. Apart from A19, we can borrow the results from Ref. [7] for these diagrams. For completeness, these results are reproduced in the present convention in Table C1 of Appendix C. Combined with the contribution of diagram A19, we have ζ (1) 13 For the momentum integration, we use the integration formulas in Appendix B. 14 Diagrams A10, A12, A17 (where the dotted line represents the ghost propagator), and A18 in these figures correspond to the conventional wave function renormalization and should be omitted in computing the operator mixing.

14/27
Downloaded from https://academic.oup.com/ptep/article-abstract/2014/6/063B02/1561362 by guest on 30 July 2018 Fig. 14. A03   Fig. 15. A04   Fig. 16. A05   Fig. 17. A06  Fig. 18. A07 Fig. 19. A08 (1) For ζ (1) 3 j (t) with j = 1 and 2, diagrams B03, B04, B08, B09, B10, B11, and B12 in Figs. 31-40 contribute. The contribution of each diagram is tabulated in Table 2. For ζ (1) 3 j (t) with j = 3, 4, and 5, diagrams B06, B07, B13, B14, B15, B16, B17, and B18 in Figs. 34-46 contribute and their contributions are tabulated in Table 3. As the sum of these contributions, we have ζ (1)   (1) 3 j in units of  Table 3. ζ (1) 3 j in units of is given by Eq. (3.27). From these, we further have and ζ (1) 55 (t) is given by the sum of the contributions of one-loop diagrams in Table 4 and the conversion factor to the ringed fields: We have now obtained all ζ (1) i j (t) in Eqs. (4.24) and (4.25). Then, since the matrix ζ i j (t) in the treelevel approximation is a unit matrix, it is straightforward to invert the matrix ζ i j (t) in the one-loop approximation; Eqs. (4.15) and (4.16) thus yield 55) Then, by using the renormalized gauge coupling in the MS scheme (A1), to the one-loop order, we have (for → 0) Since the c i (t) in Eq. (4.14) connect the finite energy-momentum tensor and UV-finite local products O iμν (t, x) constructed from (ringed) flowed fields, they should be UV finite. That our explicit oneloop calculation of c i (t) confirms this finiteness is quite reassuring. If one prefers the MS scheme instead of the MS scheme assumed in above expressions, it suffices to make the replacement where γ E is Euler's constant.

A consistency check: The trace anomaly
It is interesting to see that Eq. (4.14) with c i (t) in Eqs. ( Then the last term precisely fills the difference between 4c 2 (t) and b 0 /2. In this way, to the one-loop order, we have whose use is justified when there is no other operator in the point x, as we are assuming. We observe that our one-loop result in Eqs. (4.60)-(4.64) is consistent with the trace anomaly. In Ref. [7], for the pure Yang-Mills theory, the next-to-leading (two-loop order) term in c 2 (t) was determined as where b 0 = [1/(4π) 2 ](11/3)C 2 (G) and b 1 = [1/(4π) 4 ](34/3)C 2 (G) 2 , by imposing that the expression (4.14) reproduces the trace anomaly (2.11) to the two-loop order. For the present system with fermions, however, it seems that this requirement alone cannot fix the next-to-leading terms in c 2 (t) and in c 4 (t); so we are content with the one-loop formulas, Eqs. (4.60)-(4.64), in the present paper treating a system containing fermions.

Master formula
From Eq. (4.14), the energy-momentum tensor is given by the t → 0 limit, where operators in the right-hand side are given by Eqs. (4.1)-(4.5). One may further use the identities where m ∞ denotes the renormalization group invariant mass. For the massive fermion, m ∞ may be determined by using the method established in Ref. [43], for example. For the massless fermion, m ∞ = 0 and we can simply discard the last line of Eq. (4.70).

Equation of motion in the small flow-time limit
Let us consider the following representations for small flow-time: and δ μν m 0ψ (x)ψ(x) = e 5 (t)Õ 5μν (t, x) + O(t), (5.2) where it is understood that the vacuum expectation values are subtracted on both sides of the equations. By a renormalization group argument identical to that which led to Eqs. This implies that we can make the replacement (the subtraction of the vacuum expectation value is understood) in the master formula (4.70) [for the MS scheme, one makes the substitution (4.77)], because throughout this paper we are assuming that the energy-momentum tensor {T μν } R (x) is separated from other operators in correlation functions. This makes the expression of the energy-momentum tensor somewhat simpler. If one is interested in the trace part of the energy-momentum tensor, that is, the total divergence of the dilatation current, the above procedure leads to which is quite analogous to the trace anomaly (2.11). For the massless fermion,m(1/ √ 8t) = 0 and we end up with a quite simple expression for the trace part of the energy-momentum tensor.

Conclusion
In the present paper, on the basis of the Yang-Mills gradient flow, we constructed a formula (4.70) that provides a possible method to compute correlation functions containing the energy-momentum tensor in lattice gauge theory with fermions. This is a natural generalization of the construction in Ref. [7] for the pure Yang-Mills theory. Although the feasibility of the application in lattice Monte Carlo simulations remains to be carefully investigated, the experience in the thermodynamics of 22/27 Downloaded from https://academic.oup.com/ptep/article-abstract/2014/6/063B02/1561362 by guest on 30 July 2018 the quenched QCD [27] strongly indicates that, even with presently available lattice parameters, there exists a window (1.2) within which one can reliably carry out the extrapolation for t → 0 in Eq. (4.70). We expect various applications of the present formulation. One is the application in manyflavor gauge theories with an infrared fixed point (which are subject to recent active investigations; see contributions in the last lattice conference [44,45] for recent reviews).