QCD Thermodynamics on the Lattice from the Gradient Flow
Abstract
To obtain the precise values of the bulk quantities and transport coefficients in quarkgluonplasma phase, we propose that a direct calculation of the renormalized energymomentum tensor (EMT) on the lattice using the gradient flow. From onepoint function of EMT, authors in Ref. [1] obtained the interaction measure and thermal entropy. The results are consistent with the one obtained by the integral method. Based on the success, we try to measure the twopoint function of EMT, which is related to the transport coefficients. Advantages of our method are (1) a clear signal because of the smearing effects of the gradient flow and (2) no need to calculate the wave function renormalization of EMT. In addition, we give a short remark on a comparison of the numerical cost between the positive and adjointflow methods for fermions, needed to obtain the EMT in the (2+1) flavor QCD.
QCD Thermodynamics on the Lattice from the Gradient Flow
\FullConferenceThe 26th International Nuclear Physics Conference
1116 September, 2016
Adelaide, Australia
1 Introduction: energymomentum tensor on lattice and gradient flow
One of the most important tasks of the present heavyion physics is to determine the thermal properties of the quarkgluonplasma (QGP) phase. In this context, it is necessary to precisely calculate the thermal quantities, namely entropy, pressure and shear viscosity. In our works, we try to determine these quantities directly from the calculations of the energymomentum tensor (EMT).
Measurements of the EMT using the lattice numerical simulation have at least two difficulties: One is a conceptual difficulty. The lattice regularization manifestly breaks the general covariance, while EMT is a generator of the corresponding invariance. The other one is a numerical cost. Since the signal of twopoint function for the EMT operator becomes very noisy due to its nonzero vacuum expectation value, it is too costly to determine its renormalization factor.
In this work, we obtain the EMT using a new technique [1] based on the small flowtime expansion of the YangMills gradient flow [2, 3]. Advantages of the usage of the gradient flow are following:

At finite flowtime, we can define the “correctly renormalized EMT” from lattice data in the continuum limit

Signals become much better because of the smearing effects of the gradient flow

It is not necessary to calculate the wave function renormalization of the EMT operator thank to its UV finiteness (in quenched QCD)[4]
In this proceeding, we briefly review of the basic idea to obtain the bulk thermal quantities, namely interaction measure (trace anomaly) and thermal entropy density, from direct calculation of the onepoint function of EMT. Next, we show the twopoint function of EMT, and its flowtime dependence. In the last section, we give a short remark for the numerical costs in the case of the(2+1) flavor QCD. We compare the numerical costs between the positive and adjointflow methods proposed in Ref. [5] to solve the gradient flow for fermions.
2 Review: thermal quantities from onepoint function of EMT
In Ref. [1], one of authors (E.I.) obtained the integration measure and thermal entropy of the pure YangMills theory in finite temperature from the direct calculation of EMT on the lattice. The key relationship is given in Ref. [3] for quenched QCD in small flowtime expansion: the correctlynormalized EMT can be defined by
(1) 
where is vacuum expectation value (v.e.v.) and is the correctlynormalized conserved EMT with its v.e.v. subtracted. Here and denotes gaugeinvariant local products of dimension and they are UV finite for the positive flowtime (). Explicitly, and . Here represents the field strength constructed by the flowed gauge field (), that is a solution to the gradient flow equation as
(2) 
where denotes the original quantum gauge field variable.
The contributions from the operators of dimension or higher are suppressed for small , and the coefficients are calculated perturbatively in [3].
In paper [1], we perform the numerical simulation and obtain the bulk quantity in the quenched QCD. We utilize the Wilson plaquette gauge action and set of () for one fixed physicaltemperature to take the continuum limit. Here denotes the lattice bare coupling constant and has onetoone correspondence with a lattice spacing [6]. The number of gauge configurations for the measurements at each lattice parameter is only –. Statistical errors are estimated by the jackknife method.
In Fig. 1, and are plotted after taking the continuum limit for , , and , where denote energy density and pressure, respectively. For comparison, results of Ref. [7] obtained by the integral method are shown by blue solid lines in Fig. 1. The results of the two different approaches are consistent with each other within statistical errors.
The integral method essentially calculate the free energy of thermodynamics, and is based on the macroscopic picture in finitetemperature QCD. On the other hand, our method is based on the microscopic picture, namely the quantum field theory. It is the first numerical confirmation of the consistency between micro and macroscopic pictures of the QGP phase in (quenched) QCD.
3 Twopoint function of EMT
3.1 Transport coefficients from EMT
We now move on the calculation of the twopoint function of EMT. It is related to the shear and bulk viscosity, and here we focus on the former one, which is given by the correlation function of component. The Euclidean correlator, which can be measured on the lattice, is defined as
(3) 
which is expressed in terms of the corresponding spectral functions () as
(4) 
The shear viscosity in QGP phase is then given by
(5) 
There are several works [8, 9, 10], where the correlation function of EMT are calculated on lattice. In these works, we explained before, there exit at least two difficulties, the renormalization of the lattice bare EMT operator and the bad signal to noise ratio of the quantity. In Ref. [9], the author introduced the oneloop latticeperturbative factor to define the renormalized EMT, while in Ref. [10] they estimate factor for a diagonal component of EMT from the thermal entropy shown in Fig. 1. Our method do not need the calculation of factor if we use the renormalized coupling constant in the coefficient , thank to the UV finiteness of the flowedcomposite operators. The signal to noise ratio is also drastically improved by the gradient flow because of its smearing effects, as we will show.
3.2 Lattice setup
We consider the Wilson plaquette gauge action on and lattices with a fixed . The lattice bare coupling constant for each is tuned to realize the temperature . The corresponding for each is determined by the relation in ALPHA Collaboration [6].
Gauge configurations are generated without dynamical fermions by the pseudoheatbath algorithm with the overrelaxation. We call one pseudoheatbath update sweep plus several overrelaxation sweeps as a Sweep. To eliminate the autocorrelation, we take Sweeps between measurements. The number of gauge configurations for the measurements at finite T is , , and for and , respectively. Statistical errors are estimated by the jackknife method.
3.3 Results
Firstly, we show the improvement of the statistical uncertainty by the usage of the gradient flow. For the correlation function of operator, the similar results are reported in Ref. [10]. Here, we show the results for correlator, which is equal to the one for in the continuum limit.
In Fig. 2, we plot the correlation function of operator, which is defined by the clover leaf on the lattice, without the gradient flow (right panel) and with the gradient flow (left panel). Here the number of the measured configurations for each color in both panels is the same.
Although the data should be positive by definition and indeed so in Ref. [8] with high statistics, some unflowed data in Fig. 2 become negative due to large noises in this statistics. On the other hand, at the finite flowtime (we take with ), the correlation function is positive at all despite low statistics, demonstrating that signals are highly improved.
Figure 3 shows the correlation function for the renormalized operator at finite flowtime, which includes the dependent coefficient .
The discrepancy among data for different comes from the discretization error of the correlator in our formulation. We found that it is small and looks undercontrolled in whole range of the imaginarytime.
We also show the flowtime dependence of the shape of the correlation function.
In Fig. 4, we found that the slope of the correlation function becomes milder in the longer flowtime. This is natural, since in the large flowtime limit the operator of are smeared in whole temporal direction and then the correlation becomes a constant as a function of . We also plot the bestfit function of the unflowed data in Ref. [9] by obtaining the BreitWigner fit ansatz. Although the data in Ref. [9] is correlator with the different renormalization process, the magenta curve might be consistent with the limit of our floweddata.
3.4 Outlook and future plan
Finally, we put our future plans to obtain the shear viscosity using our method. One possible method is to extract the spectral function from the correlation function by the fit with some assumptions on the functional form of . In this case, we have to remove data at , which suffer from the (over)smearing effects. From , we obtain the at several finite flowtimes. Since the flowtime dependence of the EMT twopoint function looks stronger than that of the onepoint function (See Fig. 1 in Ref. [1]), we expect that the flowtime dependence of is also larger, so that a careful estimation is necessary.
4 (2+1)flavor QCD calculation in the positive and adjointflow
4.1 Positiveflow and adjointflow to solve the gradient flow equation for fermions
In this section, we would like to give a short remark on the calculation of the EMT components including the fermion fields for full QCD system. The basic strategy and some preliminary results before taking the continuum limit are provided in Ref. [11, 12]. In the numerical calculation in these papers, the simulation cost to solve the gradient flow for fermions using the adjointflow method are high. In this proceedings, we compere the numerical costs between two methods, namely the adjointflow and positiveflow[5], when we solve the gradient flow equation for the fermion to obtain the thermal quantities.
The small flowtime formula for the definition of the EMT in full QCD system is given in Ref. [13]. To obtain the EMT, we have to calculate following twotypes of expectation values.
(6)  
(7) 
where is the lattice volume in lattice unit, denotes a label of the quark flavor and . are solutions of the fermion flow equation, which is given by
(8) 
Note that the covariant derivative refers to the flowed gauge field at the flow time .
To solve this equation and obtain the expectation value of composite operators, twotypes of methods are proposed in Ref. [5]. The first one is the “positiveflow” method, in which we introduce the random source field at , while the second one is the “adjointflow” method where we introduce the random source at the flowtime at which we want to measure the observables and inversely solve the equation from the finite flowtime to .
The technical procedures of the adjoint flow and positive flow are summarized as follow. If we consider to obtain the vacuum expectation values of EMT at , as calculated in Ref. [11] using the adjointflow methods, we carry out the following steps.
Ad1 Solve the gauge flow and store the flowedlink variable with small flowtime interval.
Ad2 Generate the noise vector at the finite flowtime where we want to obtain the expectation value of the observables, and solve the backward flow from the finite flowtime to zero.
Ad3 Calculate the propagator at using obtained pseudofermion vector in Ad2 and calculate the expectation value of .
Here, to reduce the simulation cost of Ad2 procedure, we firstly generate and store the flowed configurations with the small interval in Ad1 (e.g. in Ref. [11] ). Steps Ad1 and Ad2 are flavor independent, so that the flavor dependent part is the calculation of the propagator in Ad3. On the other hand, the actual procedure in the positiveflow method is as follows.
Po1 Generate the noise vector and calculate the propagator for each flavor at .
Po2 Solve the flow equation for link and fermionfields simultaneously toward the positive flowtime direction and calculate the expectation value at finite flowtime.
4.2 Comparison of numerical costs between the positive and adjointflow
Now, we show the simulation costs for each procedure. All calculations have been done using MPI processes of CPU on HITACHI SR16000. We utilize a configuration generated by the Iwasakigauge action and the improved Wilson fermion. The lattice parameters of the configuration are , , and on lattice.
Table 1 shows the computation cost for each procedure. Here the length of flow time is fixed as for Ad1,Ad2 and Po2. To solve the gradient flow, we use the thirdorder RungeKutta algorithm with .
In Ad3 and Po1, the calculation of propagator is included, so that the computational cost of this part depends on the configuration. Here we compare them using the same configuration and the same convergence precision to solve the inverse of the Dirac operator.
Procedure  comp. time [sec]  Procedure  comp. time [sec] 

Ad1  14  Po1  323 
Ad2  846  Po2  1,022 
Ad3  456 
If we calculate the expectation value at flowtime , the computational cost in the positive flow for the flavor QCD is
(9) 
Here, in Po2, we can continuously solve the gradient flow equation from to , so that the overhead of I/O of configurations can be reduced rather than the times that of the time shown in Table 1. On the other hand, in the adjointflow we have to recursively carry out the procedure Ad2 due to the backward flow is unstable. To reduce the simulation costs and to know the flowtime dependence, for instance, we take data point between with the interval as shown in Ref. [11]. We repeatedly carry out the Ad2 procedure times, and then the total computational cost is
(10) 
If we need the longer flowtime simulation, the cost of the positiveflow is further cheaper than the one for the adjointflow method.
References
 M. Asakawa et al. [FlowQCD Collaboration], Phys. Rev. D 90, no. 1, 011501 (2014) [Phys. Rev. D 92, no. 5, 059902 (2015)] [arXiv:1312.7492 [heplat]].
 M. Lüscher, JHEP 1008, 071 (2010) [arXiv:1006.4518 [heplat]].
 H. Suzuki, PTEP 2013, 083B03 (2013) [PTEP 2015, 079201 (2015)] [arXiv:1304.0533 [heplat]].
 M. Lüscher and P. Weisz, JHEP 1102, 051 (2011) [arXiv:1101.0963 [hepth]].
 M. Luscher, JHEP 1304, 123 (2013) [arXiv:1302.5246 [heplat]].
 M. Guagnelli et al. [ALPHA Collaboration], Nucl. Phys. B 535, 389 (1998) [heplat/9806005].
 S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP 1207, 056 (2012) [arXiv:1204.6184 [heplat]].
 A. Nakamura and S. Sakai, Phys. Rev. Lett. 94, 072305 (2005) [heplat/0406009].
 H. B. Meyer, Phys. Rev. D 76, 101701 (2007) [arXiv:0704.1801 [heplat]].
 S. W. Mages, S. Borsányi, Z. Fodor, A. Schäfer and K. Szabó, PoS LATTICE 2014, 232 (2015).
 E. Itou, H. Suzuki, Y. Taniguchi and T. Umeda, PoS LATTICE 2015, 303 (2016) [arXiv:1511.03009 [heplat]].
 Y. Taniguchi, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, T. Umeda and N. Wakabayashi, arXiv:1609.01417 [heplat].
 H. Makino and H. Suzuki, PTEP 2014, 063B02 (2014) [PTEP 2015, 079202 (2015)] [arXiv:1403.4772 [heplat]].