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 flow. We arguethat in the t → g ( µ ( t )), the flow time dependent running coupling, where the power is determined bythe perturbation order we consider. From actual lattice data, we confirm 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 fieldtheory, 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 flow 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 flowed operators,defined from the gradient flow [4–8]. Due to UV finiteness of flowed operators, the currentrepresented with flowed operators satisfies 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 flow time, whose mass dimension is −
2; ˜ O i,µν ( t, x )’s are (dimension-four)flowed operators (whose explicit definitions are given below). The coefficients ˜ c i ( t )’s can beperturbatively calculated via the small flow time expansion [7] of the flowed 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 flowed operators ˜ O i,µν nonperturbatively and then multiply the perturbative coefficients ˜ 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 confirmed.In lattice simulation using the SF t X method, one needs to take the small flow time limit( t → This is because the expression for the EMT in the SF t X method becomes exactin the t → g ( µ ( t )), the flow timedependent running coupling, in the perturbative coefficients ˜ c i ( t ), and O ( t ) contributionsin Eq. (1.1). However one cannot directly obtain lattice data at t = 0 because lattice datasuffer from serious discretization effects when the flow 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 fixed order perturba-tive results for ˜ c i ( t ), a higher-order effect in g ( µ ( t )) should also exist. Parametrically, such aneffect is dominant in small t region. This is because g ( µ ( t )) n ≫ t Λ ∼ e − π / ( β g ( µ ( t ))) ( β is the first coefficient 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 effect in g ( µ ( t ))in order to accurately performe t → t dependence of the fixed order formulafor the EMT in the SF t X method (which means that we use fixed 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 briefly 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 defined 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 flavor f . These are bare composite operators and notfinite in general.In the SF t X method, one rewrites the EMT in terms of flowed operators. The Yang–Millsgradient flow [4–8] is defined by the following differential equations with respect to the flowtime t . For the gauge field, it is defined 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 field is defined as D µ = ∂ µ + [ B µ , · ] . (1.9) G µν is the field strength of the flowed gauge field B µ . For the fermion field, 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 define analogous operators at positive flow 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 flowed operators are finite operators [7]. We have accomplished the renormal-ization of flowed fermion fields 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 sufficient 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 flowed operators, one needs to know the relationbetween the flowed operators and unflowed operators. It can be studied through the smallflow 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 flowed operators and unflowed operators is known, one canexpress the EMT, given by the unflowed operators, in terms of the flowed 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 (finite) 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 flowed 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 coefficient k ( n ) i is a polynomial of L ( µ, t ) ≡ log (2 µ e γ E t ). Since the EMT and flowed operators are renormalization group (RG)invariant, the coefficients 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 coefficients k (2) i were obtained in Ref. [9]. For quenched QCD, thethree-loop coefficient 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 first ( 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 defined 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 unflowed 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 finite, one can determine Z from the requirement that ζ ( t ) H − Z is finite. In this way, Z has been calculated up to two-loop [9]. Now, we can express the flowed operators with therenormalized (finite) quantities,˜ O µν ( t, x ) = ζ R ( t ; g ( µ ) , µ ) ˆ O Rµν ( x ; µ ) + O ( t ) . (1.47)Here we define the renormalized matrix ζ R as ζ R ( t ; g ( µ ) , µ ) = Rζ ( t ) H − Z ( g ( µ )) , (1.48)which is finite.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 first few coefficients: β = 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 first 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 unflowed 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 confirm validityof the t → L ( µ, t ) dependence of the perturbative series forthe coefficients 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 fixed order formula for the EMT. (The NLO ζ R isused in App. G.) In App. G, we give an argument to estimate higher-order effects 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 first 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 fitted 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 fit 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 fitted 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 fit parameters. T TL µν ( x ) and T S ( x ) correspond to finalresults in which one is interested.In quenched QCD, the fit parameters k ( k +1)1 and k ( k +2)2 correspond to higher-order coeffi-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 finite 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 coefficients are known; see Sect. 5 and App. D. Ofcourse, these properties are not exactly hold due to systematic errors in fits. 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 fit parameters k ( k +1)1 and k ( k +2)2 common to all the simulated temperatures, taking the first property into account. We checkthe validity of the use of the above extrapolation functions by examining the second prop-erty (ii), behavior of the fit 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 fit 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 ) effect coming from dimension-six operators, which is a subleading effect forsufficiently 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 coefficients 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 satsfies 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 clarified 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 flowedoperators in full QCD. Hence, we consider another derivation which is applicable to fullQCD. This is Derivation II. 10e consider the difference 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 difference between the coefficients 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 flowed 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 flowed 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 flowed 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 definition 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 flowed 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 effects in g ( µ ) appear only in the coefficient ofˆ O R ,µν and the coefficient 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 specifically, 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 differ 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 flowed 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 defined 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 confirm that this result satisfies required relations for K , given in App. E.From these matrices, t -dependence of the flowed 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 effects, 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 effects 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 flowed operators ˜ O S2 and ˜ O S4 are proportional to T S . This expectation can be deniedmore explicitly by considering the difference 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 flow 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 sufficiently small t region.Nevertheless, it might be possible that the O ( t ) effect dominates over the O ( g ( µ ( t )) n ) effectin 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 specific example.We first give a general argument how to investigate the detailed t dependence comingfrom dimension-six operators. We go back to the small flow time expansion of the flowedoperators ˜ 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 coefficient 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 flow representation of theEMT [see Eq. (1.1)], T GF µν ( x ; t ) = X i ˜ c i ( t ; g ( µ )) ˜ O i,µν ( t, x ) , (4.3)differs 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 difference 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 first equality, we have used the fact that the flowed operator is finite and thus can berewritten by renormalized (finite) 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 flow 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 sufficient 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 flowed 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 flow 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 flow 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 flowed operators. See this reference forthe details of the lattice setup. We study the finite temperature effect, i.e., the difference 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 discretizationeffects under good control we need 2 e γ E t/a ≫
1. (Here we assume that the typical scaleof the flowed 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 coeffi-cients, namely we use the NLO or NNLO formula of the EMT. Finally, we extrapolate thedata to zero flow time limit using Eqs. (2.3) and (2.4) with k = 1 or 2. Then we obtain thefinal results.In the second and third steps, we rely on perturbation theory. This is valid when the flowtime satisfies 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 fit 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 fit. It indicates validity of the extrapolationfunction.We can also confirm the validity of the extrapolation function by comparing the result ofthe fit parameter k (2)1 , which corresponds to the NNLO coefficient for c ( t ), with its alreadyknown result. From the fit, 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 coefficient 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 coefficient is totally determined by the perturbative coefficients up to N k LOand the beta function; see App. D. Then we can compare the L ( µ, t ) dependence of the fitparameter 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 fit 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 different 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) , fit1 ( L = 0), which is the obtained fit 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 fit range used in the t → ig. 3 The fit parameter k (2)1 or k (3)1 as a function of log s . The data points with error barsshow the results of the fit parameter k (2)1 or k (3)1 obtained in the t → s . The blue lines show log s dependence of the perturbative coefficients k (2)1 and k (3)1 dictated from the RG equation; see Eq. (D2). We assume different log s independentconstants for three blue lines; see the main text.with k (2) , fit1 ( L = log (2 )). Here k (2) , fit1 ( L = log (2 )) is the obtained fit parameter when weset s = 2. The difference 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 fit 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) difference is given as a systematic error. Our final 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 asflat within the statistical errors. Then final 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 fitparameter 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 difference 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 fit 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 fit 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 fit range used in the t → Fig. 7
The fit 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 coefficients 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 coefficient is known to be k (3)2 ( L = 0) = − . ... [21], and wefind 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 final 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” coefficient c ( t ) meant the result up to two-loop order, weregard c ( t ) up to two-loop order as the NLO coefficient 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 coefficient in c ( t ), we notice that it is natural to call the one-loop order coefficientLO 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 fit with the range 0 . ≤ tT ≤ . t region. Note that the way to estimate this systematic error is different. 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 coefficients k (3)1 and k (4)2 in quenched QCD.Our lattice analysis is carried out at sufficiently small t region, where the coupling constantis given by ( α s =) g / (4 π ) . .
25. In full QCD, it might be difficult 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 sufficiently 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 difficult to carry out t → t dependence (e.g. higher-order effect in g ( µ ( t )) orlinear function in t ) in a fitting function, a fit would be destabilized due to the difficulty indistinguishing different functions. One conservative attitude would be giving the systematicerror associated with extrapolation functions by trying some functions. However, such adifficulty may be systematically overcome by going to finer 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 ScientificResearch 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 defined 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 define T F ≡ n f T R = n f /
2, where n f is the number of quark flavors.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 field strength is defined 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 finite. 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 finite 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 finite. Since β ( g ) / (2 Dg ) is finite as ǫ → Z (cid:18) gZ g ∂Z g ∂g (cid:19) (B9)should be finite. 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 firstequation, 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 defining a renormalization matrix in a parallel manner, we deduce that Z i − Z i + 14 Z i (B17)is finite. 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 finite operator. Each coefficient should be finitein Eq. (B21). Then, for the coefficient 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 coefficient of the gluonic operator can be obtained in this way also infull QCD. For ˆ O R, S4 , we deduce that 1 + ǫZ should be finite. (Thus, Z has only simplepoles in ǫ to all orders.) However, it seems difficult to fix its finite 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 definition, 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 definition of the renormalization factor Z is the inverse of thatof Ref. [9]. 33 . Coefficients c ( ′ ) i ( t ) The coefficients 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) satisfies all the relations (at appropriate order). F. ζ R at NLO We show the NLO result for ζ R , defined 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 confirm 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 specified 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 effect is more difficult to see than the matrix ζ R (whose higher-order effect is just given by higher power in g ( µ ( t )) ).We first 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 define γ , γ , . . . 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 difficult 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 satisfies 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 effects 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 definition [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 off 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 flowed operators (1 /g ( µ ( t )) ) ˜ O S2 ( t, x )and ˜ O S4 ( t, x ) to the unflowed 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 coefficient of ˆ O R, S2 ( x ; µ ) in Eq. (G18) and gives O ( A ), O ( A ), and O ( A, AA , A ) contributions to the coefficient 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 finite 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 flow 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 coefficientsup 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 confirm thevalidity of these results as follows. In the expression of T S ( x ) in terms of the flowed 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 flowedoperators in terms of the unflowed 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 sufficient to use the LO c , c ′ ,whose higher-order result just affects 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 define 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 flowed 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 sufficient 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 apparentdifference from the LO calculations (G33) and (G34), the O ( A ) terms are now fixed and41urned out to be zero. (The expected error is now O ( A )). One can confirm 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 difference 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 notsignificant. 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 different 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 different 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 different 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 flow,”
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 flow—inclusion of fermion fields,”
PTEP (2014) 063B02, arXiv:1403.4772 [hep-lat] . [Erratum:PTEP2015,079202(2015)].[4] R. Narayanan and H. Neuberger, “Infinite N phase transitions in continuum Wilson loop operators,”
JHEP (2006) 064, arXiv:hep-th/0601210 .[5] M. L¨uscher, “Trivializing maps, the Wilson flow 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 flow 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 flow in non-abelian gauge theories,” JHEP (2011) 051, arXiv:1101.0963 [hep-th] .[8] M. L¨uscher, “Chiral symmetry and the Yang–Mills gradient flow,” 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-flow 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-flow 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 flow 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 flow,” 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 finite 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 flow in SU(3) Yang-Mills theory at finite 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)-flavor QCD at physical point with improved Wilson fermion action usinggradient flow,”
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 finite 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 gradientflow,”
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 first order phase transition point of SU(3) gauge theory using gradientflow,”
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 coefficients in the gradient-flow 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 finite-size effects 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 flavor finite-temperature QCD using improved Wilson quarks at the physical point with thegradient flow,” 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 flow using two-loop matching coefficients,” 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 Deconfined 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 first-order deconfining phase transition ofSU(3) Yang-Mills theory using the small flow-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 flavors,” Phys. Lett. B (2014) 99–104, arXiv:1309.5258 [hep-lat] .[33]
HotQCD , A. Bazavov et al. , “Equation of state in ( 2+1 )-flavor 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 thefirst 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] ..