# t→0 extrapolation function in SF t X method for the energy-momentum tensor

aa r X i v : . [ h e p - l a t ] F e b Preprint number: KYUSHU-HET-220, KEK-TH-2294 t → extrapolation function in SF t X methodfor the energy–momentum tensor

Hiroshi Suzuki and Hiromasa Takaura , ∗ Department of Physics, Kyushu University 744 Motooka, Nishi-ku, Fukuoka, 819-0395, Japan Theory Center, High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki,305-0801, Japan ∗ E-mail: [email protected] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

We theoretically clarify the functional form to be used in t → t X) method for the energy–momentum tensor (EMT),which facilitates lattice simulation of the EMT based on the gradient ﬂow. We arguethat in the t → g ( µ ( t )), the ﬂow time dependent running coupling, where the power is determined bythe perturbation order we consider. From actual lattice data, we conﬁrm the validity ofthe extrapolation function. Using the new extrapolation function, we present updatedlattice results for thermodynamics quantities in quenched QCD; our results are con-sistent with the previous study [arXiv:1812.06444] but we obtain smaller errors due toreduction of systematic errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Subject Index B01, B31, B32, B38 typeset using PTP

TEX.cls . Introduction

The energy–momentum tensor (EMT) T µν ( x ) is a fundamental quantity in quantum ﬁeldtheory, yet its lattice simulation is not straightforward due to explicit breaking of the trans-lation invariance in lattice regularization; see Ref. [1] and references therein. In Refs. [2, 3],the so-called small ﬂow time expansion (SF t X) method was devised to solve this problem.In this method, one rewrites conserved currents in terms of the so-called ﬂowed operators,deﬁned from the gradient ﬂow [4–8]. Due to UV ﬁniteness of ﬂowed operators, the currentrepresented with ﬂowed operators satisﬁes the correct Ward-Takahashi identity in any reg-ularization. Then one can measure correctly normalized currents in lattice simulation. Forthe EMT, we can schematically write it in the form T µν ( x ) = ˜ c ( t ) ˜ O ,µν ( t, x ) + ˜ c ( t ) ˜ O ,µν ( t, x ) + ˜ c ( t ) ˜ O ,µν ( t, x ) + ˜ c ( t ) ˜ O ,µν ( t, x ) + O ( t ) . (1.1)Here t is the ﬂow time, whose mass dimension is −

2; ˜ O i,µν ( t, x )’s are (dimension-four)ﬂowed operators (whose explicit deﬁnitions are given below). The coeﬃcients ˜ c i ( t )’s can beperturbatively calculated via the small ﬂow time expansion [7] of the ﬂowed operators and thetwo-loop order results are known today [9]; see also Ref. [10]. O ( t ) represents contributionsfrom dimension-six operators. In lattice simulation, we measure the ﬂowed operators ˜ O i,µν nonperturbatively and then multiply the perturbative coeﬃcients ˜ c i ( t ) to obtain the EMT.Many lattice simulations of the EMT with the SF t X method were performed [11–27] and itsusefulness has been conﬁrmed.In lattice simulation using the SF t X method, one needs to take the small ﬂow time limit( t → This is because the expression for the EMT in the SF t X method becomes exactin the t → g ( µ ( t )), the ﬂow timedependent running coupling, in the perturbative coeﬃcients ˜ c i ( t ), and O ( t ) contributionsin Eq. (1.1). However one cannot directly obtain lattice data at t = 0 because lattice datasuﬀer from serious discretization eﬀects when the ﬂow time becomes too small compared tolattice spacing, t . a . Hence, one should take the t → t ≫ a assuming some function of t . Thus, a t → t X method, a linear function in t was mainly used in the t → O ( t ) contribution in Eq. (1.1), which is neglectedin constructing the EMT in the SF t X method. However, since we use ﬁxed order perturba-tive results for ˜ c i ( t ), a higher-order eﬀect in g ( µ ( t )) should also exist. Parametrically, such aneﬀect is dominant in small t region. This is because g ( µ ( t )) n ≫ t Λ ∼ e − π / ( β g ( µ ( t ))) ( β is the ﬁrst coeﬃcient of the beta function), where the t is an exponentially suppressed func-tion for small g ( µ ( t )). It is therefore important to identify the higher-order eﬀect in g ( µ ( t ))in order to accurately performe t → t dependence of the ﬁxed order formulafor the EMT in the SF t X method (which means that we use ﬁxed order perturbative results In actual study, we mainly use the expression of Eqs. (1.30) and (1.31). We show Eq. (1.1) justfor explanation. In principle, the continuum limit ( a → a is lattice spacing) should be taken beforethe t → c i ( t )’s). This study tells us which functional form should be used in the t → t X method a more preciseframework. Also the recent two-loop order calculation of ˜ c i ’s [9] can promote more accurateanalyses, and it is timely to discuss this issue in the SF t X method.In the rest of this this section, we brieﬂy review the SF t X method and clarify the ques-tion studied in this paper more explicitly. We also introduce quantities necessary for thesubsequent studies.The EMT in dimensional regularization is given by T µν ( x ) = 1 g (cid:20) O ,µν ( x ) − O ,µν ( x ) (cid:21) + 14 O ,µν ( x ) , (1.2)where g is the bare gauge coupling; O i,µν ’s are gauge invariant and symmetric dimension-four tensor operators deﬁned as O ,µν ( x ) ≡ F aµρ ( x ) F aνρ ( x ) , (1.3) O ,µν ( x ) ≡ δ µν F aρσ ( x ) F aρσ ( x ) , (1.4) O ,µν ( x ) ≡ X f ¯ ψ f ( x ) (cid:16) γ µ ←→ D ν + γ ν ←→ D µ (cid:17) ψ f ( x ) , (1.5) O ,µν ( x ) ≡ δ µν X f ¯ ψ f ( x ) ←→ / D ψ f ( x ) , (1.6) O ,µν ( x ) ≡ δ µν X f m f, ¯ ψ f ( x ) ψ f ( x ) , (1.7)where m f, is the bare mass of the ﬂavor f . These are bare composite operators and notﬁnite in general.In the SF t X method, one rewrites the EMT in terms of ﬂowed operators. The Yang–Millsgradient ﬂow [4–8] is deﬁned by the following diﬀerential equations with respect to the ﬂowtime t . For the gauge ﬁeld, it is deﬁned as ∂ t B µ ( t, x ) = D ν G νµ ( t, x ) + α D µ ∂ ν B ν , B µ ( t = 0 , x ) = A µ ( x ) , (1.8)with a gauge parameter α , where the covariant derivative on the gauge ﬁeld is deﬁned as D µ = ∂ µ + [ B µ , · ] . (1.9) G µν is the ﬁeld strength of the ﬂowed gauge ﬁeld B µ . For the fermion ﬁeld, it is given by ∂ t χ ( t, x ) = [∆ − α ∂ µ B µ ( t, x )] χ ( t, x ) , χ ( t = 0 , x ) = ψ ( x ) , (1.10) ∂ t ¯ χ ( t, x ) = ¯ χ ( t, x )[ ←− ∆ + α ∂ µ B µ ( t, x )] , ¯ χ ( t = 0 , x ) = ¯ ψ ( x ) , (1.11)with ∆ = D µ D µ , D µ = ∂ µ + B µ , (1.12) ←− ∆ = ←− D µ ←− D µ , ←− D µ = ←− ∂ µ − B µ . (1.13)3e deﬁne analogous operators at positive ﬂow time with tildes:˜ O ,µν ( t, x ) ≡ G aµρ ( t, x ) G aνρ ( t, x ) , (1.14)˜ O ,µν ( t, x ) ≡ δ µν G aρσ ( t, x ) G aρσ ( t, x ) , (1.15)˜ O ,µν ( t, x ) ≡ X f ˚¯ χ f ( t, x ) (cid:16) γ µ ←→ D ν + γ ν ←→ D µ (cid:17) ˚ χ f ( t, x ) , (1.16)˜ O ,µν ( t, x ) ≡ δ µν X f ˚¯ χ f ( t, x ) ←→ / D ˚ χ f ( t, x ) , (1.17)˜ O ,µν ( t, x ) ≡ δ µν X f m f ˚¯ χ f ( t, x )˚ χ f ( t, x ) . (1.18)The above ﬂowed operators are ﬁnite operators [7]. We have accomplished the renormal-ization of ﬂowed fermion ﬁelds with use of ringed variable [3]:˚ χ f ( t, x ) ≡ vuut − R )(4 π ) t D ¯ χ f ( t, x ) ←→ / D χ f ( t, x ) E χ f ( t, x ) , (1.19)˚¯ χ f ( t, x ) ≡ vuut − R )(4 π ) t D ¯ χ f ( t, x ) ←→ / D χ f ( t, x ) E ¯ χ f ( t, x ) . (1.20)Alternatively, one can also adopt the MS scheme as in Ref. [9]. We denote such an operatorset as ˜ O MS i,µν ( t, x ). The conversion from the MS scheme to the ringed variable scheme is carriedout by a matrix R by ˜ O µν ( t, x ) = R ˜ O MS µν ( t, x ) , (1.21)with R = r

00 0 0 r , (1.22)where r is given in Eq. (51) of Ref. [9] as r = ˚ Z χ /Z χ . (1.23)When we suppress the index i in the operators, it is understood as a four component vector, O µν = O ,µν O ,µν O ,µν O ,µν . (1.24)Note that we can eliminate O ,µν by using the equation of motion (EOM), (1 / O ,µν + O ,µν = 0, following Ref. [9]. Hence, it is suﬃcient to consider the four operators as a basis. In Eq. (1.18) m f denotes a renormalized mass parameter. However, in this paper, we do not needto specify it because we do not use ˜ O ,µν in the following.

4o rewrite the EMT in terms of the ﬂowed operators, one needs to know the relationbetween the ﬂowed operators and unﬂowed operators. It can be studied through the smallﬂow time expansion [7]: ˜ O MS µν ( t, x ) = ζ ( t ) O µν ( x ) + O ( t ) , (1.25)where ζ ( t ) is a 4 × O ( t )for small t . Here and hereafter, we implicitly assume that the vacuum expectation values ofthe operators are subtracted. [This is the reason why the terms proportional to the identityoperator are absent in Eq. (1.25).]Once the relation between the ﬂowed operators and unﬂowed operators is known, one canexpress the EMT, given by the unﬂowed operators, in terms of the ﬂowed operators usingthe inverse of the matrix ζ . It is conventionally expressed as T µν ( x ) = c ( t ) (cid:20) ˜ O ,µν ( t, x ) −

14 ˜ O ,µν ( t, x ) (cid:21) + c ( t ) ˜ O ,µν ( t, x )+ c ( t ) h ˜ O ,µν ( t, x ) − O ,µν ( t, x ) i + c ( t ) ˜ O ,µν ( t, x ) . (1.26)By decomposition of a (ﬁnite) symmetric tensor S µν into its traceless (TL) and scalar (S)parts as S TL µν = S µν − δ µν δ ρσ S ρσ , (1.27) S S = δ ρσ S ρσ , (1.28)the expression (1.26) is equivalent to T µν ( x ) = T TL µν ( x ) + δ µν T S ( x ) , (1.29)with T TL µν ( x ) = c ( t ) ˜ O TL1 ,µν ( t, x ) + c ( t ) ˜ O TL3 ,µν ( t, x ) , (1.30) T S ( x ) = c ( t ) ˜ O S2 ( t, x ) + c ′ ( t ) ˜ O S4 ( t, x ) , (1.31)where c ′ ( t ) = c ( t ) − c ( t ) . (1.32)The traceless ﬂowed operators are given by ˜ O TL1 ,µν ( t, x ) = ˜ O ,µν ( t, x ) − (1 /

4) ˜ O ,µν ( t, x ) and˜ O TL3 ,µν ( t, x ) = ˜ O ,µν ( t, x ) − (1 /

2) ˜ O ,µν ( t, x ). In the SF t X method, we neglect O ( t ) contributions in Eq. (1.25), and thus, strictly speaking, the O ( t ) error should be shown in the right-hand side of Eq. (1.26). c ( ′ ) i ( t ) are given by perturbative series: c ( t ) = 1 g ( µ ) ∞ X n =0 k ( n )1 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.33) c ( t ) = 1 g ( µ ) ∞ X n =1 k ( n )2 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.34) c ( t ) = ∞ X n =0 k ( n )3 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.35) c ′ ( t ) = ∞ X n =0 k ( n )4 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n . (1.36)Here µ is a renormalization scale and a perturbative coeﬃcient k ( n ) i is a polynomial of L ( µ, t ) ≡ log (2 µ e γ E t ). Since the EMT and ﬂowed operators are renormalization group (RG)invariant, the coeﬃcients c i ( t )’s are also RG invariant, namely, independent of µ . In practicalapplication, we take the t -dependent renormalization scale, µ = µ ( t ) = s/ √ e γ E t , where s isan O (1) numerical factor, so that higher-order terms in g ( µ ) vanish in the t → k (1) i ’s were calculated in Ref. [2] for quenched QCD and in Ref. [3] forfull QCD. The two-loop coeﬃcients k (2) i were obtained in Ref. [9]. For quenched QCD, thethree-loop coeﬃcient k (3)2 (only for i = 2) was obtained in Ref. [21].The purpose of this paper is to investigate the (leading) t -dependence of the next-to-...-next-to-leading order (N k LO) formula of the EMT in the SF t X method: T TL µν (N k LO) ( x ; t ) = c (N k LO)1 ( t ) ˜ O ,µν ( t, x ) + c (N k LO)3 ( t ) ˜ O ,µν ( t, x ) , (1.37) T S (N k LO) ( x ; t ) = c (N k LO)2 ( t ) ˜ O S2 ( t, x ) + c ′ (N k LO)4 ( t ) ˜ O S4 ( t, x ) . (1.38)Here c ( ′ ) i (N k LO) is given by the sum of the ﬁrst ( k + 1) terms of Eqs. (1.33)–(1.36): c (N k LO)1 ( t ) = 1 g ( µ ) k X n =0 k ( n )1 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.39) c (N k LO)2 ( t ) = 1 g ( µ ) k +1 X n =1 k ( n )2 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.40) c (N k LO)3 ( t ) = k X n =0 k ( k )3 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (1.41) c ′ k LO) ( t ) = k X n =0 k ( k )4 ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n . (1.42)Note that for c ( t ) the upper edge of the sum is set to k + 1. Our main results of this paperare summarized in Sect. 2. 6n the subsequent calculations, we mainly use hatted operators, as in Ref. [9], to discussrenormalization of the bare operators, instead of the set of Eqs. (1.3)–(1.6):ˆ O µν ( x ) = H O µν ( x ) , (1.43)where O µν ( x ) is deﬁned in Eqs. (1.3)–(1.6) and H denotes H = /g /g . (1.44)ˆ O µν are D -dimensional operators in dimensional regularization with D = 4 − ǫ . We considerthe renormalization of ˆ O µν as ˆ O µν ( x ) = Z ( g ( µ ); ǫ ) ˆ O Rµν ( x ; µ ) , (1.45)where“ R ” stands for renormalization of operators. We use the MS scheme to renormalizethe unﬂowed operators. Since the left-hand side of Eq. (1.25),˜ O MS µν ( t, x ) = ζ ( t ) O µν ( x ) + O ( t ) = ζ ( t ) H − Z ( g ( µ )) ˆ O Rµν ( x ; µ ) + O ( t ) , (1.46)is ﬁnite, one can determine Z from the requirement that ζ ( t ) H − Z is ﬁnite. In this way, Z has been calculated up to two-loop [9]. Now, we can express the ﬂowed operators with therenormalized (ﬁnite) quantities,˜ O µν ( t, x ) = ζ R ( t ; g ( µ ) , µ ) ˆ O Rµν ( x ; µ ) + O ( t ) . (1.47)Here we deﬁne the renormalized matrix ζ R as ζ R ( t ; g ( µ ) , µ ) = Rζ ( t ) H − Z ( g ( µ )) , (1.48)which is ﬁnite.We also use the following relations: g = Z g (cid:18) µe γ E / √ π (cid:19) ǫ g ( µ ) , (1.49) β ( g ) = µ dgdµ = − ∞ X i = − β i g (cid:20) g (4 π ) (cid:21) i +1 , (1.50) m = Z m m ( µ ) , (1.51) γ m = − µ d log mdµ = ∞ X i =0 γ ( i ) m (cid:20) g (4 π ) (cid:21) i +1 , (1.52)where Z g is given by Z g = 1 − β ǫ g (4 π ) + (cid:18) β ǫ − β ǫ (cid:19) (cid:20) g (4 π ) (cid:21) + · · · . (1.53) In our convention, Z is inverse of that of Ref. [9].

7n the beta function, there is β − = ǫ , which is non-zero in D -dimensional spacetime. Wealso list the ﬁrst few coeﬃcients: β = 113 C A − T F , (1.54) β = 343 C A − (cid:18) C F + 203 C A (cid:19) T F , (1.55) γ (0) m = 6 C F , (1.56) γ (1) m = 973 C A C F + 3 C F − C F T F . (1.57)An explanation of the constants C A , T F , and C F is given in App. A.The rest of the paper is organized as follows. In Sect. 2, we ﬁrst explain our main resultsand show the functional form to be used in the t → ζ R , which can be trivially obtained, and the one-loop anomalous dimensionof the unﬂowed composite operators. In Sect. 4, as an additional study, we give a generalargument how we can study in detail the t dependence caused by dimension-six operators(which is denoted roughly by O ( t ) above). An explicit result is given for the traceless partof the EMT in quenched QCD. In Sect. 5, we perform numerical analysis using lattice data.Here, we study the thermodynamic quantities, in particular the entropy density and traceanomaly, which are proportional to T TL µν ( x ) and T S ( x ), respectively. We conﬁrm validityof the t → L ( µ, t ) dependence of the perturbative series forthe coeﬃcients c ( ′ ) i ( t ) is presented. Also the one-loop results are shown. In App. E, we showthe relations that the matrix K (which is introduced in Sect. 3) should satisfy. In App. F, wepresent the NLO ζ R , which becomes necessary if one wants to study higher-order behaviorof the t dependence remaining in the ﬁxed order formula for the EMT. (The NLO ζ R isused in App. G.) In App. G, we give an argument to estimate higher-order eﬀects which areneglected in Sect. 3. In App. H, for reference, we present the results of the thermodynamicsquantities obtained with linear-type extrapolation functions. t → extrapolation functions We ﬁrst explain our results in quenched QCD. The lattice data obtained with the N k LOformulae of the EMT, T TL µν (N k LO) ( x ; t ) = c (N k LO)1 ( t ) ˜ O TL1 ,µν ( t, x ) , (2.1) T S(N k LO) ( x ; t ) = c (N k LO)2 ( t ) ˜ O S2 ( t, x ) , (2.2)where c (N k LO)1 and c (N k LO)2 are given by Eqs. (1.39) and (1.40), should be ﬁtted with thefunctions of t , T TL µν (N k LO) ( x ; t ) = − k ( k +1)1 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 ! T TL µν ( x ) , (2.3)8 S(N k LO) ( x ; t ) = − β k ( k +2)2 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 ! T S ( x ) , (2.4)in t → T TL µν ( x ), T S ( x ), k ( k +1)1 , and k ( k +2)2 are ﬁt parameters. g ( µ ( t )) is a running coupling and the renormalization scale µ ( t ) should be taken commonto that of c (N k LO)1 , ( t ). T TL µν ( x ) and T S ( x ) are the EMT we want to extract by the t → k LO formulae of the EMT, given inEqs. (1.37) and (1.38), should be ﬁtted with the functions, T TL µν (N k LO) ( x ; t ) = T TL µν ( x ) + a ( x ) (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 , (2.5) T S(N k LO) ( x ; t ) = T S ( x ) + b ( x ) (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 . (2.6)Here T TL µν ( x ), T S ( x ), a ( x ), and b ( x ) are ﬁt parameters. T TL µν ( x ) and T S ( x ) correspond to ﬁnalresults in which one is interested.In quenched QCD, the ﬁt parameters k ( k +1)1 and k ( k +2)2 correspond to higher-order coeﬃ-cients [cf. Eqs. (1.33) and (1.34)]. In this sense, (i) they do not depend on x or the typicalscale Q of a considered system (for instance, when the expectation value of the one-pointfunction of the EMT is considered at ﬁnite temperature T , the typical scale is Q = T ), and(ii) it can be predicted how they vary in response to the variation of the parameter s in µ ( t ) = s/ √ e γ E t as long as the N k LO coeﬃcients are known; see Sect. 5 and App. D. Ofcourse, these properties are not exactly hold due to systematic errors in ﬁts. Nevertheless,we can check validity of the use of the above extrapolation function by looking into theseproperties. In our analyses in Sect. 5, we will take the ﬁt parameters k ( k +1)1 and k ( k +2)2 common to all the simulated temperatures, taking the ﬁrst property into account. We checkthe validity of the use of the above extrapolation functions by examining the second prop-erty (ii), behavior of the ﬁt parameters k ( k +1)1 and k ( k +2)2 upon the variation of the choice ofthe renormalization scale.On the other hand, in full QCD, the origin of the ﬁt parameters a ( x ) and b ( x ) is not sosimple. Hence we cannot expect parallel properties to the quenched QCD case. For instance,they generally depend on x and the typical scale Q of a system. Also, it is not apparent howthey change upon the choice of the renormalization scale.

3. Derivation of t -dependence of N k LO formula

We investigate the leading t -dependence of the N k LO formula of the EMT in the SF t Xmethod and derive the results of Eqs. (2.3), (2.4), (2.5), and (2.6). Throughout this section,we neglect O ( t ) eﬀect coming from dimension-six operators, which is a subleading eﬀect forsuﬃciently small t . One should consider the ( k + 1) or higher-loop beta function for the N k LO calculation incalculating the running. The same applies to the full QCD case. Nowadays, c is known to NNLO but c is known to NLO. If one uses the coeﬃcients at thehighest order available today, it corresponds to the NLO formula of the EMT. .1. Quenched QCD In this section we consider quenched QCD. We give two derivations. In Sect. 3.1.1, wegive a simple derivation using a characteristic of quenched QCD. In Sect. 3.1.2, we give analternative derivation, which is relatively complicated but can be generalized to full QCDstraightforwardly.

We can derive Eqs. (2.3) and (2.4) in a very simple manner. Thetraceless part is given by T TL µν ( x ) = c ( t ) ˜ O TL1 ,µν ( t, x ) . (3.1)By multiplying both sides by c (N k LO)1 ( t ) c ( t ) − , we obtain from Eq. (2.1) T TL µν (N k LO) ( x ; t ) = c (N k LO)1 ( t ) ˜ O ,µν ( t, x ) = c (N k LO)1 ( t ) c ( t ) − T TL µν ( x ) . (3.2)Here, c (N k LO)1 ( t ) c ( t ) − deviates from one owing to the lack of higher-order corrections. Ifwe write c ( t ) − perturbatively as c ( t ) − = g ( µ ( t )) ∞ X i =0 p ( i )1 ( L ( µ ( t ) , t )) (cid:20) g ( µ ( t )) (4 π ) (cid:21) i , (3.3) p ( i )1 satsﬁes k (0)1 p (0)1 = 1 , X i + j = n k ( i )1 p ( j )1 = 0 for n ≥ . (3.4)Then noting that c (N k LO)1 = g ( µ ( t )) P kn =0 k ( n )1 ( L ( µ ( t ) , t )) h g ( µ ( t )) (4 π ) i n , we obtain c (N k LO)1 ( t ) c ( t ) − ≃ X i + j = k +10 ≤ i ≤ k k ( i )1 p ( j )1 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 = 1 − k ( k +1)1 ( L ( µ ( t ) , t )) (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 . (3.5)In the last line, we have used Eq. (3.4) and p (0)1 = ( k (0)1 ) − = 1. Combining this and Eq. (3.2),we obtain Eq. (2.3).The scalar part T S ( x ) = c ( t ) ˜ O S2 ( t, x ) (3.6)can be investigated in a parallel manner. In this case, one should note that there is no “ k (0)2 ”. We have already clariﬁed the t -dependence in the N k LO formulain Derivation I. However, Derivation I cannot be applied straightforwardly to full QCDbecause the traceless and scalar parts of the EMT are given by linear combinations of ﬂowedoperators in full QCD. Hence, we consider another derivation which is applicable to fullQCD. This is Derivation II. 10e consider the diﬀerence between the EMT and that of the N k LO formulae: T TL µν ( x ) − T TL µν (N k LO) ( x ; t ) = [ c ( t ) − c (N k LO)1 ( t )] ˜ O TL1 ,µν ( t, x ) , (3.7) T S ( x ) − T S(N k LO) ( x ; t ) = [ c ( t ) − c (N k LO)2 ( t )] ˜ O S2 ( t, x ) . (3.8)Because T TL µν ( x ) and T S ( x ) are t independent, the t -dependence of the N k LO formulae isexhibited by the right-hand sides of Eqs. (3.7) and (3.8). Here, the t -dependence comingfrom the diﬀerence between the coeﬃcients c i ( t ) are easily evaluated as c ( t ) − c (N k LO)1 ( t ) ≃ g ( µ ( t )) k ( k +1)1 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 , (3.9) c ( t ) − c (N k LO)2 ( t ) ≃ g ( µ ( t )) k ( k +2)2 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +2 . (3.10)Thus, the remaining task is to investigate t -dependence of the ﬂowed operators ˜ O TL1 ,µν ( t, x )and ˜ O S2 ( t, x ). From Eq. (1.47), we have ˜ O µν ( t, x ) = ζ R ( t ; g ( µ ) , µ ) ˆ O Rµν ( x ; µ ) . (3.11)As explained in the introduction, we consider a t -dependent renormalization scale, µ = µ ( t ).Accordingly, we set µ = µ ( t ):˜ O µν ( t, x ) = ζ R ( t ; g ( µ ( t )) , µ ( t )) ˆ O Rµν ( x ; µ ( t )) . (3.12)Note that t -dependence is exhibited not only by ζ R ( t ; g ( µ ( t )) , µ ( t )) but also by ˆ O Rµν ( x ; µ ( t )).We can obtain the former matrix, ζ R ( t ; g ( µ ( t )) , µ ( t )), by renormalizing the bare results ζ ( t )in Refs. [2, 3, 9] according to Eq. (1.48). (As explained below, however, for the presentpurpose, we just need the tree-level result, which can be trivially obtained.)Let us investigate the leading t -dependence of ˆ O Rµν ( x ; µ ( t )). Since the bare operatorˆ O µν ( x ) = Z ( g ( µ )) ˆ O Rµν ( x ; µ ) (3.13)is independent of the renormalization scale, we have renormalization group (RG) equationfor the renormalized operators: (cid:20) µ ddµ + γ ( g ( µ )) (cid:21) ˆ O Rµν ( x ; µ ) = 0 , (3.14)where γ ( g ( µ )) is the anomalous dimension matrix, γ ( g ) = Z − µ ddµ Z = Z − β ( g ) ddg Z. (3.15)The solution to Eq. (3.14) is given byˆ O Rµν ( x ; µ ) = K ( µ ; µ ) ˆ O Rµν ( x ; µ ) , (3.16) In quenched QCD, t dependence of these operators can actually be revealed easily. Since T TL µν ( x ) = c ( t ) ˜ O TL1 ,µν ( t, x ) is t -independent, the leading t -dependence is ˜ O TL1 ,µν ( t, x ) = c ( t ) − T TL µν ( x ) ≃ g ( µ ( t )) T TL µν ( x ). Similarly, ˜ O S ( t, x ) ≃ (4 π ) k (1)2 T S ( x ) follows. However, again for the reason that thisargument cannot straightforwardly be applied to full QCD, we develop an argument not essentiallyrelying on a characteristic of quenched QCD. K ( µ ; µ ) = P exp " − Z g ( µ ) g ( µ ) dx γ ( x ) β ( x ) . (3.17)Here, P denotes the path-ordered product. To summarize, the ﬂowed operator can be written as˜ O µν ( t, x ) = ζ R ( t ; g ( µ ( t )) , µ ( t )) K ( µ ( t ); µ ) ˆ O Rµν ( x ; µ ) . (3.18)The t -dependence of the ﬂowed operaor is exhibited by ζ R ( t ; g ( µ ( t )) , µ ( t )) and K ( µ ( t ); µ ).Then, we calculate these two matrices. Since we are interested in the leading behavior, weconsider ζ R and K at LO.At LO, ζ R [whose deﬁnition is given in Eq. (1.48)] is trivially given by ζ R ( t ; g ( µ ) , µ ) = g ( µ ) g ( µ ) ! , (3.19)and K ( µ ; µ ) is given by K ( µ ; µ ) = (cid:20)(cid:16) g ( µ ) g ( µ ) (cid:17) − (cid:21) (cid:16) g ( µ ) g ( µ ) (cid:17) , (3.20)from the anomalous dimension matrix of Eq. (C2).From these results, we obtain the t -dependence of the ﬂowed operators at LO as˜ O TL1 ,µν ( t, x ) = g ( µ ( t )) (cid:20) ˆ O R ,µν ( x ; µ ) −

14 ˆ O R ,µν ( x ; µ ) + O ( g ( µ ) ) (cid:21) + O ( g ( µ ( t )) ) , (3.21)˜ O S2 ( t, x ) = [ g ( µ ) + O ( g ( µ ) )] ˆ O R, S2 ( x ; µ ) + O ( g ( µ ( t )) ) . (3.22)Thus, from Eqs. (3.7) and (3.8), we conclude T TL µν (N k LO) ( x ; t ) = T TL µν ( x ) + O ( g ( µ ( t )) k +1) ) , (3.23) T S(N k LO) ( x ; t ) = T S ( x ) + O ( g ( µ ( t )) k +1) ) . (3.24)We have revealed the leading dependence on g ( µ ( t )), which agrees to Derivation I.However, these are not as precise as Eqs. (2.3) and (2.4). In fact, we can reproduceEqs. (2.3) and (2.4) as follows. First, let us consider the traceless part. Actually, thequantity inside the square brackets in Eq. (3.21) can be written as ˆ O R ,µν ( x ; µ ) − (1 + O ( g ( µ ) )) ˆ O R ,µν ( x ; µ ); the higher-order eﬀects in g ( µ ) appear only in the coeﬃcient ofˆ O R ,µν and the coeﬃcient of ˆ O R ,µν is exactly one. This follows from K ( µ ; µ ) = 1 at anyorder of perturbation theory, because of Z = 1 and Z = 0 to all orders [see App. B, inparticular Eq. (B12))]. [Note that g ( µ ) dependence comes only from K ( µ, µ ) and not from ζ R ( t ; g ( µ ( t )) , µ ( t )).] Now we note that ˜ O TL1 ,µν ( t, x ) is traceless. Then ˆ O R ,µν ( x ; µ ) − (1 + More speciﬁcally, the path-order product orders a product as follows: the operator whose variableis closest to g ( µ ) is brought to the most left side, and the second closest one is to the second locationfrom left, and so on. ( g ( µ ) )) ˆ O R ,µν ( x ; µ ) should be also traceless. Therefore (1 + O ( g ( µ ) )) ˆ O R ,µν ( x ; µ )should be ( δ µν /

4) ˆ O R, S1 ( x ; µ ). Using Eq. (B15), we then conclude˜ O TL1 ,µν ( t, x ) = g ( µ ( t )) (cid:20) ˆ O R ,µν ( x ; µ ) − δ µν O R, S1 ( x ; µ ) (cid:21) + O ( g ( µ ( t )) )= g ( µ ( t )) T TL µν ( x ) + O ( g ( µ ( t )) ) . (3.25)This precisely gives Eq. (2.3).Secondly, let us consider ˜ O S ( t, x ). The g ( µ ( t )) -term of ˜ O S ( t, x ), whose LO result is shownin Eq. (3.22), is given by f ( g ( µ )) O R, S2 ( x ; µ ) beyond LO with some function f . This shouldbe µ independent. Hence, f ( g ( µ )) should be f ( g ( µ )) = const. × ( − β ( g ( µ ))8 g ( µ ) ). Here we haveused the fact that the trace of the EMT, which is given by Eq. (B13), is µ independent.From the LO result (3.22), the constant is determined as const. = π ) β . Thus, we canrewrite Eq. (3.22) as˜ O S2 ( t, x ) = 8(4 π ) β (cid:18) − β ( g ( µ ))8 g ( µ ) (cid:19) ˆ O R, S2 ( x ; µ ) + O ( g ( µ ( t )) )= 8(4 π ) β T ρρ ( x ) + O ( g ( µ ( t )) ) . (3.26)This and Eq. (3.8) give Eq. (2.4); note that k (1)2 = β / We investigate the leading t -dependence of the N k LO formulae of the EMT in full QCD ina manner parallel to Derivation II. The N k LO formulae diﬀer from the exact ones by T TL µν ( x ) − T TL µν (N k LO) ( x ; t ) = [ c ( t ) − c (N k LO)1 ( t )] ˜ O TL1 ,µν ( t, x ) + [ c ( t ) − c (N k LO)3 ( t )] ˜ O TL3 ,µν ( t, x ) , (3.27) T S ( x ) − T S(N k LO) ( x ; t ) = [ c ( t ) − c (N k LO)2 ( t )] ˜ O S2 ( t, x ) + [ c ′ ( t ) − c ′ k LO) ( t )] ˜ O S4 ( t, x ) , (3.28)where c ( t ) − c (N k LO)1 ( t ) ≃ g ( µ ( t )) k ( k +1)1 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 , (3.29) c ( t ) − c (N k LO)2 ( t ) ≃ g ( µ ( t )) k ( k +2)2 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +2 , (3.30) c ( t ) − c (N k LO)3 ( t ) ≃ k ( k +1)3 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 , (3.31) c ′ ( t ) − c ′ (N k LO) ( t ) ≃ k ( k +1)4 (cid:20) g ( µ ( t )) (4 π ) (cid:21) k +1 . (3.32)Then, we investigate the leading t -dependence of ˜ O TL1 ,µν ( t, x ), ˜ O TL3 ,µν ( t, x ), ˜ O S2 ( t, x ),and ˜ O S4 ( t, x ). 13he t -dependence of the ﬂowed operators can be investigated from˜ O µν ( t, x ) = ζ R ( t ; g ( µ ( t )) , µ ( t )) K ( µ ( t ); µ ) ˆ O Rµν ( x ; µ ) , (3.33)where K ( µ ; µ ) denotes K ( µ ; µ ) = P exp " − Z g ( µ ) g ( µ ) dx γ ( x ) β ( x ) . (3.34)The anomalous dimension matrix is deﬁned in a parallel manner to Eq. (3.14).At LO, ζ R ( t ; g ( µ ) , µ ) is given by ζ R ( t ; g ( µ ) , µ ) = g ( µ ) g ( µ ) . (3.35)Now we calculate K ( µ ; µ ) at LO, where we set γ = γ g / (4 π ) and β = − β g [ g / (4 π ) ] inEq. (3.34). The g -dependence of K is determined by eigenvalues of γ , which is given inEq. (C5). They are given by λ = − β , , (2 C F + T F ). (The eigenvalue 0 is degenerate.)14hen we obtain K = 2 C F C F + T F + T F C F + T F (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = 14 (cid:18) g ( µ ) g ( µ ) (cid:19) − C F C F + T F ) − T F C F + T F ) (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = C F C F + T F ) − C F C F + T F ) (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = 3 C F β (cid:18) g ( µ ) g ( µ ) (cid:19) − (cid:20) C F C F + T F ) + 3 C F β (cid:21) + C F C F + T F ) (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = 0 ,K = (cid:18) g ( µ ) g ( µ ) (cid:19) ,K = 0 ,K = 6 C F β (cid:18) g ( µ ) g ( µ ) (cid:19) − C F β ,K = 4 T F C F + T F − T F C F + T F (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = − T F C F + T F + T F C F + T F (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = T F C F + T F + 2 C F C F + T F (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = C F C F + T F − C F C F + T F (cid:18) g ( µ ) g ( µ ) (cid:19) CF + TF )3 β ,K = K = K = 0 ,K = 1 . (3.36)One can conﬁrm that this result satisﬁes required relations for K , given in App. E.From these matrices, t -dependence of the ﬂowed operators of interest is given by˜ O TL1 ,µν ( t, x )= g ( µ ( t )) (cid:20) C F C F + T F (cid:18) ˆ O R ,µν ( x ; µ ) −

14 ˆ O R ,µν ( x ; µ ) + 14 ˆ O R ,µν ( x ; µ ) −

18 ˆ O R ,µν ( x ; µ ) (cid:19) + O ( g ( µ ) ) (cid:21) , + O ( g ( µ ( t )) ( g ( µ ( t )) /g ( µ ) ) CF + TF )3 β ) (3.37)˜ O TL3 ,µν ( t, x )= (cid:20) T F C F + T F (cid:18) ˆ O R ,µν ( x ; µ ) −

14 ˆ O R ,µν ( x ; µ ) + 14 ˆ O R ,µν ( x ; µ ) −

18 ˆ O R ,µν ( x ; µ ) (cid:19) + O ( g ( µ ) ) (cid:21) + O (( g ( µ ( t )) /g ( µ ) ) CF + TF )3 β ) , (3.38)15nd˜ O S2 ( t, x ) = g ( µ ) (cid:18) ˆ O R, S2 ( x ; µ ) + 6 C F β ˆ O R, S4 ( x ; µ ) + O ( g ( µ ) ) (cid:19) + O ( g ( µ ( t )) ) , (3.39)˜ O S4 ( t, x ) = ˆ O R, S4 ( x ; µ ) + O ( g ( µ ) ) + O ( g ( µ ( t )) ) . (3.40)Hence, we obtain T TL µν (N k LO) ( x ; t ) = T TL µν ( x ) + O ([ g ( µ ( t )) ] k +1 ) , (3.41) T S(N k LO) ( x ; t ) = T S ( x ) + O ([ g ( µ ( t )) ] k +1 ) . (3.42)We have obtained main results, Eqs. (2.5) and (2.6).For systematic calculation, for instance to study higher-order eﬀects, it would be convenientto decompose dimension-four operators into traceless parts and scalar operators. Then it isenough to treat two-by-two matrices, instead of four-by-four matrices, for each vector space.We discuss higher-order eﬀects in this treatment in Appendix G. In addition, in Appendix Gwe explain a systematic way to estimate the errors of the LO results, which are shown bythe symbol O ( . . . ) in Eqs. (3.37)–(3.40).One might wonder if the t dependent term of the N k LO formulae for the EMT is propor-tional to the EMT as in the quenched QCD case. However, this is not necessarily expectedin full QCD. The reason why we obtained such a result in quenched QCD can be explainedas follows. In quenched QCD, T TL µν ( x ) ( T S ( x )) is the only traceless (scalar) operator which isRG invariant. Then, ˜ O TL1 ,µν ( t, x ) ( ˜ O S2 ( t, x )), which is traceless (scalar) and RG invariant (ormore precisely µ independent), should be proportional to T TL µν ( T S ). Hence, the t dependentterm of the N k LO EMT expression is proportional to the exact EMT. This is the conclusionof Sect. 3.1.2. On the other hand, in the full QCD case, considering scalar operators, wehave two RG invariant operators: T S and ˆ O R . Then it is not necessary that the RG invari-ant ﬂowed operators ˜ O S2 and ˜ O S4 are proportional to T S . This expectation can be deniedmore explicitly by considering the diﬀerence between LO and NLO formula for T S . It is notproportional to [1 / (4 π ) ] (cid:0) C A − T F (cid:1) ˜ O S + ˜ O S , which is T S as t → c ( t ) and c ′ ( t ). O ( t ) correction: Contribution from dimension-six operators So far, we have neglected the contribution from dimension-six operators in the small ﬂow timeexpansion, which is suppressed for small t roughly as O ( t ). This is parameterically smallerthan O ( g ( µ ( t )) n ) terms, which we have obtained in Sect. 3, in suﬃciently small t region.Nevertheless, it might be possible that the O ( t ) eﬀect dominates over the O ( g ( µ ( t )) n ) eﬀectin the region of t where lattice simulation is practically carried out. Hence, as an additionalstudy, we discuss how we can detect the detailed t dependence coming from dimension-sixoperators, and show an explicit result of the detailed t dependence for a speciﬁc example.We ﬁrst give a general argument how to investigate the detailed t dependence comingfrom dimension-six operators. We go back to the small ﬂow time expansion of the ﬂowedoperators ˜ O i,µν ( t, x ):˜ O i,µν ( t, x ) = ζ Rij ( t ; g ( µ ) , µ ) ˆ O Rj,µν ( x ; µ ) + tη ij ( t ) O (6) j,µν ( x ) + O ( t ) , (4.1)where the second term, which was neglected in the previous section, is now focused; O (6) i,µν is abare dimension-six operator and η ij is a coeﬃcient matrix. We denote the O ( t ) contribution16y δ ˜ O i,µν ( t, x ) ≡ tη ij ( t ) O (6) j,µν ( x ) . (4.2)Hereafter, the symbol δ means an O ( t ) contribution. The gradient ﬂow representation of theEMT [see Eq. (1.1)], T GF µν ( x ; t ) = X i ˜ c i ( t ; g ( µ )) ˜ O i,µν ( t, x ) , (4.3)diﬀers from the actual EMT by the dimension-six operators δ ˜ O i,µν ( t, x ), which isroughly O ( t ). (We are now assuming that ˜ c i ’s are exactly known and this part does notinduce any error.) This is the reason why we show the superscript “GF” and t dependence.The diﬀerence is given by δT GF µν ( x ; t ) = T GF µν ( x ; t ) − T µν ( x ) = X i ˜ c i ( t ; g ( µ )) δ ˜ O i,µν ( t, x ) . (4.4)We investigate the leading t dependence of the right-hand side. For this purpose, as before,we study the t dependence of δ ˜ O i,µν ( t, x ) by rewriting it as δ ˜ O i,µν ( t, x ) = tη Rij ( t ; g ( µ ) , µ ) O (6) Rj,µν ( x ; µ )= tη Rij ( t ; g ( µ ( t )) , µ ( t )) K (6) jk ( µ ( t ) , µ ) O (6) Rk,µν ( x ; µ ) . (4.5)In the ﬁrst equality, we have used the fact that the ﬂowed operator is ﬁnite and thus can berewritten by renormalized (ﬁnite) quantities. The product of the two renormalized quantitiesis independent of the renormalization scale µ . In the second equality, we have considered theRG evolution of the dimension-six renormalized operators, (cid:20) µ ddµ + γ (6) ( g ( µ )) (cid:21) O (6) Rµν = 0 , (4.6)and K (6) is given by K (6) ( µ ; µ ) = P exp " − Z g ( µ ) g ( µ ) dx γ (6) ( x ) β ( x ) . (4.7)To reveal the t dependence, we need to know η Rij ( t ; g ( µ ( t )) , µ ( t )) and K (6) jk ( µ ( t ) , µ ).To obtain the leading t dependence concretely, we need to know the matrix η ij at leadingorder and the anomalous dimension matrix of dimension-six operators at one-loop. Theformer can be calculated without loop calculations and by using the ﬂow equation alone.(We will show this explicitly in quenched QCD below.) However, it is rather complicated toknow the anomalous dimension matrix γ (6) , whose results are not completely known untiltoday. Here, we limit ourselves to the traceless part of the EMT in quenched QCD. Wedemonstrate how we can calculate the leading t dependence coming from the dimension-sixoperators for this case.We consider T GF , TL µν ( x ; t ) in quenched QCD [ T GF , TL µν ( x ; t ) = c ( t ; g ( µ ( t ))) ˜ O TL1 ,µν ( t, x )] andinvestigate the detailed t dependence caused by dimension-six operators. We calculate η ij In this analysis, it is suﬃcient to know ˜ c i ( t ; g ( µ )) at LO.

17y studying the O ( t ) contribution of tr ( G µρ G νρ ). We will represent the O ( t ) contributionof the ﬂowed operator tr ( G µρ G νρ ) using the basis O (6)1 ,µν = 1 g tr( D ν F ρσ D µ F ρσ ) + ( µ ↔ ν ) , O (6)2 ,µν = 1 g tr( D ρ F µρ D σ F νσ ) + ( µ ↔ ν ) , O (6)3 ,µν = 1 g tr( D ν F µρ D σ F ρσ ) + ( µ ↔ ν ) , (4.8)taking into account that a similar basis is adopted in Ref. [28]. We use δG µν = δ ( ∂ µ B ν − ∂ ν B µ + [ B µ , B ν ])= ∂ µ δB ν + [ A µ , δB ν ] − ∂ ν δB µ − [ A ν , δB µ ]= D µ δB ν − D ν δB µ = t ( D µ D ρ F ρν − D ν D ρ F ρµ ) , (4.9)which follows from the ﬂow equation ∂ t B aµ ( t, x ) = D ν G aνµ ( t, x ) in Eq. (1.8) . Then we have δ tr( G µσ G νσ ) = tr( δG µσ · F νσ ) + ( µ ↔ ν )= t tr( D ρ F σρ D µ F νσ ) − t tr( D ρ F µρ D σ F νσ ) − t ∂ µ tr( D ρ F σρ · F νσ ) + t ∂ σ tr( D ρ F µρ · F νσ )+ ( µ ↔ ν ) (4.10)Hereafter, we neglect the derivative terms. This is valid when we consider the matrix elementwith the states which are translational invariant. Then, we obtain δ tr( G µσ G νσ ) = − tg O (6)2 ,µν + tg O (6)3 ,µν . (4.11)We consider its traceless part and also renormalization: δ tr([ G µσ G νσ ] TL ) = − tg ( µ ( t )) O (6) R, TL2 ,µν ( x ; µ ( t )) + tg ( µ ( t )) O (6) R, TL3 ,µν ( x ; µ ( t )) , (4.12)where we take the renormalization scale µ = µ ( t ).Now we calculate K (6) matrix. The renormalization factor for {O (6) R, TL1 ,µν , O (6) R, TL2 ,µν , O (6) R, TL3 ,µν } was calculated in Ref. [28]: O (6) TL i,µν = Z ij O (6) R, TL j,µν , (4.13)with Z = − C A ǫ g (4 π ) C A ǫ g (4 π ) − C A ǫ g (4 π ) − C A ǫ g (4 π ) − C A ǫ g (4 π ) − C A ǫ g (4 π ) − C A ǫ g (4 π ) . (4.14)From this, we obtain the anomalous dimension: (cid:18) µ ddµ + γ (6) (cid:19) O (6) R, TL µν = 0 , (4.15) In right-hand side of Eq. (1.8), we have α × (gauge non-invariant operator). However, this partdoes not contribute to the ﬂow time evolution of gauge invariant operators. γ (6) = −

23 163

83 13

43 73 C A g (4 π ) + O ( g ) . (4.16)By calculating K (6) = P exp h − R g ( µ ) g ( µ ) dx γ (6) ( x ) β ( x ) i at LO, we obtain the RG evolution of therenormalized operators: O (6) R, TL2 ,µν ( x ; µ ( t )) = 1 √ √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ + − √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ ! O (6) R, TL2 ,µν ( x ; µ )+ 1 √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ − (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ ! O (6) R, TL3 ,µν ( x ; µ ) , (4.17) O (6) R, TL3 ,µν ( x ; µ ( t )) = 4 √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ − (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ ! O (6) R, TL2 ,µν ( x ; µ )+ 12 √ ( − √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ + (1 + √ (cid:18) g ( µ ( t )) g ( µ ) (cid:19) λ ! × O (6) R, TL3 ,µν ( x ; µ ) , (4.18)where λ , , are eigenvalues of the matrix C A γ (6)0 , λ = , λ = (15 + √ , λ = (15 − √ λ ≈ . λ ≈ . λ ≈ .

49. ( λ does notappear in the result.)From Eqs. (4.12), (4.17) and (4.18), we obtain the leading behavior as δ tr([ G µρ G νρ ] TL ) ≃ Cg ( µ ( t )) λ t, (4.19)with a (dimension-six) constant C , because λ is the smallest eigenvalue. Finally fromEq. (4.4) we obtain δT GF , TL µν = c ( t ) δ ˜ O TL1 ,µν ≃ C ′ g ( µ ( t )) λ t, (4.20)with a (dimension-six) constant C ′ .

5. Numerical analysis of thermodynamics quantities in quenched QCD

In this section, we carry out lattice simulation of the EMT with the SF t X method in quenchedQCD. We study thermodynamic quantities [29–37]. We use the t → t → t was mainly used.We use the lattice data obtained in Ref. [13] for ﬂowed operators. See this reference forthe details of the lattice setup. We study the ﬁnite temperature eﬀect, i.e., the diﬀerence The authors are grateful to Takumi Iritani and Masakiyo Kitazawa for letting us use the data. s/T and trace anomaly ∆ /T , which are given by sT = ǫ + p = − T TL44 , (5.1)∆ = ǫ − p = − T S . (5.2)We measure these quantities at temperature T /T c = 0 .

93, 1 .

02, 1 .

12, 1 .

40, 1 .

68, 2 .

10, 2 . .

69, where T c denotes the critical temperature. Our lattice analysis consists of three steps. First, we carry out a → O TL1 , and ˜ O S2 . We obtain these results as a function of tT . In order to keep discretizationeﬀects under good control we need 2 e γ E t/a ≫

1. (Here we assume that the typical scaleof the ﬂowed operators is 1 / √ e γ E t .) Noting T = 1 / ( N τ a ) and N τ ≥

12 (where N τ is thenumber of the sites in Euclidean time direction) in our lattice setup, we see that the abovecondition corresponds to tT ≫ . tT ≥ . c i ( t ) [and trivial factors in Eqs. (5.1) and (5.2)]such that they correspond to the thermodynamic quantities. We use NLO or NNLO coeﬃ-cients, namely we use the NLO or NNLO formula of the EMT. Finally, we extrapolate thedata to zero ﬂow time limit using Eqs. (2.3) and (2.4) with k = 1 or 2. Then we obtain theﬁnal results.In the second and third steps, we rely on perturbation theory. This is valid when the ﬂowtime satisﬁes t ≪ Λ − . In Fig. 1, we show the size of the running coupling as a functionof tT in the relevant region, taking the renormalization scale of the running coupling as µ ( t ) = 1 / (2 e γ E t ) / . Our analysis is mainly performed for tT ≤ . tT , which is originally a function of µ ( t ) / Λ MS , we use [13] w T c = 0 . , w Λ MS = 0 . , (5.3)where w is a reference scale. (Here we only need the ratio T c / Λ MS .) Here and hereafter weuse the three-loop beta function.We summarize the setup of our central analysis:range of the used lattice data: 0 . ≤ tT ≤ . , scale setting parameters: Eq. (5.3) , renormalization scale: µ s ( t ) = s/ (2 e γ E t ) / with s = 1 . (5.4)We estimate systematic errors by varying the above conditions. In the systematic erroranalysis of the range, we extend the range to smaller t region and use the data at 0 . ≤ tT ≤ . w T c estimatedin Ref. [13], but we consider the error of Λ MS of 2 .

7% [38]. We also vary the renormalizationscale within s = 1 / √ The trace anomaly is studied only at

T /T c = 0 .

93, 1 .

02, 1 .

12, 1 .

40, and 1 .

68. This is due tolack of zero temperature simulations, which require more numerical costs. On the other hand, wedo not need zero temperature simulations for the entropy density because it is exactly zero at zerotemperature. ig. 1 The size of the running coupling constant g ( µ ( t )) / (4 π ) as a function of tT . It isshown for various T /T c .Now we present our numerical results. In Fig. 2, we show the t → k = 1. We regard the ﬁt parameter k (2)1 independent oftemperature and it is taken common to all the simulation temperatures. A characteristic ofthe extrapolation function (2.3) is that it rises sharply around t ∼

0. This stems from thesingularity of g ( µ ( t )) ∼ / ( t Λ )) at t = 0. The extrapolation function is consistent withthe (non-trivial) behavior of the lattice data at small t region, tT = 0 . . . .

015 in the ﬁt. It indicates validity of the extrapolationfunction.We can also conﬁrm the validity of the extrapolation function by comparing the result ofthe ﬁt parameter k (2)1 , which corresponds to the NNLO coeﬃcient for c ( t ), with its alreadyknown result. From the ﬁt, we obtain k (2)1 = 88(18) , (5.5)where the value inside the parentheses denotes the statistical error, and the exact NNLOresult is k (2)1 = 87 . ... [9]. Our result agrees well with the exact result.Even in the case where the higher order perturbative coeﬃcient k ( k +1)1 is not known, thefollowing analysis is possible to check the validity of the use of the extrapolation function(2.3). We focus on the property that the L ( µ, t ) = log(2 µ e γ E t ) dependence of the N k +1 LOperturbative coeﬃcient is totally determined by the perturbative coeﬃcients up to N k LOand the beta function; see App. D. Then we can compare the L ( µ, t ) dependence of the ﬁtparameter k ( k +1)1 with the predicted dependence on L ( µ, t ).In the left panel of Fig. 3, we show the log s (= L ( µ s ( t ) , t )) dependence of the ﬁt parameter k (2)1 . It indeed exhibits similar dependence to that of blue lines, which show the exact log s dependence of k (2)1 , predicted from the RG equation; see the third equation of Eq. (D2).We assume diﬀerent log independent constants k (2)1 ( L = 0) for three blue lines. For the solidline we set k (2)1 ( L = 0) in Eq. (D2) to k (2) , ﬁt1 ( L = 0), which is the obtained ﬁt parameterin our central analysis with s = 1 (or log ( s ) = 0) [i.e., the central value of Eq. (5.5)]. Forthe dashed line below, we set k (2)1 ( L = 0) in Eq. (D2) such that k (2)1 ( L = log (2 )) coincides21 ig. 2 NLO analysis of the entropy density. The blue solid lines are the t → k = 1. The error bars show statistic errors only. (The same appliesto Figs. 3-7.) The gray dashed lines show the ﬁt range used in the t → ig. 3 The ﬁt parameter k (2)1 or k (3)1 as a function of log s . The data points with error barsshow the results of the ﬁt parameter k (2)1 or k (3)1 obtained in the t → s . The blue lines show log s dependence of the perturbative coeﬃcients k (2)1 and k (3)1 dictated from the RG equation; see Eq. (D2). We assume diﬀerent log s independentconstants for three blue lines; see the main text.with k (2) , ﬁt1 ( L = log (2 )). Here k (2) , ﬁt1 ( L = log (2 )) is the obtained ﬁt parameter when weset s = 2. The diﬀerence in k (2)1 ( L = 0) between the two analyses can be regraded as anerror of the estimate of k (2)1 ( L = 0). For the dashed line above, the same size variation of k (2)1 ( L = 0) is assumed with the reversed sign.We estimate systematic errors. We vary one of the following conditions from the centralanalysis: the ﬁt range, Λ MS , and renormalization scale. We already explained how we changethese conditions above. We examine the variations from the central values caused by changinga condition, and the (largest) diﬀerence is given as a systematic error. Our ﬁnal result of theNLO analysis of the entropy density is summarized in Table 1.We move on to the NNLO analysis of the entropy density. In Fig. 4, one can see that the t dependence of the data become milder at NNLO [21]. This can be naturally understood fromour theoretical study in Sect. 3; the remaining t dependence of the NNLO formula is given by O ( g ( µ ( t )) ) at NNLO, whereas O ( g ( µ ( t )) ) at NLO. Actually, the data can be regarded asﬂat within the statistical errors. Then ﬁnal results are insensitive to extrapolation functions.We again check the validity of the extrapolation function of Eq. (2.3) (with k = 2) inthe right panel of Fig. 3. As explained above, we examine the log s dependence of the ﬁtparameter k (3)1 . Although statistical errors are large, the behavior of the central values lookshighly consistent with the expected behavior. (In drawing the blue lines in Fig. 3, we usedthe exact value of k (2)1 .) Hence, we consider the extrapolation function to be valid also inthe NNLO analysis.We summarize the NNLO result of the entropy density in Table 1 with systematic errors.We make some comments. First, the NLO and NNLO analyses are mutually consistent. It isworth noting that the variations of the central values are quite small. Here, the extrapolationfunction (2.3) plays an important role. For comparison, see App. H. Secondly, at NNLO, The largest diﬀerence in the estimate for k (2)1 ( L = 0) is caused by the s = 1 and s = 2 cases.This is why we focused on s = 2 as the error estimate. ig. 4 NNLO analysis of the entropy density. The blue solid lines are the t → k = 2. The gray dashed lines show the ﬁt range used in the t → t X method enables us to performsystematic and accurate analyses of the EMT. At NNLO, dominant uncertainties come fromthe statistical error and the error of Λ MS . 24 ig. 5 NLO analysis of the trace anomaly. The blue solid lines are the t → k = 1. The gray dashed lines show the ﬁt range used in the t → k (2)1 and k (3)1 : k (2)1 ( L = 0) = 88(18)(2)(29)(114) , (5.6) k (3)1 ( L = 0) = − . (5.7)We note that we already know k (2)1 ( L = 0) exactly and the above estimate is consistentwith it (as already mentioned above). The values inside the parentheses are statistical andsystematic errors, which are shown in the same order as in that of Table 1.We also perform the NLO and NNLO analyses for the trace anomaly; we show themin Figs. 5 and 6. The results are summarized in Table 1. We check the validity of theperturbative extrapolation function in Fig. 7. Since the situation is almost parallel to thecase of the entropy density, we do not repeat the explanation. However, we note that theerror of Λ MS does not induce a dominant error in the trace anomaly.25 ig. 6 NNLO analysis of the trace anomaly. The blue solid lines are the t → k = 2. The gray dashed lines show the ﬁt range used in the t → Fig. 7

The ﬁt parameter k (3 , as a function of log s . See the caption in Fig. 3.26 /T T /T c NLO N LO0.93 0 . . . . . . . . . . . . . . . . /T T /T c NLO N LO0.93 0 . . . . . . . . . . Table 1

NLO and NNLO results for the entropy density and trace anomaly. The val-ues inside parentheses show errors associated with (statistic)(range)(Λ MS )(renormalizationscale) in this order. The values inside square brackets show total errors, which are given bycombining all errors in quadrature.We give our estimate of the perturbative coeﬃcients k (3)2 and k (4)2 : k (3)2 ( L = 0) = − , (5.8) k (4)2 ( L = 0) = − . (5.9)The exact result of the NLO coeﬃcient is known to be k (3)2 ( L = 0) = − . ... [21], and weﬁnd a good agreement.Finally, we compare our results with recent precise studies. Our results for the entropy den-sity are consistent with Refs. [36] and [37]. For the trace anomaly, our results are consistentwith Ref. [37] but not consistent with Ref. [31] within our and their ﬁnal errors.We note that our results for the both quantities are consistent with Ref. [21] and that weobtained smaller errors in this paper due to reduction of systematic errors. One of the keyelements to the smaller errors is that in the present analysis we do not need to consider theuncertainty associated with the functional form of a t → Although in Ref. [21] the “NNLO” coeﬃcient c ( t ) meant the result up to two-loop order, weregard c ( t ) up to two-loop order as the NLO coeﬃcient in this paper. Then the “NNLO” analysisof the trace anomaly in Ref. [21] corresponds to the NLO analysis in this paper and the “NNNLO”analysis of the trace anomaly in Ref. [21] to the NNLO analysis in this paper. Since there is no tree-level perturbative coeﬃcient in c ( t ), we notice that it is natural to call the one-loop order coeﬃcientLO for the trace part as in the present paper. A merit of adopting this order counting is that the t dependence of the N k LO formulae is O ( g ( µ ( t )) k +1) ) both for the traceless and scalar (or trace)parts. This convention is common to that in a recent paper [27]. We alsonote that we did not perform a ﬁt with the range 0 . ≤ tT ≤ . t region. Note that the way to estimate this systematic error is diﬀerent. In Ref. [21], µ ( t ) ∈ [ . √ t , . √ t ]is used and in the present paper µ ( t ) ∈ [ . √ t , . √ t ] is used. . Conclusions and discussion The SF t X method is a powerful method to simulate the EMT on lattice, and in this paper werevealed the functional form to be used in t → t X method. As anadditional theoretical study, we also revealed the detailed t dependence caused by dimension-six operators for the traceless part of the EMT in quenched QCD.We carried out numerical analyses of the thermodynamics quantities in quenched QCDwith our new extrapolation function (Sect. 5). Our extrapolation function is shown to bereasonable from the lattice data and it serves to reduce systematic errors, compared withthe conventionally used linear function in t . As a by-product, we also gave estimates of theNNNLO coeﬃcients k (3)1 and k (4)2 in quenched QCD.Our lattice analysis is carried out at suﬃciently small t region, where the coupling constantis given by ( α s =) g / (4 π ) . .

25. In full QCD, it might be diﬃcult to obtain lattice dataat such small t region because lattice spacing tends to be larger. (We noted in Sect. 5 thatthe extrapolation functions do not suﬃciently look reasonable at larger t region.) In thecase where the extrapolation functions (2.5) and (2.6) do not look consistent with availablelattice data, it would be diﬃcult to carry out t → t dependence (e.g. higher-order eﬀect in g ( µ ( t )) orlinear function in t ) in a ﬁtting function, a ﬁt would be destabilized due to the diﬃculty indistinguishing diﬀerent functions. One conservative attitude would be giving the systematicerror associated with extrapolation functions by trying some functions. However, such adiﬃculty may be systematically overcome by going to ﬁner lattice and/or going to higherorder of perturbation theory.Finally we mention that it would be possible to reveal proper t → t X methods (for other conserved currents) in a similar manner.

Acknowledgments

The authors thank Takumi Iritani and Masakiyo Kitazawa for fruitful discussions and forletting us use lattice data. This work was supported by JSPS Grant-in-Aid for ScientiﬁcResearch Grant Numbers, JP16H03982 and JP20H01903 (H.S.) and JP19K14711 (H.T.).

A. Convention

We set the normalization of anti-Hermitian generators T a of the representation R of thegauge group G as tr R ( T a T b ) = − T R δ ab and T a T a = − C R . We denote tr R (1) = dim( R ).From the structure constants deﬁned by [ T a , T b ] = f abc T c , we set f acd f bcd = C A δ ab . Forexample, for the fundamental N representation of G = SU ( N ) for which dim( N ) = N , ournormalization is C A = N, T R = 12 , C F = N − N . (A1)We also deﬁne T F ≡ n f T R = n f /

2, where n f is the number of quark ﬂavors.The D -dimensional Euclidean action of the vectorial gauge theory is given by S = Z d D x g F aµν ( x ) F aµν ( x ) + Z d D x ¯ ψ ( x )( / D + m ) ψ ( x ) . (A2)29he ﬁeld strength is deﬁned by F µν ( x ) = ∂ µ A ν ( x ) − ∂ ν A µ ( x ) + [ A µ ( x ) , A ν ( x )] , (A3)for A µ ( x ) = A aµ ( x ) T a and F µν ( x ) = F aµν ( x ) T a , where g is the bare gauge coupling and m is the bare mass parameter. The covariant derivative on the fermion is D µ = ∂ µ + A µ , (A4)and / D ≡ γ µ D µ , where γ µ denotes the hermitian Dirac matrix. B. EMT and renormalization

The EMT is written by bare quantities but is ﬁnite. This property gives non-trivial infor-mation on renormalization of composite operators, as explained in Ref. [2]. We revisit thissubject with the basis ˆ O µν .First, we consider quenched QCD. The EMT is given by T µν ( x ) = ˆ O ,µν ( x ) −

14 ˆ O ,µν ( x ) . (B1)The renormalization of the composite operators is carried out by ˆ O ,µν ( x )ˆ O ,µν ( x ) ! = Z ( g ( µ )) Z ( g ( µ )) Z ( g ( µ )) Z ( g ( µ )) ! ˆ O R ,µν ( x ; µ )ˆ O R ,µν ( x ; µ ) ! , (B2)where µ denotes the renormalization scale. Since the scalar operator is not mixed with thetensor operator, Z = 0. Then, the EMT is written as T µν ( x ) = Z ˆ O R ,µν ( x ) −

14 ( Z − Z ) ˆ O R ,µν ( x ) . (B3)Noting that this quantity is ﬁnite and noting also the property of the MS scheme that arenormalization factor only has divergent terms, negative powers in ǫ , beyond LO in g ( µ )[ Z ij = δ ij + P n ≥ ,m ≥ c nm ǫ − m g ( µ ) n ], we have Z = 1 , Z − Z = 1 . (B4)This leads to the EMT in terms of the renormalized operators as T µν ( x ) = ˆ O R ,µν ( x ) −

14 ˆ O R ,µν ( x ) . (B5)Now, let us consider the trace of the EMT: δ µν T µν ( x ) = (cid:18) − D (cid:19) g F aρσ F aρσ ( x ) = ǫ D ˆ O S2 ( x ) . (B6)From Eqs. (1.49) and (1.50), we can rewrite ǫ as ǫ = − β ( g ) (cid:18) g + 1 Z g ∂Z g ∂g (cid:19) . (B7)Then, we have δ µν T µν ( x ) = − β ( g )2 Dg Z (cid:18) gZ g ∂Z g ∂g (cid:19) ˆ O R, S2 ( x ) . (B8)30his quantity should be ﬁnite. Since β ( g ) / (2 Dg ) is ﬁnite as ǫ → Z (cid:18) gZ g ∂Z g ∂g (cid:19) (B9)should be ﬁnite. From the above property of the MS scheme, we have Z (cid:18) gZ g ∂Z g ∂g (cid:19) = 1 . (B10)Hence, we obtain Z = 11 + gZ g ∂Z g ∂g (cid:18) = − β ( g ) ǫg (cid:19) . (B11)To be summarized, we obtain Z = 1 , Z = −

14 + 14 Z ,Z = 0 , Z = 11 + gZ g ∂Z g ∂g . (B12)[ Z g is given by Eq. (1.53).] We also obtain the EMT in terms of the renormalized operators: T µν ( x ) = ˆ O R ,µν ( x ; µ ) −

14 ˆ O R ,µν ( x ; µ ) ,T ρρ ( x ) = − β ( g ( µ ))8 g ( µ ) ˆ O R, S2 ( x ; µ ) . (B13)From these results, we have another non-trivial relation. By taking the trace of the ﬁrstequation, it should coincide with the second one, which impliesˆ O R, S1 ( x ) = 14 (cid:18) − β ( g )2 g (cid:19) ˆ O R, S2 ( x ) . (B14)From this, we can write the traceless part of the EMT as T TL µν ( x ) = ˆ O R ,µν ( x ) − δ µν O R, S1 ( x )= ˆ O R ,µν ( x ) − (cid:18) − β ( g )2 g (cid:19) ˆ O R ,µν ( x ) , (B15)in terms of the renormalized operators.Next, we consider full QCD. In this case, we can obtain limited results but still a few usefulrelations. The EMT is given by T µν ( x ) = ˆ O ,µν ( x ) −

14 ˆ O ,µν ( x ) + 14 ˆ O ,µν ( x ) . (B16)By deﬁning a renormalization matrix in a parallel manner, we deduce that Z i − Z i + 14 Z i (B17)is ﬁnite. Using the property of the MS scheme, we obtain Z + 14 Z = 1 ,Z − Z + 14 Z = − ,Z + 14 Z = 14 , (B18)31here we have used Z = Z = 0. Then, the EMT is expressed with the renormalizedoperators as T µν ( x ) = ˆ O R ,µν ( x ; µ ) −

14 ˆ O R ,µν ( x ; µ ) + 14 ˆ O R ,µν ( x ; µ ) . (B19)Now, let us consider the trace part: δ µν T µν ( x ) = (cid:18) − D (cid:19) g F aµν F aµν ( x ) + 12 X f ¯ ψ f ( x ) ←→ / D ψ f ( x ) = ǫ D ˆ O S ( x ) + 12 D ˆ O S ( x ) . (B20)By rewriting it by the renormalized operators, we have δ µν T µν ( x ) = − β ( g )2 Dg Z (cid:18) gZ g ∂Z g ∂g (cid:19) ˆ O R, S2 ( x ) + 12 D (1 + ǫZ ) ˆ O R, S4 ( x ) . (B21)Here, we rewrote ǫ in the same way as above [Eq. (B7)] and used Z = 1 because O ,µν ,which is proportional to O ,µν via EOM, is a ﬁnite operator. Each coeﬃcient should be ﬁnitein Eq. (B21). Then, for the coeﬃcient of ˆ O R, S2 , we obtain the same result as the quenchedcase: Z (cid:18) gZ g ∂Z g ∂g (cid:19) = 1 , (B22)which indicates δ µν T µν ( x ) = − β ( g )8 g ˆ O R, S2 ( x ) + 12 D (1 + ǫZ ) ˆ O R, S4 ( x ) . (B23)As one can see, the coeﬃcient of the gluonic operator can be obtained in this way also infull QCD. For ˆ O R, S4 , we deduce that 1 + ǫZ should be ﬁnite. (Thus, Z has only simplepoles in ǫ to all orders.) However, it seems diﬃcult to ﬁx its ﬁnite value by this argumentalone. Here, we refer to another argument giving the trace part: δ µν T µν ( x ) = − β ( g )8 g ˆ O R, S2 ( x ) − (1 + γ m ) 14 ˆ O R, S5 ( x )= − β ( g ( µ ))8 g ( µ ) ˆ O R, S2 ( x ; µ ) + 1 + γ m ( g ( µ ))8 ˆ O R, S4 ( x ; µ ) . (B24)In the last line, we used the EOM, ˆ O R ,µν + ˆ O R ,µν = 0. By this, ǫZ = γ m follows. By takingthe trace of Eq. (B19) and comparing with the above expression, we obtainˆ O R, S1 ( x ) + 14 ˆ O R, S3 ( x ) = 14 (cid:18) − β ( g )2 g (cid:19) ˆ O R, S2 ( x ) + 18 (1 + γ m ) ˆ O R, S4 ( x ) . (B25)The traceless part of the EMT is given by T TL µν ( x ) = ˆ O R ,µν ( x ; µ ) − (cid:18) − β ( g ( µ ))2 g ( µ ) (cid:19) ˆ O R ,µν ( x ; µ )+ 14 ˆ O R ,µν ( x ; µ ) − γ m ( g ( µ ))8 ˆ O R ,µν ( x ; µ ) . (B26)32 . Anomalous dimension matrix for dimension-four operators In this appendix, we give the anomalous dimension matrix for dimension-four operators. Forthe deﬁnition, see Eq. (3.15). It has a form γ ( g ) = γ g (4 π ) + γ (cid:20) g (4 π ) (cid:21) + · · · . (C1)In quenched QCD, the anomalous dimension can be obtained from Eqs. (B12) and (1.53): γ = − β − β ! , (C2) γ = − β − β ! , (C3) γ = − β − β ! . (C4)In full QCD, the leading order matrix is given by [9] γ = T F − C A − C F − C F − β − C F − T F T F C F − C F . (C5)The NLO matrix is given by [9] γ , = 427 T F (35 C A + 74 C F ) , (C6) γ , = 227 ( − C A + 56 C A T F + 5 C F T F ) , (C7) γ , = 427 C F ( − C A + 14 C F + 26 T F ) , (C8) γ , = − C F (812 C A + 85 C F − T F ) , (C9) γ , = 0 , (C10) γ , = 83 ( − C A + 10 C A T F + 6 C F T F ) , (C11) γ , = 0 , (C12) γ , = − C F (97 C A + 9 C F − T F ) , (C13) γ , = − T F (35 C A + 74 C F ) , (C14) γ , = 827 T F (34 C A + 49 C F ) , (C15) γ , = − C F ( − C A + 14 C F + 26 T F ) , (C16) γ , = 427 C F ( − C A + 4 C F + 136 T F ) , (C17) γ , = γ , = γ , = γ , = 0 . (C18)As noted in footnote 5, our deﬁnition of the renormalization factor Z is the inverse of thatof Ref. [9]. 33 . Coeﬃcients c ( ′ ) i ( t ) The coeﬃcients c ( ′ ) i ( t ) are RG invariant and the L ( µ, t ) ≡ log(2 e γ E µ t ) dependence of theperturbative series is determined by the RG equation. For i = 1 ,

2, we have c i ( t ) = 1 g ( µ ) ∞ X n =0 k ( n ) i ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n , (D1)with k (0) i ( L ( µ, t )) = k (0) i ( L = 0) ,k (1) i ( L ( µ, t )) = k (1) i ( L = 0) − β k (0) i L ( µ, t ) ,k (2) i ( L ( µ, t )) = k (2) i ( L = 0) − β k (0) i L ( µ, t ) ,k (3) i ( L ( µ, t )) = k (3) i ( L = 0) + (cid:16) β k (2) i ( L = 0) − β k (0) i (cid:17) L ( µ, t ) − β β k (0) i L ( µ, t ) ,k (4) i ( L ( µ, t )) = k (4) i ( L = 0) + (cid:16) − β k (0) i + β k (2) i ( L = 0) + 2 β k (3) i ( L = 0) (cid:17) L ( µ, t )+ (cid:26)(cid:18) − β − β β (cid:19) k (0) i + β k (2) i ( L = 0) (cid:27) L ( µ, t ) − β β k (0) i L ( µ, t ) . (D2)For i = 3 ,

4, we have c ( ′ ) i ( t ) = ∞ X n =0 k ( n ) i ( L ( µ, t )) (cid:20) g ( µ ) (4 π ) (cid:21) n (D3)with k (0) i ( L ( µ, t )) = k (0) i ( L = 0) ,k (1) i ( L ( µ, t )) = k (1) i ( L = 0) ,k (2) i ( L ( µ, t )) = k (2) i ( L = 0) + β k (1) i L ( µ, t ) ,k (3) i ( L ( µ, t )) = k (3) i ( L = 0) + (cid:16) β k (1) i + 2 β k (2) i ( L = 0) (cid:17) L ( µ, t ) + β k (1) i L ( µ, t ) . (D4)The explicit NLO result is given by [2, 3] c ( t ) = 1 g + (cid:18) − C A + 32 T F − β L ( µ, t ) (cid:19) π ) , (D5) c ( t ) = 1(4 π ) (cid:18) C A − T F (cid:19) , (D6) c ( t ) = 14 + g (4 π ) C F (cid:18)

38 + log 2 + 34 log 3 (cid:19) , (D7) c ′ ( t ) = 18 + g (4 π ) C F (cid:18) (cid:19) . (D8)34 . Properties of the matrix K ( µ ; µ ) We list some relations that the matrix K should satisfy. Since K describes evolution, itshould satisfy K ( µ = µ ; µ ) = × ,K ( µ ; µ ) K ( µ ; µ ) = × . (E1)Using the fact that the EMT [Eq. (B19)] is renormalization scale independent, we obtain K + 14 K = 1 ,K − K + 14 K = − ,K + 14 K = 14 . (E2)Here we used K = K = 0 because ˆ O R ,µν is an essentially scalar operator and is not mixedwith ˆ O R ,µν or ˆ O R ,µν . Also the trace part of the EMT [Eq. (B24)] is renormalization scaleindependent, which leads to − β ( g ( µ ))8 g ( µ ) K ( µ ; µ ) = − β ( g ( µ ))8 g ( µ ) , − β ( g ( µ ))8 g ( µ ) K ( µ ; µ ) = 18 { γ m ( g ( µ )) − γ m ( g ( µ )) } . (E3)These are exact relations, not relying on perturbation theory. One can check that K ( µ ; µ )of Eq. (3.36) satisﬁes all the relations (at appropriate order). F. ζ R at NLO We show the NLO result for ζ R , deﬁned in Eq. (1.48). For ζ R ( t ; g ( µ ) , µ ) = g g + g (4 π ) ζ R (1) + · · · , (F1)35he NLO matrix ζ R (1) is given by ζ R (1)11 = C A g (cid:18)

73 + 113 L ( µ, t ) (cid:19) ,ζ R (1)12 = − C A g (cid:18)

16 + 1112 L ( µ, t ) (cid:19) ,ζ R (1)13 = − C F g (cid:18)

718 + 23 L ( µ, t ) (cid:19) ,ζ R (1)14 = − C F g (cid:18) L ( µ, t ) (cid:19) ,ζ R (1)21 = 0 ,ζ R (1)22 = 7 C A g ,ζ R (1)23 = 0 ,ζ R (1)24 = − C F g (5 + 6 L ( µ, t )) ,ζ R (1)31 = − T F (cid:18) L ( µ, t ) (cid:19) ,ζ R (1)32 = T F (cid:18) L ( µ, t ) (cid:19) ,ζ R (1)33 = C F (cid:18) − − L ( µ, t ) (cid:19) ,ζ R (1)34 = − C F (cid:18)

49 + 43 L ( µ, t ) (cid:19) ,ζ R (1)41 = 0 ,ζ R (1)42 = 53 T F ,ζ R (1)43 = 0 ,ζ R (1)44 = C F (cid:18) − − (cid:19) . (F2)We can conﬁrm validity of this NLO result for instance as follows. By rewriting˜ O TL1 ,µν ( t, x ) = ˜ O ,µν ( t, x ) − (1 /

4) ˜ O ,µν ( t, x ) in terms of ˆ O R ( x ; µ ) with the NLO ζ R ( t ; g ( µ ) , µ )and then requiring it is traceless, we obtainˆ O R, S1 = 14 ˆ O R, S2 + g (4 π ) (cid:20) C A ˆ O R, S2 + 712 C F ˆ O R, S4 (cid:21) + O ( g ) . (F3)Similarly, from ˜ O TL3 ,µν ( t, x ) = ˜ O ,µν ( t, x ) − (1 /

2) ˜ O ,µν ( t, x ), we obtainˆ O R, S3 = 12 ˆ O R, S4 + g (4 π ) (cid:20) − T F ˆ O R, S2 + 23 C F ˆ O R, S4 (cid:21) + O ( g ) . (F4)One can check that Eq. (B25) is correctly reproduced from these results.36 . Higher order correction to the matrix K In this appendix, we clarify how we speciﬁed the parametrical errors of the LO calculations inEqs. (3.37)–(3.40). For this purpose, we need to know what kind of higher-order correctionsappears in the matrices ζ R and K in Eq. (3.33). Here, we investigate this issue particularlyfor the matrix K , because its higher-order eﬀect is more diﬃcult to see than the matrix ζ R (whose higher-order eﬀect is just given by higher power in g ( µ ( t )) ).We ﬁrst develop a general argument to detect higher-order correction to the matrix K .The RG equation for K in Eq. (3.17) is given in the form µ dKdµ = − γ ( g ) K, (G1)or dKdA = − γ ( A )2 β A ( A ) K = (cid:20) γ β A + β γ − β γ β + O ( A ) (cid:21) K, (G2)with A ≡ g ( µ ) / (4 π ) and β A ( A ) ≡ − P ∞ i =0 β i A i +2 . Here and hereafter, we deﬁne γ , γ , . . . as γ ( A ) = γ A + γ A + O ( A ) . (G3) γ and γ are not commutative generally. The formal solution to the above RG equationis given by a path ordered product, but it is diﬃcult to evaluate explicitly. In fact, we canobtain the NLO K matrix as follows. Denoting the LO solution by K LO ( A ; A ) (with A = g ( µ ) / (4 π ) ), which satisﬁes ddA K LO ( A ; A ) = γ β A K LO ( A ; A ) , (G4)and writing the NLO solution as K NLO = K LO ˜ K , we have dK NLO dA = dK LO dA ˜ K + K LO d ˜ KdA = (cid:20) γ β A + β γ − β γ β + O ( A ) (cid:21) K LO ˜ K, (G5)and hence, d ˜ KdA = (cid:18) K − β γ − β γ β K LO (cid:19) ˜ K. (G6)Since K − β γ − β γ β K LO can be regarded as a function of A and A , we can obtain ˜ K as˜ K ( A ; A ) = exp (cid:20)Z AA dx (cid:20) K − β γ − β γ β K LO (cid:21) ( x ; A ) (cid:21) . (G7)Then the NLO K matrix is given by K NLO ( A ; A ) = K LO ( A ; A ) exp (cid:20)Z AA dx (cid:20) K − β γ − β γ β K LO (cid:21) ( x ; A ) (cid:21) = K LO ( A ; A ) + Z AA dx K LO ( A ; A ) K − ( x ; A ) β γ − β γ β K LO ( x ; A ) + · · · . (G8)Now we reveal the behavior of the second term, which can be regarded as the NLO cor-rection term. For simplicity, we consider two-dimensional operator space. We denote the37igenvalues of γ / (2 β ) by λ and λ ( λ < λ ). Then K LO ( A ; A ) is given by a linearcombination of { ( A/A ) λ , ( A/A ) λ } . In Eq. (G8), we note that K LO ( A ; A ) K − ( x ; A ) = K LO ( A ; A ) K LO ( A ; x ) = K LO ( A ; x ) , (G9)because K LO is an evolution matrix. Then the integrand is given by K LO ( A ; A ) K − ( x ; A ) β γ − β γ β K LO ( x ; A )= K LO ( A ; x ) β γ − β γ β K LO ( x ; A )= (linear comb. of { ( A/x ) λ , ( A/x ) λ } ) × (linear comb. of { ( x/A ) λ , ( x/A ) λ } )= (linear comb. of { ( A/A ) λ , ( A/A ) λ , ( A λ /A λ ) x λ − λ , ( A λ /A λ ) x λ − λ } ) . (G10)After the integration, we have Z AA dx K LO ( A ; A ) K − ( x ; A ) β γ − β γ β K LO ( x ; A )= (linear comb. of { ( A/A ) λ A , ( A/A ) λ A , ( A/A ) λ A, ( A/A ) λ A } ) . (G11)Therefore either order in A or A is raised by one in the NLO correction term comparedto K LO .In the following we concretely study the NLO eﬀects in full QCD. To simplify calculations,we decompose the four operators ˆ O , . . . , ˆ O into scalar parts and traceless tensor parts.Then, anomalous dimension matrices are two-by-two matrices.Before this, for convenience we introduce ζ ′ R : ζ ′ R = /g ( µ ( t )) /g ( µ ( t )) ζ R . (G12)This matrix describes the time evolution of { /g ( µ ( t )) ˜ O , /g ( µ ( t )) ˜ O , ˜ O , ˜ O } as seenfrom the deﬁnition [cf. Eq. (1.47)]. The advantage is that the perturbation order of thismatrix is organized by the power of g .Let us begin with scalar operators. Scalar operators are given by linear combinations ofˆ O R, S2 and ˆ O R, S4 . Then, we have RG equation of [cf. Eq. (3.14)] µ ddµ ˆ O R, S2 ˆ O R, S4 ! = − γ S ˆ O R, S2 ˆ O R, S4 ! . (G13)The matrix γ S is readily read oﬀ from γ and γ given in App. C. The evolution matrix K S at LO, satisfying dK SLO dA = γ S0 β A K SLO , is obtained as K SLO = A A C F β (cid:0) A A − (cid:1) ! , (G14)where the eigenvalues of γ S0 / (2 β ) are λ = − λ = 0. From the general argument above[in particular from Eq. (G11)], we see that the NLO K matrix is given in the form K SNLO = K SLO + (linear comb. of { A /A, A , A } ) . (G15)The knowledge on the order of the NLO correction enables us to accurately estimate theparametrical error of the LO calculation. In the LO calculation, we use ζ ′ R LO = { } + O ( A )38nd K SLO = { A /A, } + O ( A /A, A , A ). Here we imply that a left-hand side is given by alinear combination of the functions inside {} and possesses the errors shown by O ( ... ). Weuse this notation hereafter. Then we have ζ ′ R, SLO K SLO = { A /A, } + O ( A /A, A , . . . ) , (G16)where ζ ′ R, S is the two-by-two matrix which relates the ﬂowed operators (1 /g ( µ ( t )) ) ˜ O S2 ( t, x )and ˜ O S4 ( t, x ) to the unﬂowed operators ˆ O R, S2 ( x ; µ ( t )) and ˜ O R, S4 ( x ; µ ( t )). Then, for(1 /g ( µ ( t )) ) ˜ O S2 ( t, x ) | LO = ( ζ ′ R, SLO K S LO ) ˆ O R, S2 ( x ; µ ) + ( ζ ′ R, SLO K S LO ) ˆ O R, S4 ( x ; µ ), from Eq. (G16)we obtain1 g ( µ ( t )) ˜ O S2 ( t, x ) = 1 g ( µ ( t )) (cid:20) g ( µ ) (cid:18) ˆ O R, S2 ( x ; µ ) + 6 C F β ˆ O R, S4 ( x ; µ ) (cid:19) + O ( g ( µ ) ) (cid:21) + O ( g ( µ ( t )) ) , (G17)where the error and neglected higher-order terms are shown by O ( ... ). (Here we set µ = µ ( t ).)This result corresponds to Eq. (3.39). For ˜ O S4 ( t, x ), we have˜ O S4 ( t, x ) = [1 + O ( A ) + O ( A ) + O ( A, AA , A )] ˆ O R, S4 ( x ; µ )+ [ O ( A ) + O ( A ) + O ( A, AA , A )] ˆ O R, S2 ( x ; µ ) , (G18)which corresponds to Eq. (3.40). In obtaining this result, we have noted the following facts.First, ˜ O S4 ( t, x ) is written in terms of ˆ O R, S2 ( x ; µ ( t )) and ˆ O R, S4 ( x ; µ ( t )) with ζ ′ R as˜ O S ( t, x ) = (1 + O ( A )) ˆ O R, S4 ( x ; µ ( t )) + O ( A ) ˆ O R, S2 ( x ; µ ( t )) . (G19)Here, the key is that ˆ O R, S2 ( x ; µ ( t )) is multiplied by a factor of O ( A ). Then in rewrit-ing ˆ O R, S2 ( x ; µ ( t )) in terms of ˆ O R, S2 ( x ; µ ) and ˆ O R, S4 ( x ; µ ) with K SLO = { A /A, } + O ( A /A, A , A ) [cf. Eqs. (G14) and (G15)], we do not have O (1 /A ) contributions from thispart. Also one can see that O ( A ) ˆ O R, S2 ( x ; µ ( t )) gives the O ( A ), O ( A ), and O ( A, AA , A )terms in the coeﬃcient of ˆ O R, S2 ( x ; µ ) in Eq. (G18) and gives O ( A ), O ( A ), and O ( A, AA , A ) contributions to the coeﬃcient of ˆ O R, S4 ( x ; µ ) in Eq. (G18). Secondly, inrewriting ˆ O R, S4 ( x ; µ ( t )) in terms of ˆ O R, S2 ( x ; µ ) and ˆ O R, S4 ( x ; µ ) with K S , the relevant compo-nents K S21 and K S22 are exactly given by ( K S ) = 0 and ( K S ) = 1, and thus ˆ O R, S4 ( x ; µ ( t )) =ˆ O R, S4 ( x ; µ ). This is because ˆ O R, S4 ( x ; µ ) is proportional to ˆ O R, S5 ( x ; µ )(:= O R, S5 ( x ; µ )) due toEOM, and ˆ O S ( x ; µ ) is a ﬁnite operator, which does not need renormalization. Hence, we donot have O (1 /A ) contributions in Eq. (G18). The O ( A ) and O ( A ) terms in Eq. (G18),respectively, can be explicitly obtained with the combination of the NLO ζ ′ R and the LO K , and that of the NLO ζ ′ R and the NLO K .Although we have obtained the main results for the scalar parts in Eqs. (G17) and (G18),we are also able to explicitly obtain the NLO K matrix with Eq. (G8). It would beinteresting to check the validity of such an explicit NLO K matrix. As a possible check,we here consider the small ﬂow time limit t →

0, where g ( µ ( t )) →

0. Then we focus on Explicitly, it is given by ζ ′ R, S = (cid:18) ζ ′ R, S22 ζ ′ R, S24 ζ ′ R, S42 ζ ′ R, S44 (cid:19) . At LO, this is the unit matrix. A rather than A . As mentioned above, we can give accurate coeﬃcientsup to A = [ g ( µ ) / (4 π ) ] in rewriting ˜ O S2 ( t, x ) and ˜ O S4 ( t, x ) in terms of ˆ O R, S2 ( x ; µ ) andˆ O R, S4 ( x ; µ ) by using the NLO K matrix and the NLO ζ ′ R . Explicitly, we obtain˜ O S2 ( t, x ) = g ( µ ) (cid:18) ˆ O R, S2 ( x ; µ ) + 6 C F β ˆ O R, S4 ( x ; µ ) (cid:19) + g ( µ ) (4 π ) β (cid:20) (cid:18) C A − C A T F − C F T F (cid:19) ˆ O R, S2 ( x ; µ )+ (cid:18) C A C F + 3 C F − C F T F (cid:19) ˆ O R, S4 ( x ; µ ) (cid:21) , (G20)˜ O S4 ( t, x ) = ˆ O R, S4 ( x ; µ ) + g ( µ ) (4 π ) T F (cid:18) ˆ O R, S2 ( x ; µ ) + 6 C F β ˆ O R, S4 ( x ; µ ) (cid:19) + (cid:20) g ( µ ) (4 π ) (cid:21) T F β (cid:20) (cid:18) C A − C A T F − C F T F (cid:19) ˆ O R, S2 ( x ; µ )+ (cid:18) C A C F + 3 C F − C F T F (cid:19) ˆ O R, S4 ( x ; µ ) (cid:21) . (G21)Here we just showed the leading order contribution for small g ( µ ( t )). We can conﬁrm thevalidity of these results as follows. In the expression of T S ( x ) in terms of the ﬂowed operators, T S ( x ) = c ( t ) ˜ O S ( t, x ) + c ′ ( t ) ˜ O S ( t, x ), when one uses the above results to rewrite the ﬂowedoperators in terms of the unﬂowed operators at µ = µ , the two-loop order expression of T S ( x ) (B24) with the renormalization scale µ = µ can be correctly reproduced with the LO c , c ′ . (Again, since we consider the small g ( µ ( t )) limit, it is suﬃcient to use the LO c , c ′ ,whose higher-order result just aﬀects higher power in g ( µ ( t )).)Now let us move on to traceless parts and do a parallel analysis. In studying tracelessparts, we have a complication that traceless operators to be considered should be changeddepending on perturbation order. At LO, the traceless operators are given byˆ O R, TL1 ,µν = ˆ O R ,µν −

14 ˆ O R ,µν , (G22)ˆ O R, TL3 ,µν = ˆ O R ,µν −

12 ˆ O R ,µν , (G23)while at NLO they are given by [cf. (F3) and (F4)]ˆ O R, TL1 ,µν = ˆ O R ,µν −

14 ˆ O R ,µν − A (cid:18) C A ˆ O R ,µν + 712 C F ˆ O R ,µν (cid:19) , (G24)ˆ O R, TL3 ,µν = ˆ O R ,µν −

12 ˆ O R ,µν − A (cid:18) − T F ˆ O R ,µν + 23 C F ˆ O R ,µν (cid:19) . (G25)When we consider traceless operators at N k LO ( k -loop), we need to know the anomalousdimension matrix for the four operators ˆ O R ,µν , . . . , ˆ O R ,µν at N k LO ( k + 1-loop). For instanceat NLO, when we consider µ ddµ ˆ O R, TL1 ,µν = µ ddµ (cid:20) ˆ O R ,µν −

14 ˆ O R ,µν − A (cid:18) C A ˆ O R ,µν + 712 C F ˆ O R ,µν (cid:19)(cid:21) , (G26)the O ( A ) term inside the square brackets gives O ( A ) term after µd/ ( dµ ) is operated. Thenfor consistency we need to know the NLO anomalous dimension matrix.40e deﬁne the anomalous dimension matrix for traceless operators as µ ddµ ˆ O R, TL1 ,µν ˆ O R, TL3 ,µν ! = − γ TL ˆ O R, TL1 ,µν ˆ O R, TL3 ,µν ! . (G27)Writing γ TL = γ TL0 A + γ TL1 A , we have γ TL0 = T F − C F − T F C F ! , (G28) γ TL1 = (35 C A T F + 74 C F T F ) − (47 C A C F − C F − C F T F ) − (35 C A T F + 74 C F T F ) (47 C A C F − C F − C F T F ) ! . (G29)The eigenvalues of γ / β are 0 and λ ≡ (2 C F + T F )2 β = C F + T F )11 C A − T F . The LO K matrix is givenby K TLLO = 12 C F + T F C F + T F (cid:16) AA (cid:17) λ C F (cid:18) − (cid:16) AA (cid:17) λ (cid:19) T F (cid:18) − (cid:16) AA (cid:17) λ (cid:19) T F + 2 C F (cid:16) AA (cid:17) λ . (G30)From Eq. (G11), the NLO K TL takes a form K TLNLO = K TLLO + (linear comb. of { A , ( A/A ) λ A , A, ( A/A ) λ A } ) . (G31)The knowledge on the order of the NLO correction enables us to accurately estimate theparametrical error of the LO calculation. We have ζ ′ R, TLLO K TLLO = ( { } + O ( A ))( { , ( A/A ) λ } + O ( A , ( A/A ) λ A , A, ( A/A ) λ A ))= { , ( A/A ) λ } + O ( A , ( A/A ) λ A , . . . ) . (G32)Then we can show the LO results of the traceless ﬂowed operators with explicit parametricerrors (and neglected higher-order terms):1 g ( µ ) ˜ O TL1 ,µν ( t, x ) = 2 C F C F + T F (cid:18) ˆ O R, TL1 ,µν ( µ ) + 14 ˆ O R, TL3 ,µν ( µ ) + O ( A ) (cid:19) + O (( A/A ) λ ) , (G33)˜ O TL3 ,µν ( t, x ) = 4 T F C F + T F (cid:18) ˆ O R, TL1 ,µν ( µ ) + 14 ˆ O R, TL3 ,µν ( µ ) + O ( A ) (cid:19) + O (( A/A ) λ ) . (G34)These correspond to Eqs. (3.37) and (3.38). Here we understand ˆ O R, TL1 , µν ( µ ) as the LO onesgiven by Eqs. (G22) and (G23). They have O ( A ) errors and this is consistent with theabove errors.Calculating the NLO K matrix, we can explicitly show higher-order results. Again we areinterested in higher order in g ( µ ) and just show the leading contribution for small g ( µ ( t )).With the NLO K matrix (and LO ζ ′ R matrix, which is suﬃcient for the present purpose),we have 1 g ( µ ( t )) ˜ O TL1 ,µν ( t, x ) = 2 C F C F + T F (cid:18) ˆ O R, TL1 ,µν ( x ; µ ) + 14 ˆ O R, TL3 ,µν ( x ; µ ) (cid:19) , (G35)˜ O TL3 ,µν ( t, x ) = 4 T F C F + T F (cid:18) ˆ O R, TL1 ,µν ( x ; µ ) + 14 ˆ O R, TL3 ,µν ( x ; µ ) (cid:19) . (G36)Now ˆ O R, TL1 , µν ( µ ) are the ones at NLO [Eqs. (G24) and (G25)]. Although there is no apparentdiﬀerence from the LO calculations (G33) and (G34), the O ( A ) terms are now ﬁxed and41urned out to be zero. (The expected error is now O ( A )). One can conﬁrm that the one-loopexpression of T TL (B26) with the renormalization scale µ = µ can be correctly reproducedfrom T TL µν ( x ) = c ( t ) ˜ O TL1 ,µν ( t, x ) + c ( t ) ˜ O TL3 ,µν ( t, x ) with the LO c , c and the above results. H. Results of thermodynamics quantities with other t → extrapolationfunctions In this appendix, for reference, we present the results obtained with other t → t for the entropy density and trace anomaly, anduse a linear function with the anomalous dimension, g ( µ ( t )) λ t , given in Eq. (4.20), for theentropy density.The results are summarized in Tables H1 and H2. One can see that the perturbativeextrapolation functions, Eqs. (2.3) and (2.4), gives smaller diﬀerence in the NLO and N LOresults than the other linear-type functions. We can also see from Table H1 that the numer-ical impact of the inclusion of the anomalous dimension of dimension-six operators is notsigniﬁcant. This can be seen also in Fig. H1, where we show the t → s/T (NLO) T /T c Perturbative function Linear Linear with anomalous dim.0.93 0 . . . . . . . . . . . . . . . . . . . . . . . . s/T (N LO)

T /T c Perturbative function Linear Linear with anomalous dim.0.93 0 . . . . . . . . . . . . . . . . . . . . . . . . Table H1

Comparison of the results for the entropy density obtained with diﬀerent t → k = 1 or k = 2,“Linear” means a linear function in t , and “Linear with anomalous dim.” means Eq. (4.20).We only show statistical errors. 42 ig. H1 NLO analysis of the entropy density with diﬀerent t → t → /T (NLO) T /T c Perturbative function Linear0.93 0 . . . . . . . . . . /T (N LO)

T /T c Perturbative function Linear0.93 0 . . . . . . . . . . Table H2

Comparison of the results for the trace anomaly obtained with diﬀerent t → k = 1 or k = 2, and“Linear” means a linear function in t . We only show statistical errors.44 eferences [1] H. Suzuki, “Energy–momentum tensor on the lattice: recent developments,” PoS

LATTICE2016 (2017) 002, arXiv:1612.00210 [hep-lat] .[2] H. Suzuki, “Energy–momentum tensor from the Yang–Mills gradient ﬂow,”

PTEP (2013) 083B03, arXiv:1304.0533 [hep-lat] . [Erratum: PTEP2015,079201(2015)].[3] H. Makino and H. Suzuki, “Lattice energy–momentum tensor from the Yang–Mills gradient ﬂow—inclusion of fermion ﬁelds,”

PTEP (2014) 063B02, arXiv:1403.4772 [hep-lat] . [Erratum:PTEP2015,079202(2015)].[4] R. Narayanan and H. Neuberger, “Inﬁnite N phase transitions in continuum Wilson loop operators,”

JHEP (2006) 064, arXiv:hep-th/0601210 .[5] M. L¨uscher, “Trivializing maps, the Wilson ﬂow and the HMC algorithm,” Commun. Math. Phys. (2010) 899–919, arXiv:0907.5491 [hep-lat] .[6] M. L¨uscher, “Properties and uses of the Wilson ﬂow in lattice QCD,”

JHEP (2010) 071, arXiv:1006.4518 [hep-lat] . [Erratum: JHEP03,092(2014)].[7] M. L¨uscher and P. Weisz, “Perturbative analysis of the gradient ﬂow in non-abelian gauge theories,” JHEP (2011) 051, arXiv:1101.0963 [hep-th] .[8] M. L¨uscher, “Chiral symmetry and the Yang–Mills gradient ﬂow,” JHEP (2013) 123, arXiv:1302.5246 [hep-lat] .[9] R. V. Harlander, Y. Kluth, and F. Lange, “The two-loop energy–momentum tensor within the gradient-ﬂow formalism,” Eur. Phys. J.

C78 no. 11, (2018) 944, arXiv:1808.09837 [hep-lat] .[10] J. Artz, R. V. Harlander, F. Lange, T. Neumann, and M. Prausa, “Results and tech-niques for higher order calculations within the gradient-ﬂow formalism,”

JHEP (2019) 121, arXiv:1905.00882 [hep-lat] . [Erratum: JHEP 10, 032 (2019)].[11] FlowQCD , M. Asakawa, T. Hatsuda, E. Itou, M. Kitazawa, and H. Suzuki, “Thermodynamicsof SU(3) gauge theory from gradient ﬂow on the lattice,”

Phys. Rev. D no. 1, (2014) 011501, arXiv:1312.7492 [hep-lat] . [Erratum: Phys.Rev.D 92, 059902 (2015)].[12] Y. Taniguchi, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, T. Umeda, andN. Wakabayashi, “Exploring N f = 2+1 QCD thermodynamics from the gradient ﬂow,” Phys. Rev. D no. 1, (2017) 014509, arXiv:1609.01417 [hep-lat] . [Erratum: Phys.Rev.D 99,059904 (2019)].[13] M. Kitazawa, T. Iritani, M. Asakawa, T. Hatsuda, and H. Suzuki, “Equation of Statefor SU(3) Gauge Theory via the Energy–Momentum Tensor under Gradient Flow,” Phys. Rev. D no. 11, (2016) 114512, arXiv:1610.07810 [hep-lat] .[14] S. Ejiri, R. Iwami, M. Shirogane, N. Wakabayashi, K. Kanaya, M. Kitazawa, H. Suzuki, Y. Taniguchi,and T. Umeda, “Determination of latent heat at the ﬁnite temperature phase transition of SU(3) gaugetheory,” PoS

LATTICE2016 (2017) 058, arXiv:1701.08570 [hep-lat] .[15] M. Kitazawa, T. Iritani, M. Asakawa, and T. Hatsuda, “Correlations of the energy–momentum tensor via gradient ﬂow in SU(3) Yang-Mills theory at ﬁnite temperature,”

Phys. Rev. D no. 11, (2017) 111502, arXiv:1708.01415 [hep-lat] .[16] WHOT-QCD , K. Kanaya, S. Ejiri, R. Iwami, M. Kitazawa, H. Suzuki, Y. Taniguchi, and T. Umeda,“Equation of state in (2+1)-ﬂavor QCD at physical point with improved Wilson fermion action usinggradient ﬂow,”

EPJ Web Conf. (2018) 07023, arXiv:1710.10015 [hep-lat] .[17]

WHOT-QCD , Y. Taniguchi, S. Ejiri, K. Kanaya, M. Kitazawa, A. Suzuki, H. Suzuki, and T. Umeda,“Energy–momentum tensor correlation function in N f = 2 + 1 full QCD at ﬁnite temperature,” EPJ Web Conf. (2018) 07013, arXiv:1711.02262 [hep-lat] .[18] R. Yanagihara, T. Iritani, M. Kitazawa, M. Asakawa, and T. Hatsuda, “Distribution of Stress Tensoraround Static Quark–Anti-Quark from Yang-Mills Gradient Flow,”

Phys. Lett. B (2019) 210–214, arXiv:1803.05656 [hep-lat] .[19] T. Hirakida, E. Itou, and H. Kouno, “Thermodynamics for pure SU(2) gauge theory using gradientﬂow,”

PTEP no. 3, (2019) 033B01, arXiv:1805.07106 [hep-lat] .[20] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, Y. Taniguchi, and T. Umeda,“Equation of state near the ﬁrst order phase transition point of SU(3) gauge theory using gradientﬂow,”

PoS

LATTICE2018 (2018) 164, arXiv:1811.04220 [hep-lat] .[21] T. Iritani, M. Kitazawa, H. Suzuki, and H. Takaura, “Thermodynamics in quenched QCD:energy–momentum tensor with two-loop order coeﬃcients in the gradient-ﬂow formalism,”

PTEP no. 2, (2019) 023B02, arXiv:1812.06444 [hep-lat] .[22] Y. Taniguchi, A. Baba, A. Suzuki, S. Ejiri, K. Kanaya, M. Kitazawa, T. Shimojo, H. Suzuki, andT. Umeda, “Study of energy–momentum tensor correlation function in N f = 2 + 1 full QCD for QGPviscosities,” PoS

LATTICE2018 (2019) 166, arXiv:1901.01666 [hep-lat] .

23] M. Kitazawa, S. Mogliacci, I. Kolb´e, and W. A. Horowitz, “Anisotropic pressure inducedby ﬁnite-size eﬀects in SU(3) Yang-Mills theory,”

Phys. Rev. D no. 9, (2019) 094507, arXiv:1904.00241 [hep-lat] .[24] K. Kanaya, A. Baba, A. Suzuki, S. Ejiri, M. Kitazawa, H. Suzuki, Y. Taniguchi, and T. Umeda, “Studyof 2+1 ﬂavor ﬁnite-temperature QCD using improved Wilson quarks at the physical point with thegradient ﬂow,” PoS

LATTICE2019 (2019) 088, arXiv:1910.13036 [hep-lat] .[25]

WHOT-QCD , Y. Taniguchi, S. Ejiri, K. Kanaya, M. Kitazawa, H. Suzuki, and T. Umeda,“ N f = 2+1 QCD thermodynamics with gradient ﬂow using two-loop matching coeﬃcients,” Phys. Rev. D no. 1, (2020) 014510, arXiv:2005.00251 [hep-lat] . [Erratum: Phys.Rev.D 102,059903 (2020)].[26] R. Yanagihara, M. Kitazawa, M. Asakawa, and T. Hatsuda, “Distribution of Energy–MomentumTensor around a Static Quark in the Deconﬁned Phase of SU(3) Yang-Mills Theory,”

Phys. Rev. D no. 11, (2020) 114522, arXiv:2010.13465 [hep-lat] .[27] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, Y. Taniguchi andT. Umeda, “Latent heat and pressure gap at the ﬁrst-order deconﬁning phase transition ofSU(3) Yang-Mills theory using the small ﬂow-time expansion method,”

PTEP (2021) 013B08, arXiv:2011.10292 [hep-lat] .[28] H. Kim and S. H. Lee, “Renormalization of dimension 6 gluon operators,”

Phys. Lett. B (2015) 352–355, arXiv:1503.02280 [hep-ph] .[29] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier, and B. Petersson, “Thermo-dynamics of SU(3) lattice gauge theory,”

Nucl. Phys. B (1996) 419–444, arXiv:hep-lat/9602007 .[30]

CP-PACS , M. Okamoto et al. , “Equation of state for pure SU(3) gauge theory with renormalizationgroup improved action,”

Phys. Rev. D (1999) 094510, arXiv:hep-lat/9905005 .[31] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, “Precision SU(3) lattice thermody-namics for a large temperature range,” JHEP (2012) 056, arXiv:1204.6184 [hep-lat] .[32] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, “Full result for the QCDequation of state with 2+1 ﬂavors,” Phys. Lett. B (2014) 99–104, arXiv:1309.5258 [hep-lat] .[33]

HotQCD , A. Bazavov et al. , “Equation of state in ( 2+1 )-ﬂavor QCD,”

Phys. Rev. D (2014) 094503, arXiv:1407.6387 [hep-lat] .[34] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, and M. Kitazawa, “Latent heat at theﬁrst order phase transition point of SU(3) gauge theory,” Phys. Rev. D no. 1, (2016) 014506, arXiv:1605.02997 [hep-lat] .[35] L. Giusti and M. Pepe, “Equation of state of a relativistic theory from a moving frame,” Phys. Rev. Lett. (2014) 031601, arXiv:1403.0360 [hep-lat] .[36] L. Giusti and M. Pepe, “Equation of state of the SU(3) Yang–Mills theory: A precise determinationfrom a moving frame,”

Phys. Lett. B (2017) 385–390, arXiv:1612.00265 [hep-lat] .[37] M. Caselle, A. Nada, and M. Panero, “QCD thermodynamics from lattice calculations withnonequilibrium methods: The SU(3) equation of state,”

Phys. Rev. D no. 5, (2018) 054513, arXiv:1801.03110 [hep-lat] .[38] Flavour Lattice Averaging Group , S. Aoki et al. , “FLAG Review 2019: Flavour Lattice AveragingGroup (FLAG),”

Eur. Phys. J. C no. 2, (2020) 113, arXiv:1902.08191 [hep-lat] ..