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 space-time translations, from which also those associated to the larger Poincaré group and scale invariance can be built. Apart the theoretical one, the great interest in the energy-momentum 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 two-point 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 energy-momentum 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 10-dimensional 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.
An important step forward was made a few years ago by noticing that useful WIs to fix the renormalization constants are obtained by considering the theory in a finite box, where the Euclidean Lorentz symmetry is also softly broken by its shape [9][10][11]. In particular, if the length in one (temporal) direction L 0 is chosen to be finite (thermal theory), interesting WIs follows. They become particularly simple and of practical use if the periodicity axes are tilted with respect to the lattice grid (moving reference frame), i.e. if the hard breaking of the Poincaré symmetry due to the lattice discretization and the soft one due to the finite temporal direction are not aligned. This set-up has a natural implementation in the Euclidean path-integral formulation in terms of shifted boundary conditions [9][10][11][12]. These ideas led for the first time to a non-perturbative definition of T µν in the SU(3) Yang-Mills theory [13].
The purpose of this paper is to present the WIs for defining non-perturbatively the energy-momentum tensor in QCD. By generalizing what has been done in the SU(3) Yang-Mills theory, we work in the framework of shifted boundary conditions supplemented by the presence of an imaginary chemical potential. It is the latter that gives us the handle to solve the problem of the mixing between the gluonic and fermionic parts of the tensor since the chemical potential couples differently to quark and gluons. Most interestingly, the derived relations can be used in practice to carry out the non-perturbative numerical computation of the renormalization constants of T µν [14].
Although the strategy is general, for definiteness we consider the Wilson formulation of quarks on the lattice with and without the Sheikholeslami-Wohlert improving term [15]. In the presence of the latter, we also carry out the analysis for the onshell O(a)-improvement of the energy-momentum tensor field. The implementation of the renormalization and of the improvement program turn out to be greatly simplified because the thermal theory enjoys de-facto automatic O(a)-improvement at finite temperature.
The validity of the proposal is scrutinized analytically to 1-loop order in lattice perturbation theory with shifted and twisted (for quarks only) boundary conditions. The results obtained this way serve also to give a 1-loop improved definition of the nonperturbative renormalization constants with the aim of reducing discretization effects. For completeness, we notice that over the last few years alternative methods, based on the Yang-Mills gradient flow [16], have also been explored for renormalizing nonperturbatively the energy-momentum tensor [17][18][19][20][21].
The paper is organized as follows. In section 2 we derive the basic identities which are later enforced on the lattice to fix the renormalization constants of the energymomentum tensor, and we construct the renormalization-group-invariant (RGI) definitions of its gluonic and fermionic components. In the following section we show how these relations define non-perturbatively the energy-momentum tensor on the lattice, while in section 4 we discuss the O(a) improvement. In the following section we investigate the renormalization conditions in perturbation theory, and we compute the renormalization constants and the improvement coefficients of the gluonic and fermionic parts of the energy-momentum tensor to 1-loop order. Section 6 contains our conclusions. Notations, conventions, and technical details are reported in several appendices.

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 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 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 [22]. 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 The free-energy density is given by 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 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] 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 [11], e.g. (no summation over repeated indices) Such relations can be imposed on the lattice to fix their relative normalization. Finally for the trace it holds [13] 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 minimal subtraction (MS) scheme. The anomalous-dimension matrix takes the form where the β-function is defined in Eq. (B.14). At leading order in g 2 , the coefficient matrix is where [6,23] It is then straightforward to verify that By following Refs. [24][25][26], we can define the RGI gluonic and fermionic components of the energy-momentum tensor as  where Θ τ (g) satisfies the differential equation 31) and the initial condition 2 The solution has the same structure as Z τ , and it can be written as From these equations it follows that are two independent RGI fields. Once defined, they allow for an unambiguous scaleand 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  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. By inserting these findings in Eq. (2.21), and taking the limit ǫ → 0, one obtains the result [27][28][29] which is valid to all orders in perturbation theory. Also for the singlet, we can define the RGI gluonic and fermionic operators as  where the 2 × 2 matrix Θτ (g) satisfies the differential equation and with the anomalous-dimension matrix given by By inserting the solution into Eq (2.43), we can finally rewrite Eq. (2.42) for the trace anomaly in terms of the RGI operators as 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 [30] ψ 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] 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 in direction k is null while the other components are unchanged, the condition (2.54) and L k / 1 + ξ 2 k = L j have to be met together with In the simplified case in which the only non-null component of the shift is ξ k , the conditions in Eqs. (2.54)-(2.56) 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 Sect. 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 where 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 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 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 {i} G and Z {i} F 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.53) 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.54)-(2.56) 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 [31,32]. 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.

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 [33], 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 [30]. 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 [34].
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 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 [35] 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. [35]. 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.
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 [36], 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 µν 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 1loop 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 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.
where the derivatives on the r.h.s. of these equations can be found in Eqs. (F.2), (H.17) and (I.6) respectively.

Renormalization constants of the sextet
By imposing Eq. (3.16) for two different values 9 of θ 0 = θ A 0 , θ B 0 , the tree-level values can be defined as 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 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 [38] after having used the Fast Fourier Transform algorithm for computing the gluon and quark propagators. A hierarchical procedure for sums has been implemented in order to preserve a high numerical accuracy. For some of them, however, we needed to run in quadruple precision due to large cancellations taking place at large volumes.
The final results for the various contributions to Z {6} G and Z

{6}
F 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 the triplet
where the renormalization constant of the fermion component is 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.54). 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.57). 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 [39] 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) 0 and by replacing   . In Table 2  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 energy-momentum 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 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.

Acknowledgments
We thank Martin Lüscher for many illuminating discussions, especially on topics in sections 2.2 and 4.4, and for comments on a preliminary version of this paper. The numerical integrals needed in lattice perturbation theory have been computed on the PC clusters Marconi at CINECA (CINECA-INFN and CINECA-Bicocca agreements) and Wilson at Milano-Bicocca. We thank these institutions for the computer resources and the technical support. We also acknowledge PRACE and ISCRA for awarding us access to MareNostrum at Barcelona Supercomputing Center (BSC), Spain (n. 2018194651) and to Marconi at CINECA (EoSQCD) respectively, where the non-perturbative computations are being performed. We acknowledge partial support by the INFN project "High performance data network".

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 1-loop 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 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 while, thanks to Eqs. (A.8), traces of products of an odd number of γ-matrices vanish. 10 Summation over repeated indices is always understood.

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 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. [24], see also Ref. [40] 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 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 D = D w + aD sw . (C.5) 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 The second term is the Sheikholeslami-Wohlert operator defined by [15] is the clover discretization of the field strength tensor 12 (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,30]. 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 [41,42]. This occurs, for instance, in a finite volume without boundaries [42] or in the thermal theory at high temperature [43]. 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. [41,42]. In the massless limit, the action of the Symanzik's effective continuum theory reads where S 0 is defined as in Eq. (B.2) but with the bare coupling replaced by the renormalized one, and The very same conclusion can be reached by using the R 5 discrete symmetry in Ref. [41] 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 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.53), 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 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 where S F P,0 is the tree-level term, while

E.6 Quark-gluon interaction
By inserting Eq. (E.1) into Eq. (C.4) and by expanding in g 0 , the fermionic action reads 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 free-energy 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 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.

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 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 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.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 (H.31) The contribution due to the improvement term to the fermionic part of the energymomentum tensor is made up of 3 terms

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.