Thermodynamics in quenched QCD: energy--momentum tensor with two-loop order coefficients in the gradient flow formalism

Recently, Harlander et al.\ [Eur.\ Phys.\ J.\ C {\bf 78}, 944 (2018)] have computed the two-loop order (i.e., NNLO) coefficients in the gradient-flow representation of the energy--momentum tensor (EMT) in vector-like gauge theories. In this paper, we study the effect of the two-loop order corrections (and the three-loop order correction for the trace part of the EMT, which is available through the trace anomaly) on the lattice computation of thermodynamic quantities in quenched QCD. The use of the two-loop order coefficients generally reduces the $t$~dependence of the expectation values of the EMT in the gradient-flow representation, where $t$~is the flow time. With the use of the two-loop order coefficients, therefore, the $t\to0$ extrapolation becomes less sensitive to the fit function, the fit range, and the choice of the renormalization scale; the systematic error associated with these factors is considerably reduced.


Introduction
The energy-momentum tensor (EMT) T µν (x) is a fundamental physical observable in quantum field theory. It has been pointed out in Refs. [1,2] that a "universal" representation of the EMT can be written down by utilizing the so-called gradient flow [3][4][5][6][7] and its small flow-time expansion [6]. This representation of the EMT is universal in the sense that it is independent of the adopted regularization. The representation can thus be employed in particular with the lattice regularization that makes nonperturbative computations possible. An advantage of this approach to the lattice EMT is that the expression of the EMT is known a priori and it is not necessary to compute the renormalization constants involved in the lattice EMT [8]. 1 This approach instead requires the limit t → 0, where t is the flow time (see below), because the representation is obtained in the small flow-time limit. In actual lattice simulations, however, since t is limited as t a 2 by the lattice spacing a, the t → 0 limit has to be obtained by the extrapolation from the range of t satisfying t a 2 ; this t → 0 extrapolation can be a source of systematic error. By employing this gradient-flow representation of the EMT, expectation values and correlation functions of the EMT have been computed to study various physical questions [12][13][14][15][16][17][18][19][20][21].
One important application of the lattice EMT is the thermodynamics of gauge theory at finite temperature; see Refs. [22][23][24][25][26] and more recent works in Refs. [27][28][29][30][31][32] on this problem. Two independent thermodynamic quantities, such as the energy density ε and the pressure p, can be computed from the finite-temperature expectation value of the traceless part and the trace part of the EMT, respectively, as 2 In the gradient-flow approach, moreover, the computation of isotropic/anisotropic Karsch coefficients (i.e., the lattice β function) is not necessary [33], because the expression of the EMT is a priori known.
In this paper, we investigate the thermodynamics in quenched QCD (quantum chromodynamics), i.e., the pure Yang-Mills theory, in the gradient-flow approach. The EMT in the gradient-flow representation is obtained as follows. Assuming dimensional regularization, the EMT in the pure Yang-Mills theory is given by where g 0 is the bare gauge coupling and F a µν (x) is the field strength. 3 Note that this is an expression in D ≡ 4 − 2 -dimensional spacetime and is not generally traceless.
One can express any composite operator in gauge theory such as the EMT (1.3) as a series of flowed composite operators through the small flow-time expansion [6]. That is, one can . In these expressions, the "flowed" gauge field B a µ (t, x) is defined by the gradient flow [3][4][5], i.e., a one-parameter evolution of the gauge field by The parameter t ≥ 0, which possesses the mass dimension −2, is termed the flow time. Since Eq. (1.4) is finite [6], one can set D = 4 and the first term on the right-hand side (that is proportional to c 1 (t)) is traceless. The coefficients in this small flow-time expansion, which are analogous to the Wilson coefficients in OPE, can be calculated by perturbation theory [6] as where g denotes the renormalized gauge coupling. Throughout this paper, we assume the MS scheme, in which Here, µ is the renormalization scale, γ E is the Euler constant, and Z g is the renormalization factor. In Eq. (1.6), On the other hand, there is no "k (0) 2 " in Eq. (1.6) because the EMT is traceless in the tree-level approximation (the trace anomaly emerges from the one-loop order).
In Eq. (1.6), the one-loop order (i.e., NLO) coefficients k (1) i (t) (i = 1, 2) were computed in Refs. [1,2] (see also Ref. [35]). Recently, in Ref. [34], Harlander et al. have computed the two-loop order (i.e., NNLO) coefficients k (2) i (t) for general vector-like gauge theories; see also Ref. [36]. The purpose of the present paper is to study the effect of the two-loop corrections given in Ref. [34] by performing the lattice computation of thermodynamic quantities in quenched QCD. For the trace part of the EMT, we also examine the use of the threeloop order coefficient, k 2 , which is presented in this paper; this higher-order coefficient can be obtained for quenched QCD by combining a two-loop result in Ref. [34] and the trace anomaly [37][38][39], as we will explain below. From analyses using lattice data obtained in Ref. [14], we find that the use of the two-loop order coefficients generally reduces the t dependence of the expectation values of the EMT in the gradient-flow representation. With the use of the two-loop order coefficients, therefore, the t → 0 extrapolation becomes less sensitive to the fit function, the fit range, and the choice of the renormalization scale; the systematic error associated with these factors is considerably reduced. We expect that this improvement brought about by the two-loop order coefficients also persists in wider applications of the gradient-flow representation of the EMT, such as the thermodynamics of full QCD. This paper is organized as follows. In Sect. 2, we explain our treatment of perturbative coefficients c 1 (t) and c 2 (t) of Eq. (1.6). We list the known expansion coefficients k ( ) i , and present the three-loop coefficient for c 2 (t), k 2 . In Sect. 3, we perform numerical analyses of the thermodynamic quantities, which are mainly based on FlowQCD 2016 [14]. We give conclusions in Sect. 4.

Expansion coefficients 2.1. β function and the running gauge coupling constant
The β function corresponding to Eq. (1.7) is given by with coefficients [40][41][42][43][44][45][46] where C A is the quadratic Casimir for the adjoint representation defined by To compute expectation values of the EMT by employing the representation (1.4), we take the limit t → 0 [1,2]. First of all this limit removes the last O(t) term in Eq. (1.4), the contribution of operators of higher (≥ 6) mass dimensions. It also justifies finite-order truncation of perturbative expansions of the coefficients c i (t); we treat c i (t) in Eq. (1.4) as follows. We apply the renormalization group improvement [1,2], i.e., we set µ ∝ 1/ √ t in k ( ) i and concurrently replace the coupling constant g(µ) with the running gauge coupling satisfying where the coupling is now a function of µ ∝ 1/ √ t. Then the t → 0 limit allows us to neglect the higher-order terms in g 2 because the running gauge coupling g(µ ∝ 1/ √ t) goes to zero due to the asymptotic freedom. We note that this renormalization group improvement is legitimate since the coefficients c i (t) (i = 1, 2) are independent of the renormalization scale µ (when the bare coupling g 0 is kept fixed). This can be seen from the fact that the EMT (1.3) and the operator G a µρ (t, x)G a νρ (t, x) are bare quantities. Although the above argument shows that in principle the coefficients c i (t) are independent of the choice of the relation between the renormalization scale µ and the flow time t, i.e., of the parameter c in µ = c/ √ t, this independence does not exactly hold in practical calculations based on fixed-order perturbation theory. In other words, the difference caused by different choices of c implies the remaining perturbative uncertainty. Following Ref. [34], we introduce the combination L(µ, t) ≡ ln(2µ 2 t) + γ E . (2.8) A conventional choice of µ is given by All the numerical experiments on the basis of the representation (1.4) so far [12][13][14][15][16][17][18][19][20][21] have adopted this choice. On the other hand, in Ref. [34], it is argued that would be an optimal choice on the basis of the two-loop order coefficients. 5 In the following numerical analyses, we will examine both choices µ = µ 0 (t) and µ = µ d (t). The difference in the results with these two choices gives an estimate of higher-order uncertainty, where Let us now list the known coefficients in Eq. (1.6).

One-loop order (NLO) coefficients
The number L is defined by Eq. (2.8).

Two-loop order (NNLO) coefficients
The two-loop order coefficients in Ref. [34] specialized to the pure Yang-Mills theory are In the pure Yang-Mills theory, if one has the small flow-time expansion of the renormalized operator {F a µν F a µν } R (x) in the th-loop order, it is possible to further obtain the coefficient of c 2 (t) one loop higher, k ( +1) 2 , by using information on the trace anomaly [1]. The twoloop order (NNLO) coefficient (2.14) can also be obtained in this way from a one-loop order calculation and has already been used in numerical experiments in quenched QCD. Repeating this argument, we can now obtain the three-loop order coefficient, k We recall the trace anomaly [37][38][39] where the β function is given by Eq. (2.1). According to Eq. (64) of Ref. [34], we now have the small flow-time expansion of {F a µν F a µν } R (x) to the two-loop order: Plugging this into Eq. (2.15) and using Eq. (2.1), we have Comparing this with the trace of Eq. (1.4), we obtain One can also confirm that Eqs. (2.12) and (2.14) are correctly reproduced. We will examine the use of this N 3 LO coefficient for the trace anomaly in the numerical analyses below.

Numerical analyses
In what follows, we use the lattice data obtained in Ref. [14] for the G = SU (3) pure Yang-Mills theory and study the effects of the higher-order coefficients quantitatively. We do not repeat the explanation of the lattice setups; see Ref. [14] for details. We measure the entropy density and trace anomaly (which are normalized by the temperature). To obtain these thermodynamic quantities, the double limit, a → 0 and t → 0, is required because Eq. (1.4) is the relation in continuum spacetime and in the small flow-time limit. We first take the continuum limit a → 0 and then take the small flow-time limit t → 0 [14]. In continuum extrapolation while keeping t in physical units fixed, we adopt only the range where the flow time satisfies t a 2 . This is because the flow time t is meaningful only for t a 2 with finite lattice spacing a; the lattice data actually exhibit violently diverging behavior for t a 2 due to finite lattice spacing effects. Hence, at this stage, we cannot obtain the continuum limit for t = 0. Thus, we carry out the t → 0 extrapolation by assuming a certain functional form of Eq. (3.1) with respect to t. Let us start with the entropy density, ε + p. From Eq. (1.1), (ε + p)/T 4 is obtained by taking the t → 0 limit of the thermal expectation value: In Fig. 1, we plot Eq. (3.1) as a function of tT 2 ; the temperature is T /T c = 1.68. The plots for the other temperatures listed in the left-most column of Table 2 are deferred to Appendix A.
In each panel of Fig. 1, we show lattice results for Eq.  23) and (A2) in Ref. [14]. We use the four-loop running gauge coupling in the MS scheme, adopting the approximate formula (9.5) in Ref. [47]. 6 We then take the continuum limit a → 0 at each fixed value of tT 2 . (See Ref. [14] for details of this procedure.) The continuum limit results are shown by gray bands. In Fig. 1 and in the corresponding figures in Appendix A, Figs. A1-A7, we observe that the use of the twoloop order coefficient generally reduces the t dependence of the continuum limit (it becomes flatter in t). 7 This behavior generally allows us to perform stable t → 0 extrapolation.
We carry out the t → 0 extrapolation with a linear function in t as Ref. [14]. The following three reasonable fit ranges are adopted for the t → 0 extrapolation [14] In Table 1, the coefficients of the linear term in t, determined in t → 0 extrapolation using Range 1, are shown. The tendency for the slopes to get smaller at N 2 LO than at NLO is quantitatively observed. In addition to the linear function in t, we also use the linear function of [g(µ) 2 /(4π) 2 ] +1 , where = 1 for the NLO approximation and = 2 for the N 2 LO approximation. 9 This functional form is suggested from a detailed study of the asymptotic behavior of Eq. (3.1) for t → 0 (H. Suzuki and H. Takaura, manuscript in preparation). 6 In Ref. [14], the Λ parameter is obtained at the three-loop level, while we use the four-loop running coupling. The effect of the discrepancy in the perturbation order would be able to be taken into account by varying Λ MS in our error analysis below. 7 An explanation for this flatter behavior will be provided in another paper (H. Suzuki and H. Takaura, manuscript in preparation). 8 The finite lattice spacing and volume effects are controlled by tT 2 [14]. 9 Recall that the renormalization scale µ is a function of t through Eq. (2.8).
We estimate the systematic error associated with the t → 0 extrapolation by examining the variation obtained from the different extrapolation function. As mentioned, the use of the two-loop coefficient leads to a flatter behavior with respect to t. Hence, t → 0 extrapolation becomes less sensitive to the fit function, the fit range, and the choice of the renormalization scale.
In Table 2 Table 1: Linear coefficients in the t → 0 extrapolation. The renormalization scale µ 0 (t) and Range 1 are used.
function are greatly reduced in the N 2 LO approximation. This clearly shows an advantage of the two-loop order coefficient. In Fig. 2, ( + p)/T 4 is plotted as a function of T /T c . The error bar represents the total error, obtained by combining all the errors in quadrature. Our N 2 LO results are consistent with the results of Refs. [14,22,24,30,32], especially with Refs. [22,32].   Table 2: Summary of the entropy density and the trace anomaly obtained by using coefficients in different orders of perturbation theory. The statistical errors are shown in the first parentheses. The numbers in the other parentheses show systematic errors: the error associated with the fit range, the 3% uncertainty of Λ MS , the renormalization scale, and the t → 0 extrapolation function are shown in the second, third, fourth, and fifth parentheses, respectively. The results of Ref. [14] (FlowQCD 2016) are also tabulated in the last column: the numbers in the first, second, and third parentheses are the statistical error, the systematic error associated with the choice of the fit range (only the linear t fit is adopted in Ref. [14]), and the systematic error associated with Λ MS (varying it by 1% in Ref. [14]). (The difference in the central values between FlowQCD 2016 and our analyses at the lower orders stems from the choice of the renormalization scale; µ 0 (t) is adopted in the present work rather than µ d (t)).
Reference [14] is the preceding analysis performed with the same method as this work but with the NLO coefficient. In our NLO analysis, the systematic errors are investigated in more detail, including some additional error sources. This leads to larger errors than in Ref. [14], while the central values of Ref. [14] are consistent with the present work. Figure 2 clearly shows that the use of the two-loop coefficient generally reduces the systematic errors.
We now turn to the trace anomaly, ε − 3p, which is investigated in a parallel manner to the entropy density. The expectation value as a function of tT 2 is plotted in Fig. 3 for T /T c = 1.68. (Results for other temperatures are deferred to Appendix A.) As we noted, the two-loop order (N 2 LO) and the three-loop order (N 3 LO) coefficients are available for the trace anomaly. We observe that already with the N 2 LO coefficient the continuum limit (the gray band) is almost constant in t within our fit ranges. Thus, naturally, the extrapolation of the continuum limit to t = 0 is quite insensitive to the choices N 2 LO or N 3 LO, and µ = µ 0 (t) or µ = µ d (t). Similarly to the entropy density above, we use the linear function of t for the t → 0 extrapolation. We also use the linear function of [g(µ) 2 /(4π) 2 ] , where = 2 for the N 2 LO approximation and = 3 for the N 3 LO approximation, as suggested from the study of the asymptotic t → 0 behavior of Eq. (3.5) (H. Suzuki and H. Takaura, manuscript in preparation). The difference caused by this is treated as a systematic error. The results are summarized in Table 2 and Fig. 4. All the results are almost degenerate as seen from Fig. 4. However, it is worth noting that the use of the N 3 LO coefficient certainly reduces the dependences on the choice of the renormalization scale and the fit function, as seen from the fourth and the fifth parentheses in Table 2.    Fig. 4: Summary of the trace anomaly (ε − 3p)/T 4 as a function of T /T c . In the righthand panel, the region 1.00 (ε − 3p)/T 4 2.75 is magnified. The results from the present paper are the red circles (N 2 LO) and the blue squares (N 3 LO). The error bars include the systematic error as well as the statistical error; see Table 2 for details. For comparison, we also show the results of Refs. [14,22,24,30,32].

Conclusions
We investigated the thermodynamics in quenched QCD using the gradient-flow representation of the EMT. In particular, we studied the effect of the N 2 LO coefficients in the gradient-flow formalism, which have become available recently. For the trace anomaly, we used the N 3 LO coefficient, which was obtained in this paper for quenched QCD. It turned out that the use of the N 2 LO (or N 3 LO) coefficients considerably reduces the systematic errors, especially concerning the choice of the renormalization scale and the t → 0 extrapolation. We expect that the use of the N 2 LO coefficients will also make precise studies possible in full QCD, which has been investigated with the NLO coefficients so far.