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

We study latent heat and pressure gap between the hot and cold phases at the first order deconfining phase transition temperature of the SU(3) Yang-Mills theory. Performing simulations on lattices with various spatial volumes and lattice spacings, we calculate the gaps of the energy density and pressure using the small flow time expansion (SFtX) method. We find that the latent heat $\Delta \epsilon$ in the continuum limit is $\Delta \epsilon /T^4 = 1.117 \pm 0.040$ for the aspect ratio $N_s/N_t=8$ and $1.349 \pm 0.038$ for $N_s/N_t=6$ at the transition temperature $T=T_c$. We also confirm that the pressure gap is consistent with zero, as expected from the dynamical balance of two phases at $T_c$. From hysteresis curves of the energy density near $T_c$, we show that the energy density in the (metastable) deconfined phase is sensitive to the spatial volume, while that in the confined phase is insensitive. Furthermore, we examine the effect of alternative procedures in the SFtX method - the order of the continuum and the vanishing flow time extrapolations, and also the renormalization scale and higher order corrections in the matching coefficients. We confirm that the final results are all well consistent with each other for these alternatives.

We study latent heat and pressure gap between the hot and cold phases at the first order deconfining phase transition temperature of the SU(3) Yang-Mills theory. Performing simulations on lattices with various spatial volumes and lattice spacings, we calculate the gaps of the energy density and pressure using the small flow time expansion (SFtX) method. We find that the latent heat ∆ǫ in the continuum limit is ∆ǫ/T 4 = 1.117 ± 0.040 for the aspect ratio N s /N t = 8 and 1.349 ± 0.038 for N s /N t = 6 at the transition temperature T = T c . We also confirm that the pressure gap is consistent with zero, as expected from the dynamical balance of two phases at T c . From hysteresis curves of the energy density near T c , we show that the energy density in the (metastable) deconfined phase is sensitive to the spatial volume, while that in the confined phase is insensitive. Furthermore, we examine the effect of alternative procedures in the SFtX method -the order of the continuum and the vanishing flow time extrapolations, and also the renormalization scale and higher order corrections in the matching coefficients. We confirm that the final results are all well consistent with each other for these alternatives.

Introduction
First order phase transition appears in various important thermodynamic systems including the high density QCD and many-flavor QCD, and is associated with various particular phenomena such as the phase coexistence and the latent heat. It is worth developing numerical techniques to investigate first order phase transitions. The SU(3) Yang-Mills theory, i.e., the quenched approximation of QCD, at finite temperature has a first order deconfining phase transition, and provides us with a good testing ground for that.
At a first order phase transition point, metastable states corresponding to two phases coexist at the same time. To keep a dynamical balance among them, their pressures must be the same. Therefore, one can check the validity of computational methods by measuring the pressure in each metastable state. In fact, the problem of the non-vanishing pressure gap has played an important role in early development of lattice methods to study thermodynamic quantities [1][2][3][4][5][6]. On the other hand, the energy density takes different values in each phase. The difference of them is the latent heat, which is one of the fundamental observables characterizing the first order phase transition.
In Refs. [4,7], we studied the first order transition of the SU(3) Yang-Mills theory using the derivative method [1]. We have numerically confirmed that the pressure gap vanishes when we adopt 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 adopting a new technique to calculate thermodynamic observables: the small flow time expansion (SFtX) method based on the gradient flow [8,9]. The gradient flow is an evolution of the fields in terms of a fictitious "time" t [10][11][12][13][14]. Fields at positive flow-time t > 0 can be viewed as smeared fields averaged over a mean-square physical radius of √ 8t in four dimensions, and the operators constructed by flowed fields at t > 0 are shown to be free from the ultraviolet divergences and short-distance singularities. Making use of this strict finiteness, the SFtX method provides us with a general method to correctly calculate any renormalized observables on the lattice. The SFtX method has been successfully applied to the evaluation of thermodynamic quantities in the Yang-Mills gauge theory [15][16][17][18] and in QCD with 2 + 1 flavors of dynamical quarks [19][20][21][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 with the energy-momentum tensor [23][24][25][26][27].
We apply the SFtX method to calculate the energy density and pressure, separately in hot and cold phases at the first order deconfining phase transition temperature of the SU(3) Yang-Mills theory. We perform simulations on lattices with several lattice spacings and spatial volumes to carry out the continuum extrapolation and also to study the finite volume effect. With confirming that the pressure gap is consistent with zero, we evaluate the latent heat.
We also test several alternative procedures in the SFtX method: To obtain a physical result by the SFtX method, a double extrapolation to the continuum and vanishing flow time limits, a → 0 and t → 0, is required. When a proper fitting range avoiding small t singularities at a > 0 is chosen, the final results should be insensitive to the order of these limits. We confirm this for the latent heat and pressure gap in the SU(3) Yang-Mills theory. The second point to be addressed is the quality of the matching coefficients which relate the physical observable and flowed operators in the SFtX 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 matching coefficients. 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 the final results are affected by the choice of them. We show that the final results for the latent heat and pressure gap are insensitive to the choice of these alternatives. We also show that the NNLO matching coefficients help reducing lattice artifacts in the calculation of the latent heat.
In the next section, we introduce the SFtX method and several details of the method to be discussed in this paper. Then, our simulations at the first order transition of SU(3) Yang-Mills theory 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 SFtX method. We then study the hysteresis of the energy density near the phase transition temperature and its spatial volume dependence in Sec. 5. Our conclusions are given in Sec. 6. Appendix A is devoted to show the results of the latent heat and pressure gap with the NLO matching coefficients, with which we study the effect due to the truncation of the perturbative series for the matching coefficients. In Appendix B, we discuss on the choice of the lattice operator for the field strength.

Small flow-time expansion method
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 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 √ 8t 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 non-perturbatively and thus can be calculated directly on the lattice.
The small flow-time expansion (SFtX) method provides us with a general method to evaluate any renormalized observables in terms of flowed operators at small and positive flow time t, which are 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 the coefficients ("matching coefficients") relating the renormalized observables to the flowed operators at small flow time t [13]. Contamination of O(t) terms can be removed by a vanishing flow time extrapolation t → 0.

Energy-momentum tensor
To calculate the energy-momentum tensor by the SFtX method, we consider the following flowed operators at flow time t which are gauge-invariant and local, With these dimension-four operators, the correctly renormalized energy-momentum tensor T R µν is given by [8] where E(t, x) 0 is the zero temperature subtraction. Here, the matching coefficients c 1 (t) and c 2 (t) are calculated by the perturbation theory [13] as 1 where g = g(µ) is the renormalized gauge coupling in the MS scheme at the renormalization scale µ.
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 1 (t) contains k (0) 1 and k (1) 1 , and that for c 2 (t) contains k (1) 2 and k (2) 2 , while the next to next to leading order (NNLO) expression for c 1 (t) contains terms up to k (2) 1 , and that for c 2 (t) contains terms up to k 2 . We mainly adopt NNLO matching coefficients in this study. To study the effect due to the truncation of perturbative series, we also perform calculation using NLO matching coefficients in Appendix A.
The coefficients k (1) i in the one-loop level are calculated in Refs. [8,9,29], and k (2) i in the twoloop level in Refs. [28]. (See also Ref. [30].) As pointed out in Ref. [8], in pure gauge Yang-Mills theories, k can be deduced by ℓ-loop coefficients using the trace anomaly. A concrete form of k 2 is given in Ref. [18]. Collecting these results for the case of pure gauge Yang-Mills theories, we where are the first three coefficients of the beta function in Yang-Mills theories, and was introduced in Ref. [28] with γ E the Euler-Mascheroni constant. The factor C A is the quadratic Casimir for the adjoint representation defined by and C A = N c for the gauge group SU(N c ).

Latent heat and pressure gap
The energy density and the pressure are obtained from the diagonal elements of the energymomentum tensor, 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 the critical point, we calculate the energy density and the pressure in each phase just at the transition temperature, to estimate the latent heat and the pressure gap defined by where ǫ (hot/cold) and p (hot/cold) indicate ǫ and p in the hot and cold phase, respectively.
In this paper, we calculate these quantities in the following conventional combinations of the trace anomaly and the entropy density, When ∆p = 0 as expected, ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 should coincide with each other. We note that ∆(ǫ − 3p)/T 4 is computed by the operator E(t, x) with the matching coefficient c 2 (t), while ∆(ǫ + p)/T 4 is computed by the operator U µν (t, x) with c 1 (t).

Renormalization scale
The matching coefficients c i (t) given in the previous subsection are written in terms of the renormalization scale µ. However, because the matching coefficients relate the energy-momentum tensor and flowed operators, both physical and thus both finite and independent of the renormalization scale µ, the matching coefficients themselves should also be finite and independent of µ -the explicit µ-dependence from L(µ, t) should be cancelled by that of the running coupling g, in principle. In practice, however, because the perturbative series for the matching coefficients are truncated at a finite order, the choice of µ may affect the magnitude of systematic errors due to the neglected higher order terms.
In early studies with the SFtX method [15,16,19,23], the µ d -scale has been adopted because µ d is a natural scale of local flowed operators which have physical smearing extent of √ 8t. On the other hand, in Ref. [28], it is argued that would be a good choice of µ because it keeps the magnitude of two-loop contributions small. The µ 0 -scale was tested in quenched QCD [18] and in 2 + 1 flavor QCD [21,22], and it was reported that the µ 0 -scale may improve the signal when the lattice is not quite fine [21,22]. In this paper, we examine the both choices µ = µ 0 (t) and µ d (t). From the difference in the results between these two choices, we learn the magnitude of the effect due to truncation of higher order terms in the perturbation theory.

Continuum and vanishing flow time extrapolations
To obtain a physical result with the SFtX method, the continuum extrapolation a → 0 and the vanishing flow time extrapolation t → 0 are both mandatory to remove contamination of unwanted higher dimension operators in Eq. (5) as well as errors from the lattice regularization. In Refs. [16-18, 23, 24, 26], this double extrapolation is carried out by first taking a → 0 at each flow time t, and then take t → 0 using the continuum-extrapolated results. In Refs. [19][20][21]25], t → 0 is taken first at fixed finite a, reserving the continuum extrapolation for a later stage. In the following, we call the way to first take t → 0 at fixed finite a and then take a → 0 as "method 1," and the way to first take a → 0 at fixed finite t and then take t → 0 as "method 2." On finite lattices, additional unwanted operators may contaminate in 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 lattice operator in the curly bracket of Eq. (5) is expanded as where T R µν is the physical energy-momentum tensor, S µν and S ′ µν are contaminations of dimensionsix operators with the same quantum number, and A µν , C µν , and D µν are those from dimension-four operators. In higher orders of a, we also have terms proportional to (a 2 /t) 2 etc. [31]. In Eq. (20), the term proportional to a 2 /t is singular at t = 0 and thus is problematic when we want to take t → 0 at finite a in the method 1. When we take the a → 0 limit first, the term proportional to a 2 /t is removed in principle. However, in practice, the a 2 /t term is problematic also in the method 2 when we want to take a → 0 at small t where the a 2 /t term dominates the data -we cannot perform the a → 0 extrapolation reliably at small t. Therefore, in the both methods, we need to find a range of t at each a where the a 2 /t and more singular terms are negligible. Let us call the range of t in which the a 2 /t and more singular terms as well as the O(t 2 ) terms look negligible as "linear window" [19,21].
On the other hand, when we could find such range of t and perform t → 0 and a → 0 extrapolations using data in this range only, because the flowed operator on the lattice will be well approximated as in this range, the final results of the method 1 and the method 2 for T R µν should be identical. The consistency of the both methods is thus a good test of the SFtX method. In the following, we adopt both method 1 and method 2 using data in linear windows and perform extrapolations linear in t and a 2 . When the both methods are consistent with each other, we get a strong support in our identification of the linear windows.
Another issue in the study of thermodynamic properties is the possible dependence on the physical volume of the system. We have to keep the physical volume fixed in the a → 0 and t → 0 extrapolations to identify finite volume effects. In our study, T is adjusted to the transition temperature T c whose value is a physical constant. On a lattice with the size N 3 s × N t , the lattice spacing defined as a = 1/(N t T c ) is thus changed by changing N t , and the spatial volume in physical 3 , is fixed when the aspect ratio N s /N t is fixed. In this paper, we take the double extrapolation for the cases of N s /N t = 6 and 8 to study the finite volume effect.

Simulation parameters
We perform simulations of the SU(3) Yang-Mills theory with the standard Wilson action at several inverse gauge couplings β = 6/g 2 around the deconfining transition point β c . We have studied the lattices with temporal lattice size of N t = 8, 12 and 16 with several different spatial lattice sizes N s . Some of the configurations are taken from Ref. [7]. Although we have studied N t = 6 lattices too, we do not use them in this study because the linear window turned out to be not clear enough for N t = 6. Our simulation parameters are summarized in Table 1. The configurations are generated by a pseudo heat bath algorithm followed by over-relaxation updates. Each measurement is separated by N sep heat-bath sweeps. Data are taken at 3 to 6 β values for each (N s , N t ), and are combined using the multipoint reweighting method [32,33]. The statistical errors are estimated by the jack-knife 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 3 s ( Ω 2 − Ω 2 ). 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 gauge action [12]. Numerical solution of the flow equation is obtained by the third order Runge-Kutta method. For the field strength G µν (t, x) in the observables, several candidates are available for the Table 1 Simulation parameters: the spatial lattice size N s , the temporal lattice size N t , coupling parameter β = 6/g 2 0 and the critical coupling β c . The measurements have been performed N conf times at intervals of N sep sweeps. The data for 48 3 × 8, 64 3 × 8 and 64 3 × 12 are taken from our previous study [7]. lattice operator. From a test of the plaquette and the clover-type representations for G µν (t, x) given in Appendix B, we adopt the clover-type representation.

Phase separation around the first order transition point
To evaluate the latent heat and the pressure gap, we need to separate the configurations around the first order transition point into the hot and cold phases. In the left panel of Fig. 1, we show the histogram of the absolute value of the Polyakov loop |Ω| obtained on the 96 3 × 12 lattice, where Ω is averaged over the spatial coordinates. The coupling parameter β is adjusted to the transition point β c using the multipoint reweighting method. The two peaks in Fig. 1 correspond to the hot and cold phases. As discussed in Ref. [7], though the two peaks of the Polyakov loop histogram are well separated, the peaks of the plaquette histogram are overlapping. It is thus meaningful to use the 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 =2.0 , which is constructed with the flowed gauge field at t/a 2 = 2.0. We find that, though the scale of the horizontal axis is different, its shape is quite similar to that shown in the left panel. This suggests a strong correlation between the original Polyakov loop Ω and the flowed Polyakov loopΩ. In Fig. 2, we show the two-dimensional distribution of (|Ω|, |Ω| t/a 2 =2.0 ) measured on the 96 3 × 12 lattice. Darker tone means higher probability. We find that Ω andΩ have a strong linear correlation with each other. These results show that the gradient flow does not affect the separation, and thus there is no merit to use the flowed 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 (101 configurations, except for the 128 3 × 16 lattice for which 45 configurations) around the current configuration. We then classify the configurations into the hot, cold, and mixed phases by the value of this time-smeared Polyakov loop. With this classification, we observe many flip-flops between the hot and cold phases during our Monte-Carlo steps, while mixed configurations in which the two phases coexist are found to be rare on our lattices. We discard the configurations in the mixed phase in the following.
After the phase separation, we carry out the gradient flow on each of the configurations to measure flowed operators given in Eqs. (3) and (4). We then combine their expectation values in each phase by the multipoint reweighting method to obtain their values just at the transition point β c .

Latent heat and pressure gap on finite lattices
In Figs. 3-5, we show the results of ∆(ǫ − 3p)/T 4 , ∆(ǫ + p)/T 4 and ∆ǫ/T 4 as functions of t/a 2 , using the NNLO matching coefficients. The blue, green and red symbols are the results of ∆(ǫ − 3p)/T 4 , ∆(ǫ + p)/T 4 and ∆ǫ/T 4 , respectively. The plots in the left panels are the results adopting the renormalization scale µ = µ 0 and those in the right panels adopting µ = µ d . The rapid decrease of these observables at small t signals the appearance of the O(a 2 /t) lattice discretization effect discussed in Sec. 2.4. Similar behavior at similar values of t has been reported for (ǫ − 3p)/T 4 and (ǫ + p)/T 4 themselves in previous studies of SU(3) Yang-Mills theory [15,16,18]. We thus disregard the data in this region in the t → 0 and a → 0 extrapolations. We note that the results behave almost linearly for t/a 2 1.0, while the O(t) error sometimes starts to contaminate at large t/a 2 . We also note that, though the dependence on the renormalization scale µ is small, there is a general tendency that the µ 0 -scale (left panels) leads to a slightly smaller slope towards the t → 0 limit than the µ d -scale (right panels).
It should be noted that, within the statistical errors, the results of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 are consistent with each other in a wide range of t. This means that the pressure gap between two phases ∆p vanishes, which must be satisfied in the final results after the double extrapolation of t → 0 and a → 0. This suggests that, for these observables, the errors from the lattice discretization and 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 procedure of Ref. [8]. The results are summarized in Appendix A. We find that the NNLO matching coefficients drastically 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 discussed in the next section, we also confirm that the results of ∆(ǫ − 3p)/T 4 , ∆(ǫ + p)/T 4 and ∆ǫ/T 4 , extrapolated to the t → 0 limit are consistent between the NNLO and NLO analyses. The same as Fig. 3, but on the 48 3 × 12, 64 3 × 12, 72 3 × 12, 96 3 × 12 lattices (from top to bottom).  Let us first analyze the data adopting the method 1, i.e., we first take the small flow time limit t → 0 on each lattice, and then take the continuum limit a → 0. We note that the data shown in Figs. 3-5 are well linear in a wide range of t within the statistical error. In this study, we choose the fit range as follows: First, we require the fit range to satisfy tT 2 c ≤ 0.025. This value of the upper bound is determined by consulting the linearity of the data and also from the consistency with the fit range adopted in the method 2 to be discussed in Sec. 4.3. We confirm that this upper bound satisfies the requirement that the smearing by the gradient flow is not overlapping around the periodic lattice (t/a 2 < t 1/2 defined in Ref. [19]). We also require t/a 2 ≥ √ 2 so that the smearing radius by the gradient flow well covers the nearest neighbor lattice sites. For the case of N t = 8, however, because this makes the available range too narrow, we relaxed it to the minimum requirement of t/a 2 ≥ 1. The fit ranges we adopt as well as the results of linear t → 0 extrapolations are given in Table 2. Corresponding linear fit lines are shown in Figs. 3-5. In Table 2, we also show the results with the NLO matching coefficients discussed in Appendix A. The errors in Table 2 are statistical only. We have also tested different fit ranges, but the variation of the results due to the fit range turned out to be within the statistical errors given in Table 2.
From Table 2, we note that results of different renormalization scales as well as with NNLO and NLO matching coefficients are all well consistent with each other. We also find that the differences between ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 are less than about one sigma for all lattices, in accordance with the expectation that ∆p should vanish.
In Fig. 6, we plot the results of ∆(ǫ − 3p)/T 4 (red) and ∆(ǫ + p)/T 4 (blue) at t = 0 measured on N 3 s × 12 lattices with N s = 48, 64, 72 and 96, as functions of the inverse spatial volume in lattice units 1/N 3 s . The lattice spacing a = 1/(N t T c ) is the same for all data shown in Fig. 6. We find that ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 decreases as the spatial volume increases. Therefore, we have to take care the finite volume effect.
We now perform the continuum extrapolation a → 0 using the results shown in Table 2. To study the finite volume effect, we perform this on lattices with fixed spatial volume V . As discussed in Sec. 2.4, this is achieved by fixing the aspect ratio N s /N t in our study. In Fig. 7 Fig. 8 The same as Fig. 7, but for the aspect ratio N s /N t = 6. Table 2 The results of t → 0 extrapolation of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 on each lattice. Errors are statistical only. Systematic errors due to the choice of the fit range are smaller than the statistical errors. ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 at t = 0 as functions of 1/N 2 t = (T c a) 2 for N s /N t = 8. The left and middle panels are the results adopting µ = µ 0 and µ d , respectively. We also show the results with NLO matching coefficients adopting µ = µ d , discussed in Appendix A. Corresponding results for N s /N t = 6 are shown in Fig. 8.
We make the a → 0 extrapolation by fitting these results by a linear function of (T c a) 2 = 1/N 2 t . We find that the lattice spacing dependence is small in these observables. However, we also note that the slight slope toward the continuum limit changes its sign between N s /N t = 8 and 6. This may be suggesting that, when we consider various systematic errors, the lattice spacing dependence is approximately absent within errors.  Table 3. We find that the results for ∆ǫ/T 4 , ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 , obtained using the NNLO or NLO matching coefficients and with the µ 0 or µ d scale, are all well consistent with each other at each aspect ratio. From the results of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 with the NNLO matching coefficients and the µ 0 -scale, we obtain the pressure gap to be ∆p/T 4 = 0.015(17) (N s /N t = 8), ∆p/T 4 = 0.015 (14) which are only about 1% of the latent heat and are consistent with zero.  Next, we adopt the method 2 to take the double limit (t, a) → (0, 0), i.e., we first take the a → 0 limit at each flow time t in physical units, and then take the t → 0 limit. We do them with fixing the aspect ratio N s /N t to take care the finite size effect. In Fig. 9, we plot the results of ∆(ǫ − 3p)/T 4 (left) and ∆(ǫ + p)/T 4 (right) as functions of the flow time in a physical unit, tT 2 c = t/(aN t ) 2 , 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 µ = µ 0 and µ d are given in the top and bottom panels of each figure, respectively. The red, green and blue symbols are 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 small as the flow time increases.
We take the a → 0 limit by fitting the data with a linear function of 1/N 2 t at each flow time tT 2 c . The orange curves in Figs. 9 and 10 are the results in the continuum limit. The center curve represents the central value, and the upper and lower curves represent the statistical error. For information, we also show the results obtained outside the linear window in these figures. In Fig. 11, we show the results of ∆(ǫ − 3p)/T 4 (circle) and ∆(ǫ + p)/T 4 (square) in the continuum limit as  Fig. 10 The same as Fig. 9, but for N s /N t = 6. Table 3 ∆ǫ/T 4 , ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 in the continuum limit with the NNLO and NLO matching coefficients adopting the renormalization scale µ 0 or µ d . The left and right three columns are results of the method 1 (t → 0 followed by a → 0) and the method 2 (a → 0 followed by t → 0), respectively. Errors are statistical only. Systematic errors due to variation of the fit ranges are smaller than the statistical errors.

method 1 method 2
Ns Nt functions of tT 2 c for N s /N t = 8 (filled symbols) and N s /N t = 6 (open symbols). The results adopting µ = µ 0 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 results  using data with NLO matching coefficients with the µ d -scale given in Appendix A. Comparing upper and bottom panels, we find that, by adopting the NNLO matching coefficients, the slopes in t are much reduced and the pressure gap ∆p/T 4 at finite t is much suppressed.
We then perform t → 0 extrapolation using the continuum extrapolated results. We avoid the data in the small t region where contamination of the O(a 2 /t) terms is suspected. To estimate a systematic error due to the t → 0 extrapolation, we test the following five fit ranges: Range 1: 0.010 < tT c < 0.020, Range 2: 0.010 < tT c < 0.025, Range 3: 0.005 < tT c < 0.020, Range 4: 0.015 < tT c < 0.025, and Range 5: 0.005 < tT c < 0.025. In Fig. 12, we show the results of these liner extrapolations in tT 2 c using data with the NNLO matching coefficients and the µ 0 -scale. The triangle, square and circle symbols are for ∆ǫ/T 4 , ∆(ǫ + p)/T 4 and ∆(ǫ − 3p) at t = 0, and the filled and open symbols are the results of N s /N t = 8 and 6, respectively. We find that the results adopting different fit ranges are well consistent with each other. In Fig. 12, we also show the results of the method 1 obtained in the previous subsection at the left end of the plot, which are also consistent with the results of the method 2 within statistical errors. Our final results of the 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 µ 0 -scale using the fit range 2. Errors are statistical only. As shown in Fig. 12, systematic errors due to the choice of the fit range are smaller than the statistical errors.

Results of the SFtX method for the latent heat
Finally, we summarize our results of the SFtX method for the latent heat. From Table 3 and Fig. 12, we see that the results of the method 1 and method 2 agree well with each other. We take our central values from the method 2 with the NNLO matching coefficients and the µ 0 -scale. We obtain for N s /N t = 8 ∆(ǫ + p)/T 4 = 1.135(43), ∆ǫ/T 4 = 1.117(40), ∆(ǫ + p)/T 4 = 1.372(43), ∆ǫ/T 4 = 1.349(38).
Systematic errors estimated from the differences between the method 1 and method 2 as well as among different fit ranges are smaller than the statistical errors quoted in these equations. Performing heuristic extrapolations of these two data points to the thermodynamic limit, we find ∆ǫ/T 4 ∼ 0.95 ± 0.07 by a 1/V linear fit and ∆ǫ/T 4 ∼ 1.23 ± 0.03 by a constant fit.

Comparison to the derivative method
Finally, we compare our results of the SFtX method with those obtained by the derivative method [7]. The computational cost for the SFtX method is much smaller than that for the derivative method, since the nonperturbative calculation of the Karsch coefficients needed in the derivative method requires large systematic study with high statistics. We note that the statistical errors are reduced with the SFtX method thanks the smearing nature of the gradient flow and also the avoidance of the nonperturbative Karsch coefficients.
In Ref. [7], using the derivative method, a rather low value of ∆ǫ/T 4 = 0.75 ± 0.17 was obtained for the latent heat in the continuum limit using data at N t = 6, 8 and 12. In Ref. [7], because the spatial volume dependence in the latent heat was found to be small in comparison with the statistic errors, which are much larger than the errors in the current study mainly due to the error from the nonperturbative Karsch coefficients, spatial volume dependence was assumed to be absent.
To avoid uncertainties due to the continuum and infinite spatial volume extrapolations, we directly compare the results of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 at finite lattice spacings and spatial volumes in Fig. 13. The horizontal axis is N s /N t = 3 √ V T c . The filled symbols are the results of the SFtX method in the t → 0 limit with the NNLO matching coefficients and the µ 0 -scale, obtained on 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 = 6 and 8 by the method 2. The open symbols are the results of the derivative method obtained on the lattices with N t = 6 (black pentagon), 8 (red diamond) and 12 (green circle), respectively [7]. The large statistical errors for ∆(ǫ + p)/T 4 with the derivative method are caused by the large error of nonperturbative Karsch coefficients which increase rapidly with N t . From these plots, we find that the SFtX 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 SFtX method is small for N t = 8, 12, 16. The lattice spacing dependence in the results by the derivative method is also small for N t = 8 and 12. The left and right panels of Fig. 14 are the results of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 as functions of 1/N 2 t = a 2 T 2 c , i.e., the lattice spacing squared normalized by T c . The filled and open symbols are the results by the SFtX 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 data at N t = 6 with the derivative method, shown at the right end of these plots, show deviation from general tendency of data on finer lattices, suggesting a large discretization error at N t = 6. As suggested from Fig. 9 of Ref. [7], by removing the data at N t = 6, the derivative method may also lead to a larger latent heat in the continuum, though a reliable extrapolation is not possible with only two data points.

Hysteresis of the energy density near 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 the origin of the spatial volume dependence in our latent heat, we study (ǫ + p)/T 4 separately in the hot and cold phases at temperatures near T c using the multipoint reweighting method.  The results obtained on 48 3 × 8 and 64 3 × 8 lattices are plotted by the cross and circle symbols in Fig. 15. The flow time is t/a 2 = 1.4 -as shown in Sec. 4.2, the result is about the same as that in the t → 0 limit at this value of t. The NNLO matching coefficients with the µ 0 -scale are used. The red symbols are for the results obtained in the hot phase and the blue symbols for the cold phase, while the green symbols show the results without the phase separation. The horizontal axis is the temperature T = 1/(N t a) normalized by the critical temperature T c , where the relation between the lattice spacing a and β is determined by the critical point β c as function of N t [7]. The results for the hot phase below T c and for the cold phase above T c are those obtained in corresponding metastable 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) cold phase. We thus conclude that the spatial volume dependence of the latent heat is due to that in the contribution of the hot phase. From this figure, we also note that latent heat is sensitive to the value of the critical point β c . A careful determination is required for β c .
In Fig. 16, we show the corresponding plot by the derivative method, obtained on the 48 3 × 8 (cross) and 64 3 × 8 (circle) lattices [4,7]. In this calculation, because our nonperturbative determination of the Karsch coefficients is possible only at β c [4,7], we use the Karsch coefficients obtained at β c also at other values of β around β c . We see similar spatial volume dependence of the results in the (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 separation are not large in the reweighting calculation when T is not sufficiently close to T c . Further study with higher statistics is needed to precisely determine the volume dependence of the hot phase, and thus to carry out a reliable extrapolation of the latent heat to the thermodynamic limit.

Conclusions
To study the latent heat and the pressure gap at the first order deconfining phase transition temperature in the SU(3) Yang-Mills theory, we performed simulations on lattices with various spatial volumes and lattice spacings around the phase transition temperature T c . Separating generated configurations into hot and cold phases using the value of the Polyakov loop, and adjusting the temperature to T c by the multipoint reweighting method, we calculated the energy density and the pressure gap in each phase at T c , using the small flow time expansion (SFtX) method based on the 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 the continuum limit are summarized in Sec. 4.4. We obtained ∆ǫ/T 4 = 1.117 ± 0.040 for the spatial volume corresponding to the aspect ratio N s /N t = 8, and 1.349 ± 0.038 for N s /N t = 6. From a study of hysteresis curves for (ǫ + p)/T 4 around T c , we note that (ǫ + p)/T 4 in the (metastable) deconfined phase is sensitive to the spatial volume, while that in the confined phase is insensitive. The value of (ǫ + p)/T 4 decreases as the volume increases in the (metastable) deconfined phase. Study with higher statistics is needed to precisely determine the latent heat in the thermodynamic limit.
Furthermore, using the systematic data at various lattice spacings, we examined the effect of alternative procedures in the SFtX method. We compared two orders of the double extrapolation (a, t) → (0, 0) -the method 1 (first t → 0 and then a → 0) and the method 2 (first a → 0 and then t → 0). We also studied the effects of the renormalization scale and the truncation of higher order terms in the matching coefficients. We confirmed that the final results with alternative procedures are all consistent with each other. We also found that the use of NNLO matching coefficients improve the signal of the latent heat with the SFtX method.

A Energy-momentum tensor with NLO matching coefficients
In this appendix, we calculate the energy-momentum tensor using the NLO matching coefficients, following the original paper of the SFtX method [8]: where  For the NLO calculation in this appendix, following Ref. [8], we keep terms up to the next to leading order only in α U (t) and α E (t) and adopt the µ d -scale. We then havē  Method 1: t → 0 followed by a → 0 The results of the t → 0 extrapolation on each finite lattice are shown in Figs. A1, A2, and A3 by the lines, and the results in the t → 0 limit are summarized in Table 2. Continuum extrapolation of the results given in Table 2 is discussed in Sec. 4.2.
Method 2: a → 0 followed by t → 0 To carry out the continuum extrapolation at each t in physical units, we replot the data shown in Figs. A1, A2, and A3 as function of tT 2 c . For the aspect ratio N s /N t = 8 and 6, we obtain Figs. A4 and A5. The red, green and blue symbols are the results of N t = 8, 12 and 16. These figures corresponds to Figs. 9 and 10 with the NNLO matching coefficients. We find that the difference among the results with different lattice spacing becomes smaller as the flow time t increases. Repeating the analyses of Sec. 4.3, we obtain the ∆(ǫ − 3p)/T 4 using clover-shaped operator (red) and plaquette (blue) in the continuum limit for N s /N t = 8 (upper) and 6 (lower) by the NNLO calculation with µ = µ 0 , respectively. continuum limit shown by the orange curves in Figs. A4 and A5. These results of ∆(ǫ + p)/T 4 and ∆(ǫ − 3p)/T 4 in the continuum limit are summarized in Fig. 11 as functions of t, and the results of the t → 0 extrapolation are given in Table 3.

B Choice of lattice field strength operator
In this appendix, we test the influence of lattice operators for the field strength G µν (t, x) to define the operators U µν (t, x) and E(t, x) in Eqs. (3) and (4). Two major choices for G µν (t, x) are the plaquette and clover-shaped lattice operators.
In the left and right panels of Fig. B1, we compare the results of ∆(ǫ − 3p)/T 4 and ∆(ǫ + p)/T 4 in the continuum limit obtained by the clover-shaped and plaquette operators adopting the method 2 discussed in Sec. 2.4. The upper and lower panels are the results for N s /N t = 8 and 6, respectively. The blue and red symbols are for results with clover-shaped and by plaquette operators, respectively. The two results agree well with each other at sufficiently large t. The disagreement at tT 2 c < ∼ 0.01 is caused by the lattice discretization error, and the data there should be removed in the t → 0 extrapolation needed in the SFtX method. Because we have a wider range in t for the linear t → 0 extrapolation with the clover-shaped operator, we use them for the calculation of U µν (t, x) and E(t, x) in this study. With the clover-shaped operator, we may use data down to tT 2 c ∼ 0.005.