Measurement of thermodynamics using gradient flow
We analyze bulk thermodynamics and correlation functions of the energy-momentum tensor in pure Yang-Mills gauge theory using the energy-momentum 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 energy-momentum tensor can be successfully defined and observed on the lattice with moderate numerical costs with the gradient flow.
Measurement of thermodynamics using gradient flow
Department of Physics, Kyushu University, 6-10-1 Hakozaki, Higashi-ku, Fukuoka, 812-8581, Japan
In this manuscript, we consider the energy-momentum tensor
As is well known, the energy-momentum tensor (EMT) is one of the most fundamental quantities in physics. It is fundamental because this operator is related to the symmetry of space-time; 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,
with a choice of discretized definitions for the field strength on the lattice, does not give the correct EMT even for off-diagonal components in the sense that it does not correspond to the EMT in the continuum limit . 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 non-perturbatively, 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
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  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 
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 
where is the Kubo’s canonical correlation and . The transport coefficients are important quantities to understand dynamical time evolution of relativistic heavy ion collisions . The measurement of transport coefficients with Kubo formula, however, requires real-time correlation function of the EMT, while the lattice simulations can obtain only imaginary time correlators directly [3, 2]. The analytic continuation from imaginary-time 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  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.  on the basis of the Yang-Mills gradient flow  and the so-called small flow time expansion . 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 .
2 Gradient flow
Because we use the Yang-Mills gradient flow to define the EMT, let us first give a brief review on this concept. The gradient flow for the Yang-Mills gauge field is the continuous transformation of the field defined by the differential equation 
with the Yang-Mills 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
Neglecting the gauge dependent terms, Eq. (2.0) is the diffusion equation in four-dimensional space. For positive , therefore, the gradient flow acts as the cooling of the gauge field with the smearing radius . In Ref. , 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 . 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
where on the right-hand 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 space-time.
The relation Eq. (3.0) is reminiscent of the operator product expansion (OPE). While a product of operators at different space-time 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 . For this purpose, we first consider dimension-four gauge-invariant operators for operators of the left-hand side in Eq. (3.0). In pure gauge theory, there are two such operators;
Although Eqs. (3.0) and (3.0) are quite similar to the traceless-part 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 right-hand 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 gauge-invariant operators are the dimension-four EMTs. Up to this order, therefore, the small flow time expansions of Eqs. (3.0) and (3.0) are given by
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 .
The coefficients and are calculated perturbatively in Ref.  as
Here denotes the running gauge coupling in the scheme with the choice, , and
There are two important observations: (i) The right-hand 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 space-time region of radius , the statistical noise in calculating the right-hand 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 gauge-invariant 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 space-time 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 non-perturbative 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 four-dimensional 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.
The numerical analyses reported in the following have been performed with the following two settings.
Setting 1: The first analysis is performed on lattices . Gauge configurations are generated by the pseudo-heatbath algorithm with the over-relaxation, mixed in the ratio of . We call one pseudo-heatbath update sweep plus five over-relaxation 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 . The resultant values of are then converted to by using the result at in Ref. . 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  by the measurements of  and  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 four-loop running coupling with the scale parameter determined by the ALPHA Collaboration, . 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 . 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 . 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.  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  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 . 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.   and  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 , 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,
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 . 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 above-mentioned 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
The constant of gives the thermal fluctuation of the total energy in grand canonical ensemble, i.e.
Combining Eq. (8.0) with the linear-response 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
In the analysis of the specific heat of the gauge theory in Ref.  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. . This result opens a possibility to measure the specific heat in a novel method using the linear-response relation Eq. (1.0).
9 Summary and outlook
In this study, we analyzed the expectation values and correlation functions of the energy-momentum 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 energy-energy correlator is consistent with the conservation law and the linear-response 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 energy-energy 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.
Numerical simulation for this study was carried out on NEC SX-8R and SX-9 at RCNP, Osaka University, and Hitachi SR16000 and IBM System Blue Gene Solution at KEK under its Large-Scale Simulation Program (Nos. T12-04 and 13/14-20). The work of M. A., M. K., and H. S. are supported in part by a Grant-in-Aid 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.
-  S. Caracciolo, G. Curci, P. Menotti, and A. Pelissetto, The Energy Momentum Tensor for Lattice Gauge Theories, Annals Phys. 197 (1990) 119.
-  H. B. Meyer, Transport Properties of the Quark-Gluon Plasma: A Lattice QCD Perspective, Eur.Phys.J. A47 (2011) 86, [arXiv:1104.3708].
-  A. Nakamura and S. Sakai, Transport coefficients of gluon plasma, Phys.Rev.Lett. 94 (2005) 072305, [hep-lat/0406009].
-  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, [hep-lat/9602007].
-  L. D. Landau and E. M. Lifshitz, Statistical Physics. Butterworth-Heinemann, 1980.
-  R. Kubo, Statistical Mechanics. North Holland-Amsterdam, 1965.
-  P. Romatschke, New Developments in Relativistic Viscous Hydrodynamics, Int.J.Mod.Phys. E19 (2010) 1–53, [arXiv:0902.3663].
-  G. Bali, K. Schilling, and C. Schlichter, Observing long color flux tubes in SU(2) lattice gauge theory, Phys.Rev. D51 (1995) 5165–5198, [hep-lat/9409005].
-  H. Suzuki, Energy-momentum tensor from the Yang-Mills gradient flow, PTEP 2013 (2013), no. 8 083B03, [arXiv:1304.0533].
-  M. Luscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 1008 (2010) 071, [arXiv:1006.4518].
-  M. Luscher and P. Weisz, Perturbative analysis of the gradient flow in non-abelian gauge theories, JHEP 1102 (2011) 051, [arXiv:1101.0963].
-  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].
-  M. Luscher, Future applications of the Yang-Mills gradient flow in lattice QCD, arXiv:1308.5598.
-  A. Ramos, “Wilson flow and renormalization.” this proceedings.
-  L. Del Debbio, A. Patella, and A. Rago, Space-time symmetries and the Yang-Mills gradient flow, JHEP 1311 (2013) 212, [arXiv:1306.1173].
-  M. Luscher, Chiral symmetry and the Yang–Mills gradient flow, JHEP 1304 (2013) 123, [arXiv:1302.5246].
-  H. Makino and H. Suzuki, Lattice energy-momentum tensor from the Yang-Mills gradient flow – inclusion of fermion fields, PTEP 2014 (2014), no. 6 063B02, [arXiv:1403.4772].
-  Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, et al., The lattice gradient flow at tree-level and its improvement, JHEP 1409 (2014) 018, [arXiv:1406.0827].
-  A. Patella, L. Del Debbio, and A. Rago, Poincare symmetries and the Yang-Mills gradient flow, PoS LATTICE2013 (2014) 324.
-  CP-PACS 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, [hep-lat/9905005].
-  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].
-  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].
-  ALPHA collaboration Collaboration, M. Guagnelli, R. Sommer, and H. Wittig, Precision computation of a low-energy reference scale in quenched lattice QCD, Nucl.Phys. B535 (1998) 389–402, [hep-lat/9806005].
-  FlowQCD Collaboration Collaboration, M. Asakawa, T. Hatsuda, T. Iritani, E. Itou, M. Kitazawa, and H. Suzuki. in preparation.
-  S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al., High-precision scale setting in lattice QCD, JHEP 1209 (2012) 010, [arXiv:1203.4469].
-  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, [hep-lat/9810063].
-  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, [hep-lat/0412036].