Non-perturbative definition of the QCD energy-momentum tensor on the lattice

We present a strategy to define non-perturbatively the energy-momentum tensor in Quantum Chromodynamics (QCD) which satisfies the appropriate Ward identities and has the right trace anomaly. The tensor is defined by regularizing the theory on a lattice, and by fixing its renormalization constants non-perturbatively by suitable Ward identities associated to the Poincaré invariance of the continuum theory. The latter are derived in thermal QCD with a non-zero imaginary chemical potential formulated in a moving reference frame. A renormalization group analysis leads to simple renormalization- group-invariant definitions of the gluonic and fermionic contributions to either the singlet or the non-singlet components of the tensor, and therefore of their form factors among physical states. The lattice discussion focuses on the Wilson discretization of quark fields but the strategy is general. Specific to that case, we also carry out the analysis for the on-shell O(a)-improvement of the energy-momentum tensor. The renormalization and improvement programs profit from the fact that, as shown here, the thermal theory enjoys de-facto automatic O(a)-improvement at finite temperature. The validity of the proposal is scrutinized analytically by a study to 1-loop order in lattice perturbation theory with shifted and twisted (for quarks only) boundary conditions. The latter provides also additional useful insight for a precise non-perturbative calculation of the renormalization constants. The strategy proposed here is accessible to Monte Carlo computations, and in this sense it provides a practical way to define non-perturbatively the energy-momentum tensor in QCD.


Introduction
The energy-momentum tensor, T µν , is a central quantity in a quantum field theory since it groups together the currents associated to the invariance of the theory under spacetime translations, from which also those associated to the larger Poincaré group and scale invariance can be built. Apart from the theoretical one, the great interest in the energymomentum tensor is manifold. For instance, in general relativity it enters the Einstein field equations acting as the source of space-time curvature generated by the fields. In thermal field theories its expectation values provide the equation of state of the theory, while its twopoint correlation functions allow one to measure the transport coefficients of the plasma. The only known non-perturbative regularization of QCD is the lattice where, however, the Poincaré group is explicitly broken into discrete subgroups, and the full symmetry is recovered in the continuum limit. As a consequence, a given definition of the energymomentum tensor on the lattice needs to be properly renormalized to guarantee that the associated charges are the generators of the Poincaré group in the continuum limit, and that the trace anomaly is correctly reproduced.
In order to construct the renormalized energy-momentum tensor, the way to proceed is to impose suitable Ward Identities (WIs) at fixed lattice spacing that hold up to cutoff effects which vanish in the continuum limit. Indeed that problem was addressed in refs. [1,2] for the Yang-Mills theory and QCD, where it was shown that on the lattice the 10dimensional symmetric tensor T µν breaks into the sum of a sextet, a triplet and a singlet representation of the hypercubic group. Each one of those three parts picks up finite renormalization constants which were computed to 1-loop order in perturbation theory [3][4][5], see also [6][7][8]. However, it was not clear how to define the renormalization constants so that they could be computed non-perturbatively. JHEP04(2020)043 JHEP04(2020)043 2 The energy-momentum tensor in the continuum In this section we consider QCD in the continuum, for definitions and conventions see appendices A and B. We are interested in correlation functions of the energy-momentum tensor T µν with gauge-invariant operators inserted at a physical distance. As reviewed in [4,13], it is appropriate in those cases to consider the symmetric gauge-invariant definition of the energy-momentum tensor given by where the first term is the gluonic component 1 while the second is the fermionic one given by where ← → D µ is defined in eq. (B.10).

Thermal QCD in the presence of shift and imaginary chemical potential
For our strategy it is instrumental to consider the thermal theory with a non-zero imaginary chemical potential formulated in a moving reference frame. Its grand canonical partition function reads Z(L 0 , ξ, µ I ) = Tr{e −L 0 ( H−iξ· P −iµ I N ) } , (2.4) where the trace is over all the states of the Hilbert space, L 0 is the finite length of the temporal direction, and N is the quark number operator (three times the baryon number). The Hamiltonian H, the total momentum operator P , and the imaginary chemical potential µ I are expressed in a reference frame where the system is moving at a velocity v which, by analytic continuation, we fix to the imaginary value v = iξ.
In the Euclidean path-integral formalism the partition function (2.4) is given by Z(L 0 , ξ, µ I ) = Z(L 0 , ξ, θ 0 ) θ 0 =−L 0 µ I , (2.5) where Z(L 0 , ξ, θ 0 ) is the ordinary QCD path integral as defined in eq. (B.1) with fields which, in the time direction, satisfy periodic boundary conditions up to a shift L 0 ξ [9][10][11] and a twist of the fermion fields [23]. Specifically, the gauge field A µ obeys while the fermion fields, on top of the usual minus sign, pick up also a non-trivial twist phase at the boundaries so that JHEP04(2020)043 The free-energy density is given by ln Z(L 0 , ξ, θ 0 ) , (2.8) where V = L 1 L 2 L 3 is the spatial volume of the box. In the thermodynamic limit, which is always assumed in this section, the invariance of the dynamics under the SO(4) group implies [11] f (L 0 , ξ, θ 0 ) = f L 0 1 + ξ 2 , 0, θ 0 , (2.9) where the parameter θ 0 remains the same on the two sides of the equation since the conserved quark number charge is a relativistic invariant. The effect of the shift is therefore to lower the physical temperature and the chemical potential of the system of the same factor, i.e. from T = L −1 0 to T = (L 0 1 + ξ 2 ) −1 and from µ I to µ I / 1 + ξ 2 respectively, with respect to the system with the same values of L 0 and µ I but no shift. The entropy density reads where V 0 =ψγ 0 ψ. This formula generalizes the one in ref. [11] to the case of a non-zero chemical potential. The expectation values of the space-time components of the energymomentum tensor dictate the dependence on ξ of the free-energy density to be [9][10][11] T 0k ξ,θ 0 = − ∂ ∂ξ k f (L 0 , ξ, θ 0 ) , (2.11) while its θ 0 -dependence is determined by the average value of the temporal component of the vector current so to satisfy As a result, the significant dependence of T 0k ξ,θ 0 on θ 0 can be written as a relation which turns out to be useful on the lattice.
The expectation values of the traceless diagonal components of the energy-momentum tensor are related to those of the space-time components via the Lorentz transformation discussed in sections 2 and 4 of ref. [11], e.g. (no summation over repeated indices)
(2. 16) JHEP04(2020)043 By combining eq. (2.13) with eqs. (2.14)-(2.16), the θ 0 -dependence of the expectation values of the diagonal components of T µν can be related to the one of the vector current as well.

Finiteness of T µν and trace anomaly
In dimensional regularization the energy-momentum tensor (2.1) is decomposed as The singlet operator isτ =τ G +τ F , and up to terms that vanish by the equation of motion of the fermion fields with D = 4 − 2 .

The non-singlet
The dimension-four gauge-invariant fields τ G µν and τ F µν are parity and charge conjugation invariant, and transform as a two-index traceless symmetric irreducible representation of the SO(D) group. No other gauge-invariant field with the same quantum numbers and dimension ≤ 4 can be constructed. Since the derivative of the free-energy density with respect to the shift is finite once the bare parameters of the theory have been renormalized [9][10][11], i.e. eq. (2.11), we may choose the renormalization pattern to be where the superscript R labels renormalized quantities in the desired scheme. The anomalous-dimension matrix takes the form

JHEP04(2020)043
where at leading order in g 2 the coefficient matrix is with [6,24] It is then straightforward to verify that By following refs. [25][26][27], we can define the RGI gluonic and fermionic components of the energy-momentum tensor as where Θ τ (g) satisfies the differential equation and with b 0 being the first coefficient of the β-function defined in eq. (B.16). The solution has the same structure as Z τ , and it can be written as where are two independent RGI fields. Once defined, they allow for an unambiguous scale-and scheme-independent separation of the energy-momentum tensor in two parts, and therefore of its form factors among physical states, which tend to the free gluonic and fermionic contributions when g → 0. At finite temperature, for instance, this allows for an unambiguous split of the entropy density in two parts which tend, in the infinite temperature limit, to the Stefan-Boltzmann values for gluons and fermions respectively.

The singlet
The gluonic and fermionic components ofτ are singlets under SO(D), and therefore they mix with the identity. A natural prescription for subtracting this contribution leads to the renormalization pattern where . . . 0 indicates the expectation value for L 0 → ∞ (zero temperature), and Also for the singlet, we can define the RGI gluonic and fermionic operators as where the 2 × 2 matrix Θτ (g) satisfies the differential equation with the initial condition The anomalous-dimension matrix is given by where at leading order in g 2 , see below, the coefficient matrix is In dimensional regularization, by deriving eq. (2.11) with respect to the renormalized coupling at fixed renormalized quark mass M , see eqs. (2.8), (B.11) and (B.12), we obtain Analogously by deriving with respect to the renormalized quark masses, it follows that we can take (Zτ ) 22 = 1. These values of the renormalization constants correspond exactly to the MS prescription, and therefore the superscript R in eq. (2.38) labels renormalized fields defined in the MS scheme in this case. The anomalous dimension matrix reads By inserting these findings in eq. (2.21), and taking the limit → 0, one obtains the result [28][29][30]τ we can finally rewrite eq. (2.51) for the trace anomaly in terms of RGI operators as

JHEP04(2020)043
Analogously to the non-singlet case, we can writē Considerations analogous to those made for the non-singlet components of the energymomentum tensor apply also in this case.

Finite volume
In a finite volume the theory needs to be further specified by the boundary conditions imposed on the fields in the spatial directions. For the gauge field we supplement eq. (2.6) by standard periodic boundary conditions in the spatial directions, wherek is the unit vector in the direction k. Note that this implies −L k /2 ≤ L 0 ξ k < L k /2. For quark and anti-quark fields the boundary conditions in eqs. (2.7) are supplemented by periodic boundary conditions in the spatial directions up to a non-trivial twist phase θ k [31] ψ(x 0 , x +kL k ) = e iθ k ψ(x 0 , x) , The finite length of the spatial directions and the angles θ k softly break the SO(4) group, and θ k breaks charge conjugation as well. As a consequence eqs. (2.14)-(2.16) are modified in a finite volume. For instance, eq. (2.14) needs to be replaced by eq. (4.25) of ref. [11]. In the latter there is an extra term proportional to the expectation value of T 0k computed in the reference frame where the velocity in direction k is null while the other velocity components are unchanged. For this extra term to be null, the space-time geometry must satisfy the condition [11] L k ξ k where b µ = θ µ /L µ . Similarly eq. (2.15) picks up two extra terms: one proportional to the expectation value of T 0k as before, and a second one proportional to the expectation value of (T kk − T jj ). For both these terms to be null in the reference frame where the velocity JHEP04(2020)043 in direction k is null while the other components are unchanged, the condition (2.59) and In the simplified case in which the only non-null component of the shift is ξ k , the conditions in eqs. (2.59)-(2.61) can be summarized as It must be said, however, that in the presence of a mass gap the extra terms acquired by eqs. (2.14)-(2.16) vanish exponentially with the length of the spatial directions, and therefore they are expected to be negligible in large enough volumes.

The energy-momentum tensor on the lattice
A non-perturbative definition of the theory presented in section 2 is achieved by introducing a four-dimensional Euclidean hypercubic lattice of spacing a which acts as an ultraviolet regulator. The gauge potential A µ is replaced as usual by the SU(N c )-valued gauge fields U µ residing on the links of the lattice, while the quark and anti-quark fields ψ, ψ are defined on the sites. 3 Although the non-perturbative renormalization strategy presented in this paper is general, for definiteness we consider the Wilson formulation for gluons and quarks, see appendix C for the details. A discretization of the energy-momentum tensor is obtained by replacing the fields and the derivatives appearing in the continuum expressions (2.1)-(2.3) with their lattice counterparts. In particular we define [1][2][3][4] with the gluonic component given by and with F µν being the clover discretization of the field strength tensor in eq. (C.9) which, at variance with F µν , is not traceless. For the fermionic part we take

JHEP04(2020)043
where ← → ∇ µ and ← → ∇ * µ are defined in eq. (C.12). In the naive continuum limit a → 0, these expressions tend to the continuum ones in eqs. (2.2) and (2.3). As any other discretization of T µν , however, they need to be properly renormalized to guarantee that their correlation functions satisfy the continuum Ward identities, up to discretization errors which vanish in the continuum limit.
The target energy-momentum tensor in the continuum is a gauge-invariant operator of dimension 4, which is a combination of a traceless two-index symmetric and a singlet irreducible representations of SO(4) even under parity and charge conjugation. Since on the lattice the SO(4) symmetry reduces to the hypercubic group SW 4 , the traceless twoindex symmetric representation splits into a sextet (non-diagonal components) and a triplet (diagonal traceless components). At finite lattice spacing, the energy-momentum tensor is thus a combination of gauge-invariant operators of dimension d ≤ 4 which, under the hypercubic group, transform as one of those two representations and the singlet [1]. In QCD there are seven such operators plus the identity. We can take as a basis (no summation over the repeated indices µ and ν in eqs. (3.6) and (3.9)) for the purely gluonic fields and for the fermionic ones. 4 Since the hypercubic group is an exact symmetry of the lattice theory, a given field in eqs. (3.5)-(3.11) can mix under renormalization only with those in the same irreducible representation.
In this paper we focus on the definition of the sextet and triplet whose renormalization pattern is given by where the renormalization constants Z {i} G and Z {i} F are finite and depend on g 2 0 only since in the continuum the nonet component of the energy-momentum tensor has vanishing anomalous dimension. In a regularization which breaks chiral symmetry explicitly, the JHEP04(2020)043 renormalization pattern of the singlet component requires a rather different and more involved analysis which fills a different publication by itself. Departure from scale invariance can, however, be extracted directly from the Callan-Symanzik renormalization group equations or, for instance, from the r.h.s. of eq. (2.16). Finally we note that, as anticipated in section 2.2.1, the determination of the renormalized fields 13) and in particular of their RGI counterparts, gives access to the RGI gluonic and fermionic components of the energy-momentum tensor. While the renormalization of T {i} µν is fixed by WIs, that of (3.13) can be obtained, for instance, by imposing suitable conditions on its expectation values in the presence of shifted and twisted boundary conditions. A detail investigation of this problem is in progress.

Non-perturbative renormalization conditions
The renormalization constants Z can be determined non-perturbatively by enforcing on the lattice the relations (2.11) and (2.14) up to discretization effects which vanish in the continuum limit. To this aim, shifted boundary conditions for the links are while eqs. (2.7) and (2.58) fix the boundary conditions for quark and anti-quark fields on the lattice too, where the components of ξ are discretized in integer units of a/L 0 . As very commonly adopted in lattice QCD, we opt for a mass-independent renormalization scheme. The renormalized coupling g R and the renormalized quark mass m R for each given flavour are thus related to the bare parameters as where µ is a renormalization scale, m q = m 0 − m c (g 2 0 ) is the subtracted mass, and m c (g 2 0 ) is the critical mass.
The two renormalization constants of the sextet are fixed by requiring that, for two different set of values θ A and θ B , where θ = (θ 0 , θ), it holds is a symmetric discrete approximation of the derivative of the free energy with respect to the k-th component of the shift. From a practical point of view, it is useful to combine eq. (3.16) for one set of values of the angles, e.g. θ = θ A , with the lattice realization of the identity (2.13)  20) and the expectation value of the renormalized current on the r.h.s. of eq. (3.18) corresponds to the bare one. 5 Once the sextet renormalization constants have been fixed, those of the triplet can be determined by enforcing on the lattice the analogous of eqs. (2.14) or (2.15) for two different set of values θ = θ A , θ B of the twist angles. In a finite volume these relations are satisfied only up to exponentially small finite-size effects, unless lattice sizes and twist phases according to the constraints given in eqs. (2.59)-(2.61) are considered. The above renormalization conditions are imposed at zero quark masses, i.e. all bare masses are set to the critical value m 0 = m cr (g 0 ). In practice this is possible thanks to the presence of a spectral gap in the lattice Dirac operator at finite temperature. Finally, it is important to emphasize that the above renormalization conditions are valid non-perturbatively, and are designed to be accessible to numerical Monte Carlo computations. The eq. (3.16) can be studied numerically with a strategy analogous to the one already successfully implemented for the Yang-Mills theory [13], while eqs. (3.18) and (3.21) require standard numerical techniques. In this sense these conditions provide a practical strategy to define non-perturbatively the energy-momentum tensor in QCD.

O(a)-improvement
The Symanzik improvement programme has the purpose of accelerating the approach to the continuum limit of field correlators. It is achieved by adding suitable counterterms to the lattice action and to the fields multiplied by numerical coefficients which are properly adjusted so to cancel discretization errors order by order in the lattice spacing [32,33]. In the following we discuss how to implement this programme to O(a)-improve on-shell matrix elements of the sextet and triplet components of the energy-momentum tensor. For the clarity of the presentation, we discuss separately the massless, mass-degenerate, and mass non-degenerate cases. 5 It is interesting to notice that the renormalization constant and the improvement coefficients of the flavour-singlet lattice vector currents can be determined by comparing V c 0 ξ,θ with the analogous one for the local current and similarly for higher-point correlation functions.

Massless quarks
The first step consists in identifying complete bases of dimension-5 gauge-invariant fields which, in the Symanzik effective continuum theory, are parity and charge conjugation invariant and transform as sextets and triplets under the SO(4) hypercubic subgroup. When all quarks are massless, i.e. m 0 = m c (g 0 ), there are no operators made of the gauge field only. Complete bases of O(a)-counterterms are built by projecting the fields on their sextet and triplet components. Terms proportional to the second and third fields in (4.1) do not contribute to matrix elements between initial and final states with the same four-momentum. Since we restrict our analysis to those cases, we can discard them. The improved sextet and triplet fields are then obtained by replacing in eqs. (3.8) and (3.9) respectively, with (no summation over µ, ν in (4.4)) By performing a classical expansion of T

F,{i} µν
in the lattice spacing [34], it turns out that the tree-level values of the coefficients c {i} F are null as well as those multiplying the other counteterms in eq. (4.1).

Mass-degenerate quarks
When all quarks have the same non-vanishing mass, the renormalization of the coupling and of the mass can be obtained from eqs. (3.15) by replacing where m q is the subtracted bare mass common to all flavours. The two improvement coefficients b g and b m must be chosen so to have a mass-independent renormalization scheme where the renormalized coupling and the mass are free from O(a)-effects [31]. Defined this way, both coefficients do not depend on the renormalization conditions chosen to set Z g and Z m . The perturbative expansion of b g starts at O(g 2 0 ) since it arises from sea-quark loop contributions to a gluonic quantity [35].
To improve the sextet and triplet parts of T µν , two more O(a)-counterterms, made of the original gluon and fermion components multiplied by the quark mass, have to be taken into account. As a result, the O(a)-improved fields read 6

JHEP04(2020)043
Notice that the renormalization constants appearing in eq. (4.6) must be evaluated at the value g 2 0 . On the other hand, it is consistent to evaluate c

Mass non-degenerate quarks
When all quark masses are different, the pattern of improvement is further complicated. An O(a)-improved renormalized coupling can be defined from eqs. (3.15) by replacing [36] while the expressions for the renormalized improved quark masses are quite involved and, since they are not needed in the following, we refer the interested reader directly to eq. (26) of ref. [36]. The O(a)-improved fields turn out to be and analogously for T

Improvement conditions
As we have seen in the previous subsections, improving the energy-momentum tensor may require the tuning of several parameters. A decoupling of the equations that fix them, however, occurs naturally within our strategy. Indeed, in the massless limit chiral symmetry is expected to be either exact or effectively restored when the temperature T is much larger than the typical scale of the strong interactions, a few hundred MeV or so. As a consequence, the expectation values of the chiral non-singlet counterterms in eq. (4.1) either vanish or become quickly negligible at high temperature. The origin of this result can be traced back to the more general fact that the thermal theory with massless quarks enjoys de-facto automatic O(a)-improvement at high temperature, see appendix D for a detailed discussion in the presence of a generic number of flavours.

JHEP04(2020)043
For the particular case we are concerned here, namely the expectation values of the energy-momentum tensor, stronger results can be proven. 7 In the Symanzik effective continuum theory with two or more massless flavours, it holds where ∂R is the union of the top and bottom lids of an hyper-cylinder R containing the origin, A a k =ψγ k γ 5 T a ψ (k = 1, 2, 3) with T a (a = 1, . . . , N 2 f −1) being the generators of the group SU(N f ), and δT a,F,{i} µν is defined as δT F,{i} µν but with the replacement ψ → γ 5 T a ψ. At large temperature, where the theory has a mass gap, the integrand in eq. (4.10) decreases exponentially with the distance |x|. Its integral is therefore null exactly since the lids can be sent to infinity, and the conclusions of appendix D are recovered. At smaller temperature, where the theory develops Goldstone bosons [37], the integrand decreases power-like in |x|. At large distances, the leading behaviour is dictated by the single Goldstone-boson contribution so that where ∆(x) is the free propagator of a massless boson while the dots indicate sub-leading corrections. At distances |x| much larger than the inverse temperature, ∆(x) ∝ 1/|x|, the correlator (4.11) decreases as |x| −4 , and the surface integral in (4.10) is again null.
The final important outcome of this analysis is that, if the lattice action is O(a)-improved and quarks are massless, the expectation values T The improvement coefficients in eq. (4.8) can be determined non-perturbatively by imposing the very same equations (3.16), (3.21) or (3.22) to be valid up to O(a) terms for several values of the quark masses, and by remembering that the free-energy density is already improved once the Sheikholeslami-Wohlert term has been included in the action. A detailed implementation of this strategy to 1-loop order in perturbation theory is discussed in section 5.2.3.

Perturbative analysis
In order to verify analytically the validity of the strategy proposed in this paper, we have computed the free-energy density and the expectation values of the energy-momentum tensor components to 1-loop order in lattice perturbation theory in the presence of shifted and twisted boundary conditions in the infinite spatial volume limit. The calculation has been carried out by regularizing gluons with the Wilson plaquette action and quarks with the O(a)-improved Wilson operator, see appendices C and E for the definitions of the actions, free propagators, and lattice vertices. This computation serves also to determine, for the first time, the renormalization constants of the sextet and triplet components of T µν 7 We thank Martin Lüscher for suggesting to us this line of argumentation.

JHEP04(2020)043
in the O(a)-improved theory to 1-loop order in perturbation theory, as well as their O(a)improvement coefficients. As a byproduct we have confirmed the 1-loop expressions of the renormalization constants in the unimproved theory which were determined in refs. [4][5][6][7][8].
Finally the combination of the results in this section with eqs. (3.16), (3.21) or (3.22) leads to perturbatively improved versions of these non-perturbative renormalization conditions. As a consequence, the numerical non-perturbative determinations of the renormalization constants are free from discretization effects up to order g 2 0 .

Free-energy density
The 1-loop expansion of the bare free-energy density defined in eq. (2.8) is where and are the tree-level and 1-loop contributions respectively. 8 The functions f G(0) and f F (0) are the tree-level gluonic and fermionic contributions, f G(1,Nc) and f G(1, 1 Nc ) are the 1-loop gluonic parts, and f F (1,N f ) collects the 1-loop fermionic contributions. All these functions are reported in appendix F. In the perturbative computations presented in this paper we assume always to have N f quarks with equal masses and twist angles. The formulas for the generic case, however, can be easily obtained by summing the contributions of each individual flavour rather than multiplying the single-fermion contribution by N f , e.g for the free-energy density the terms N f f F (0) and N f f F (1,N f ) must be replaced by the sums over the flavours of the f F (0) and f F (1,N f ) functions computed for the mass and the twist angles of each single flavour respectively.
Once the bare free-energy density has been calculated, the 1-loop renormalized expression is obtained by re-writing the bare parameters g 0 and the common bare quark mass m 0 through the renormalized counterparts defined in eq. (3.15). By properly combining eq. (2.8), the renormalized 1-loop expression at finite lattice spacing and its continuum limit, a 1-loop perturbative improved definition of the free-energy density can also be obtained.

Energy-momentum tensor
The bare expectation values of the sextet and triplet components of T µν are where the tree-level values are

JHEP04(2020)043
while the 1-loop contributions are All functions on the r.h.s. of eqs. (5.5)-(5.6) are given in appendices H and I.
Once the bare parameters of the theory have been renormalized, the definition of the sextet and the triplet components of T µν require the calculation of the renormalization constants defined in eq. (3.12). At one loop they can be written as where we can define To impose the renormalization conditions in eqs. (3.16), (3.21) or (3.22) in the massless limit m R = 0, we remind that at 1-loop the critical mass is where m where the discrete derivative ∆ with respect to the shift is defined as in eq. (3.17). As expected, in the limit L 0 /a → ∞ it holds Z , depend on the interaction between quarks and gluons. If we define the combination where the renormalization constant of the fermion component is have been computed numerically on lattices with temporal extension L 0 /a ranging from 4 to 32 in steps of 2 and spatial size L 1 = L 2 = L 3 = RL 0 with R = 6, 8, 10 and 12. The calculations have been carried out for two values of the shift, ξ = (1, 0, 0) and (1/2, 1/2, 0), and for three values of the fermionic phase in the temporal direction, θ 0 = 0, π/16 and π/4. At fixed value of the shift, data for the three values of θ 0 have been analyzed by considering two independent differences. This large amount of data allowed us to extrapolate the results to infinite spatial volume and to a/L 0 → 0 limit with confidence. The latter extrapolation has been performed by fitting the results in powers of (a/L 0 ) 2 supplemented by terms multiplied by ln (a/L 0 ) for the 1-loop coefficients. The required lattice sums have been computed in coordinate space [39] after having used the Fast Fourier Transform algorithm for computing the gluon and quark propagators. A hierarchical procedure for sums has been implemented JHEP04(2020)043 This work Ref. [4][5][6][7][8] This work Ref. [4][5][6][7][8]   in the limit a/L 0 → 0 are listed in table 1 together with the analogous values present in the literature [4][5][6][7][8]. Our error bars have been estimated by changing the fit range in (a/L 0 ), by considering the spread over the two differences in θ 0 , and by analyzing the dependence on R. For all values that can be compared with refs. [4][5][6][7][8], the agreement is excellent.

Renormalization constants of the triplet
where the renormalization constant of the fermion component is F have been computed on lattices with temporal extension L 0 /a ranging from 4 to 32 in steps of 2 and spatial size L 1 = L 2 = L 3 = RL 0 with R = 10 and 15. The calculations have been carried out for ξ = (2, 0, 0) so to satisfy the constraint in eq. (2.59). Three values of the fermionic phase, θ 0 = 0, π/16 and π/4, have been considered and a non vanishing phase along the direction1 has been chosen according to the third constraint in eq. (2.62). Analogously to the sextet case, data for three values of θ 0 have been analyzed by considering two possible independent differences. The results show no relevant dependence from the spatial volume. The a/L 0 → 0 extrapolation has been performed again by fitting the results in powers of (a/L 0 ) 2 supplemented by terms multiplied by ln (a/L 0 ) for the 1-loop coefficients. As for the sextet, the lattice sums have been computed in coordinate space, and also in this case a hierarchical procedure was implemented.
The final results for the various contributions to Z {3} G and Z {3} F in the limit a/L 0 → 0 are listed in table 1 together with the analogous ones present in the literature [4][5][6][7][8]. The error bars have been estimated by changing the fit range in (a/L 0 ), and by considering the spread over the two differences in θ 0 . Whenever a comparison with refs. [4][5][6][7][8] is possible, the agreement is excellent.

Improvement coefficients
We conclude the perturbative analysis of the strategy proposed in this paper by computing to 1-loop order the improvement coefficients of the sextet and triplet components of T µν introduced in eq. (4.8). The terms δT  To this aim, by using eqs. (3.15) and (4.5), we remind that the O(a)-improved quark mass is given by where at 1-loop in perturbation theory [40] Z m (g 2 0 , aµ) = 1 + (10) . (1,Nc) m The free-energy density and the expectation values of the energy-momentum tensor at bare mass m 0 (g 0 ) are obtained by evaluating the expressions (5.1)-(5.6) at m 0 = m 0 and by replacing  The solutions of the eqs. (3.16), (3.21) or (3.22) at finite quark mass are obtained by replacing that we have obtained in the limit a/L 0 → 0 by numerical computations and analyses analogous to those carried out in the previous two subsections.

Conclusions
Shifted boundary conditions in the presence of an imaginary chemical potential offer an extremely powerful tool to non-perturbatively renormalize composite operators on the lattice. In this work we have applied this framework to the case of the energy-momentum tensor. The strategy proposed here is the natural extension of the one already applied successfully to the SU(3) Yang-Mills theory [13]. The inclusion of quarks, however, complicates the problem because the gluonic and fermionic parts of the tensor mix together. Introducing a non-zero imaginary chemical potential gives the handle to solve that problem since, via a conserved charge, it couples differently to quark and gluons, in particular directly to quarks but only indirectly to gluons through their interaction with quarks. Ward identities can thus be written, both for the sextet and the triplet, which are different enough to resolve the mixing between the gluonic and the fermionic parts so that the computation becomes feasible non-perturbatively.
In view of that application and in order to check the whole construction, we have applied the method in lattice perturbation theory by computing the renormalization constants and the O(a)-improvement coefficients of the sextet and triplet components of the energymomentum tensor to 1-loop order. The agreement with the results in the literature for the unimproved theory represents a very non-trivial test of the entire strategy proposed in this paper. A further confirmation of theory expectations is the ξ and θ independence of the JHEP04(2020)043 renormalization constants once extrapolated to the a/L 0 → 0 limit. An important byproduct of these computations is the possibility of defining 1-loop perturbative improved estimators of the renormalization constants and of the expectation values of the energy-momentum tensor, so to reduce discretization effects in their non-perturbative determinations.
Once fixed so to satisfy the WIs, the renormalization constants are part of the definition of the energy-momentum field itself. Up to discretization errors, they do not depend on the particular correlator or kinematic conditions they were employed to fix them. They can be directly used to renormalize the energy-momentum tensor inserted at a physical distance from other fields in any correlator of QCD, e.g at zero or non-zero temperature with or without chemical potential.

A Conventions and useful identities
Here we summarize the conventions for the generators of the SU(N c ) group and for the Dirac matrices, γ µ , together with some standard identities that we have used in the 1loop perturbative computation. Let T a , a = 1, . . . , (N 2 c − 1), be the hermitean traceless generators of the group SU(N c ) normalized as Their commutation and anti-commutation relations are 10

JHEP04(2020)043
where f abc is the completely antisymmetric tensor of the structure constants while d abc is completely symmetric. It then holds Other useful identities are as well as d aae = 0 and f ace d bce = 0. By defining γ 5 = γ 0 γ 1 γ 2 γ 3 , the Euclidean anti-commutation relations Useful trace identities are and

B Continuum theory
In the Euclidean space-time, the path integral of QCD is defined as where the integration measures on the various fields are defined as usual. The action is defined as 11 11 Throughout the paper we assume the strong CP violation term to be absent.
where g 0 is the bare coupling constant, λ 0 is the gauge-fixing parameter, the trace is over the color index and The quark and anti-quark fields, ψ and ψ have N f -flavour components ψ f , ψ f , f = 1, . . . , N f and, accordingly, the mass matrix M 0 = diag(m 0,1 , m 0,2 , m 0,3 , . . .) is a N f × N f matrix, whose entries on the diagonal are the bare quark masses. It turns out to be useful also to define

B.1 Dimensional regularization
Here we report the essential formulas of dimensional regularization which are needed in this paper. We follow the conventions of ref. [25], see also ref. [41] for a recent review. By replacing d 4 x → d D x, one defines the renormalized coupling g and quark mass matrix M as where D = 4 − 2 and µ is a mass parameter. A generic renormalization constant, including those of composite operators, is expanded in g 2 . In the MS scheme is then implicitly fixed by requiring it to be a polynomial in 1/ with no constant term The β-function of the theory is defined as where

JHEP04(2020)043
with the first coefficient given by The anomalous dimension of the quark mass is defined as where the first coefficient is

C Lattice theory
The action of the lattice theory reads where S G and S F are the gluonic and fermionic contributions respectively. For the gluonic one we consider the Wilson plaquette action where g 0 is the bare gauge coupling. The plaquette field is defined by whereμ,ν are unit vectors oriented along the directions µ, ν respectively. The fermionic part reads where M 0 is the bare quark mass matrix, and for D we choose the O(a)-improved Wilson-Dirac operator The first operator on the r.h.s. is the massless Wilson-Dirac operator defined by where ∇ * µ , ∇ µ are covariant lattice derivatives acting on the quark fields as follows

JHEP04(2020)043
The second term is the Sheikholeslami-Wohlert operator defined by [15] D sw ψ(x) = c sw (g 0 ) is the clover discretization of the field strength tensor 12 and (C.10) The coefficient c sw is tuned in order to remove O(a) lattice artifacts generated by the action in on-shell correlation functions [15,31]. The gauge-invariant path integral is It is also useful to define with ∇ µ , ∇ * µ being the lattice covariant derivatives in eq. (C.7), and

D Automatic O(a)-improvement with massless quarks
In the absence of spontaneous chiral symmetry breaking, chirally symmetric correlators of Wilson fermions in the presence of an even number of massless quarks are proven to be automatically O(a)-improved [42,43]. This occurs, for instance, in a finite volume without boundaries [43] or in the thermal theory at high temperature [44]. To extend this result to a generic number of flavours N f > 1 but still small enough to have asymptotic freedom, we consider the discrete axial symmetry S 5 defined as This non-anomalous element of the U(1) A group indeed allows for a simple generalization of the line of argumentation given in refs. [42,43]. In the massless limit, the action of the Symanzik's effective continuum theory reads

JHEP04(2020)043
where S 0 is defined as in eq. (B.2) but with the bare coupling replaced by the renormalized one, and The leading discretization effects in the connected correlation function of a multi-local renormalized field O are then given by the corresponding continuum correlation functions with the insertion either of S 1 or of the O(a)-counterterm δO for the field O, If we restrict ourselves to fields O inv invariant under parity and the chiral symmetry S 5 , such as (the fermionic component of) T µν in eq. (2.1) in the massless limit, the invariance of the measure and of the action S 0 implies when spontaneous chiral symmetry breaking is absent. Moreover δO con = 0 as well because fields with the same quantum numbers as O inv but of one dimension higher are not invariant under the symmetry S 5 since they must have and extra derivative and therefore an extra γ-matrix with respect to O inv , a fact that was essential to prove the automatic O(a)-improvement in the twisted mass QCD regularization [42]. The eq. (D.4) then reads The very same conclusion can be reached by using the R 5 discrete symmetry in ref. [42] for N f = 2, and an element of the axial subgroup Z N f of the non-Abelian chiral symmetry group for larger N f .

E Propagators and vertices for perturbation theory
On the lattice, perturbation theory is normally set-up in terms of algebra-valued fields A µ (x) defined as where we opt for a minus sign in the exponential so to recover, in the naive continuum limit, the widely used conventions in appendix B. By inserting (E.1) in the expressions in appendix C and by expanding to the appropriate order, the analogous continuum formulas given in appendix B are recovered after the usual field rescaling A µ (x) → A µ (x)/g 0 . The free theory in recovered in the limit g 0 → 0.
In the presence of the boundary conditions (3.14), the Fourier transform can be written as 13 A

JHEP04(2020)043
where, for a generic function f (p), the finite-volume integration is defined to be and n µ = 0, . . . , L µ /a − 1. For lattice fermions satisfying the boundary conditions in eqs. (2.7) and (2.58), the Fourier transform can be written as where the integration p ξ,θ is defined as in eq. (E.3) but for the set of lattice momenta The term π/L 0 in p 0 is due to the anti-periodicity of the boundary conditions along the temporal direction. When L k → ∞ for all three spatial directions k, the integration in eq. (E.3) becomes where BZ stands for the Brillouin zone. The analogous holds for the fermion integrals (E.5) which become independent on θ k . If also L 0 → ∞, the sum over n 0 in eq. (E.7) is replaced by the integral as well, which becomes independent on the shift ξ and the twist θ 0 . As expected by general quantum field theory arguments, the expressions of the propagators and of the vertices of the theory with shifted and twisted boundary conditions are equal to those valid for periodic boundary conditions provided the definition of the momenta are replaced by those in eqs. (E.4) and (E.6), i.e. ξ and θ enter the values of the allowed momenta only. For consistency and to define our conventions, however, we report their definitions in the rest of this appendix.

E.1 Gluon propagator
By adding the gauge-fixing contribution 14 to the gluonic action (C.2), the free-gluon propagator in Feynman gauge reads (E.10) 14 As usual the backward lattice derivative ∂ * µ is defined as eq. (C.7) but with the link omitted.

E.2 Ghost propagator
By following the usual Faddeev-Popov (FP) procedure, we add to the gluonic action the contribution coming from the FP determinant and Φ(A µ ) is a matrix in the adjoint representation of the algebra of SU(N c ) whose matrix elements are [Φ(A µ )] ab = iag 0 f abc A c µ . The ghost propagator then reads Notice that the ghost and anti-ghost fields satisfy the same boundary conditions of the gauge field.

E.3 Fermion propagator
The free-fermion propagator for a single flavour is given by

E.4 Gluonic interaction
The perturbative expansion of the Wilson action (C.2) can be written as where S G,0 is the tree-level gluonic action, while the contributions from the three-and the four-gluon vertices are · X abcd µνλρ (k, q, r, s) + Y abcd µνλρ (k, q, r, s) , The Jacobian from the Haar integration measure due to the change of variables (E.1) can be recast in the form of an extra contribution to the action which reads

E.5 Ghost-gluon interaction
The expansion of the FP action (E.11) reads S F P = S F P,0 + g 0 S F P,1 + g 2 0 S F P,2 + O(g 3 0 ) (E. 24) where S F P,0 is the tree-level term, while

JHEP04(2020)043
where S F,0 is the tree-level contribution, while and The S SW ,2 has two quark and two gluonic lines, but it does not contribute to the quantities we are interested in this paper due to its color structure.

F The free energy density
In this appendix we report the coefficients of the 1-loop perturbative expansion of the freeenergy density defined in eqs.
where F (8) is defined in appendix J. The 1-loop contributions result from connected diagrams with no external legs. The gluonic ones are from which we have

JHEP04(2020)043
where The fermionic contribution is and r = p + k.

F.1 O(a)-improved action
The insertion of the improvement term modifies only f F (1,N f ) so that The expressions of F F 1 and F F 2 in terms of the integrals defined in appendix J are and . (F.14) As an explicit check of the whole computation, we compared the derivative of the free-energy density with respect to θ µ to the expectation value of the µ-component of the conserved vector current (3.19), e.g. for µ = 0 this corresponds to verify the lattice analog of eq. (2.12). The required 1-loop computation of the expectation value of the vector current is reported in the following appendix.

JHEP04(2020)043 G Expectation value of the vector current
By expanding the conserved current in eq. (3.19) to order g 2 its expectation value at one loop can be written as The tree-level value is given by (no summation over µ) and its derivative with respect to the bare mass is (7) µ . (G.4) The 1-loop contribution is The expressions of V 1 µ , V 2 µ , and V 3 µ in terms of the integrals defined in appendix J are × ar µ m 0 (k)m 0 (p) −k µpµ + c µ (r) m 0 (p)k µ + m 0 (k)p µ , (G. 10) and where we have defined r = p + k.

G.1 O(a)-improved action
The Sheikholeslami-Wohlert terms leads to 3 additional terms to the 1-loop coefficient where The expressions of V 4 µ , V 5 µ , and V 6 µ in terms of the integrals defined in appendix J are

H Expectation values of sextet components of T µν
In this appendix we report the coefficients of the 1-loop perturbative expansion of the expectation values of the sextet components of the energy-momentum tensor defined in eqs. (5.4)-(5.6). By expanding the field strength F µν to order g 2 the tree-level value of the gluonic part is given by The 1-loop contributions to the gluonic component are (H.6) The expressions of T G1 µν , . . . , T G7 µν in terms of the integrals defined in appendix J are µµα JHEP04(2020)043 (H.14) By expanding the fermion part of the energy-momentum tensor in eq.

H.1 O(a)-improved action
The 1-loop gluonic contribution to the off-diagonal components of the energy-momentum due to the improvement term is The contribution due to the improvement term to the fermionic part of the energymomentum tensor is made up of 3 terms respectively.

J Integrals
In this appendix we report the definitions of the tree-level integrals which appear in the appendices F, G, H, and I. The functions c µ (p), s µ (p), D G (p), D F (p) and p ξ,θ are given in appendix E. Repeated indices are not summed over. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.