Energy-momentum tensor correlation function in Nf = 2 + 1 full QCD at finite temperature

We measure correlation functions of the nonperturbatively renormalized energy-momentum tensor in Nf = 2 + 1 full QCD at finite temperature by applying the gradient flow method both to the gauge and quark fields. Our main interest is to study the conservation law of the energy-momentum tensor and to test whether the linear response relation is properly realized for the entropy density. By using the linear response relation we calculate the specific heat from the correlation function. We adopt the nonperturbatively improved Wilson fermion and Iwasaki gauge action at a fine lattice spacing = 0.07 fm. In this paper the temperature is limited to a single value T ≃ 232 MeV. The u, d quark mass is rather heavy with mπ/mρ ≃ 0.63 while the s quark mass is set to approximately its physical value.


Introduction
The nucleons are expected to undergo a phase transition at high temperature and density and turn to be the quark gluon plasma (QGP). Several heavy ion collision experiments are going on and many phenomena have been observed which support the transition to QGP. One of the most fascinating observation is that of strongly coupled hydrodynamical collective motion, which is well described with the ideal fluid model [1].
At the same time there appear several challenges to calculate the viscosity of QGP by using the lattice QCD without quarks [2][3][4]. These pioneering works calculate the viscosity from correlation functions of the energy-momentum tensor via the spectral function by making use of Kubo formula.
However there are two difficulties in this strategy. One is the ill defined inverse problem to extract the continuum spectral function from the correlation function on lattice with finite temporal extent. The other is the renormalization of the energy-momentum tensor. The energy-momentum tensor is the conserved current associated with the translational invariance and is not uniquely given on lattice. Exception is diagonal components of the tensor for the pure gluon system, for which nonperturbative renormalization factors are given by using the thermodynamical relation.
In this paper we shall avoid the latter difficulty by applying the gradient flow [5][6][7] both to the gluon and quark fields and use it as a nonperturbative renormalization scheme [8,9]. The gradient flow renormalization scheme can be applied every component of the energy-momentum tensor even with quarks. The method has already been applied to the finite temperature QCD with N f = 2 + 1 flavors [10,11] and successful results are given for the equation of state of QCD, the chiral condensate and the topological susceptibility. The formulation is also used for a calculation of the energy-momentum tensor correlation function in pure Yang-Mills theory at finite temperature [12]. We apply the strategy to N f = 2 + 1 flavors QCD and calculate the energy-momentum tensor correlation function. The purpose in this paper is to investigate (i) conservation law, (ii) restoration of rotational symmetry, (iii) consistency between several derivation of the entropy density and (iv) specific heat.

Nonperturbative renormalization using gradient flow
We adopt the flow equations given in Refs. [5,7] for the gauge and quark fields where the field strength and the covariant derivative are given in terms of the flowed gauge field f = u, d, s, denotes the flavor index. The statement of the gradient flow is that any composite operator made of flowed fields is already renormalized and is free from the UV divergences if multiplied with the wave function renormalization factor for quark fields [5,7]. The renormalization scale is given in terms of the flow time µ = 1/ √ 8t. The gradient flow can be used as a nonperturbative renormalization scheme when applied to the lattice QCD Monte Carlo simulation. An important manipulation in a nonperturbative renormalization is a conversion to commonly used schemes like the MS or MS scheme, which are defined in perturbation theory. This can be accomplished by two steps; (i) in the nonperturbative scheme visit a perturbative high energy region nonperturbatively for example by using the step scaling function, (ii) convert to the MS or MS scheme with a matching factor calculated perturbatively.
For the gradient flow scheme the first step is accomplished by following the flow towards the t → 0 limit, which is realized easily in numerical simulation. The matching factor needed for the second step is calculated in Refs. [8,9] for the energy-momentum tensor at one loop level. According to the strategy the properly normalized energy-momentum tensor which satisfies the Ward-Takahashi identity in the continuum limit is given by taking a limit of the flowed tensor operator where ⟨· · · ⟩ 0 stands for the vacuum expectation value (VEV), i.e., the expectation value at zero temperature. The flowed tensor operator is given by a linear combination of five operatorsÕ iµν (t, x) with matching coefficients c i (t) given in Refs. [8,9] using MS scheme 1 .
We notice that five operators are expected to be evaluated on lattice and the continuum limit is taken before the t → 0 limit in (8). The t → 0 limit is necessary in order to resolve a mixing with irrelevant dimension six operators which is proportional to the flow time t.

Observables
In this paper we shall calculate two point correlation functions of the fluctuation of the energymomentum tensor where and τ is the Euclidean time. ⟨· · · ⟩ T is the one point function at the same temperature. The renormalized energy-momentum tensor consists of gluon contribution given by the first line of (9) and that from quarks given by the second and third lines in (9). In two point correlation functions, they lead to contributions categorized into two types; (i) two point functions. The first terms are calculated by using the same technique used in Ref. [10], where the noise estimator is adopted for the quark one point functions by inserting the noise vector at flow time t. In this paper we put the noise vector at flow time t = 0 in order to save the numerical cost.
For the second contribution from connected quark diagrams we need the flowed quark propagator [7].
where c fl is the O(a) improvement factor and we adopt c fl = 1/2 at tree level. � S and � C are given by where v and w stands for a apace-time point at t = 0. S f (v, w) is the bare quark propagator at t = 0. K and K † is the flow kernel which satisfies the ordinary flow equation and the adjoint flow equation [7] ( In an evaluation of the propagator (12)

Numerical results
Measurements of the energy-momentum tensor are performed on N f = 2+1 gauge configurations generated for Ref. [13]. As given in (9) we need to subtract the zero-temperature values of the operators. The zero temperature gauge configurations are also prepared which were generated for Ref. [14].
The nonperturbatively O(a)-improved Wilson quark action and the renormalization-group improved Iwasaki gauge action are adopted. The bare coupling constant is set to β = 2.05, which corresponds to a = 0.0701(29) fm (1/a ≃ 2.79 GeV). The hopping parameters are set to κ u = κ d ≡ κ ud = 0.1356 and κ s = 0.1351, which correspond to heavy u and d quarks, m π /m ρ ≃ 0.63, and almost physical s quark, m η ss /m ϕ ≃ 0.74. See Ref. [10] for a detailed explanation of numerical parameters. The temperature is limited to a single value T ≃ 232 MeV.
In this section we perform three tests for the conservation law, the rotational symmetry and the linear response relations. We calculate the specific heat by using the linear response relation.

Conservation law
Spatial integral of temporal component of the energy-momentum tensor is a conserved charge in the continuum limit and satisfies the conservation law Three corresponding correlation functions C 00;00 , C 20;20 and C 00;22 are plotted as a function of the Euclidean time in Fig. 1 at flow time t/a 2 = 0.5, 1.0, 1.5, 2.0. As the gradient flow goes ahead statistical quality of the correlation function is improved. In the middle of the Euclidean time 5 ≤ τ/a ≤ 7 the function is flat within one sigma with high statistical precision for t/a 2 ≥ 1.5, which we consider a realization of the conservation law on lattice. On the other hand the correlation function is far from being flat near the boundary where we set the point source. This is considered to be due to the a 2 /τ 2 lattice artifact.

Rotational symmetry
Under the three dimension rotational symmetry the spatial component of the correlation function should take the form In Fig. 2 we plot the left hand side of (19) as a function of the Euclidean time. The signal becomes better as we proceed the gradient flow. The linear combination is consistent with zero for all region within 1σ -1.5σ.

Linear response relations
In this subsection we test linear response relations using the entropy density. The entropy density is expressed in three different ways. First it is written as (ϵ + p)/T by using the Maxwell's relation and the integrable condition for the entropy which is given by the expectation value of one point function of the energy momentum tensor. On the other hand starting from a linear response of the pressure against a variation of the temperature the entropy density is given by a derivative of the energy momentum tensor Substituting the statistical physics relation for the expectation value we have the second expression in terms of the two point function of the energy momentum tensor The last representation is given by using a linear response to the infinitesimal Lorentz boost The correlation functions (23) and (24) are plotted as a function of the Euclidean time in the middle and right panel of Fig. 1 for i = 2. According to the discussion given in subsection 4.1 three data at 5 ≤ τ/a ≤ 7 are employed and fitted by a constant. The results are plotted in the left panel of Fig. 3 as a function of the flow time by blue down triangles (23) and red up triangles (24) for the entropy density divided by the temperature (ϵ + p)/T 4 . The contribution from the one point function (20) is averaged over space-time and is plotted by black circles.  After taking the t → 0 limit (8) three contributions give the entropy density (ϵ + p)/T 4 . We adopt the same strategy used in Refs. [10,11] to take the limit. The flowed operator (9) evaluated on the lattice behaves as follows as a function of the flow time where A µν appears due to the lattice artifact before taking the continuum limit. S µν and R µν are higher dimensional irrelevant operator. What we find in our data is the linear window where the nonlinear contributions a 2 /t and t 2 are negligible and data behaves linearly in the flow time. In the figure the linear window is indicated by the black and red vertical dotted lines for the one point function (20) and two point correlators (23), (24). The data are fitted linearly in t within the window to take the t → 0 limit. The entropy density given by the linear fit is plotted by the corresponding filled symbols near the origin. Three different representations of the entropy density are consistent with each other within the statistical error. This indicates a success of the evaluation the energy-momentum tensor correlation functions.

Specific heat
The specific heat is given by a response of the energy against a variation of the temperature Inserting the statistical physics relation it is given in terms of the two point function of the energy momentum tensor The correlation function is plotted in the left panel of Fig. 1. Here again data at 5 ≤ τ/a ≤ 7 are taken and fitted by a constant. The result are plotted in the right panel of Fig. 3 as a function of the flow time. As in the previous subsection we find the linear window indicated by the red vertical lines and perform the linear fit to take the t → 0 limit. The resultant specific heat is plotted by the filled black circle at the origin, whose value is c V = 38(27).

Conclusion
The two point correlation function of the renormalized energy-momentum tensor is calculated in N f = 2 + 1 QCD on lattice. On the lattice we confirmed the conservation law is realized apart from the point source. The spatial rotational symmetry is well shown within the statistical error. The nonperturbative renormalization on lattice is done by applying the gradient flow to the gluon and quark fields and taking a → 0 and then t → 0 limits. Both limits are realized by making use of the linear window and applying the linear fit. This procedure is tested using the entropy density. Three difference representations are consistent with each other in the t → 0 limit. Finally we calculate the specific heat by applying the method.
Future application shall be derivation of the bulk and shear viscosity using the corresponding two point correlation functions shown in Fig. 4. We are also planning to calculate the correlation function  on gauge configurations with physical quark mass, which is used for a derivation of equation of state in Ref. [15].