Latent heat and pressure gap at the first-order deconfining phase transition of SU(3) Yang-Mills theory using the small flow-time expansion method
Mizuki Shirogane, Shinji Ejiri, Ryo Iwami, Kazuyuki Kanaya, Masakiyo Kitazawa, Hiroshi Suzuki, Yusuke Taniguchi, Takashi Umeda
aa r X i v : . [ h e p - l a t ] J a n Preprint number: UTHEP-753, UTCCS-P-135, KYUSHU-HET-218, J-PARC-TH-0232
Latent heat and pressure gap at the first-order deconfining phasetransition of SU(3) Yang-Mills theory using the small flow-timeexpansion method
Mizuki Shirogane , Shinji Ejiri , Ryo Iwami , Kazuyuki Kanaya , Masakiyo Kitazawa ,Hiroshi Suzuki , Yusuke Taniguchi , and Takashi Umeda Graduate School of Science and Technology, Niigata University, Niigata 950-2181, Japan Department of Physics, Niigata University, Niigata 950-2181, Japan Track Maintenance of Shinkansen, Rail Maintenance 1st Department, East Japan RailwayCompany Niigata Branch, Niigata 950-0086, Japan Tomonaga Center for the History of the Universe, University of Tsukuba, Tsukuba,Ibaraki 305-8571, Japan Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies, KEK,203-1, Shirakata, Tokai, Ibaraki, 319-1106, Japan Department of Physics, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395,Japan Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577,Japan Graduate School of Humanities and Social Sciences, Hiroshima University,Higashihiroshima, Hiroshima 739-8524, Japan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
We study latent heat and pressure gap between the hot and cold phases at the first-order decon-fining phase transition temperature of the SU(3) Yang-Mills theory. Performing simulations onlattices with various spatial volumes and lattice spacings, we calculate the gaps of the energydensity and pressure using the small flow time expansion (SF t X) method. We find that thelatent heat ∆ ǫ in the continuum limit is ∆ ǫ/T = 1 . ± .
040 for the aspect ratio N s /N t = 8and 1 . ± .
038 for N s /N t = 6 at the transition temperature T = T c . We also confirm thatthe pressure gap is consistent with zero, as expected from the dynamical balance of two phasesat T c . From hysteresis curves of the energy density near T c , we show that the energy density inthe (metastable) deconfined phase is sensitive to the spatial volume, while that in the confinedphase is insensitive. Furthermore, we examine the effect of alternative procedures in the SF t Xmethod — the order of the continuum and the vanishing flow time extrapolations, and also therenormalization scale and higher order corrections in the matching coefficients. We confirm thatthe final results are all well consistent with each other for these alternatives. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Subject Index B01, B38, B64
PTP
TEX.cls
Introduction
First-order phase transitions appear in various important thermodynamic systems including thehigh density QCD and many-flavor QCD, and are associated with various phenomena such as phasecoexistence and latent heat. It is worth developing numerical techniques to investigate first-orderphase transitions. The SU(3) Yang-Mills theory, i.e., the quenched approximation of QCD, at finitetemperature has a first-order deconfining phase transition, and provides us with a good testingground for that.At a first-order phase transition point, metastable states corresponding to two phases coexistat the same time. To keep a dynamical balance between them, their pressures must be the same.Therefore, one can check the validity of computational methods by measuring the pressure in eachmetastable state. In fact, the problem of the nonvanishing pressure gap has played an importantrole in early development of lattice methods to study thermodynamic quantities [1–6]. On theother hand, the energy density takes different values in each phase. The difference between themis the latent heat, which is one of the fundamental observables characterizing the first-order phasetransition.In Refs. [4, 7], we studied the first-order transition of the SU(3) Yang-Mills theory using thederivative method [1]. We have numerically confirmed that the pressure gap vanishes when weadopt the nonperturbatively evaluated values of the Karsch coefficients (anisotropy coefficients),and then computed the latent heat using the nonperturbative Karsch coefficients.In this paper, we study the first-order phase transition of the SU(3) Yang-Mills theory adoptinga new technique to calculate thermodynamic observables: the small flow time expansion (SF t X)method based on the gradient flow [8, 9]. The gradient flow is an evolution of the fields in terms of afictitious “time” t [10–14]. Fields at positive flow-time t > √ t in four dimensions, and the operators constructed byflowed fields at t > t X method provides us with a generalmethod to correctly calculate any renormalized observables on the lattice. The SF t X method hasbeen successfully applied to the evaluation of thermodynamic quantities in the Yang-Mills gaugetheory [15–18] and in QCD with 2 + 1 flavors of dynamical quarks [19–22]. In the case of SU(3)Yang-Mills theory we study, agreement with other methods has been established within a 1%level [18]. The method has also been applied recently to study various observables associated withthe energy-momentum tensor [23–28].We apply the SF t X method to calculate the energy density and pressure separately in hot andcold phases at the first-order deconfining phase transition temperature of the SU(3) Yang-Millstheory. We perform simulations on lattices with several lattice spacings and spatial volumes tocarry out the continuum extrapolation and also to study the finite-volume effect. On confirmingthat the pressure gap is consistent with zero, we evaluate the latent heat.We also test several alternative procedures in the SF t X method: To obtain a physical result bythe SF t X method, a double extrapolation to the continuum and vanishing flow-time limits, a → t →
0, is required. When a proper fitting range avoiding small t singularities at a > t X method. We test the next-to-leading order (NLO) and next-to-next-to-leading order (NNLO)expressions as well as the dependence on the choice of the renormalization scale for the matchingcoefficients. Though the final results should be insensitive to the details of the matching coefficients,the width of the available fitting range for the extrapolations and thus the numerical quality of thefinal results are affected by the choice of them. We show that the final results for the latent heatand pressure gap are insensitive to the choice of these alternatives. We also show that the NNLOmatching coefficients help to reduce lattice artifacts in the calculation of the latent heat.In the next section, we introduce the SF t X method and several details of the method to bediscussed in this paper. Then, our simulations at the first-order transition of SU(3) Yang-Millstheory are explained in Sec. 3. In Sec. 4, we show the results of the latent heat and pressure gap,mainly using the NNLO matching coefficients, and test several alternative procedures in the SF t Xmethod. We then study the hysteresis of the energy density near the phase transition temperatureand its spatial volume dependence in Sec. 5. Our conclusions are given in Sec. 6. Appendix Ais devoted to showing the results of the latent heat and pressure gap with the NLO matchingcoefficients, with which we study the effect due to the truncation of the perturbative series for thematching coefficients. In Appendix B, we discuss the choice of lattice operator for the field strength.
Our flow equation for the gradient flow is the simplest one proposed in Ref. [12]. The “flowed”gauge field B aµ ( t, x ) at flow time t is obtained by solving the flow equation ∂ t B aµ ( t, x ) = D ν G aνµ ( t, x ) ≡ ∂ ν G aνµ ( t, x ) + f abc B bν ( t, x ) G cνµ ( t, x ) (1)with the initial condition B aµ (0 , x ) = A aµ ( x ), where G aµν ( t, x ) ≡ ∂ µ B aν ( t, x ) − ∂ ν B aµ ( t, x ) + f abc B bµ ( t, x ) B cν ( t, x ) (2)is the flowed field strength constructed from B aµ ( t, x ). Because Eq. (1) is a kind of diffusion equation,we can regard B aµ ( t, x ) as a smeared field of the original gauge field A aµ ( x ) over a physical range of √ t in four dimensions. It was shown that operators constructed from B aµ ( t, x ) (“flowed operators”)have no ultraviolet divergences nor short-distance singularities at finite and positive t . Therefore,the gradient flow defines a kind of renormalization scheme, which is formulated nonperturbativelyand thus can be calculated directly on the lattice.The small flow-time expansion (SF t X) method provides us with a general method to evaluateany renormalized observables in terms of flowed operators at small and positive flow time t , whichare strictly finite and thus can be calculated on the lattice without further renormalization [8].In asymptotic free theories such as QCD, we can apply the perturbation theory to calculate thecoefficients (“matching coefficients”) relating the renormalized observables to the flowed operatorsat small flow time t [13]. Contamination of O ( t ) terms can be removed by a vanishing flow-timeextrapolation t →
0. 3 .1 Energy-momentum tensor
To calculate the energy-momentum tensor by the SF t X method, we consider the following flowedoperators at flow time t which are gauge-invariant and local, U µν ( t, x ) ≡ G aµρ ( t, x ) G aνρ ( t, x ) − δ µν G aρσ ( t, x ) G aρσ ( t, x ) , (3) E ( t, x ) ≡ G aµν ( t, x ) G aµν ( t, x ) . (4)With these dimension-four operators, the correctly renormalized energy-momentum tensor T Rµν isgiven by [8] T Rµν ( x ) = lim t → { c ( t ) U µν ( t, x ) + 4 c ( t ) δ µν [ E ( t, x ) − h E ( t, x ) i ] } , (5)where h E ( t, x ) i is the zero-temperature subtraction. Here, the matching coefficients c ( t ) and c ( t )are calculated by the perturbation theory [13] as c ( t ) = 1 g ∞ X ℓ =0 k ( ℓ )1 ( µ, t ) (cid:20) g (4 π ) (cid:21) ℓ , c ( t ) = 1 g ∞ X ℓ =1 k ( ℓ )2 ( µ, t ) (cid:20) g (4 π ) (cid:21) ℓ , (6)where g = g ( µ ) is the renormalized gauge coupling in the MS scheme at the renormalization scale µ .Because F aµρ ( x ) F aνρ ( x ) = G aµρ ( t, x ) G aνρ ( t, x ) + O ( t ) in the tree-level approximation, we have k (0)1 = 1.On the other hand, there is no “ k (0)2 ” because the tree-level energy-momentum tensor is traceless—the trace anomaly emerges from the one-loop order. Therefore, the next-to-leading order ( NLO )expression for c ( t ) contains k (0)1 and k (1)1 , and that for c ( t ) contains k (1)2 and k (2)2 , while the next-to-next-to-leading order ( NNLO ) expression for c ( t ) contains terms up to k (2)1 , and that for c ( t )contains terms up to k (3)2 . We mainly adopt NNLO matching coefficients in this study. To study theeffect due to the truncation of perturbative series, we also perform calculation using NLO matchingcoefficients in Appendix A.The coefficients k (1) i in the one-loop level are calculated in Refs. [8, 9, 30], and k (2) i in the two-loop level in Refs. [29]. (See also Ref. [31].) As pointed out in Ref. [8], in pure gauge Yang-Millstheories, k ( ℓ +1)2 can be deduced by ℓ -loop coefficients using the trace anomaly. A concrete form of k (3)2 is given in Ref. [18]. Collecting these results for the case of pure gauge Yang-Mills theories, we Note that our convention for c ( t ) differs from that of Ref. [29]. Our c ( t ) corresponds to c ( t ) + (1 / c ( t )in Ref. [29]. k (1)1 ( µ, t ) = − β L ( µ, t ) − C A = C A (cid:18) − L ( µ, t ) − (cid:19) , (7) k (1)2 ( µ, t ) = 18 β = 1124 C A , (8) k (2)1 ( µ, t ) = − β L ( µ, t ) + C A (cid:18) −
14 482405 −
16 546135 ln 2 + 118710 ln 3 (cid:19) = C A (cid:18) − L ( µ, t ) −
14 482405 −
16 546135 ln 2 + 118710 ln 3 (cid:19) , (9) k (2)2 ( µ, t ) = 18 β − β C A = C A (cid:18) − (cid:19) , (10) k (3)2 ( µ, t ) = 18 β − β C A + β C A (cid:18) − L ( µ, t ) − − (cid:19) = C A (cid:18) − L ( µ, t ) − − (cid:19) , (11)where β = 113 C A , β = 343 C A , and β = 285754 C A (12)are the first three coefficients of the beta function in Yang-Mills theories, and L ( µ, t ) ≡ ln(2 µ t ) + γ E (13)was introduced in Ref. [29] with γ E the Euler-Mascheroni constant. The factor C A is the quadraticCasimir for the adjoint representation defined by f acd f bcd = C A δ ab , (14)and C A = N c for the gauge group SU( N c ). The energy density and the pressure are obtained from the diagonal elements of the energy-momentum tensor, ǫ = − (cid:10) T R ( x ) (cid:11) , p = 13 X i =1 , , h T Rii ( x ) i . (15)Separating configurations in the Monte Carlo time history into the hot and cold phases (see Sec. 3.2),and adopting the multipoint reweighting method [32, 33] to fine-tune the coupling parameter to thecritical point, we calculate the energy density and the pressure in each phase just at the transitiontemperature, to estimate the latent heat and the pressure gap defined by∆ ǫ = ǫ (hot) − ǫ (cold) , ∆ p = p (hot) − p (cold) , (16)where ǫ (hot / cold) and p (hot / cold) indicate ǫ and p in the hot and cold phases, respectively.5n this paper, we calculate these quantities in the following conventional combinations of thetrace anomaly and the entropy density:∆( ǫ − p ) T and ∆( ǫ + p ) T . (17)When ∆ p = 0 as expected, ∆( ǫ − p ) /T and ∆( ǫ + p ) /T should coincide with each other. Wenote that ∆( ǫ − p ) /T is computed by the operator E ( t, x ) with the matching coefficient c ( t ),while ∆( ǫ + p ) /T is computed by the operator U µν ( t, x ) with c ( t ). The matching coefficients c i ( t ) given in the previous subsection are written in terms of therenormalization scale µ . However, because the matching coefficients relate the energy-momentumtensor and flowed operators, both physical and thus both finite and independent of the renor-malization scale µ , the matching coefficients themselves should also be finite and independent of µ —the explicit µ dependence from L ( µ, t ) should be canceled by that of the running coupling g ,in principle. In practice, however, because the perturbative series for the matching coefficients aretruncated at a finite order, the choice of µ may affect the magnitude of systematic errors due tothe neglected higher order terms.In early studies with the SF t X method [15, 16, 19, 23], the µ d scale µ = µ d ( t ) ≡ √ t , L ( µ, t ) = − γ E (18)has been adopted because µ d is a natural scale of local flowed operators that have a physicalsmearing extent of √ t . On the other hand, in Ref. [29], it is argued that µ = µ ( t ) ≡ √ t e γ E ≃ . µ d ( t ) , L ( µ, t ) = 0 , (19)would be a good choice of µ because it keeps the magnitude of two-loop contributions small. The µ scale was tested in quenched QCD [18] and in 2 + 1-flavor QCD [21, 22], and it was reportedthat the µ scale may improve the signal when the lattice is not quite fine [21, 22].In this paper, we examine both choices, µ = µ ( t ) and µ d ( t ). From the difference in the resultsbetween these two choices, we learn the magnitude of the effect due to truncation of higher-orderterms in the perturbation theory. To obtain a physical result with the SF t X method, the continuum extrapolation a → t → a → t ,and then take t → t → a , reserving the continuum extrapolation for a later stage. In the following, wecall the way to first take t → a and then take a → method 1 ,” and the way tofirst take a → t and then take t → method 2 .”6n finite lattices, additional unwanted operators may contaminate the right-hand side of Eq. (5)due to the lattice artifacts. To the leading order of the lattice spacing a , we expect that the latticeoperator in the curly bracket of Eq. (5) is expanded as T µν ( t, x, a ) = T Rµν ( x ) + tS µν ( x ) + A µν ( x ) a t + C µν ( x )( aT ) + D µν ( x )( a Λ QCD ) + a S ′ µν ( x ) + O ( a , t ) , (20)where T Rµν is the physical energy-momentum tensor, S µν and S ′ µν are contaminations of dimension-six operators with the same quantum number, and A µν , C µν , and D µν are those from dimension-fouroperators. In higher orders of a , we also have terms proportional to ( a /t ) etc. [34]. In Eq. (20),the term proportional to a /t is singular at t = 0 and thus is problematic when we want to take t → a in method 1. When we take the a → a /t is removed in principle. However, in practice, the a /t term is also problematic in method 2 whenwe want to take a → t where the a /t term dominates the data—we cannot performthe a → t . Therefore, in both methods, we need to find a rangeof t at each a where the a /t and more singular terms are negligible. Let us call the range of t in which the a /t and more singular terms as well as the O ( t ) terms look negligible the “linearwindow” [19, 21].On the other hand, when we can find such a range of t and perform t → a → T µν ( t, a ) ≈ T Rµν + tS µν + a C ′ µν + O ( a , t ) in this range, the final results ofmethod 1 and method 2 for T Rµν should be identical. The consistency of both methods is thus agood test of the SF t X method. In the following, we adopt both method 1 and method 2 using datain linear windows and perform extrapolations linear in t and a . When both methods are consistentwith each other, we get strong support in our identification of the linear windows.Another issue in the study of thermodynamic properties is the possible dependence on thephysical volume of the system. We have to keep the physical volume fixed in the a → t → T is adjusted to the transitiontemperature T c whose value is a physical constant. On a lattice with the size N s × N t , the latticespacing defined as a = 1 / ( N t T c ) is thus changed by changing N t , and the spatial volume in physicalunits, V = ( N s a ) = N s / ( N t T c ) , is fixed when the aspect ratio N s /N t is fixed. In this paper, wetake the double extrapolation for the cases of N s /N t = 6 and 8 to study the finite-volume effect. We perform simulations of the SU(3) Yang-Mills theory with the standard Wilson action atseveral inverse gauge couplings β = 6 /g around the deconfining transition point β c . We havestudied lattices with temporal lattice sizes of N t = 8, 12 and 16 with several different spatial latticesizes N s . Some of the configurations are taken from Ref. [7]. Although we have studied N t = 6lattices too, we do not use them in this study because the linear window turned out not to be clearenough for N t = 6. Our simulation parameters are summarized in Table 1. The configurations aregenerated by a pseudo-heat-bath algorithm followed by over-relaxation updates. Each measurement7 Ω |Polyakov loop histogram t/a = 0.00 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35~| Ω |Polyakov loop histogram t/a = 2.00 Fig. 1
Histograms of the original Polyakov loop (left) and the flowed Polyakov loop after gradientflow to t/a = 2 . T = T c on the 96 ×
12 lattice.
0 0.005 0.01| Ω | 0 0.1 0.2 0.3 ~ | Ω | ( t/ a = . ) Fig. 2
Double distribution of the original Polyakov loop | Ω | and the flowed Polyakov loop | ˜Ω | t/a =2 . measured on the 96 ×
12 lattice. The probability observed at each simulation point wastranslated to that at the transition point by the multipoint reweighting method.is separated by N sep heat-bath sweeps. Data are taken at 3 to 6 β values for each ( N s , N t ), and arecombined using the multipoint reweighting method [32, 33]. The statistical errors are estimated bythe jackknife method with the bin size chosen such that the errors are saturated.We define the transition point as the peak position of the Polyakov loop susceptibility χ Ω = N s ( h Ω i − h Ω i ). The results of the critical point β c determined in this way are also listed in Table 1.We employ the discretized representation of the flow equation (1) defined from the Wilson gaugeaction [12]. The numerical solution of the flow equation is obtained by the third-order Runge-Kuttamethod. For the field strength G µν ( t, x ) in the observables, several candidates are available for the8 able 1 Simulation parameters: the spatial lattice size N s , the temporal lattice size N t , couplingparameter β = 6 /g and the critical coupling β c . The measurements have been performed N conf times at intervals of N sep sweeps. The data for 48 ×
8, 64 ×
8, and 64 ×
12 are taken from ourprevious study [7]. N s N t β β c N sep N conf
48 8 6.056 6.06160(18) 20 100006.058 100006.060 100006.062 100006.065 110006.067 1000064 8 6.0585 6.06247(14) 20 47506.061 1030006.063 125006.065 255006.068 8100048 12 6.333 6.33472(13) 50 525006.335 400006.337 5250064 12 6.332 6.33493(17) 50 60006.3335 175006.335 58006.3375 96006.339 880072 12 6.334 6.33531(10) 50 80006.335 80006.337 800096 12 6.334 6.33532(11) 50 205506.335 120006.336 1440096 16 6.543 6.54667(20) 50 43506.545 48006.547 4250128 16 6.544 6.54616(23) 200 27506.545 30006.546 30006.547 30006.548 3000 lattice operator. From a test of the plaquette and the clover-type representations for G µν ( t, x ) givenin Appendix B, we adopt the clover-type representation. To evaluate the latent heat and the pressure gap, we need to separate the configurations aroundthe first-order transition point into the hot and cold phases. In the left panel of Fig. 1, we show thehistogram of the absolute value of the Polyakov loop | Ω | obtained on the 96 ×
12 lattice, where Ωis averaged over the spatial coordinates. The coupling parameter β is adjusted to the transitionpoint β c using the multipoint reweighting method. The two peaks in Fig. 1 correspond to the hotand cold phases. As discussed in Ref. [7], though the two peaks of the Polyakov loop histogram are9ell separated, the peaks of the plaquette histogram are overlapping. It is thus meaningful to usethe Polyakov loop to classify configurations into the hot and cold phases.In the right panel of Fig. 1, we show the histogram of the flowed Polyakov loop | ˜Ω | t/a =2 . , whichis constructed with the flowed gauge field at t/a = 2 .
0, as a typical example. We find that, thoughthe scale of the horizontal axis is different, its shape is quite similar to that shown in the left panel.The same similarity is also observed at other values of t/a . This suggests a strong correlationbetween the original Polyakov loop Ω and the flowed Polyakov loop ˜Ω. In Fig. 2, we show the 2Ddistribution of ( | Ω | , | ˜Ω | t/a =2 . ) measured on the 96 ×
12 lattice. A darker tone means a higherprobability. We find that Ω and ˜Ω have a strong linear correlation with each other. These resultsshow that the gradient flow does not affect the separation, and thus there is no merit in using theflowed Polyakov loop for the phase separation.We thus follow the method of Ref. [7] to classify configurations into the hot and cold phases:First, we remove short-time-range fluctuations by averaging | Ω | over several configurations (101configurations, except for the 128 ×
16 lattice, for which we use 45 configurations) around thecurrent configuration. We then classify the configurations into hot, cold, and mixed phases bythe value of this time-smeared Polyakov loop. With this classification, we observe many flip-flopsbetween the hot and cold phases during our Monte Carlo steps, while mixed configurations in whichthe two phases coexist are found to be rare on our lattices. We discard the configurations in themixed phase in the following.After the phase separation, we carry out the gradient flow on each of the configurations tomeasure flowed operators given in Eqs. (3) and (4). We then combine their expectation valuesin each phase by the multipoint reweighting method to obtain their values just at the transitionpoint β c . In Figs. 3–5, we show the results of ∆( ǫ − p ) /T , ∆( ǫ + p ) /T and ∆ ǫ/T as functions of t/a , using the NNLO matching coefficients. The blue, green and red symbols are the results of∆( ǫ − p ) /T , ∆( ǫ + p ) /T and ∆ ǫ/T , respectively. The plots in the left panels are the resultsadopting the renormalization scale µ = µ and those in the right panels adopting µ = µ d . The rapiddecrease of these observables at small t signals the appearance of the O ( a /t ) lattice discretizationeffect discussed in Sec. 2.4. Similar behavior at similar values of t has been reported for ( ǫ − p ) /T and ( ǫ + p ) /T themselves in previous studies of SU(3) Yang-Mills theory [15, 16, 18]. We thusdisregard the data in this region in the t → a → t/a & .
0, while the O ( t ) error sometimes starts to contaminate at large t/a . We also note that, though the dependence on the renormalization scale µ is small, there is ageneral tendency that the µ scale (left panels) leads to a slightly smaller slope towards the t → µ d scale (right panels).It should be noted that, within the statistical errors, the results of ∆( ǫ − p ) /T and ∆( ǫ + p ) /T are consistent with each other in a wide range of t . This means that the pressure gap betweentwo phases ∆ p vanishes, which must be satisfied in the final results after the double extrapolation of10 x8 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x8 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x8 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x8 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. 3 ∆( ǫ − p ) /T (blue), ∆( ǫ + p ) /T (green), and ∆ ǫ/T (red) calculated on the 48 × × µ = µ and µ = µ d , respectively. The horizontal axis is the flow time in lattice units t/a . The symbols at t/a = 0 and the straight lines are the results of t → t → a →
0. This suggests that, for these observables, the errors from the lattice discretizationand the small flow-time expansion are well under control in these ranges of t .To study the influence of the truncation of perturbative series for the matching coefficients,we repeat the calculation with the NLO matching coefficients following the original procedureof Ref. [8]. The results are summarized in Appendix A. We find that the NNLO matching coefficientsdrastically improve the signal over the NLO matching coefficients in the sense that the ∆ p at finite t and a becomes much smaller and also the observables show smaller slope in t . As discussedin the next section, we also confirm that the results of ∆( ǫ − p ) /T , ∆( ǫ + p ) /T , and ∆ ǫ/T extrapolated to the t → x12 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. 4
The same as Fig. 3, but on the 48 ×
12, 64 ×
12, 72 ×
12, 96 ×
12 lattices (from topto bottom). 12 x16 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x16 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x16 NNLO µ ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x16 NNLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. 5
The same as Fig. 3, but on the 96 ×
16 (top) and 128 ×
16 (bottom) lattices. s3 NNLO µ ∆(ε +p)/T ∆(ε -3p)/T Fig. 6
Spatial volume dependence of ∆( ǫ + p ) /T (red) and ∆( ǫ − p ) /T (blue) with the NNLOmatching coefficients adopting µ = µ on N s ×
12 lattices with N s = 48 ,
64, 72, and 96.13 t2 N s /N t =8 NNLO µ ∆(ε +p)/T ∆(ε -3p)/T t2 N s /N t =8 NNLO µ d ∆(ε +p)/T ∆(ε -3p)/T t2 N s /N t =8 NLO µ d ∆(ε +p)/T ∆(ε -3p)/T Fig. 7 ∆( ǫ + p ) /T (square) and ∆( ǫ − p ) /T (circle) at t = 0 for the aspect ratio N s /N t = 8.The horizontal axis is 1 /N t = ( T c a ) . The open square and open circle at 1 /N t = 0 are the resultsof method 2, first taking a → t →
0. The left and middle panels are the results withNNLO matching coefficients adopting µ = µ and µ d , respectively. The right panel is the resultswith NLO matching coefficients adopting µ = µ d . t → followed by a → t → a →
0. We note that the data shown in Figs. 3–5are very linear in a wide range of t within the statistical error. In this study, we choose the fit rangeas follows: First, we require the fit range to satisfy tT c ≤ . t/a < t / defined in Ref. [19]). We also require t/a ≥ √ N t = 8, however, becausethis makes the available range too narrow, we relaxed it to the minimum requirement of t/a ≥ t → ǫ − p ) /T and ∆( ǫ + p ) /T are less than about one sigma for all lattices, in accordancewith the expectation that ∆ p should vanish.In Fig. 6, we plot the results of ∆( ǫ − p ) /T (red) and ∆( ǫ + p ) /T (blue) at t = 0 measuredon N s ×
12 lattices with N s = 48 ,
64, 72, and 96 as functions of the inverse spatial volume in latticeunits 1 /N s . The lattice spacing a = 1 / ( N t T c ) is the same for all data shown in Fig. 6. We find that∆( ǫ − p ) /T and ∆( ǫ + p ) /T decrease as the spatial volume increases. Therefore, we have to takethe finite-volume effect into account.We now perform the continuum extrapolation a → V . As discussedin Sec. 2.4, this is achieved by fixing the aspect ratio N s /N t in our study. In Fig. 7, we plot14 t2 N s /N t =6 NNLO µ ∆(ε +p)/T ∆(ε -3p)/T t2 N s /N t =6 NNLO µ d ∆(ε +p)/T ∆(ε -3p)/T t2 N s /N t =6 NLO µ d ∆(ε +p)/T ∆(ε -3p)/T Fig. 8
The same as Fig. 7, but for the aspect ratio N s /N t = 6. Table 2
The results of t → ǫ − p ) /T and ∆( ǫ + p ) /T on each lattice.Errors are statistical only. Systematic errors due to the choice of the fit range are smaller than thestatistical errors. N s N t Fit range NNLO with µ NNLO with µ d NLO with µ d t/a ǫ − p ) T ∆( ǫ + p ) T ∆( ǫ − p ) T ∆( ǫ + p ) T ∆( ǫ − p ) T ∆( ǫ + p ) T
48 8 1.0–1.6 1.256(34) 1.202(26) 1.270(34) 1.205(26) 1.266(34) 1.188(25)64 8 1.0–1.6 1.171(28) 1.132(20) 1.181(28) 1.135(20) 1.175(28) 1.111(19)48 12 1.4–3.6 1.449(17) 1.476(19) 1.461(17) 1.480(19) 1.451(18) 1.441(18)64 12 1.4–2.0 1.347(35) 1.375(40) 1.349(35) 1.379(40) 1.325(35) 1.307(40)72 12 1.4–3.6 1.222(45) 1.257(46) 1.232(46) 1.260(46) 1.230(46) 1.237(44)96 12 1.4–3.6 1.083(34) 1.120(34) 1.094(34) 1.123(34) 1.093(34) 1.101(34)96 16 1.4–6.4 1.283(36) 1.309(45) 1.290(36) 1.311(45) 1.291(37) 1.281(42)128 16 1.4–6.4 1.092(31) 1.118(41) 1.097(31) 1.121(41) 1.067(33) 1.065(41) ∆( ǫ − p ) /T and ∆( ǫ + p ) /T at t = 0 as functions of 1 /N t = ( T c a ) for N s /N t = 8. The left andmiddle panels are the results adopting µ = µ and µ d , respectively. We also show the results withNLO matching coefficients adopting µ = µ d , discussed in Appendix A. Corresponding results for N s /N t = 6 are shown in Fig. 8.We do the a → T c a ) = 1 /N t . Wefind that the lattice spacing dependence is small in these observables. However, we also note thatthe slight slope toward the continuum limit changes its sign between N s /N t = 8 and 6. This maybe suggesting that, when we consider various systematic errors, the lattice spacing dependence isapproximately absent within errors. The continuum-extrapolated values by the linear fit are plottedby the filled square and filled circle at 1 /N t = 0 in Figs. 7 and 8. The open square and open circleat 1 /N t ≈ ǫ/T , ∆( ǫ − p ) /T and ∆( ǫ + p ) /T , obtained using the NNLO orNLO matching coefficients and with the µ or µ d scale, are all very consistent with each other ateach aspect ratio. From the results of ∆( ǫ − p ) /T and ∆( ǫ + p ) /T with the NNLO matchingcoefficients and the µ scale, we find the pressure gaps to be∆ p/T = 0 . N s /N t = 8) , ∆ p/T = 0 . N s /N t = 6) , (21)which are only about 1% of the latent heat and are consistent with zero.15 ∆ ( ε - ) / T tT c2 N s /N t =8 NNLO µ continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =8 NNLO µ continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε - ) / T tT c2 N s /N t =8 NNLO µ d continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =8 NNLO µ d continuum limitN t = 8N t = 12N t = 16 Fig. 9 ∆( ǫ − p ) /T (left) and ∆( ǫ + p ) /T (right) with the NNLO matching coefficients adopt-ing µ = µ (top) and µ = µ d (bottom) on N s /N t = 8 lattices. Three orange curves represent theresult in the continuum limit. (Thick orange curves represent the central values and thin curvesrepresent the range of their errors.) a → followed by t → t, a ) → (0 , a → t in physical units, and then take the t → N s /N t to take the finite size effect into account. In Fig. 9, we plot the resultsof ∆( ǫ − p ) /T (left) and ∆( ǫ + p ) /T (right) as functions of the flow time in a physical unit, tT c = t/ ( aN t ) , on lattices with N s /N t = 8. Figure 10 shows corresponding results for N s /N t = 6.For these figures, the NNLO matching coefficients are used, and the results adopting µ = µ and µ d are given in the top and bottom panels of each figure, respectively. The red, green and blue symbolsare for the results of N t = 8, 12 and 16. Recall that the lattice spacing a = 1 / ( N t T c ) is varied by N t . From these figures, we find that the lattice spacing dependence becomes rapidly smaller as theflow time increases.We take the a → /N t at each flow time tT c . The orange curves in Figs. 9 and 10 are the results in the continuum limit. The center curverepresents the central value, and the upper and lower curves represent the statistical error. Forinformation, we also show the results obtained outside the linear window in these figures. In Fig. 11,we show the results of ∆( ǫ − p ) /T (circle) and ∆( ǫ + p ) /T (square) in the continuum limit as16 ∆ ( ε - ) / T tT c2 N s /N t =6 NNLO µ continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =6 NNLO µ continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε - ) / T tT c2 N s /N t =6 NNLO µ d continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =6 NNLO µ d continuum limitN t = 8N t = 12N t = 16 Fig. 10
The same as Fig. 9, but for N s /N t = 6. Table 3 ∆ ǫ/T , ∆( ǫ − p ) /T and ∆( ǫ + p ) /T in the continuum limit with the NNLO andNLO matching coefficients adopting the renormalization scale µ or µ d . The left and right threecolumns are results of method 1 ( t → a →
0) and method 2 ( a → t → N s N t ∆ ǫT ∆( ǫ − p ) T ∆( ǫ + p ) T ∆ ǫT ∆( ǫ − p ) T ∆( ǫ + p ) T NNLO µ µ d µ d µ µ d µ d tT c for N s /N t = 8 (filled symbols) and N s /N t = 6 (open symbols). The results adopting µ = µ and µ d are given in the upper left and upper right panels. We find that the dependence on µ is almost negligible in these quantities. In the bottom panel of Fig. 11, we also show the results17 c2 N s /N t = 8 NNLO µ ∆(ε +p)/T ∆(ε -3p)/T c2 N s /N t = 8 NNLO µ d ∆(ε +p)/T ∆(ε -3p)/T c2 N s /N t = 8 NLO µ d ∆(ε +p)/T ∆(ε -3p)/T Fig. 11 ∆( ǫ + p ) /T (red square) and ∆( ǫ − p ) /T (blue circle) in the continuum limit with theNNLO matching coefficients adopting µ = µ (top left) and µ = µ d (top right), and with the NLOmatching coefficients adopting µ = µ d (bottom). The filled symbols are the results of N s /N t = 8and the open symbols are those of N s /N t = 6. The symbols on the vertical axes are the results inthe continuum limit.using data with NLO matching coefficients with the µ d scale given in Appendix A. Comparing theupper and lower panels, we find that, by adopting the NNLO matching coefficients, the slopes in t are much reduced and the pressure gap ∆ p/T at finite t is much suppressed.We then perform t → t region where contamination of the O ( a /t ) terms is suspected. To esti-mate a systematic error due to the t → . < tT c < . . < tT c < . . < tT c < . . < tT c < . . < tT c < . tT c using data with the NNLO matching coefficients and the µ scale. Thetriangle, square and circle symbols are for ∆ ǫ/T , ∆( ǫ + p ) /T and ∆( ǫ − p ) at t = 0, and thefilled and open symbols are the results of N s /N t = 8 and 6, respectively. We find that the resultsadopting different fit ranges are very consistent with each other. In Fig. 12, we also show the resultsof method 1 obtained in the previous subsection at the left end of the plot, which are also consistentwith the results of method 2 within statistical errors.18 .80.911.11.21.31.41.5 ∆ ( ε + p )/ T ∆ ( ε -3 p )/ T ∆ ε / T ∆ p / T Fig. 12
Fit range dependence of ∆ ǫ/T (triangle), ∆( ǫ + p ) /T (square), ∆( ǫ − p ) /T (circle),and ∆ p/T (diamond), determined in the continuum and t = 0 limits by method 2, with the NNLOmatching coefficients adopting the µ scale. The filled symbols are the results of N s /N t = 8 andthe open symbols are those of N s /N t = 6. The numbers 1–5 on the horizontal axis are for the fitranges 1–5 explained in the text. For comparison, results of method 1 are also shown at the leftend of the plot.We note that ∆( ǫ + p ) /T and ∆( ǫ − p ) /T at t = 0 are very consistent with each other withinthe statistical errors. Similar results are obtained also with the µ d scale and the NLO matchingcoefficients.The results of the direct calculation of ∆ p are plotted by the diamond symbols in the bottompanel of Fig. 12. We find that the values of ∆ p are only about 1% of the latent heat and areconsistent with the results of method 1 given in Eq. (21). We also note that the values of ∆ p are even smaller than the errors of the latent heat. Due to correlation between ∆( ǫ + p ) /T and∆( ǫ − p ) /T , the jackknife statistical errors for ∆ p turned out to be quite small in comparisonto the errors by the error propagation formula using the errors of ∆( ǫ + p ) /T and ∆( ǫ − p ) /T .With the small errors, the mean values of ∆ p deviate from zero by about 2-3 σ statistical errorsdepending on the fit range. We find that these nonvanishing values in the t → N t = 8 used in the continuum extrapolation at each fixed flow time t : Asseen from Fig. 3, the values of ∆ p for N t = 8 clearly deviate from zero, while those for N t = 12 and16 are consistent with zero in the fit range of the t → N t = 8, we obtain ∆ p consistent with zero. Thus,taking account of the systematic error due to the continuum extrapolation which is larger than thestatistical error for the case of ∆ p , we think that ∆ p is consistent with zero.Our final results of method 2 for the latent heat are summarized in the right columns of Table 3,for which we adopt the results with the NNLO matching coefficients and the µ scale using the fitrange 2. Errors are statistical only. As shown in Fig. 12, systematic errors due to the choice of thefit range are smaller than the statistical errors. 19 ∆ ( ε - ) / T N s /N t NNLO µ N t = 681216continuum limit 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 4 5 6 7 8 9 10 11 ∆ ( ε + p ) / T N s /N t NNLO µ N t = 681216continuum limit Fig. 13
Comparison of ∆( ǫ − p ) /T (left) and ∆( ǫ + p ) /T (right) between the derivativemethod and SF t X method as functions of N s /N t = √ V T c computed by the NNLO calculationwith µ = µ . The open symbols are the result of the derivative method [7]. We shifted the value of N s /N t to avoid overlapping symbols. The original values of N s /N t are 4, 16 / ,
6, 8, and 32 / . t X method for the latent heat
Finally, we summarize our results of the SF t X method for the latent heat. From Table 3and Fig. 12, we see that the results of method 1 and method 2 agree well with each other. Wetake our central values from method 2 with the NNLO matching coefficients and the µ scale. Weobtain for N s /N t = 8 ∆( ǫ − p ) /T = 1 . , (22)∆( ǫ + p ) /T = 1 . , (23)∆ ǫ/T = 1 . , (24)and for N s /N t = 6 ∆( ǫ − p ) /T = 1 . , (25)∆( ǫ + p ) /T = 1 . , (26)∆ ǫ/T = 1 . . (27)Systematic errors estimated from the differences between method 1 and method 2 as well as amongdifferent fit ranges are smaller than the statistical errors quoted in these equations. Performingheuristic extrapolations of these two data points to the thermodynamic limit, we find ∆ ǫ/T ∼ . ± .
07 by a 1 /V linear fit and ∆ ǫ/T ∼ . ± .
03 by a constant fit.
Finally, we compare our results of the SF t X method with those obtained by the derivativemethod [7]. The computational cost for the SF t X method is much smaller than that for the deriva-tive method, since the nonperturbative calculation of the Karsch coefficients needed in the derivativemethod requires large systematic study with high statistics. We note that the statistical errors are20 ∆ ( ε - ) / T t2 NNLO µ N s /N t = 416/36832/3 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 0 0.005 0.01 0.015 0.02 0.025 0.03 ∆ ( ε + p ) / T t2 NNLO µ N s /N t = 416/36832/3 Fig. 14
Comparison of ∆( ǫ − p ) /T (left) and ∆( ǫ + p ) /T (right) between the SF t X methodand the derivative method as functions of 1 /N t = a T c . Filled symbols are the results of the SF t Xmethod with the NNLO matching coefficients and the µ scale, and open symbols are the resultsof the derivative method given in Ref. [7].reduced with the SF t X method thanks the smearing nature of the gradient flow and also theavoidance of the nonperturbative Karsch coefficients.In Ref. [7], using the derivative method, a rather low value of ∆ ǫ/T = 0 . ± .
17 was obtainedfor the latent heat in the continuum limit using data at N t = 6, 8 and 12. In Ref. [7], because thespatial volume dependence in the latent heat was found to be small in comparison with the statisticalerrors, which are much larger than the errors in the current study mainly due to the error from thenonperturbative Karsch coefficients, spatial volume dependence was assumed to be absent.To avoid uncertainties due to the continuum and infinite spatial volume extrapolations, wedirectly compare the results of ∆( ǫ − p ) /T and ∆( ǫ + p ) /T at finite lattice spacings and spatialvolumes in Fig. 13. The horizontal axis is N s /N t = √ V T c . The filled symbols are the results of theSF t X method in the t → µ scale, obtainedon the lattices with N t = 8 (red diamond), 12 (green circle) and 16 (blue triangle), respectively(Table 2). The magenta square symbols are for the results in the continuum limit at N s /N t = 6and 8 by method 2. The open symbols are the results of the derivative method obtained on thelattices with N t = 6 (black pentagon), 8 (red diamond) and 12 (green circle), respectively [7]. Thelarge statistical errors for ∆( ǫ + p ) /T with the derivative method are caused by the large errors ofnonperturbative Karsch coefficients which increase rapidly with N t . From these plots, we find thatthe SF t X method and the derivative method lead to quite similar results at each N s /N t .As discussed in Sec. 4.3, the lattice spacing dependence of the latent heat by the SF t X methodis small for N t = 8, 12, 16. The lattice spacing dependence in the results by the derivative method isalso small for N t = 8 and 12. The left and right panels of Fig. 14 are the results of ∆( ǫ − p ) /T and∆( ǫ + p ) /T as functions of 1 /N t = a T c , i.e., the lattice spacing squared normalized by T c . Thefilled and open symbols are the results by the SF t X method and the derivative method, respectively.The red diamonds and green circles are for N s /N t = 6 and 8, respectively. We note that the dataat N t = 6 with the derivative method, shown at the right end of these plots, show deviation fromthe general tendency of data on finer lattices, suggesting a large discretization error at N t = 6. As21uggested from Fig. 9 of Ref. [7], by removing the data at N t = 6, the derivative method may alsolead to a larger latent heat in the continuum, though a reliable extrapolation is not possible withonly two data points. c hotcold (ε +p)/T Fig. 15 ( ǫ + p ) /T by the SF t X method calculated using configurations in the hot phase (red),the cold phase (blue) and all configurations (green) with the NNLO matching coefficients and the µ scale, obtained on the 48 × × T c In Sec. 4, we found that the latent heat is clearly dependent on the spatial volume of the system.At a first-order phase transition point, however, because the correlation length does not diverge,the spatial volume dependence should be mild when the volume is sufficiently large. To find theorigin of the spatial volume dependence in our latent heat, we study ( ǫ + p ) /T separately in thehot and cold phases at temperatures near T c using the multipoint reweighting method.The results obtained on 48 × × t/a = 1 . t → t . The NNLO matching coefficients with the µ scale are used.The red symbols are for the results obtained in the hot phase and the blue symbols are for the coldphase, while the green symbols show the results without the phase separation. The horizontal axis isthe temperature T = 1 / ( N t a ) normalized by the critical temperature T c , where the relation betweenthe lattice spacing a and β is determined by the critical point β c as function of N t [7]. The resultsfor the hot phase below T c and for the cold phase above T c are those obtained in correspondingmetastable states. We find that the spatial volume dependence appears only in the (metastable)hot phase around T c , while no apparent volume dependence is visible in the (metastable) coldphase. We thus conclude that the spatial volume dependence of the latent heat is due to that inthe contribution of the hot phase. From this figure, we also note that latent heat is sensitive to thevalue of the critical point β c . A careful determination is required for β c .22 c hotcold (ε +p)/T Fig. 16
The same as Fig. 15, but by the derivative method.In Fig. 16, we show the corresponding plot by the derivative method, obtained on the 48 × × β c [4, 7], we use the Karsch coefficients obtainedat β c at other values of β around β c . We see similar spatial volume dependence of the results inthe (metastable) hot phase, though the statistical errors are large.To draw a definite statement about the nature of the metastable state, we need more statistics,because the effective number of configurations for the metastable states after the phase separationis not large in the reweighting calculation when T is not sufficiently close to T c . Further study withhigher statistics is needed to precisely determine the volume dependence of the hot phase, and thusto carry out a reliable extrapolation of the latent heat to the thermodynamic limit. To study the latent heat and the pressure gap at the first-order deconfining phase transitiontemperature in the SU(3) Yang-Mills theory, we performed simulations on lattices with variousspatial volumes and lattice spacings around the phase transition temperature T c . Separating gen-erated configurations into hot and cold phases using the value of the Polyakov loop and adjustingthe temperature to T c by the multipoint reweighting method, we calculated the energy density andthe pressure gap in each phase at T c , using the small flow-time expansion (SF t X) method based onthe gradient flow.We confirmed that the pressure gap between the hot and cold phases is consistent with zero,as expected from the dynamical balance of two phases at T c . Our results for the latent heat in thecontinuum limit are summarized in Sec. 4.4. We obtained ∆ ǫ/T = 1 . ± .
040 for the spatialvolume corresponding to the aspect ratio N s /N t = 8, and 1 . ± .
038 for N s /N t = 6. From astudy of hysteresis curves for ( ǫ + p ) /T around T c , we note that ( ǫ + p ) /T in the (metastable)deconfined phase is sensitive to the spatial volume, while that in the confined phase is insensitive.The value of ( ǫ + p ) /T decreases as the volume increases in the (metastable) deconfined phase.Study with higher statistics is needed to precisely determine the latent heat in the thermodynamiclimit. 23 x8 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x8 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. A1 ∆( ǫ − p ) /T (blue), ∆( ǫ + p ) /T (green) and ∆ ǫ/T (red) calculated on 48 × × t/a . The straight lines are the results of t → t X method. We compared two orders of the double extrapolation( a, t ) → (0 ,
0) — method 1 (first t → a →
0) and method 2 (first a → t → t X method.
Acknowledgements
We thank other members of the WHOT-QCD Collaboration for discussions and comments.This work was in part supported by JSPS KAKENHI Grant Numbers JP20H01903, JP19H05146,JP19H05598, JP19K03819, JP18K03607, JP17K05442 and JP16H03982, and the Uchida EnergyScience Promotion Foundation. This research used computational resources of COMA, Oakforest-PACS, and Cygnus provided by the Interdisciplinary Computational Science Program of Center forComputational Sciences, University of Tsukuba, K and other computers of JHPCN through theHPCI System Research Projects (Project ID:hp17208, hp190028, hp190036, hp200089) and JHPCNprojects (jh190003, jh190063, jh200049), OCTOPUS at Cybermedia Center, Osaka University, ITOat Research Institute for Information Technology, Kyushu University, and Grand Chariot at Infor-mation Initiative Center, Hokkaido University, SR16000 and BG/Q by the Large Scale SimulationProgram of High Energy Accelerator Research Organization (KEK) (Nos. 14/15-23, 15/16-T06,15/16-T-07, 15/16-25, 16/17-05). The authors also thank the Yukawa Institute for TheoreticalPhysics at Kyoto University for the workshop YITP-W-19-09.24 x12 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x12 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. A2
The same as Fig. A1, but on the 48 ×
12 (top left), 64 ×
12 (top right), 72 × ×
12 (bottom right) lattices. x16 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T x16 NLO µ d ∆ε /T ∆(ε +p)/T ∆(ε -3p)/T Fig. A3
The same as Fig. A1, but on the 96 ×
16 (left) and 128 ×
16 (right) lattices.25 ∆ ( ε - ) / T tT c2 N s /N t =8 NLO µ d continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =8 NLO µ d continuum limitN t = 8N t = 12N t = 16 Fig. A4 ∆( ǫ − p ) /T (left) and ∆( ǫ + p ) /T (right) with the NLO matching coefficients,obtained on the N s /N t = 8 lattices. Three orange curves represent the result in the continuumlimit. (Thick orange curves represent the central values and thin curves represent the range of theirerrors.) A Energy-momentum tensor with NLO matching coefficients
In this appendix, we calculate the energy-momentum tensor using the NLO matchingcoefficients, following the original paper of the SF t X method [8]: T Rµν ( x ) = lim t → (cid:26) α U ( t ) U µν ( t, x ) + δ µν α E ( t ) [ E ( t, x ) − h E ( t, x ) i ] (cid:27) , (A1)where α U ( t ) = 1 c ( t ) = g (cid:20) β (4 π ) ¯ s g + O ( g ) (cid:21) (A2) α E ( t ) = 1 c ( t ) = (4 π ) β (cid:20) β (4 π ) ¯ s g + O ( g ) (cid:21) . (A3)For the NLO calculation in this appendix, following Ref. [8], we keep terms up to the next-to-leadingorder only in α U ( t ) and α E ( t ) and adopt the µ d scale. We then have¯ s = 722 + 12 γ E − ln 2 ≃ − . , (A4)¯ s = 2144 − β β = 27484 ≃ . . (A5)The details of the simulations are the same as the NNLO calculations given in Sec. 3.1.The results of ∆( ǫ − p ) /T , ∆( ǫ + p ) /T and ∆ ǫ/T are plotted as functions of flow time t/a in Figs. A1, A2, and A3. The blue, green, and red symbols are the results of ∆( ǫ − p ) /T ,∆( ǫ + p ) /T , and ∆ ǫ/T , respectively. Method 1: t → followed by a → t → t → Method 2: a → followed by t → t inphysical units, we replot the data shown in Figs. A1, A2, and A3 as a function of tT c . For the26 ∆ ( ε - ) / T tT c2 N s /N t =6 NLO µ d continuum limitN t = 8N t = 12N t = 16 0.6 0.8 1 1.2 1.4 1.6 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =6 NLO µ d continuum limitN t = 8N t = 12N t = 16 Fig. A5
The same as Fig. A4, but on the N s /N t = 6 lattices.aspect ratio N s /N t = 8 and 6, we obtain Figs. A4 and A5. The red, green and blue symbols arethe results of N t = 8, 12 and 16. These figures correspond to Figs. 9 and 10 with the NNLOmatching coefficients. We find that the difference between the results with different lattice spacingbecomes smaller as the flow time t increases. Repeating the analyses of Sec. 4.3, we obtain thecontinuum limit shown by the orange curves in Figs. A4 and A5. These results of ∆( ǫ + p ) /T and∆( ǫ − p ) /T in the continuum limit are summarized in Fig. 11 as functions of t , and the resultsof the t → B Choice of lattice field strength operator
In this appendix, we test the influence of lattice operators for the field strength G µν ( t, x ) todefine the operators U µν ( t, x ) and E ( t, x ) in Eqs. (3) and (4). Two major choices for G µν ( t, x ) arethe plaquette and clover-shaped lattice operators.In the left and right panels of Fig. B1, we compare the results of ∆( ǫ − p ) /T and ∆( ǫ + p ) /T in the continuum limit obtained by the clover-shaped and plaquette operators adoptingmethod 2 discussed in Sec. 2.4. The upper and lower panels are the results for N s /N t = 8 and6, respectively. The blue and red symbols are for results with the clover-shaped and plaquetteoperators, respectively. The two results agree well with each other at sufficiently large t . Thedisagreement at tT c < ∼ .
01 is caused by the lattice discretization error, and the data there shouldbe removed in the t → t X method. Because we have a widerrange of t for the linear t → U µν ( t, x ) and E ( t, x ) in this study. With the clover-shaped operator, we may use datadown to tT c ∼ . References [1] F. Karsch, “SU( N ) Gauge Theory Couplings on Asymmetric Lattices,” Nucl. Phys. B , 285 (1982).[2] J. Engels, F. Karsch, H. Satz and I. Montvay, “Gauge Field Thermodynamics for the SU(2) Yang-Mills System,”Nucl. Phys. B , 545 (1982).[3] G. Burgers, F. Karsch, A. Nakamura and I.O. Stamatescu, “QCD on anisotropic lattices,” Nucl. Phys. B ,587 (1988).[4] S. Ejiri, Y. Iwasaki and K. Kanaya, “Nonperturbative determination of anisotropy coefficients in lattice gaugetheories,” Phys. Rev. D , 094505 (1998). ∆ ( ε - ) / T tT c2 N s /N t =8 NNLO µ plaqclov 0.8 0.9 1 1.1 1.2 1.3 1.4 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =8 NNLO µ plaqclov 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε - ) / T tT c2 N s /N t =6 NNLO µ plaqclov 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 0 0.005 0.01 0.015 0.02 0.025 ∆ ( ε + p ) / T tT c2 N s /N t =6 NNLO µ plaqclov Fig. B1 ∆( ǫ − p ) /T using clover-shaped operator (red) and plaquette (blue) in the continuumlimit for N s /N t = 8 (upper) and 6 (lower) by the NNLO calculation with µ = µ , respectively. [5] J. Engels, F. Karsch and T. Scheideler, “Determination of anisotropy coefficients for SU(3) gauge actions fromthe integral and matching methods,” Nucl. Phys. B , 303 (2000).[6] S. Ejiri, “Remarks on the multiparameter reweighting method for the study of lattice QCD at nonzerotemperature and density,” Phys. Rev. D , 094506 (2004).[7] M. Shirogane, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa [WHOT-QCD Collaboration] “Latent heat at thefirst order phase transition point of the SU(3) gauge theory,” Phys. Rev. D , 014506 (2016).[8] H. Suzuki, “Energy-momentum tensor from the Yang-Mills gradient flow,” PTEP , 083B03 (2013)Erratum: [PTEP , 079201 (2015)].[9] H. Makino and H. Suzuki, “Lattice energy-momentum tensor from the Yang-Mills gradient flow—inclusion offermion fields,” PTEP , 063B02 (2014) Erratum: [PTEP , 079202 (2015)][10] R. Narayanan and H. Neuberger, “Infinite N phase transitions in continuum Wilson loop operators,” JHEP , 064 (2006)[11] M. L¨uscher, “Trivializing maps, the Wilson flow and the HMC algorithm,” Commun. Math. Phys. , 899(2010)[12] M. L¨uscher, “Properties and uses of the Wilson flow in lattice QCD,” JHEP , 071 (2010) Erratum: [ ,092 (2014)][13] M. L¨uscher and P. Weisz, “Perturbative analysis of the gradient flow in nonabelian gauge theories,” JHEP , 051 (2011)[14] M. L¨uscher, “Chiral symmetry and the Yang-Mills gradient flow,” JHEP , 123 (2013)[15] M. Asakawa, T. Hatsuda, E. Itou, M. Kitazawa, H. Suzuki [FlowQCD Collaboration], “Thermodynamics ofSU(3) gauge theory from gradient flow on the lattice,” Phys. Rev. D , 011501 (2014) Erratum: , 059902(2015)[16] M. Kitazawa, T. Iritani, M. Asakawa, T. Hatsuda, and H. Suzuki , “Equation of State for SU(3) Gauge Theoryvia the Energy-Momentum Tensor under Gradient Flow,” Phys. Rev. D , 114512 (2016).[17] T. Hirakida, E. Itou and H. Kouno, “Themodynamics for pure SU(2) gauge theory using gradient flow,” PTEP , 033B01 (2019).[18] T. Iritani, M. Kitazawa, H. Suzuki, and H. Takaura, “Thermodynamics in quenched QCD: energy-momentumtensor with two-loop order coefficients in the gradient-flow formalism,” Prog. Theor. Exp. Phys. , 023B02(2019).[19] Y. Taniguchi, S. Ejiri, R. Iwami, K. Kanaya, M. Kitazawa, H. Suzuki, T. Umeda, and N. Wakabayashi [WHOT-QCD Collaboration], “Exploring N f = 2 + 1 QCD thermodynamics from gradient flow,” Phys. Rev. D ,014509 (2017); Erratum: , 059904(E) (2019).[20] Y. Taniguchi, K. Kanaya, H. Suzuki, and T. Umeda [WHOT-QCD Collaboration], “Topological susceptibilityin finite temperature (2 + 1)-flavor QCD using gradient flow,” Phys. Rev. D , 054502 (2017)[21] Y. Taniguchi, S. Ejiri, K. Kanaya, M. Kitazawa, H. Suzuki, T. Umeda [WHOT-QCDCollaboration], “ N f =2+1 QCD thermodynamics with gradient flow using two-loop matching coefficients,” Phys. Rev. D , 014510(2020); Erratum: , 059903 (2020).[22] K. Kanaya, A. Baba, A. Suzuki, S. Ejiri, M. Kitazawa, H. Suzuki, Y. Taniguchi, and T. Umeda, “Study of 2 + 1flavor finite-temperature QCD using improved Wilson quarks at the physical point with the gradient flow,”Proc. Sci., LATTICE 2019, 088 (2020).[23] M. Kitazawa, T. Iritani, M. Asakawa and T. Hatsuda, “Correlations of the energy-momentum tensor via gradientflow in SU(3) Yang-Mills theory at finite temperature,” Phys. Rev. D , no. 11, 111502 (2017),[24] R. Yanagihara, T. Iritani, M. Kitazawa, M. Asakawa and T. Hatsuda, “Distribution of Stress Tensor aroundStatic Quark–Anti-Quark from Yang-Mills Gradient Flow,” Phys. Lett. B , 210 (2019).[25] M. Kitazawa, S. Mogliacci, I. Kolb´e and W. Horowitz, “Anisotropic pressure induced by finite-size effects inSU(3) Yang-Mills theory,” Phys. Rev. D , 094507 (2019).[26] R. Yanagihara, M. Kitazawa, M. Asakawa and T. Hatsuda, “Distribution of Energy-Momentum Tensor arounda Static Quark in the Deconfined Phase of SU(3) Yang-Mills Theory,” arXiv:2010.13465 [hep-lat].[27] Y. Taniguchi, S. Ejiri, K. Kanaya, M. Kitazawa, A. Suzuki, H. Suzuki, and T. Umeda, “Energy-momentumtensor correlation function in N f = 2 + 1 full QCD at finite temperature,” EPJ Web Conf. , 07013 (2018);[28] Y. Taniguchi, A. Baba, A. Suzuki, S. Ejiri, K. Kanaya, M. Kitazawa, T. Shimojo, H. Suzuki, and T. Umeda,“Study of energy-momentum tensor correlation function in N f = 2 + 1 full QCD for QGP viscosities,” Proc.Sci. , LATTICE 2018, 166 (2019).[29] R.V. Harlander, Y. Kluth, and F. Lange, “The two-loop energy-momentum tensor within the gradient-flowformalism,” Eur. Phys. J. C 78:944 (2018).[30] H. Suzuki, “Background field method in the gradient flow,” PTEP , no. 10, 103B03 (2015).[31] J. Artz, R. V. Harlander, F. Lange, T. Neumann and M. Prausa, “Results and techniques for higher ordercalculations within the gradient-flow formalism,” J. High Energy Phys. , 121 (2019).[32] A.M. Ferrenberg and R. H. Swendsen, “Optimized Monte Carlo analysis,” Phys. Rev. Lett. , 1195 (1989).[33] R. Iwami, S. Ejiri, K. Kanaya, Y. Nakagawa, D. Yamamoto, T. Umeda [WHOT-QCD Collaboration],“Multipoint reweighting method and its applications to lattice QCD,” Phys. Rev. D , 094507 (2015).[34] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, and C. H. Wong, “The lattice gradient flow at tree-leveland its improvement,” J. High Energy Phys. , 018 (2014)., 018 (2014).