Themodynamics for pure SU($2$) gauge theory using gradient flow

We study the equation of state of pure SU($2$) gauge theory using Monte Carlo simulations. The scale-setting of lattice parameters has been carried by using the gradient flow. We propose a reference scale $t_0$ for the SU($2$) gauge theory satisfying $t^2 \langle E \rangle|_{t=t_0} =0.1$, which is fixed by a natural scaling-down of the standard $t_0$-scale for the SU($3$) case based on perturbative analyses. We also show the thermodynamic quantities as a function of $T/T_c$, which are derived by the energy-momentum tensor using the small flow-time expansion of the gradient flow.


Introduction and motivation
Thermodynamic quantities (energy, pressure, thermal entropy, etc.) and transport coefficients (heat capacity, shear and bulk viscosities, etc.) play an important role to understand features of Quantum ChromoDynamics (QCD). In lattice ab initio calculations, the thermodynamic quantities [1,2] for the quenched QCD have been precisely obtained. These precise thermodynamic data are well explained by the ideal gluon gas and the Hard-Thermal-Loop (HTL) model [3] in the high temperature regime (5T c T ) and by the massive free glueball model in T T c , where T c denotes the critical temperature of the confined/deconfined phase transition. The consistency between Lattice QCD and the model studies gives us an intuitive picture of QCD state in these two limits.
As for the thermodynamic quantities, recently, Suzuki proposed a new method to calculate the energy-momentum tensor (EMT) using the gradient flow [4]. In the quenched QCD, the method works well for calculating the pressure and the energy density, and the statistical errors of the thermodynamic quantities reduce in the gradient flow method [5]. Furthermore, the method has been also extended to the full QCD [6] and the numerical simulations have been successfully performed [7][8][9].
Other important quantities are the transport coefficients. In the intermediate temperature, experimental data show the small shear viscosity-to-thermal entropy ratio (η/s), which is c The Author(s) 2012. Published by Oxford University Press on behalf of the Physical Society of Japan. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The trace anomaly (∆/T 4 ) normalized by the pressure in the Stefan-Boltzmann (SB) limit (P SB /T 4 ) as a function of (T c /T ) 2 . The circle (red) and the triangle (cyan) symbols, which are given in this work and in Ref. [20], denote the results for the SU (2) gauge theory, respectively. The dashed lines represent the interpolating functions of the results for several SU(N c ≥ 3) gauge theories, which are given by Eq. (16) in Ref. [21]. The data of the SU(N c ≥ 3) cases in the range (T c /T ) 2 ≤ 0.9 are well-fitted by the interpolating functions.
scale-setting function, the thermodynamic quantities in T c T 2T c are directly calculated from the EMT by the gradient flow method. We carefully take the continuum limit of the lattice thermodynamic data and then compare the results to the other calculations; the numerical integral-method and the analytical HTL calculations. Several independent methods are proposed to determine the equation of state by using the numerical simulations; the integral method [18], the moving-frame method [25], the gradient flow method [5], and the non-equilibrium methods [26], and it can extend these methods to the QCD including dynamical fermions [7,[27][28][29]. Among them, the gradient flow method has two advantages; the statistical uncertainties are smaller and the wave-function renormalization of the EMT is not necessary for pure gauge theories. This is a promising approach to overcome difficulties (i) and (ii). The structure of the paper is following: In §. 2, we briefly review the gradient flow equation. We propose a reference scale for the scale setting using the flowed action density in the SU (2) gauge theory, which is obtained by a natural scale-down of the standard t 0 scale for the SU(3) theory. Furthermore, we give a brief review of the calculation for the EMT using the flowed field variable. In §. 3, the simulation setup for the configuration generation and how to solve the gradient flow equation are explained. In §. 4, we show our results for zero-temperature simulations. We obtain the scale-setting equation, which gives the relationship between the lattice bare coupling constant and the lattice spacing. The ratios between our referencescale and the other quantities are calculated. In §. 5, the results for the finite-temperature simulations are shown. We obtain the thermodynamic quantities as a function of temperature 3/25 and compare the results with the other analyses. Section 6 gives the summary and future directions of this work.

Definition of the observables 2.1. Reference scale for the scale-setting equation using the gradient flow
The Yang-Mills gradient flow equation, which is a key equation in this work, is defined by where t denotes a fictitious time, which is called "flow-time", and A µ and B µ are quantum and deformed gauge fields by the flow in the SU(N c ) gauge theory, respectively. Note that the flow-time has mass dimension −2. The operator G µν is the field strength consisting of B µ . Thus the right hand side in Eq. (1) shows the same form as the equation of motion for the deformed gauge field B µ . The solution of the equation parametrized by the flow-time defines a transformation of the gauge field toward the stationary points of the gauge action. At the time t, the high frequency mode, whose momentum is larger than 1/ √ t, is suppressed (see Eq. (2.18) in Ref. [30]). The deformed field can be considered to the renormalized field by the nonperturbative transformation and the flow-time can be identified as a typical energy scale of the renormalization. Now, we introduce a length scale of the order of √ t to the local operators through the flow equation. This scale gives the renormalization condition to the local operators in spite of the fact that the operator itself has no characteristic length. Moreover, the flows do not separate the modes of the gauge fields in contrast to the case of the ordinary cutoff scale, and the gauge invariance is preserved. The one-loop computation [30] shows that the energy density E = (1/4)G a µν G a µν can be renormalized at the momentum scale of 1/ √ 8t and can be represented by the renormalized coupling constant; Here, γ E denotes the Euler constant and α s is the running coupling constant, α s = g 2 /(4π), at the scale µ R with the renormalized coupling constant g in the MS scheme. The next-nextto-leading order calculation of the E is given in Ref. [31]. The β-function for the SU(N c ) gauge theory is given by the following form; where the coefficients b 0 -b 3 are calculated in Ref. [32]; Here, ζ is the Riemann zeta-function and ζ 3 = 1.2020569 · · · . It is well-known that the βfunction does not depend on the renormalization scheme up to two-loop order.

4/25
The solution of the β-function gives the running coupling constant, which is explicitly written by where we putt ≡ ln µ 2 R Λ 2 , and µ 2 R = 1/(8t) at the gradient flow-time t. Once we fix the Λ MS , we can estimate the expectation value of E(t) as a function of flow-time. Note that E(t) in Eq. (2) takes a finite value in non-zero flow-time even in the continuum theory.
For the SU(3) gauge theory, a new reference scale, namely t 0 -scale, has been proposed in Ref. [30], where t 2 E(t) is a dimensionless quantity. One advantage of the usage of this reference scale is less statistical uncertainties since the operator E is a local operator. In general, we can take any reference value of A(≡ t 2 E(t) ). The discretization error becomes larger for smaller A.
The statistical error and the numerical cost grow up for the larger A. In A = 0.3 for the SU(3) gauge theory, the t 0 -scale gives almost the same scale as the Sommer scale (r 0 ) [30], and it is a suitable reference scale to investigate the thermal phase transition in the SU(3) gauge theory. Similar scales for the SU(2) gauge theory with the definition t 2 E = 0.2 and 0.3 [20] and the others [33] have been discussed. In this work, we would like to introduce the t 0 -scale in the SU(2) gauge theory as The perturbative analysis (Eq. (2)) shows that the coefficient of α s in the leading order is proportional to the number of gauge bosons. The reference value in Eq. (8) is chosen as an approximate scaling-down of the standard t 0 -scale in the SU(3) theory. We will show that our t 0 -scale is almost the same scale with the Necco-Sommer scale (r c ) and is good for the study of thermodynamics in the temperature region of T c T 2T c .

Thermodynamic quantities from the energy-momentum tensor
In this work, we will obtain the thermodynamic quantities from the energy-momentum tensor (EMT). In general, measurements of the EMT using the lattice numerical simulation are essentially difficult, since the lattice regularization manifestly breaks the general covariance, while EMT is a generator of the corresponding invariance [34,35]. Here, we calculate the EMT by using a new technique [5] based on the small flow-time expansion of the Yang-Mills gradient flow [4,30].
The key relation is given in Ref. [4] for the quenched QCD based on the small flow-time expansion.
where E 0 is the expectation value of the trace-part of the EMT. Here, we define the "correctly-renormalized conserved EMT" (T R µν (x)) by the subtraction of vacuum expectation values (v.e.v.) of the EMT.
Let us briefly review how to obtain the relation. There are two gauge-invariant local products of dimension-4, U µν (t, x) and E(t, x), in finite flow-time. These composite operators are UV finite in the positive flow-time (t > 0) and are manifestly given by As an advantage of the usage of the gradient flow, it is not necessary to calculate the wave function renormalization of operators thanks to their UV finiteness in pure gauge theories [36]. They can expand the dimension-4 and gauge covariant operators in terms of small flow-time in continuum theory; The EMT (T R µν ) is a conserved quantity, so that it must be a scheme-independent. Thus, the coefficients (α U (t) and α E (t)) should be also scheme-independent. However, it is hard to determine the coefficients nonperturbatively, so that we practically utilize the calculated values in the one-loop order in Ref. [4], where α s is the renormalized coupling constant in the MS scheme (Eq. (6)) at the scale µ R = 1/ √ 8t. The coefficientss 1 ands 2 depend on the renormalization scheme, which is the same as the scheme for α s . In the MS scheme, Although the coefficients are derived perturbatively, we assume that the relation is available in the nonperturbative regime. Then, we can obtain the renormalized EMT from Eq. (9), if we can calculate the values of U µν and E in a nonperturbative method, for instance, the lattice simulation. In practice, the operators U µν and E can be calculated by the clover-leaf operators on the lattice. In the finite-temperature system, the following components of the EMT correspond to a combination of the energy density (ε) and the pressure (P ) at the temperature T , 6/25 Here, s denotes the thermal entropy density, and the EMT at the temperature T is given by

Configuration generation
We utilize the standard Wilson-Plaquette gauge action defined on a four-dimensional Euclidean lattice, Here, g 0 and P µν denote the lattice bare coupling constant and the plaquette, respectively. Here U µ (x) represents the link-variable from a site x to its neighbor in µdirection and takes the value with the SU(N c ) group elements. Hereafter, we introduce the lattice parameter β = 2N c /g 2 0 . In this work, we adopt periodic boundary conditions for all directions. Gauge configurations are generated by using the pseudo-heatbath algorithm with over-relaxation. We use the word "a sweep" to refer to the combination of one pseudo-heatbath update step followed by these multiple over-relaxation steps. The mixed ratio of the combination is 1:20 for zero-temperature configuration, and for finite temperature, the ratio is 1:N τ . In order to eliminate the influence of autocorrelation, we take large enough number of sweeps (100 sweeps) between measurements. To see that, we observe the topology changing between each measured configuration (see Appendix A). The number of configuration is 100-600 and the lattice extent is N 4 s = 32 4 for the zero-temperature simulation. For the finite-temperature simulation, the number of configurations is 200 on N 3 s × N τ with the fixed aspect ratio N s /N τ = 4 for N τ = 6, 8, 10, and 12.

Gradient flow equation on the lattice
In §. 2.1, we briefly explained how the gradient flow works in the continuum theory. The method is also applicable to lattice formulations. The flow equation on the lattice is given in Refs. [30,37], Here, U µ (x) is the SU(N c ) link variable and V t is the deformed one. The S W in the right hand side denotes the Wilson-Plaquette action. The uniqueness and smoothness of the gradient flow on the lattice are guaranteed. Furthermore, the expectation value of a local operator after the flow by the lattice simulation has a well-defined continuum limit. The equation is an ordinary first-order differential equation so that we numerically integrate it using the Runge-Kutta (RK) method. The third-order RK formula with the initial 7/25 configuration V 0 = U µ (x) was given in Ref. [30]. Each step of the integration is obtained as where The error per step is of order ǫ 3 in the third-order RK method. We take ǫ = 0.01.

Scale setting 4.1. Lattice data of t 2 E(t) in perturbative regime
Before carrying out the scale setting simulation, we compare the operator E as a function of the flow-time between the perturbative result in the continuum theory and the numerical data on the lattice. To set the lattice parameter in the perturbative regime, we take β = 2.85.
The operator E = G a µν G a µν /4 can be constructed by the clover-leaf operator on the lattice. Furthermore, there is the other simpler definition; whereP t µν (x) represents a flowed-plaquette and P x is the set of unoriented plaquette in which the site of lower-left corner is x. In the continuum limit, the results should be independent of the definitions on the lattice.
It is worth to summarize several values for the SU(2) gauge theory with the Wilson-Plaquette action on the lattice.
The running coupling constant in the Schrödinger function (SF) scheme for the SU(2) gauge theory has been investigated in Ref. [38]. The nonperturbative β-function is estimated by the following polynomial b ef f 2 = 0.35 (12).
Here, L denotes the scale in the finite volume scaling and b ef f 2 is determined by fitting the nonperturbative running coupling constant with the functional form given in Eq. (26). On the large volume, where g 2 SF (L SF 0 ) = 4.765, we can calculate L SF 0 Λ SF by solving the β function; On the tree-level SF action with β = 2.85, the value of renormalized coupling g 2 SF = 4.765 is realized on a/L SF 0 = 0.0834 (5). Using the ratio of Λ parameters between SF and MS schemes, Λ SF /Λ MS = 0.44567 [38], we obtain the Λ MS in lattice units at β = 2.85;  2)) and the lattice result at β = 2.85. Figure 2 shows t 2 E as a function of t/a 2 , where the operator E on lattice is calculated by the clover-leaf. The thick (brown) bound denotes the lattice numerical result with 1-σ statistical error bar, and the black solid curve is the perturbative calculation using Eqs. (2) and (28). Here, the flow-time corresponds to the renormalization scale, so that the small flow-time regime must deal with the perturbative behavior. We confirm that the lattice and the perturbative calculations are consistent with each other in the small flow-time regime, as expected.

t 0 -scale in the SU(2) gauge theory
Now, we measure the t 2 E(t) in the range of β = 2.400-2.900 on 32 4 lattice. Figure 3 is the plot for the expectation values of t 2 E(t) at β = 2.500, 2.700 and 2.850. Here, we show two types of the data for each β; one is given by the plaquette definition (Eq. (25)) as the lattice operator E, while the other is obtained from the clover definition. It is known that there is a linear-scaling region of t 2 E(t) as a function of t for the SU(3) gauge theory. We found that the data in the SU(2) gauge theory also present the same property.
The values of t 0 in lattice units are summarized in Table 1. At β = 2.400, the value of t 0 in lattice units is smaller than unity. To avoid a strong lattice artifact, we drop the data at β = 2.400 in the following analysis. On the other hand, to avoid a finite volume effect we concentrate on 8t/a 2 ≤ 32/2 regime, since beyond the regime all link variables on the lattice are averaged under the periodic boundary condition. We also calculate the data at β = 2.900, but it does not grow up larger than the reference value in the regime so that we do not include the data at β = 2.900 in our analysis.
The data of ln(t 0 /a 2 ) are well interpolated using a quadratic function of β. The parametrization of ln(t 0 /a 2 ) is determined by the best fit-function;  (17) 600 Table 1 t 0 /a 2 and the number of configurations for each β. Figure 4 shows the ratio of the lattice spacing ln(a 2 /a 2 0 ), where we take the reference lattice spacing a 0 with the value at β = 2.50. As a comparison, the previous data, which are obtained from the scale-setting equation in Ref. [17], are also shown; for 2.27 ≤ β ≤ 2.60. Here, we estimate the error bar by using the covariance matrix of the chi-square fit. There is a precise agreement between the scalings given by two scale-setting functions within 1-σ error bar in 2.42 ≤ β ≤ 2.60, where both functions are available. Our data cover the higher β regime, which enables us to investigate the physics in T c T . We also do a similar calculation for the other reference scales. The results are summarized in Appendix B.

Relationship between t 0 and the other reference scales
The relationship between the following typical scales in the theory is useful to understand the dynamics. Left (right) panel in Fig. 5 shows the continuum extrapolation of the ratio between 10 Fig. 4 Ratio of the lattice spacing (ln(a 2 /a 2 0 )) as a function of β. Here, a 0 is the lattice spacing at β = 2.50. The black data are calculated by the scale-setting function Eq.(3.1) in Ref. [17], and the red data are given by Eq. (29). The error bars of our result are the same with the symbol size.
√ 8t 0 and Sommer scale (r 0 ) [39] (Necco-Sommer scale (r c ) [40]). Here, r 0 and r c are defined via the static quark-antiquark force; r 2 F (r)| r=r0 = 1.65 and r 2 F (r) r=rc = 0.65, respectively. We estimate these values in lattice units using the data of a 2 F (r) at β = 2.50, 2.60, and 2.70 shown in Table 1 of Ref. [39].  Fig. 5 Left (right) panel shows the continuum extrapolation for the ratio between √ 8t 0 and r 0 (r c ). Filled-square (blue) data are obtained from the clover representation in the calculation of E, while the open-circle (red) data are calculated using the plaquette representation in both panels.
The filled-square (blue) and open-circle (red) symbols in each panel are obtained as a function of t 0 /a 2 by using the clover and plaquette definitions for the operator E, respectively. 11/25 It is shown that the ratio of t 0 to r 0 or r c in the continuum limit takes the universal value in spite of the lattice definition of the operator. The gentle slope of the data obtained by the clover-definition indicates a small discretization error. In fact, the size of order a effects depends on the combination of the lattice action in the configuration generation, the lattice action in the gradient flow, and the lattice definition of the operators. The discretization effects for these choices in the tree-level are systematically discussed in Refs. [41,42].
The values in the continuum limit with the clover-definition of the operator E are given by √ 8t 0 r 0 = 0.6020(86) (40), where the first bracket denotes the statistical 1-σ error, and the second one shows the systematic error estimated by the deviation from the value calculated in the plaquette representation. The discretization effect is well under control since the systematic uncertainty is smaller than the statistical error. Although in the SU (3) Table 4 in Ref. [43], while the ones of the square (red) symbol are given in Table 1 of Ref. [17] .

12/25
It shows the 6-σ consistency, and the difference mainly comes from the discrepancy of the data of a 2 σ at β = 2.50 between Refs. [17] and [43]. Using these 8t 0 σ and Eq.(31), we obtain σr 2 0 = 1.51(8) by the former data, while σr 2 0 = 1.34(7) by the latter data. In Ref. [39], the value is independently obtained as σr 2 0 = 1.39 (50) via the lattice-size as a reference scale, so that it is consistent with our value within 1-σ error bar, which is obtained via the t 0 -scale.

Simulation parameters
The simulation parameters (β and N τ ) for several T /T c are determined by using Eq.(29) as shown in Table 2. The value of T c = [N τ a] −1 is fixed at β = 2.4265, N τ = 6, which is given in Ref. [44]. The estimated values of β at T c on the other N τ are obtained by Eq.(29) as β c = 2.514, 2.583, and 2.650 for N τ = 8, 10, and 12, respectively. On the other hand, the direct measurement of the critical values of β, which is determined by the Polyakov-loop susceptibility for each N τ , is summarized in Ref. [33]; β c = 2.510363(71), 2.57826 (14), and 2.63625 (35) for N τ = 8, 10, and 12, respectively. Our estimated values are consistent with the direct determinations in the two-or three-digit accuracy. In this work, we determine the value of β for each temperature in the three-digit accuracy in the finite-temperature simulation.  (29)). We did not use these parameters in the main analyses.
The thermodynamic quantities have been obtained using 200 configuration with 100sweep separations for each lattice parameter. The interval of the gradient flow-time for the measurement of the EMT is ∆t/a 2 = 0.01. 13/25

Simulation results: Thermodynamic quantities
The procedure to calculate the EMT on the lattice is summarized as following four steps in Ref. [5]; Step 1: Generate gauge configurations at t = 0 on a space-time lattice with the lattice spacing a and the lattice size N 3 s × N τ .
Step 2: Solve the gradient flow for each configuration to obtain the flowed link variables in the fiducial window, a ≪ √ 8t ≪ R, to suppress the discretization and the finite volume effects. Here, R is an infrared cutoff scale such as Λ −1 QCD or T −1 = N τ a.
Step 3: Construct U µν (t, x) and E(t, x) in Eq. (10) using the flowed link variables and average them over the gauge configurations at each t.
We have to carefully estimate the propagation of errors, in particular, taking the doublelimits in Step 4, since each flow-time data after taking the continuum extrapolation is correlated with each other. In this work, we use the jackknife method. Thus, firstly we generate the jackknife samples for the lattice raw data of U µν and E in each lattice parameter, and independently obtain the renormalized EMT by taking the double limits for each jackknife sample. Finally, we calculate the standard error from the deviation of the obtained results.
The left two panels in Fig. 7 show the dimensionless trace anomaly (∆/T 4 = (ε − 3P )/T 4 ) and the dimensionless entropy density (s/T 3 = (ε + P )/T 4 ) at T = 1.12T c as a function of the dimensionless flow parameter √ 8tT . Here, the error bar denotes the statistical errors. The lower limit of the fiducial window in Step 2 is indicated by the dashed lines in Fig. 7. This lower limit, which is related to the discretization error, is set to be √ 8t min = 2a, since we consider the clover-leaf operator with a size 2a. The upper limit of the fiducial window is shown in the dotted line. It is set to be √ 8t max = 0.35N τ a [45], so that the smearing by the gradient flow does not exceed the temporal lattice size. The data with statistical errors in Fig. 7 show the plateau inside the fiducial window (2/N τ ≤ √ 8tT ≤ 0.35) for each N τ . The right two panels in Fig. 7 present the temperature dependence of raw data in N τ = 12, and we can find that the plateaus inside the fiducial window are shown for all temperature regimes. It suggests that the contributions of the higher dimensional operator, which is proportional to t, are small in Eq. (11).
In Step 2, we carry out both constant and linear extrapolations for 160 datasets, namely 16 fixed flow-time in increments of 0.01 from 0.25 to 0.35 for each temperature listed in Table 2. Figure 8 shows the example plots of these continuum extrapolations at the fixed flow-time for T = 0.98T c , 1.08T c , and 1.76T c . The solid and dashed lines denote the constant and linear extrapolations of (T a) 2 = 1/N τ 2 , respectively. As central analyses, we utilize the constant extrapolation of three data, N τ = 8, 10, and 12, for 0.95 ≤ T /T c ≤ 1.76. Only for T = 2.07T c , where the number of the data points is few and we expect that the discretization error in N τ = 6 would be acceptable, the constant extrapolations by two data (N τ = 8 and 10) in the short flow-time regime ( √ 8tT < 1/3) and by three data (N τ = 6, 8, and 10) in the long flow-time regime ( √ 8tT ≥ 1/3), are carried as central analyses. The linear-fit results are used to estimate the systematic error from the scaling violation. In the continuum limit, 14 Fig. 7 Flow-time dependence for the expectation values of the trace anomaly (top two panels) and entropy density (bottom two panels) in each lattice parameter. Left two panels show the N τ (lattice spacing) dependence at T = 1.12T c . Right two panels show the temperature dependence at the fixed N τ = 12.
two extrapolated values of the trace anomaly (the entropy density) are consistent each other within 2-σ (3-σ), where σ denotes the statistical standard error. Figure 9 shows the t → 0 extrapolation. The solid and dashed lines denote the constant and linear extrapolations of t, respectively. We take the result, which has the better reduced χ 2 , as central values, and estimate a systematic uncertainty from the discrepancy depending on the extrapolation function. The systematic error for the trace anomaly is almost the same as the statistical error, while that of the entropy density is at most three times as large as the statistical one.
We also estimate the other systematic error coming from the uncertainty of the lattice determination of Λ MS . It changes the values of the coefficients α U (t) and α E (t) via the renormalized coupling constant, and the difference becomes larger in the lower temperature regime since the running coupling constant rapidly grows up in that regime. The corresponding systematic error for the trace anomaly is negligible in comparison with the statistical error, while the one for the entropy density changes the central value at most by 6%.
We plot the trace anomaly (red-circle symbols) and entropy density (black-square symbols) as a function of T /T c after taking the double (a, t) → (0, 0) limits in the left panel of Fig. 10. We also summarize the data in Appendix C. Here, the error bars in the figure and the table denote only statistical errors. As a comparison, the trace anomaly (cyan-triangle 15 symbols) and the entropy density (brown-diamond symbols), which are obtained by the integral method using the improved action [20] at the fixed temporal lattice extent N τ = 5, are also shown in the figure. In T /T c ≥ 1.12, our results including the above three types of systematic errors agree with the results given by the integral method. Note that our statistical errors obtained by a fewer number of configurations are smaller than the ones in the integration methods. In low temperature regime, there are small discrepancies among ours and the results given in Refs. [17,20]. Our data of the trace anomaly are larger than the results in Ref. [17], which utilizes the same lattice action and takes the continuum extrapolation. Furthermore, the result in Ref. [17] is larger than the one in Ref. [20], which utilizes the improved action and the fixed N τ . There are several possible reasons for the discrepancy, which appears especially in the low temperature. The one would come from the discretization errors, that exist in both methods, and it may be needed to take the extrapolations more carefully. The other reason for the gradient flow method would come from the usage of the one-loop approximation. The right panel of Fig. 10 shows the equation of state in T ≥ T c , namely the relationship between the energy density and the pressure. The linear function, P = ε/3, presents the 16 Fig. 10 (Left) Results of the trace anomaly (circle-red) and the entropy density (squareblack) as a function of temperature. The triangle (cyan) and diamond (brown) symbols denote the trace anomaly and the entropy density at N τ = 5 given in Ref. [20], respectively. All error bars denote only the statistical errors. (Right) The equation of state, namely the relationship between the energy density and the pressure. The dashed line, P = ε/3, denotes the case with vanishing trace anomaly. The diamond (blue) symbol represents the SB limit (ε/T 4 , P/T 4 ) = (π 2 /5, π 2 /15). case with vanishing trace anomaly, and the diamond (blue) symbol denotes the values in the ideal gas (Stefan-Boltzmann (SB)) limit. The nonperturbative interaction reduces the values of pressure in the whole energy-density regime. In high temperature regime, our result heads toward the point in the SB limit, but still, the lattice data with the highest ε, which correspond to T ≃ 2T c , are almost 70-85% of the value in the limit. It is evidence that the state of the two-color "QGP" phase around T ≤ 2T c cannot be described by the ideal gas model yet. 17 results in the next-next-to-leading order (NNLO) are also shown [3]. In both panels, the solid (black) curve denotes the result of HTL with µ HT L /(2πT ) = 1.0, and the two dashed curves present the ones with µ HT L /(2πT ) = 0.5 and 2.0, respectively. The brown bounds show the HTL results with the error bar coming from T c /Λ MS = 1.23 (11). The triangle (cyan) symbols in the right panel denote the results of the trace anomaly at N τ = 5 given in Ref. [20]. Now, let us compare our results with the analytic prediction, namely the results of the Hard-Thermal-Loop (HTL) model. The left panel of Fig. 11 shows the comparison of the energy density normalized by the value in the SB limit between our lattice result and the HTL calculation in N c = 2 case [3]. In HTL analyses, we use the next-next-to-leading order (NNLO) formula with the three-loop running coupling constant shown in Eq. (6) with the renormalization scheme µ HT L = 2πT . In the calculation, we use T c /Λ MS = 1.23 (11) obtained by the lattice data. The (brown) bound in Fig. 11 shows the uncertainty coming from the error in T c /Λ MS , where the renormalization scale is fixed as µ HT L /(2πT ) = 1.0. The systematic uncertainty of the choice of µ HT L is shown as two dashed curves obtained by using µ HT L /(2πT ) = 0.5, and 2.0. Our results are nicely consistent with the HTL results until near T c .
Finally, to see the scaling law of the trace anomaly more precisely, we also compare the results between the lattice data and the HTL analyses. The trace anomaly has a leading correction term of 1/T 2 in the high temperature regime, as we have shown in Fig. 1. The data can be well-fitted by the linear function of 1/T 2 in (T c /T ) 2 ≤ 0.6 regime. On the other hand, the nonperturbative logarithmic correction term for ∆(T )/T 2 is predicted by the HTL analyses. To see these correction terms, we take both axes to a logarithmic scale in the right panel of Fig. 11. In 1.3T c T , the lattice data exhibit almost plateau and approaches to the HTL results. In the SU(3) gauge theory [2], the similar behavior also appears around T = 2T c , and in the further high temperature the lattice data become consistent with the HTL and the perturbative analyses. We may consider that it also occurs in the SU(2) gauge theory.

Summary and future directions
In this work, we numerically investigate the thermodynamics of the pure SU(2) gauge theory. The theory is a good testing ground for the studies on the methodology and the N cdependence for the pure SU(N c ) and QCD theory. We determine the scale-setting function and thermodynamic quantities by using the gradient flow method in the N c = 2 case.
For the precise scale-setting of lattice parameters, we propose a reference scale t 0 for the SU(2) gauge theory, which satisfies t 2 E | t=t0 = 0.1. This value is determined by a natural scaling-down of the standard t 0 -scale for the SU(3) gauge theory. We have shown that the reference scale is suitable to study thermodynamics in the temperature region of T c T 2T c .
We also obtain the thermodynamic quantities, which are directly calculated from the EMT by the small flow-time expansion of the gradient flow method. This work is the first application of the gradient flow method to the thermodynamic quantities for the SU(2) theory. We take the double limits (a, t) → (0, 0), first a → 0 and then t → 0, to remove the artifact in the gradient flow method. The statistical error is smaller than the one given by the integration method, nevertheless, we utilize only a few hundred configurations. We precisely show that the trace anomaly in the pure SU(2) gauge theory has a different scaling property from the N c ≥ 3 cases. It shows the gently curved behavior near T c , and the linear behavior of (1/T 2 ) appears in (T c /T ) 2 ≤ 0.6 regimes. We also find a strong tendency toward the consistency with the HTL results in the high-temperature regime.
For future works, we address the following points.

Universality of the critical phenomena
The finite-temperature phase transition in the pure SU(2) gauge theory is characterized by the Polyakov loop, and it is Z 2 -symmetric/ broken phase transition. It is believed that the phase transition belongs to the same universality class as the Ising model in three dimensions since it has the same symmetry and spatial dimensions [46]. However, it is hard tasks to obtain precisely the same critical exponents with the one for the Ising model. Actually, we also observed the critical exponent of the scaling of the Polyakov loop, but still it is the similar value with the one in Ref. [47] and is not consistent with the value of the Ising model, β = 0.3265 (3). On the other hand, as for trace anomaly, the recent lattice Monte Carlo study on the 3d Ising model shows the evidence of the traceless of the EMT [48]. Meanwhile, the trace anomaly of the pure SU(2) gauge theory at the critical temperature with one compact-dimension has a finite value. However, this is not so strange, since the EMT can be affected by the all scale fluctuations near the transition point and is not expected to behave universally as in the 3d Ising model. Exact determination of α U , α E coefficients In our analysis, we utilize the one-loop calculation for the coefficients; α U and α E . Our results are consistent with the ones given by the integral methods in the high temperature regime, but there is a small discrepancy in T < T c even if the systematic errors are included. The discrepancy may come from the usage of the one-loop approximation in these coefficients. The analytical calculation of these coefficients in the higher order is future work. Furthermore, the nonperturbative determination of these coefficients has been proposed [49,50]. That is also a valuable direction.

19/25
Calculation of viscosities One of the motivations for this work is to prepare the determination of η/s, where η is the shear viscosity. The SU(N c ) gauge theory with large-N c has a lower bound of η/s = 1/(4π) based on the AdS/CFT correspondence. It is also suggested that the η/s takes the minimum value at the phase transition in a wide class of the systems [51]. However, the 1/N c corrections, even their signs of the correction terms, are unclear [12]. The lattice results on the pure SU(3) approach the bound near the critical temperature [52,53], and η/s on the pure SU(2) at T ≃ 1.2T c is also very close to the bound [54]. The systematic study on the temperature-dependence and the N c -dependence could reveal properties of the vacuum in the QCD(-like) theories via these viscosities. Constructing an effective model of two color QCD Once the temperature dependences of the energy density (or pressure) and the expectation value of the Polyakov loop are determined by the pure gauge lattice simulations, one can construct the effective Polyakov-loop potential used in effective model such as the PNJL model [55][56][57][58][59][60]. Using the effective model, we can analyze the physics in full QCD which contains dynamical quark contributions. Since two-color full QCD simulation is easier to be treated than three color one, it is also easier to check the efficiency of the constructed effective model in two color case. In particular, it is able to investigate the efficiency even at finite quark number density [61], since the lattice simulation of two-color full QCD is also feasible at finite quark chemical potential.

A. Distribution of the topological charge at finer lattice
To investigate the autocorrelation of the generated configuration, we measure the topological charge using the gradient flow. The topological charge is related with the vacuum structure, and it has a long autocorrelation among the observables in Yang-Mills theory. The gluonic definition of the topological charge in Euclidean space-time is given by where ǫ µνρσ denotes the totally antisymmetric tensor. In this gluonic definition, the topological charge of quantum configurations in lattice simulations generally does not take an integer-value because of UV fluctuations. The application of the gradient flow suppresses the UV fluctuations and recovers almost integer-valued quantity [30,62]. Figure A1 shows the topological charge for each 600 configuration in β = 2.85, in which the lattice spacing is the smallest and being most likely to occur a topological freezing in our simulations. Here, we observe the topological charge at the gradient flow-time t/a 2 = 32, where the effective smeared-regime is the same with the half of the lattice extent. We found 20/25 that the topological charge takes an almost integer-value as expected for each configuration and it frequently changes within one configuration separation. We conclude that the autocorrelation can be negligible in our data sets.

B. Scale-setting function for several reference scales
In our analysis in §.4.2, the scale-setting function is valid for 2.420 ≤ β ≤ 2.850. The range is determined by the choice of the reference scale (t 0 ) to avoid a strong lattice artifact (t 0 /a 2 ≥ 1) and a finite-volume effect ( 8t 0 /a 2 ≤ N τ /2) with the reference value A ≡ t 2 E | t=t0 = 0.1. However, if we take the other reference values, then we can extend the available regime of β.

Table B1
Results of t 0 /a 2 for A = 0.08, 0.09, 0.10 and 0.12. Table B1 shows the values of t 0 in lattice units for A = 0.08, 0.09, 0.10 and 0.12. In β = 2.400 and 2.420, the proper reference values should be taken in A ≥ 0.12 and A ≥ 0.10, respectively.
To see the scale violation coming from the choice of the reference values, we show the ratio of the lattice spacing (a/a 0 ) for each β in Table B2. Here, we rescale the lattice spacing with 21/25 the value at β = 2.600 for each A.  (18) Table B2 Ratios of the lattice spacings a/a 0 for A = 0.08, 0.09, 0.10 and 0.12, where we take the reference lattice spacing a 0 with the value at β = 2.600. almost independent of the choice of the reference values.
Finally, we obtain the scale-setting function with the following interpolating function; ln(t 0 /a 2 ) = α 0 + α 1 (β − 2.600) + α 2 (β − 2.600) 2 , where α 0 , α 1 and α 2 are the fitting parameters. The obtained parameters and the fit range for each reference value are summarized in Table B3. Although the constant term α 0 , which corresponds to ln(t 0 /a 2 ) at β = 2.600, depends on the choice of A, the other coefficients α 1 and α 2 do not change so much. It suggests that we can connect several scale-setting functions obtained by the different choice of A, and can obtain the extended scale-setting function, which is available in the wider regime.

Table B3
The coefficients of the scale-setting function for A = 0.08, 0.09, 0.10 and 0.12.