Measurement of thermodynamics using gradient flow
Abstract
We analyze bulk thermodynamics and correlation functions of the energymomentum tensor in pure YangMills gauge theory using the energymomentum tensor defined by the gradient flow and small flow time expansion. Our results on thermodynamic observables are consistent with those obtained by the conventional integral method. The analysis of the correlation function of total energy supports the energy conservation. It is also addressed that these analyses with gradient flow require less statistics compared with the previous methods. All these results suggest that the energymomentum tensor can be successfully defined and observed on the lattice with moderate numerical costs with the gradient flow.
Measurement of thermodynamics using gradient flow
Hiroshi Suzuki
Department of Physics, Kyushu University, 6101 Hakozaki, Higashiku, Fukuoka, 8128581, Japan
Email: hsuzuki@phys.kyushuu.ac.jp
\abstract@cs
1 Preliminary
In this manuscript, we consider the energymomentum tensor
As is well known, the energymomentum tensor (EMT) is one of the most fundamental quantities in physics. It is fundamental because this operator is related to the symmetry of spacetime; the EMT can be defined as the Noether current of translational symmetry, and it also constitutes generators of Poincaré transformation. Moreover, elements of the EMT, energy and momentum densities and stress tensor, are basic quantities in physics. The EMT also appears in various fundamental equations such as Einstein equation and hydrodynamic equations.
However, if one wants to define and analyze the EMT in lattice QCD Monte Carlo simulations, there are highly nontrivial problems. First, the definition of the EMT itself is nontrivial on the lattice, because the lattice regularization explicitly breaks translational symmetry. For example, a naïve definition,
(1.0) 
with a choice of discretized definitions for the field strength on the lattice, does not give the correct EMT even for offdiagonal components in the sense that it does not correspond to the EMT in the continuum limit [1]. Although the EMT on the lattice can be constructed by an adequate linear combination of operators, the renormalization factor of each operators has to be determined nonperturbatively, which requires additional numerical simulations for each lattice parameter [1, 2]. Second, even if one employs a definition of the EMT on the lattice, its measurement is noisy. For example, the calculation of the correlation function of the EMT requires extremely large statistics [3, 2]. The analysis becomes more noisy as the lattice spacing becomes finer. Because of these reasons, there are not many studies which have analyzed the EMT directly on the lattice.
On the other hand, if one were able to deal with the EMT on the lattice in a controllable manner, various interesting analyses would become possible, because the EMT is one of the most fundamental quantities in physics. For example, if we have the correctly normalized EMT operator on the lattice one can calculate thermodynamic quantities, energy density and pressure , at nonzero temperature by directly taking the thermal expectation values of the diagonal components of the EMT as
(1.0) 
The measurement of bulk thermodynamics on the lattice itself can be done without the EMT; the method based on the thermodynamic relations called the integral method [4] is well established. The measurement with simple relations Eq. (1.0), however, will provide us with more intuitive way to investigate these quantities. Next, the fluctuations and correlation functions of the EMT are also interesting observables. The fluctuation of the total energy of the system , for example, is related to the specific heat per unit volume as [5]
(1.0) 
This equation suggests a novel way to analyze through the measurement of the energy fluctuation, which becomes possible with the correctly normalized EMT operator. More interesting quantities related to the EMT correlators are transport coefficients. Shear viscosity, for example, is related to the correlation function of via the Kubo formula as [6]
(1.0) 
where is the Kubo’s canonical correlation and . The transport coefficients are important quantities to understand dynamical time evolution of relativistic heavy ion collisions [7]. The measurement of transport coefficients with Kubo formula, however, requires realtime correlation function of the EMT, while the lattice simulations can obtain only imaginary time correlators directly [3, 2]. The analytic continuation from imaginarytime functions to real time is needed, which, however, is a nontrivial procedure requiring an accurate measurement of the EMT correlator.
The application of the EMT operator is not limited to the physics of nonzero temperature. For example, with the correctly normalized EMT operator one would be able to measure the energy and momentum distributions in hadrons and flux tube [8] directly on the lattice. The correlation functions of the EMT in vacuum are also a useful probe to understand vacuum structure of the gauge theory.
Recently, a construction of the EMT in gauge theory is proposed in Ref. [9] on the basis of the YangMills gradient flow [10] and the socalled small flow time expansion [11]. In this method, the constructed EMT is renormalized and does not have divergences, and this property is kept even for the lattice regularization. Moreover, the statistical error is expected to be reduced thanks to the cooling nature of the gradient flow. Therefore, this method can become a solution for the longstanding problem on the measurement of the EMT on the lattice. In this manuscript, we report the first numerical analysis of the EMT in this method in lattice simulations [12].
2 Gradient flow
Because we use the YangMills gradient flow to define the EMT, let us first give a brief review on this concept. The gradient flow for the YangMills gauge field is the continuous transformation of the field defined by the differential equation [10]
(2.0) 
with the YangMills action composed of . Color indices are suppressed for simplicity. The initial condition at is taken for the field in the conventional gauge theory; . The flow time , which has a dimension of inverse mass squared, is a parameter which controls the transformation. The gauge field is transformed along the steepest descent direction as increases. At the tree level, Eq. (2.0) is rewritten as
(2.0) 
Neglecting the gauge dependent terms, Eq. (2.0) is the diffusion equation in fourdimensional space. For positive , therefore, the gradient flow acts as the cooling of the gauge field with the smearing radius . In Ref. [11], it is rigorously proved that all composite operators composed of take finite values for . This property ensures that observables at are automatically renormalized.
Recently, the gradient flow has been applied for various purposes, such as the scale setting, analysis of the running coupling, measurement of topology, as well as the studies on their properties on the lattice; see reviews in Refs. [13, 14] and references [11, 9, 15, 16, 17, 18]. In the present study, we use the gradient flow to define the EMT operator using the small flow time expansion [11, 9, 15].
3 Small flow time expansion and the EMT
In the present study, we consider the EMT defined by the small flow time expansion [9]. In order to illustrate the concept of small flow time expansion, let us consider a composite operator composed of the field at positive flow time. The small flow time expansion asserts that in the small limit this operator can be written by a superposition of operators of the original gauge theory at as
(3.0) 
where on the righthand side represents renormalized operators in some regularization scheme in the original gauge theory at with the subscript denoting different operators, and represents the coordinate in four dimensional spacetime.
The relation Eq. (3.0) is reminiscent of the operator product expansion (OPE). While a product of operators at different spacetime points are expanded by local operators in the OPE, Eq. (3.0) relates an operator at with those of . The validity of this expansion is expected from the cooling nature of the gradient flow; since the operator depends on the fundamental gauge theory only inside the smearing radius around , should be described by the local property of the fundamental theory at in limit. Similarly to the Wilson coefficients in the OPE, the coefficients in Eq. (3.0) can be calculated perturbatively [11, 9].
Using Eq. (3.0), one can define the renormalized EMT [9]. For this purpose, we first consider dimensionfour gaugeinvariant operators for operators of the lefthand side in Eq. (3.0). In pure gauge theory, there are two such operators;
(3.0)  
(3.0) 
Although Eqs. (3.0) and (3.0) are quite similar to the tracelesspart and trace of the EMT, they are not the EMT since is defined at nonzero flow time . Because these operators are gauge invariant, when they are expanded as in Eq. (3.0) only gauge invariant operators can appear in the righthand side. Such operator with the lowest dimension is an identity operator. In the expansion of the traceless operator Eq. (3.0), however, the constant term cannot appear. The next gaugeinvariant operators are the dimensionfour EMTs. Up to this order, therefore, the small flow time expansions of Eqs. (3.0) and (3.0) are given by
(3.0)  
(3.0) 
where is vacuum expectation value and is the correctly normalized conserved EMT with its vacuum expectation value subtracted. Abbreviated are the contributions from the operators of dimension or higher, which are proportional to powers of because of dimensional reasons, and thus suppressed for small .
Combining relations Eqs. (3.0) and (3.0), we have
(3.0) 
The coefficients and are calculated perturbatively in Ref. [9] as
(3.0)  
(3.0) 
Here denotes the running gauge coupling in the scheme with the choice, , and
(3.0)  
(3.0) 
with , , and . Note that a nonperturbative determination of is also proposed recently [15, 19].
There are two important observations: (i) The righthand side of Eq. (3.0) is independent of the regularization because of its UV finiteness, so that one can take, e.g., the lattice regularization scheme; (ii) since flowed fields at depend on the fundamental fields at in the spacetime region of radius , the statistical noise in calculating the righthand side of Eq. (3.0) is suppressed for finite .
4 Numerical procedure on the lattice
The formula Eq. (3.0) indicates that can be obtained by the small limit of the gaugeinvariant local operators defined through the gradient flow. To perform the numerical analysis on the lattice, we take the following steps:

Generate gauge configurations at on a spacetime lattice with a standard algorithm with the lattice spacing and the lattice size .

Solve the gradient flow for each configuration to obtain the flowed link variables for .

Construct and in terms of the flowed link variables at each .

Carry out an extrapolation toward .
Note that this analysis has to be performed in the fiducial window, . Here, is an infrared cutoff scale such as or at which the perturbative analysis used in Eq. (3.0) is violated. The first inequality is necessary to suppress finite corrections, while the second one is required to suppress nonperturbative corrections and finite volume effects. We will later see that the upper limit of the fiducial window is given by by for temperatures analyzed in this study.
5 Numerical setup
To demonstrate that Eq. (3.0) can successfully be used in the measurement of the EMT on the lattice, we consider the gauge theory defined on a fourdimensional Euclidean lattice, whose thermodynamics has been extensively studied by the integral method [4, 20, 21, 22]. We consider the Wilson plaquette gauge action under the periodic boundary condition with several different with being the bare coupling constant. The gradient flow in the direction is numerically solved by the Runge–Kutta method. We have checked that the accumulation errors due to the Runge–Kutta method is more than two order smaller than the statistical errors in all analyses.
6  8  10  

6.20  6.40  6.56  1.65  
6.02  6.20  6.36  1.24  
5.89  6.06  6.20  0.99 
The numerical analyses reported in the following have been performed with the following two settings.

Setting 1: The first analysis is performed on lattices [12]. Gauge configurations are generated by the pseudoheatbath algorithm with the overrelaxation, mixed in the ratio of . We call one pseudoheatbath update sweep plus five overrelaxation sweeps as a “Sweep”. To eliminate the autocorrelation, we take – Sweeps between measurements. The number of gauge configurations for the measurements at finite is . Statistical errors are estimated by the jackknife method. To relate and corresponding for each , we first use the relation between ( is the Sommer scale) and given by the ALPHA Collaboration [23]. The resultant values of are then converted to by using the result at in Ref. [4]. Nine combinations of and corresponding obtained by this procedure are shown in Table 1.

Setting 2: The second analysis is performed on finer lattices with and with several values of . The lattice spacing of each are determined originally [24] by the measurements of [10] and [25] with the gradient flow, and the simulation is performed for several values of . In this manuscript, we present the result only for . The parameters , and for this are shown in Table 2. To eliminate the autocorrelation, we take Sweeps between configurations. The estimate of the statistical error by the jackknife method suggests that the autocorrelation between each configuration is well suppressed. The typical number of gauge configurations for each measurement at finite is .
To analyze the EMT with Eq. (3.0), we measure written in terms of the clover leaf representation on the lattice. To subtract out the contribution, , we carry out simulations on a lattice for each . Note that this vacuum subtraction is required only for the measurement of the trace anomaly , while no subtraction is needed for the measurement of the entropy density . For in and given by Eqs. (3.0) and (3.0), we use the fourloop running coupling with the scale parameter determined by the ALPHA Collaboration, [26]. In the following, we present the numerical results of the thermodynamics for two Settings in Secs. 6 and 7, respectively, and then discuss the correlation function of the EMT in Setting 2 in Sec. 8.
6 Thermodynamics: Setting 1
In this section, we first present the numerical results on bulk thermodynamics obtained by the Setting 1 [12]. In Fig. 1, we show the results for the dimensionless trace anomaly and the dimensionless entropy density at as a function of the dimensionless flow parameter [12]. The bold bars denote the statistical errors, while the thin (light color) bars show the statistical and systematic errors including the uncertainty of . We found that the statistical error is dominant in the small region for both and , while the systematic error originating from the scale parameter becomes significant for in the large region.
The fiducial window discussed in Sec. 4 is indicated by the dashed lines in Fig. 1. The lower limit of the window, beyond which the lattice discretization error grows, is set to be , where we take into account the fact that our clover leaf operator extends the size . The upper limit of the window, beyond which the smearing by the gradient flow exceeds the temporal lattice size, is set to be .
The data in Fig. 1 show, within the error bars, that (i) the plateau appears inside the preset fiducial window () for each , and (ii) the plateau extends to the smaller region as increases or equivalently as decreases. It should be remarked that only gauge configurations are necessary to obtain these results. Similar plateaues as in Fig. 1 also appear inside the fiducial window for other temperatures, and , with comparable error bars. These features imply that the double extrapolation is indeed possible.
Our lattice results at fixed with three different lattice spacings allow us to take the continuum limit. First, we pick up a flow time which is in the middle of the fiducial window. Then we extract and for each set of and . In the left panel of Fig. 2, resultant values taking into account the statistical errors (bold error bars) and the statistical plus systematic errors (thin error bars) are shown. The lattice data for with the same lattice setup at and in Ref. [4] are also shown by the cross (green) symbols in the top panel; our results with gauge configurations have substantially smaller error bars at these points.
We note that the continuum extrapolation in this analysis [12] is taken with fixed ; this procedure is graphically shown in Fig. 3 (a). Strictly speaking, however, the double limit has to be taken as mentioned above. Although we have checked that different choices of do not change the final results within the error bar as long as it is in the plateau region in the analysis in Setting 1, in the next section we will see that the limit with fixed shows a deviation compared with the double limit when the statistics is improved.
The horizontal axis of the left panel of Fig. 2, , is a variable suited for making continuum extrapolation of the thermodynamic quantities [4]. We consider two ways of extrapolation: A linear fit with the data at , , and (the solid lines in the left panel of Fig. 2), and a constant fit with the data at and (the dashed lines in the left panel of Fig. 2). In both fits, the correlation between the errors due to the common systematic error from is taken into account. The former fit is used to determine the central value in the continuum limit whose error is within even at our lowest temperature. The latter is used to estimate the systematic error from the scaling violation whose typical size is at high temperature and at low temperature.
We have analyzed various systematic errors; the perturbative expansion of , the running coupling , the scale parameter, and the continuum extrapolation. We found that the dominant errors in the present lattice setup are those from and the continuum extrapolation, which are included in the left panel of Fig. 2. To reduce these systematic errors, finer lattices are quite helpful: They make the plateau in wider by reducing the lower limit of the fiducial window, so that the continuum extrapolation becomes easier. Also, larger aspect ratio would be helpful to guarantee the thermodynamic limit.
Finally, we plot, in the right panel of Fig. 2, the continuum limit of and obtained by the linear fit of the , , and data (the solid lines) in the left panel of Fig. 2 for , , and . For comparison to the existing data of the gauge theory after the continuum extrapolation, the results of Refs. [4] [20] and [22] obtained by the integral method are also plotted in the right panel of Fig. 2. The results of the two different approaches are consistent with each other within sigma level. The figure, however, indicates that our results underestimate the observables compared with the previous ones. This discrepancy in part come from the continuum extrapolation with fixed . The correct double limit will be discussed in the next section.
7 Thermodynamics: Setting 2
Next, let us see the numerical results on thermodynamics with finer lattices with Setting 2. In the left panel of Fig. 4, we show the dimensionless entropy density as a function of for with and for several values of . In this result, one can take a closer look at the dependence of thanks to the improved statistics and finer lattices. In the previous section, we discussed that a plateau is observed in the fiducial window, , in Setting 1. Because the statistics is significantly improved in Setting 2, in Fig. 4 one sees a dependent structure in the fiducial window with . Next, while the dependence at small becomes more stable as the lattice spacing becomes finer, the dependence does not show a plateau even on the finest lattice with . In the right panel of Fig. 4, the same result is shown as a function of instead of . The figure indicates that behaves as a linear function of as the lattice spacing becomes finer.
The origin of the linear term can be understood as the contribution of higher dimensional operators in Eq. (3.0) giving rise to contribution to this formula. In fact, the contribution of dimension six operators are proportional to in limit. Figure 4 shows that the effect of these terms has nonnegligible contribution. Such a contribution has to be subtracted to measure thermodynamic quantities. A similar result is obtained for , while the dependence is weaker than the one in .
The effect discussed here can be eliminated by taking the double limit . In order to take this double extrapolation, in this analysis we take limit with fixed ; this procedure is graphically shown in Fig. 3 (b). Here, the ratio has to be sufficiently large so that and satisfy the condition , or . Once this condition is satisfied, however, the double limit can be taken safely irrespective of the value of . In this extrapolation, it is expected that the discretization effect is of order . This is because the leading discretization error on thermodynamic observables is order [4], while the error arising from nonzero is of order .
In Fig. 5, we show the result of the continuum extrapolation for and with fixed for and . We choose for the horizontal axis and take a linear extrapolation because of the dependence discussed above. The result shows that the continuum extrapolation is stable and the extrapolations with different are consistent. Because of the negative slope in Fig. 4, the continuum extrapolation of obtained in this method gives slightly larger value compared with the result in the previous section. The final results after the double extrapolation are consistent with the existing results in the integral method [4, 22].
8 Correlation function
Once the operator of the EMT on the lattice is obtained, the application of this operator is not limited to the bulk thermodynamics which is given by expectation values of one point functions. In this section, as an example of other applications we consider the correlation function of the EMT. As discussed in Sec. 1, the correlation functions of the EMT are important observables related to transport properties of the medium and fluctuations.
In this manuscript, as an example of such analyses we present the numerical result on the Euclidean correlation function of the total energy,
(8.0) 
with and
When one investigates correlation functions of the EMT with Eq. (3.0), one has to take care of the cooling nature of the gradient flow. When correlation functions are observed on this field, the smearing nature gives rise to artificial correlations in two point correlators when the two operators are not separated twice the smearing length. Because of Eq. (2.0), the smearing radius defined by the mean square distance in four dimensional space is [10]. When one focuses on a smearing along one definite direction, on the other hand, the smearing length is . Therefore, the correlation function Eq. (8.0) receives an artificial effect due to the smearing for .
Keeping this precaution in mind, let us see the numerical result on Eq. (8.0). In Fig. 6, we show dependence of the correlator Eq. (8.0) for several values of the flow time for and . From the figure, one can see several interesting features. First, takes a large value for in which the abovementioned smearing problem affects the correlation function. The enhancement of at has a clear dependence, which indicates that this enhancement is indeed the artificial behavior arising from the gradient flow. On the other hand, the correlator for does not have a statistically significant dependence. This result implies that the correlation function can be obtained in this range. Second, the errorbars of the result are more suppressed for larger . For , the error is less than with configurations. Compared with the previous studies on the EMT correlators [3, 2], this number of statistics is significantly small.
Figure 6 also shows that the correlator Eq. (8.0) takes a constant within statistics in the range at which the smearing effect is well suppressed. This result can be understood as a consequence of the energy conservation, because the energy conservation requires
(8.0) 
for . The emergence of a plateau in Fig. 6, therefore, shows that the EMT defined in Eq. (3.0) satisfies the energy conservation.
The constant of gives the thermal fluctuation of the total energy in grand canonical ensemble, i.e.
(8.0) 
where the value of in this formula must satisfy in order to avoid the smearing problem. In the continuum theory the value of in Eq. (8.0) is arbitrary except for because of Eq. (8.0).
Combining Eq. (8.0) with the linearresponse relation Eq. (1.0), one can make an estimate of the specific heat by reading the value of in the plateau region in Fig. 6. Here, we use the value of at midpoint . The specific heat is then obtained as
(8.0) 
In the analysis of the specific heat of the gauge theory in Ref. [27] using the differential method, the values are estimated as for and for . Although we have not taken the continuum extrapolation, it is notable that our result in Eq. (8.0) is consistent with the one in Ref. [27]. This result opens a possibility to measure the specific heat in a novel method using the linearresponse relation Eq. (1.0).
9 Summary and outlook
In this study, we analyzed the expectation values and correlation functions of the energymomentum tensor (EMT) on the lattice using the EMT operator Eq. (3.0) defined by the gradient flow and the small flow time expansion. Our results show that the numerical analyses of the EMT on the lattice with Eq. (3.0) work quite successfully. In fact, the thermodynamics observables analyzed in our approach agree with the previous results using integral method. Besides, the result on the energyenergy correlator is consistent with the conservation law and the linearresponse relation Eq. (1.0). Moreover, our approach can perform these measurements with a significantly small statistics compared with the previous methods.
As discussed at the beginning of this manuscript, the measurement of the EMT on the lattice has been a difficult subject, although the EMT is one of the most fundamental quantities in physics. Since the EMT defined by Eq. (3.0) can become a solution of this longstanding problem, there are a lot of applications of the present study. Among them, the accurate measurement of the correlation functions of the EMT is one of the most interesting subject. While we have only presented the numerical result on the energyenergy correlator in this manuscript, the same analysis can be performed for various channels. In particular, the correlation functions of the spatial components of the EMT, and , are interesting because they are related with shear and bulk viscosity, respectively, via Kubo formulas. The measurements of these correlators with Eq. (3.0) thus would enable us a stable analysis of these quantities on the lattice. Furthermore, the measurements of the energy and momentum distributions in hadrons and flux tube will also be possible using the EMT operator in Eq. (3.0). These analyses will reveal various aspects of the hadron structure and confinement. The application of this method to full QCD with fermion fields is another important future study [16, 17]. These analyses will be reported elsewhere.
Acknowledgments
Numerical simulation for this study was carried out on NEC SX8R and SX9 at RCNP, Osaka University, and Hitachi SR16000 and IBM System Blue Gene Solution at KEK under its LargeScale Simulation Program (Nos. T1204 and 13/1420). The work of M. A., M. K., and H. S. are supported in part by a GrantinAid for Scientific Researches 23540307 and 26400272, 25800148, and 23540330, respectively. E. I. is supported in part by Strategic Programs for Innovative Research (SPIRE) Field 5. T. H. is partially supported by RIKEN iTHES Project.
References
 [1] S. Caracciolo, G. Curci, P. Menotti, and A. Pelissetto, The Energy Momentum Tensor for Lattice Gauge Theories, Annals Phys. 197 (1990) 119.
 [2] H. B. Meyer, Transport Properties of the QuarkGluon Plasma: A Lattice QCD Perspective, Eur.Phys.J. A47 (2011) 86, [arXiv:1104.3708].
 [3] A. Nakamura and S. Sakai, Transport coefficients of gluon plasma, Phys.Rev.Lett. 94 (2005) 072305, [heplat/0406009].
 [4] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, et al., Thermodynamics of SU(3) lattice gauge theory, Nucl.Phys. B469 (1996) 419–444, [heplat/9602007].
 [5] L. D. Landau and E. M. Lifshitz, Statistical Physics. ButterworthHeinemann, 1980.
 [6] R. Kubo, Statistical Mechanics. North HollandAmsterdam, 1965.
 [7] P. Romatschke, New Developments in Relativistic Viscous Hydrodynamics, Int.J.Mod.Phys. E19 (2010) 1–53, [arXiv:0902.3663].
 [8] G. Bali, K. Schilling, and C. Schlichter, Observing long color flux tubes in SU(2) lattice gauge theory, Phys.Rev. D51 (1995) 5165–5198, [heplat/9409005].
 [9] H. Suzuki, Energymomentum tensor from the YangMills gradient flow, PTEP 2013 (2013), no. 8 083B03, [arXiv:1304.0533].
 [10] M. Luscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 1008 (2010) 071, [arXiv:1006.4518].
 [11] M. Luscher and P. Weisz, Perturbative analysis of the gradient flow in nonabelian gauge theories, JHEP 1102 (2011) 051, [arXiv:1101.0963].
 [12] FlowQCD Collaboration Collaboration, M. Asakawa, T. Hatsuda, E. Itou, M. Kitazawa, and H. Suzuki, Thermodynamics of SU(3) gauge theory from gradient flow on the lattice, Phys.Rev. D90 (2014), no. 1 011501, [arXiv:1312.7492].
 [13] M. Luscher, Future applications of the YangMills gradient flow in lattice QCD, arXiv:1308.5598.
 [14] A. Ramos, “Wilson flow and renormalization.” this proceedings.
 [15] L. Del Debbio, A. Patella, and A. Rago, Spacetime symmetries and the YangMills gradient flow, JHEP 1311 (2013) 212, [arXiv:1306.1173].
 [16] M. Luscher, Chiral symmetry and the Yang–Mills gradient flow, JHEP 1304 (2013) 123, [arXiv:1302.5246].
 [17] H. Makino and H. Suzuki, Lattice energymomentum tensor from the YangMills gradient flow – inclusion of fermion fields, PTEP 2014 (2014), no. 6 063B02, [arXiv:1403.4772].
 [18] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, et al., The lattice gradient flow at treelevel and its improvement, JHEP 1409 (2014) 018, [arXiv:1406.0827].
 [19] A. Patella, L. Del Debbio, and A. Rago, Poincare symmetries and the YangMills gradient flow, PoS LATTICE2013 (2014) 324.
 [20] CPPACS Collaboration Collaboration, M. Okamoto et al., Equation of state for pure SU(3) gauge theory with renormalization group improved action, Phys.Rev. D60 (1999) 094510, [heplat/9905005].
 [21] T. Umeda, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya, et al., Fixed Scale Approach to Equation of State in Lattice QCD, Phys.Rev. D79 (2009) 051501, [arXiv:0809.2842].
 [22] S. Borsanyi, G. Endrodi, Z. Fodor, S. Katz, and K. Szabo, Precision SU(3) lattice thermodynamics for a large temperature range, JHEP 1207 (2012) 056, [arXiv:1204.6184].
 [23] ALPHA collaboration Collaboration, M. Guagnelli, R. Sommer, and H. Wittig, Precision computation of a lowenergy reference scale in quenched lattice QCD, Nucl.Phys. B535 (1998) 389–402, [heplat/9806005].
 [24] FlowQCD Collaboration Collaboration, M. Asakawa, T. Hatsuda, T. Iritani, E. Itou, M. Kitazawa, and H. Suzuki. in preparation.
 [25] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al., Highprecision scale setting in lattice QCD, JHEP 1209 (2012) 010, [arXiv:1203.4469].
 [26] ALPHA Collaboration Collaboration, S. Capitani, M. Luscher, R. Sommer, and H. Wittig, Nonperturbative quark mass renormalization in quenched lattice QCD, Nucl.Phys. B544 (1999) 669–698, [heplat/9810063].
 [27] R. V. Gavai, S. Gupta, and S. Mukherjee, The Speed of sound and specific heat in the QCD plasma: Hydrodynamics, fluctuations and conformal symmetry, Phys.Rev. D71 (2005) 074013, [heplat/0412036].