Equation of state in (2+1)-flavor QCD at physical point with improved Wilson fermion action using gradient flow

We study the energy-momentum tensor and the equation of state as well as the chiral condensate in (2+1)-flavor QCD at the physical point applying the method of Makino and Suzuki based on the gradient flow. We adopt a nonperturbatively O(a)-improved Wilson quark action and the renormalization group-improved Iwasaki gauge action. At Lattice 2016, we have presented our preliminary results of our study in (2+1)-flavor QCD at a heavy u, d quark mass point. We now extend the study to the physical point and perform finite-temperature simulations in the range T \simeq 155--544 MeV (Nt = 4--14 including odd Nt's) at a \simeq 0.09 fm. We show our final results of the heavy QCD study and present some preliminary results obtained at the physical point so far.


Introduction
Ab initio calculations of the QCD Equation of State (EoS) at vanishing baryon density have been in the focus of research for more than a decade, see [1,2,3,4,5,6,7] and references therein. These calculations were performed with different varieties of improved staggered fermions, nearly physical quark masses and 2+1 or 2+1+1 flavors of sea quarks. As a result, the EoS in the continuum limit has been determined for 2+1 flavors using two different staggered discretizations (HISQ [8] and stout) up to temperatures of about T = 400 MeV [4,5] and for 2+1+1 flavors using only the stout formulation up to temperatures of about T = 1 GeV [6].
Overall the different 2+1 flavor calculations are in good agreement, although there is a slight tension at the highest temperatures. By comparing the 2+1 and 2+1+1 flavor calculations using the stout formulation it was concluded that the contribution from the charm sea is already significant at T ≈ 400 MeV. However, as the presence of the charm quark in the sea may lead to a suppression of the other contributions to the EoS, the overall effect may be subject to subtle compensations. Moreover, the 2+1 and 2+1+1 flavor results are compared in the region where the tension between both formulations has been found. From this discussion it is clear that dynamical charm quarks should be included in the lattice determination of the EoS at higher temperatures for use in phenomenology and that an independent calculation is desirable. Efforts in this direction with the HISQ formulation are still on-going [9]. In particular, the large discretization errors associated with the charm quark mass strongly suggest that such a calculation should be performed with the HISQ action. Yet the 2+1 flavor EoS is needed with higher precision for a robust, quantitative determination of the correct charm contribution to the pressure including all of the subtle effects.
For a comparison to weak-coupling calculations a lattice determination without the charm quark may actually be advantageous. This is evident from the absence of massive quark degrees of freedom in the weak-coupling calculations [10,11] at the highest available order, NNNLO 1 . At temperatures in the range between T ≈ 400 MeV and some multiple of the charm quark mass, m c ≡ m c (m c ) = 1.277(10) GeV [12], corrections due to the finite charm quark mass are almost certainly relevant and large. Unless the applicability of weak coupling is firmly established for the respective observables at these temperatures, there is a risk of conflating the uncertainties associated with the temperature scale with those of the charm mass corrections. Precise results with 2+1 flavors at high temperatures from the lattice alleviate the concerns regarding applicability of weak coupling considerably and establish a robust baseline for determining the actual charm quark contribution.
For these reasons we extend the previous calculations of the 2+1 flavor EoS using the HISQ formulation to higher temperatures. We use lattices with larger light sea quark mass, m l = m s /5 (compared to m l = m s /20 in [5]), and discuss the light quark mass dependence of the trace anomaly. We discuss the cutoff dependence of the trace anomaly and the pressure in detail and determine the continuum limit or continuum estimates using different strategies. The rest of this paper is organized as follows. In Sec. 2, we discuss the details of the lattice calculations. In Sec. 3, we discuss the trace anomaly and its cutoff effects. In Sec. 4, we obtain results for the pressure and scrutinize its cutoff effects. We compare to weak-coupling results in Sec. 5

Details of the lattice calculations
The objective of our research is extending the calculation of the 2+1 flavor QCD Equation of State in [5] to higher temperatures and compare the EoS to results obtained in the weak-coupling limit. Following [5] we use the tree-level improved Symanzik gauge action and the Highly Improved Staggered Quark (HISQ) action for quarks. We use gauge ensembles at finite temperature with N τ = 12, 10, 8, 6 and 4, which have been generated for studies of the pseudocritical temperature [13] and the EoS [5] by the HotQCD collaboration, and for studies of the entropy shift due to a heavy quark in the thermal medium [14] and color screening [15] by the TUMQCD collaboration.
In order to cancel the UV divergences in the trace anomaly at T > 0, we subtract the trace anomaly at T = 0 using the same bare parameters. We use gauge ensembles at T = 0 from the HotQCD collaboration [5,13] with a pion mass of m π = 160 MeV in the continuum limit. Since the corresponding lattice spacings are insufficient for reaching higher temperatures with small discretization errors, i.e. large N τ , we generated new ensembles with finer lattice spacing but larger pion mass, m π = 320 MeV [7]. We use the rational hybrid Monte Carlo at five values of the coupling, β = 10/g 2 , see Tab. 1 for the parameters. For the three finest ensembles, β ≥ 8, topological tunneling only takes place during pre-thermalized molecular dynamics evolution, but is too suppressed after thermalization. We have generated multiple streams of MD trajectories with the same lattice parameters that randomly stall in different topological sectors, Q ∈ {0, 1, 2}. We find among the observables contributing to the EoS only for the light quark condensate ψψ l a statistically significant dependence on the topological charge Q. The respective changes of the light quark contribution to the trace anomaly are smaller than the statistical error of the much larger contribution from the gauge action. We show in the left panel of Fig. 1 that the static energy is consistent for different values of Q in the finest ensemble β = 8.4. For the other ensembles the picture is similar. We conclude that systematic errors due to frozen topology are numerically irrelevant for our study.
We fix the lattice scale with the static energy, i.e. we calculate r 1 and r 2 defined via  0.3% lower for β = 7.825 we use r 2 to set the scale for the finer ensembles, which is consistent with the non-perturbative beta function of [5]. The observed variation of the static energy due to the frozen topology charge is smaller than the quark mass dependence and safely covered within the statistical errors. Cominbing the new and old results for the lattice scale we obtain a consistent beta function with smaller uncertainty [7], which we parameterize with an Allton-type Ansatz where b 0 = 9/(16π 2 ) and b 1 = 1/(4π 4 ). The other coefficients c 0 = 43.12 (18), c 2 = 347008(32131) and d 2 = 5584(599) are determined with a fit, see [7]. In our analysis we combine ensembles with different quark masses, using the ensembles with larger pion mass for the UV subtraction only at high temperatures, T 470 MeV, where numerical effects of m l are known to be small for m l < 0.4m s [3]. Due to T = 1/(aN τ ) these ensembles provide access up to T = 2 GeV with N τ = 4.

Trace anomaly
In order to obtain the 2+1 flavor QCD Equation of State we first calculate the trace anomaly, i.e. the trace of the energy-momentum tensor, Θ µ µ = ε − 3p on the lattice. Θ µ µ can be expressed in terms of the expectation values of the gauge action and the scalar condensates of the quark flavors in the sea, i.e. light and strange. Since these densities diverge linearly in the continuum limit, we consider differences between the results at zero temperature and at finite temperature. Finally, we account for multiplicative renormalization bu rescaling these differences with the non-perturbative beta function R β and mass renormalization function R m . For the HISQ action we may write using the same notation as in [5]. R β and R m are defined via We show the two contributions to the trace anaomaly in Fig. 2. The gauge field contribution, which we show in the left panel, is systematically larger for the ensembles with larger sea quark mass. Although we cannot correct for this implicit quark mass dependence, the difference is covered by the statistical errors for temperatures above T 400 MeV, i.e. in the range where we use these data. For N τ ≥ 8 discretization errors are quite visible, but mostly covered by the rather large statistical errors in this temperature window. The fermion contribution, which we show in the right panel, has an explicit dependence on the sea quark mass. As we use a value m l = m s /20 instead of m l = m s /5 in (3.3) to account for this, we cannot resolve a residual quark mass dependence in Θ µ µ F anymore for T 400 MeV. Moreover, statistical errors and cutoff effects are very small. We separately discuss the trace anomaly for low and high temperatures.   We revisit the low temperature results for the trace anomaly since we require higher precision at low temperatures in order to take the continuum limit in terms of the pressure. In this study we use gauge ensembles at low temperatures with high statistics that have been generated for a study of the Polyakov loop in the crossover region by the TUMQCD collaboration [14]. These ensembles provide access to T = 123 MeV with N τ = 10 and to T = 133 MeV as well as T = 140 MeV with N τ = 12. We show our low temperature results for the trace anomaly with N τ = 12 and 10 in the left panel of Fig. 3 together with bands obtained from interpolating splines whose errors have been estimated from the bootstrap method. The difference between both in the window 160 MeV < T < 180 MeV is indicative of rather large cutoff effects. We compare our results directly to predictions of hadron resonance gas (HRG) models. Namely, we consider a model with only the states listed by the PDG -labeled HRG-PDG -and a model including also missing states, i.e. states predicted by the quark model that have not been confirmed experimentally -labeled HRG-QM. We find good agreement of both HRG models with the data up to T < 140 MeV. For higher temperatures the difference between both models is significantly and they cannot describe the data well.  At high temperatures T > 400 MeV, discretization errors and quark mass dependendence are very much suppressed. We show the HISQ results together with older results obtained with the p4 action and a quark mass of m l = m s /10. Whereas the results with N τ = 4 and N τ = 6 are systematically lower, the results for N τ ≥ 8 do not depend on N τ or the light quark mass within their uncertainties. The smallness of the cutoff effects can be understood from the weak-coupling picture, where the trace anomaly has discretization errors starting at three-loop order, i.e. O(α 2 s ). Since the running of α s is controlled by the temperature, discretization errors are strongly suppressed at high temperatures and their respective temperature dependence is mild, i.e. logarithmic.

Pressure
Eventually we obtain the Equation of State and the pressure from Θ µ µ via the integral method, Here, we assume that we know the pressure p(T 0 ) at the reference temperature T 0 . If T 0 is sufficiently small, the pressure may be set to zero or taken from a calculation using an HRG model. For the temperature range in our simulations, we have to adhere to the latter case. Moreover, we have to account in the HRG for the cutoff effects due to distortion of the hadron spectrum with the HISQ formulation. In short, for each species of hadrons there are multiple tastes of the same hadron with masses that are enlarged by discretization errors to varying degree. This distortion is most severe in the sector of pseudo-Goldstone bosons. Hadrons with larger masses contribute less to the pressure. The effect of this distortion is suppressed in the trace anomaly since hadrons with larger masses contribute more strongly to the trace anomaly and partly compensate for the distortion. We use the HRG-QM where only the pseudo-Goldstone bosons are modified, and estimate a systematic uncertainty from the difference to the HRG-QM where all ground states or all ground  and excited states are modified, see [7] for a more detailed discussion. In the left panel of Fig. 4 we show the lattice result for the pressure at low temperatures together with a prediction of the HRG with a distorted spectrum. In the right panel, we show the lattice QCD results obtained via the integral method together with the HRG prediction. We summarize the reference temperatures  and pressures for different N τ in Tab. 2. Clearly the dominant cutoff effects in the pressure at low temperatures are due to the taste-symmetry violation in the sector of pseudo-Goldstone bosons, i.e. they are of the form O(α s a 2 ) for the HISQ action. At fixed temperature T = 1/(aN τ ), we may neglect the running of α s with T and parameterize cutoff effects by powers of 1/N 2 τ . For high temperatures the picture is rather different and considerations assuming weak coupling apply. For the HISQ or p4 actions the pressure and the quark number susceptibilities (QNS), are quite close to the result for the ideal gas limit. The dominant discretization errors in the ideal gas limit are starting with a one-loop contribution, i.e. O(α 0 s ), and can be parameterized as starting with the power 1/N 4 τ = (aT ) 4 at fixed temperature. Inspection of the lattice results for the pressure and second order QNS shown in Fig. 5 indicates that the lattice results for both are about 15% below the ideal gas limit for T 400 MeV and follow the same pattern of cutoff effects.
These considerations open up four different approaches to obtain a continuum result for the pressure. First, we may follow [5] and obtain the continuum limit of the trace anomaly and then   Figure 5: Left: The pressure in the entire temperature range. The horizontal lines correspond to the free theory result. Also shown are the results for the pressure obtained with p4 action and N τ = 6 or 8 [1,2]. Right: The discretization errors of the second order QNS follow the same pattern. [16].
obtain the continuum pressure via the integral method. Second, we may use the integral method at finite N τ to obtain the pressure and extrapolate the pressure at fixed temperature to the continuum limit using the appropriate parameterizations of the cutoff effects as discussed in the preceding paragraphs. For an intermediate temperature window, 200 MeV < T < 400 MeV, the picture is less clear and cutoff effects can be described by 1/N 2 τ , 1/N 4 τ or combinations thereof. Continuum results for either parameterization are consistent within errors in this window, although the continuum limit for 1/N 2 τ form is slightly higher and the continuum limit for the fits with combined form has large uncertainties due to having only one degree of freedom. We consider the differences as a measure of the systematic uncertainty and use the continuum result for 1/N 4 τ form as our central value. We double the corresponding statistical error in this window, since this is numerically even more conservative than adding the systematic error in quadrature. For temperatures above T > 670 MeV only data with N τ ≤ 10 is available and we may consistently extrapolate with controlled uncertainties using the same form. Above T > 800 MeV we have only N τ ≤ 8 and cannot control the uncertainty of an extrapolation, but provide only a continuum estimate. For T ≈ 800 MeV a tentative extrapolation using only N τ ≤ 8 yields within errors the same coefficients as an extrapolation using N τ ≤ 10. Thus, the continuum estimate seems to be on rather solid ground. This is not completely surprising, since we expect that the coefficients have only a mild temperature dependence due to the logarithmic running of the coupling. In the same spirit, we use the coefficient of 1/N 4 τ at T = 1 GeV to estimate the correction for the pressure with N τ = 6 at higher temperatures, up to T = 1.33 GeV. Third, we may use our knowledge of the cutoff dependence of the QNS to correct for the discretization errors in the pressure. In the weak-coupling picture, the pressure can be written as a sum of quark and gluon pressures p q and p g , where the latter has only negligible cutoff effects for an improved gauge action. Then we may assume that the cutoff dependence is completely carried by p q , which is dominated by the light quarks, and that it ought to be similar to the cutoff dependence of the second order light QNS χ l 2 . Finally, we may assume that the continuum limit of p q at high temperatures can be estimated using the rescaled ideal quark gas pressure p q id . On the grounds of these considerations, we may develop a scheme for removing the discretization errors of the pressure by using the known results for QNS, , p q (T ) ≈ 0.85p q id (T ). (4. 3) The advantage of determining corrections from using QNS is that these calculations do not rely on the availability of ensembles at zero temperature with the same parameters. As such, it is considerably cheaper to determine the discretization errors of χ q 2n than for the pressure. We use the second and third approaches to determine the continuum limit and results that are corrected for cutoff effects up to T = 1.33 GeV. In Fig. 6 we show these results together with results obtained with p4 action [1] that have been subjected to the same correction approach using corresponding QNS results with p4 action [17]. We interpret this as a validation of the correction procedure that the corrected results with different actions approach the same continuum limit from opposite sides, cf. Fig. 5, and are consistent for T 400 MeV with N τ ≥ 8. We note that our continuum result for T = 500 MeV is one and a half standard errors larger than the 2+1 flavor result for stout action [4] and agrees well with the previous result for HISQ action [5], although our result has significantly smaller errors. Eventually this is still not sufficient for reaching the highest   temperatures T ≈ 2 GeV and we need a fourth approach. Hence, we use the first approach to obtain a continuum estimate for T > 1.33 GeV. Recalling that discretization errors in the trace anomaly with N τ ≥ 8 are small for T > 300 MeV, we obtain a continuum estimate for Θ µ µ by performing a spline interpolation of these data for T ≥ 300 MeV. Comparison of this continuum estimate and the data with N τ = 6 or N τ = 4 reveals that these are in the temperature window 800 MeV < T < 1 GeV smaller by a factor 1.4 or 1.2 respectively. Assuming a mild temperature dependence for the trace anomaly at high temperatures as suggested by the weak-coupling limit, we use the same rescaling factors 1.4 or 1.2 respectively to estimate corrected results for the trace anomaly from the data with N τ = 6 or N τ = 4. We then perform a combined spline interpolation of the trace anomaly with N τ ≥ 8 for T ≥ 400 MeV and the corrected data with N τ = 6 or N τ = 4 for T > 1 GeV, where we assign a 40% or 20% systematic error to these corrected data to be conservative. Finally we apply the integral method to this continuum estimate and fix the reference point through the continuum limit of the pressure at T = 660 MeV. In order to be sufficiently conservative we double the errors of this continuum estimate for the pressure and show it in the right panel of Fig. 6. We stress that none of these error enlargements are necessary in order to obtain full consistency with the continuum limit or the QNS corrected lattice data over the whole common temperature range.

Comparison to weak coupling
Armed with the continuum result developed throughout the last section we embark on a comparison to the weak-coupling calculations. In Fig. 7 we show the lattice data for the trace anomaly, where data with N τ = 6 or 4 have been corrected, see the discussion in Sec. 4. We show the data together with analytic results obtained using hard thermal loop (HTL) perturbation theory at three loop [11] and dimensionally reduced effective field theory, namely the electrostatic QCD (EQCD) at order O(g 6 ) [10]. The uncertainty band of the weak-coupling results are due to variation of the resummation scale µ. We see fair agreement between the different calculations. In Fig. 8, we compare the lattice result for the pressure and entropy density s = ∂ p/∂ T to the weak-coupling calculations. Our result is higher than the central value of the HTL result by about one sigma, but still covered by its scale uncertainty, and a few percent lower than the EQCD prediction. We also compare the entropy to the resummed calculation in next-to-leading log approximation (NLA) [18], which is higher than our result, but consistent within uncertainties for T > 1.3 GeV. From our comparison it is evident that the presented lattice data are sufficiently precise in order to determine through a comparison with weak-coupling results, whether the thermodynamics of the quark-gluon plasma can be understood in terms of the weak-coupling picture in the given temperature range.

Conclusions
We . We show the HTL [11] (red) and EQCD [10] (green) results each with a line (µ = 2πT ) and a band (πT to 4πT ). For the entropy density the blue band corresponds to the resummed calculation in next-to-leading log approximation (NLA) [18].
action to higher temperatures. In order to extend the previous result to higher temperatures we have used gauge ensembles with a light sea quark mass m l = m s /5 and demonstrated that the quark mass dependence of the trace anomaly is covered by the statistical uncertainty for T > 400 MeV. We have determined the pressure from the trace anomaly using the integral method. At low temperatures, we have used a hadron resonance gas with distorted spectrum as reference and have demonstrated that it can describe the trace anomaly for T < 140 MeV. We have studied the cutoff dependence of the pressure at low and high temperatures and extrapolated the pressure to the continuum limit. Moreover, we have developed a correction scheme for the cutoff effects in the pressure based on the known cutoff dependence of the quark number susceptibilities and the known weak-coupling limit at high temperatures. We demonstrated the validity of this correction scheme by applying it to results with the HISQ and p4 actions at finite N τ , which become numerically consistent with the continuum limit after correction. Finally, we have developed another correction scheme for the trace anomaly with N τ < 8 at high temperatures, have demonstrated that it agrees with the other results and have extracted a continuum estimate up to T ≈ 2 GeV. We have compared this estimate including conservative errors to weak-coupling predictions for the trace anomaly, the pressure and the entropy density and find fair agreement at high temperatures. The lattice results for pressure and entropy density are in the middle between the predictions of different weak-coupling calculations.