QCD Thermodynamics on the Lattice from the Gradient Flow

TL;DR

This paper introduces a lattice QCD method using gradient flow to accurately compute the energy-momentum tensor, enabling precise determination of bulk and transport properties of quark-gluon plasma with reduced computational complexity.

Contribution

It presents a novel approach for calculating the renormalized energy-momentum tensor on the lattice via gradient flow, improving signal clarity and eliminating the need for wave function renormalization.

Findings

Consistent results for interaction measure and thermal entropy with existing methods

Successful measurement of two-point functions related to transport coefficients

Comparison of numerical costs between different fermion flow methods

Abstract

To obtain the precise values of the bulk quantities and transport coefficients in quark-gluon-plasma phase, we propose that a direct calculation of the renormalized energy-momentum tensor (EMT) on the lattice using the gradient flow. From one-point 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 two-point 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 adjoint-flow methods for fermions, needed to obtain the EMT in the (2+1) flavor QCD.

Tables

Table 1. Table 1: Simulation cost for each procedure. The length of flow-time is Δ ​ t / a 2 = 0.1 Δ 𝑡 superscript 𝑎 2 0.1 \Delta t/a^{2}=0.1 in Ad-1 , Ad-2 and Po-2 .

Procedure	comp. time [sec]	Procedure	comp. time [sec]
Ad-1	14	Po-1	323
Ad-2	846	Po-2	1,022
Ad-3	456

Equations

$$T_{\mu \nu }^R(x)=\underset{t\to 0}{\lim }\left\{\frac{1}{{\alpha }_U(t)}{U}_{\mu \nu }(t,x)+\frac{{\delta }_{\mu \nu }}{4{\alpha }_E(t)}\left[E(t,x)-{⟨E(t,x)⟩}_0\right]\right\},$$

$${\partial }_t{B}_{\mu }=D_{\nu }{G}_{\nu \mu },B_{\mu }(t=0,x)=A_{\mu }(x),$$

$$C(\tau )=\frac{1}{T^5}\int d\vec{x}⟨T_{12}(0,\vec{0})T_{12}(\tau ,\vec{x})⟩,$$

$$C(\tau )=\frac{1}{T^5}{\int }_0^{\infty }d\omega \rho (\omega )\frac{\cosh \omega \left(\frac{1}{2T}-\tau \right)}{\sinh (\frac{\omega }{2T})}.$$

$$\eta (T)=\pi {\frac{d\rho }{d\omega }|}_{\omega =0}.$$

$$t_{\mu \nu }^f(t)$$

$$s^f(t)$$

$${\partial }_t{\chi }^f(t,x)$$

$${\partial }_t{\bar{\chi }}^f(t,x)$$

$$\text{mbox}positive-flow:21,104[\sec ].$$

$$\text{mbox}adjoint-flow:196,180[\sec ].$$

Full text

1 Introduction: energy-momentum tensor on lattice and gradient flow

One of the most important tasks of the present heavy-ion physics is to determine the thermal properties of the quark-gluon-plasma (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 energy-momentum 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 two-point 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 flow-time expansion of the Yang-Mills gradient flow [2, 3]. Advantages of the usage of the gradient flow are following:

• •

At finite flow-time, 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 one-point function of EMT. Next, we show the two-point function of EMT, and its flow-time 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 adjoint-flow methods proposed in Ref. [5] to solve the gradient flow for fermions.

2 Review: thermal quantities from one-point function of EMT

In Ref. [1], one of authors (E.I.) obtained the integration measure and thermal entropy of the pure Yang-Mills 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 flow-time expansion: the correctly-normalized EMT can be defined by

[TABLE]

where $\langle\cdot\rangle_{0}$ is vacuum expectation value (v.e.v.) and $T_{\mu\nu}^{R}(x)$ is the correctly-normalized conserved EMT with its v.e.v. subtracted. Here $U_{\mu\nu}(t,x)$ and $E(t,x)$ denotes gauge-invariant local products of dimension $4$ and they are UV finite for the positive flow-time ($t>0$). Explicitly, $U_{\mu\nu}(t,x)\equiv G_{\mu\rho}(t,x)G_{\nu\rho}(t,x)-\delta_{\mu\nu}E(t,x)$ and $E(t,x)\equiv\frac{1}{4}G_{\mu\nu}(t,x)G_{\mu\nu}(t,x)$. Here $G_{\mu\nu}$ represents the field strength constructed by the flowed gauge field ($B_{\mu}(t,x)$), that is a solution to the gradient flow equation as

[TABLE]

where $A_{\mu}(x)$ denotes the original quantum gauge field variable.

The contributions from the operators of dimension $6$ or higher are suppressed for small $t$, and the coefficients $\alpha_{U},\alpha_{E}$ 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 $3$-set of ($\beta,N_{\tau}$) for one fixed physical-temperature to take the continuum limit. Here $\beta=6/g_{0}^{2}$ denotes the lattice bare coupling constant and has one-to-one correspondence with a lattice spacing [6]. The number of gauge configurations for the measurements at each lattice parameter is only $100$–$300$. Statistical errors are estimated by the jackknife method.

In Fig. 1, $(\epsilon-3P)/T^{4}=\langle\sum_{\mu=1}^{4}T_{\mu\mu}\rangle/T^{4}$ and $(\epsilon+P)/T^{4}=\langle T_{11}-T_{44}\rangle/T^{4}$ are plotted after taking the continuum limit for $T/T_{c}=1.65$, $1.24$, and $0.99$, where $\epsilon,P$ 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 finite-temperature 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 macro-scopic pictures of the QGP phase in (quenched) QCD.

3 Two-point function of EMT

3.1 Transport coefficients from EMT

We now move on the calculation of the two-point 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 $T_{12}$ component. The Euclidean correlator, which can be measured on the lattice, is defined as

[TABLE]

which is expressed in terms of the corresponding spectral functions ($\rho(\omega)$) as

[TABLE]

The shear viscosity in QGP phase is then given by

[TABLE]

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 one-loop lattice-perturbative $Z(g_{0})$ factor to define the renormalized EMT, while in Ref. [10] they estimate $Z$-factor for a diagonal component of EMT from the thermal entropy shown in Fig. 1. Our method do not need the calculation of $Z$-factor if we use the renormalized coupling constant in the coefficient $\alpha_{U},\alpha_{E}$, thank to the UV finiteness of the flowed-composite 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 $N_{\tau}=8,10,12$ and $16$ lattices with a fixed $N_{s}/N_{\tau}=4$. The lattice bare coupling constant $\beta$ for each $N_{\tau}$ is tuned to realize the temperature $T=1.65T_{c}$. The corresponding $\beta$ for each $N_{\tau}$ is determined by the relation in ALPHA Collaboration [6].

Gauge configurations are generated without dynamical fermions by the pseudo-heatbath algorithm with the over-relaxation. We call one pseudo-heatbath update sweep plus several over-relaxation sweeps as a Sweep. To eliminate the autocorrelation, we take $200$ Sweeps between measurements. The number of gauge configurations for the measurements at finite T is $3000$, $5000$, $2500$ and $800$ for $N_{\tau}=8,10,12$ and $16$, 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 $(T_{11}-T_{22})/2$ operator, the similar results are reported in Ref. [10]. Here, we show the results for $T_{12}$ correlator, which is equal to the one for $(T_{11}-T_{22})/2$ in the continuum limit.

In Fig. 2, we plot the correlation function of $U_{12}$ 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 un-flowed data in Fig. 2 become negative due to large noises in this statistics. On the other hand, at the finite flow-time (we take $r_{\mathrm{smear}}T=0.25$ with $r_{\mathrm{smear}}=\sqrt{8t}$), the correlation function is positive at all $\tau$ despite low statistics, demonstrating that signals are highly improved.

Figure 3 shows the correlation function for the renormalized $T_{12}$ operator at finite flow-time, which includes the $\beta$ dependent coefficient $1/\alpha_{U}^{2}$.

The discrepancy among data for different $N_{\tau}$ comes from the discretization error of the $T_{12}$ correlator in our formulation. We found that it is small and looks under-controlled in whole range of the imaginary-time.

We also show the flow-time 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 flow-time. This is natural, since in the large flow-time limit the operator of $T_{12}$ are smeared in whole temporal direction and then the correlation becomes a constant as a function of $\tau$. We also plot the best-fit function of the un-flowed data in Ref. [9] by obtaining the Breit-Wigner fit ansatz. Although the data in Ref. [9] is $(T_{11}-T_{22})/2$ correlator with the different renormalization process, the magenta curve might be consistent with the $t\rightarrow 0$ limit of our flowed-data.

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 $\rho(\omega)$ from the correlation function by the fit with some assumptions on the functional form of $\rho(\omega)$. In this case, we have to remove data at $\tau/N_{\tau}\leq r_{\mathrm{smear}}T$, which suffer from the (over-)smearing effects. From $\rho(\omega)$, we obtain the $\eta/T^{3}$ at several finite flow-times. Since the flow-time dependence of the EMT two-point function looks stronger than that of the one-point function (See Fig. 1 in Ref. [1]), we expect that the flow-time dependence of $\eta/T^{3}$ is also larger, so that a careful estimation is necessary.

4 (2+1)-flavor QCD calculation in the positive- and adjoint-flow

4.1 Positive-flow and adjoint-flow 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 adjoint-flow method are high. In this proceedings, we compere the numerical costs between two methods, namely the adjoint-flow and positive-flow[5], when we solve the gradient flow equation for the fermion to obtain the thermal quantities.

The small flow-time 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 two-types of expectation values.

[TABLE]

where $N_{\Gamma}$ is the lattice volume in lattice unit, $f$ denotes a label of the quark flavor and $\Delta\chi^{f}(t,a)=D_{\mu}D_{\mu}\chi^{f}(t,x)$. $\chi(t,x),\bar{\chi}(t,x)$ are solutions of the fermion flow equation, which is given by

[TABLE]

Note that the covariant derivative refers to the flowed gauge field at the flow time $t$.

To solve this equation and obtain the expectation value of composite operators, two-types of methods are proposed in Ref. [5]. The first one is the “positive-flow” method, in which we introduce the random source field at $t=0$, while the second one is the “adjoint-flow” method where we introduce the random source at the flow-time at which we want to measure the observables and inversely solve the equation from the finite flow-time to $t=0$.

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 $t/a^{2}=2.0$, as calculated in Ref. [11] using the adjoint-flow methods, we carry out the following steps.

Ad-1 Solve the gauge flow and store the flowed-link variable with small flow-time interval.

Ad-2 Generate the noise vector at the finite flow-time where we want to obtain the expectation value of the observables, and solve the backward flow from the finite flow-time to zero.

Ad-3 Calculate the propagator at $t=0$ using obtained pseudo-fermion vector in Ad-2 and calculate the expectation value of $t_{\mu\nu}(t),s(t)$.

Here, to reduce the simulation cost of Ad-2 procedure, we firstly generate and store the flowed configurations with the small interval in Ad-1 (e.g. $\Delta t/a^{2}=0.1$ in Ref. [11] ). Steps Ad-1 and Ad-2 are flavor independent, so that the flavor dependent part is the calculation of the propagator in Ad-3. On the other hand, the actual procedure in the positive-flow method is as follows.

Po-1 Generate the noise vector and calculate the propagator for each flavor at $t=0$.

Po-2 Solve the flow equation for link- and fermion-fields simultaneously toward the positive flow-time direction and calculate the expectation value at finite flow-time.

4.2 Comparison of numerical costs between the positive- and adjoint-flow

Now, we show the simulation costs for each procedure. All calculations have been done using $32$-MPI processes of $1$ CPU on HITACHI SR16000. We utilize a configuration generated by the Iwasaki-gauge action and the $N_{f}=2+1$ $O(a)$-improved Wilson fermion. The lattice parameters of the configuration are $\beta=1.900$, $c_{SW}=1.715$, $\kappa_{u,d}=0.1368$ and $\kappa_{s}=0.1364$ on $32^{3}\times 8$ lattice.

Table 1 shows the computation cost for each procedure. Here the length of flow time is fixed as $\Delta t/a^{2}=0.1$ for Ad-1,Ad-2 and Po-2. To solve the gradient flow, we use the third-order Runge-Kutta algorithm with $\epsilon=0.01$.

In Ad-3 and Po-1, 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.

If we calculate the expectation value at flow-time $t/a^{2}=2.0$, the computational cost in the positive flow for the $(2+1)$-flavor QCD is

[TABLE]

Here, in Po-2, we can continuously solve the gradient flow equation from $t=0$ to $t/a^{2}=2.0$, so that the overhead of I/O of configurations can be reduced rather than the $20$ times that of the time shown in Table 1. On the other hand, in the adjoint-flow we have to recursively carry out the procedure Ad-2 due to the backward flow is unstable. To reduce the simulation costs and to know the flow-time dependence, for instance, we take $20$ data point between $0\leq t/a^{2}\leq 2.0$ with the interval $\Delta t/a^{2}=0.1$ as shown in Ref. [11]. We repeatedly carry out the Ad-2 procedure $210$ times, and then the total computational cost is

[TABLE]

If we need the longer flow-time simulation, the cost of the positive-flow is further cheaper than the one for the adjoint-flow method.