Equation of state of the SU($3$) Yang-Mills theory: a precise determination from a moving frame

The equation of state of the SU($3$) Yang-Mills theory is determined in the deconfined phase with a precision of about 0.5%. The calculation is carried out by numerical simulations of lattice gauge theory with shifted boundary conditions in the time direction. At each given temperature, up to $230\, T_c$ with $T_c$ being the critical temperature, the entropy density is computed at several lattice spacings so to be able to extrapolate the results to the continuum limit with confidence. Taken at face value, above a few $T_c$ the results exhibit a striking linear behaviour in $\ln(T/T_c)^{-1}$ over almost 2 orders of magnitude. Within errors, data point straight to the Stefan-Boltzmann value but with a slope grossly different from the leading-order perturbative prediction. The pressure is determined by integrating the entropy in the temperature, while the energy density is extracted from $T s=(\epsilon + p )$. The continuum values of the potentials are well represented by Pad\'e interpolating formulas, which also connect them well to the Stefan-Boltzmann values in the infinite temperature limit. The pressure, the energy and the entropy densities are compared with results in the literature. The discrepancy among previous computations near $T_c$ is analyzed and resolved thanks to the high precision achieved.


Introduction
The equation of state (EoS) of strongly-interacting matter is of absolute interest in particle and nuclear physics, and in cosmology. Apart the obvious theoretical interest in such a basic property, the collective behaviour of strongly-interacting particles has determined the evolution of the Universe in its early stages. Those extreme conditions are now being reproduced and investigated at heavy-ion colliders, where the EoS is a crucial input in the analysis of data.
Lattice gauge theory is the theoretical framework where the EoS is determined from first principles. A first measurement in the SU(3) Yang-Mills theory was performed in Ref. [1]. The pressure p, the entropy density s and the energy density were computed with a numerical accuracy of about 1-2% up to temperatures of ∼ 5T c . The strategy that was used, the "integral" method 1 , has become the most popular technique in numerical investigations of the EoS. It is based on the direct measurement of the trace anomaly by Monte Carlo simulations. The pressure is then obtained by integrating in the temperature, while the entropy and the energy densities are calculated using the thermodynamic relation T s = ( + p). In Ref. [3] a refined version of the integral method was used to determine the EoS in a much broader range of temperatures, from 0 up to ∼ 1000 T c , with a target accuracy at the permille level. In the region near T c , the results of the two computations show significant discrepancies.
Severe limitations hinder the integral method. The need for 1 For a variant see Ref. [2]. the (zero temperature) subtraction of the ultraviolet power divergence, and the complicated procedure for determining the lines of constant physics make the computation very demanding numerically and technically involved, see Refs. [4,5] for recent reviews. Despite the impressive progress over the last few years [6,7], uncertainties in the EoS of Quantum Chromodynamics (QCD) are still rather large, and temperatures higher than a few hundreds MeV are still unreachable with staggered fermions. The computation remains prohibitive with Wilson fermions. These obstacles are not rooted in the physics of the EoS, but in the method adopted for its computation. Recently there has been an intense activity to design new numerical strategies for simulating thermal field theories on the lattice and in particular to compute the EoS [8,9,10,11,12,13,14,15].
Among the proposed new methods, we have been focusing on the formulation of a thermal theory in a moving reference frame [8,9,10]. In this approach the entropy density is the primary observable computed by Monte Carlo simulations. It is extracted from the expectation value of the space-time component of the energy-momentum tensor T 0k in presence of shifted boundary conditions [10,12]. This quantity has no ultraviolet power divergences, and the finite multiplicative renormalization constant can be computed by imposing suitable Ward identities [10,16]. The pressure is then computed by integrating the entropy in the temperature, while the energy density is determined by integrating the temperature in the entropy or equivalently by using the relation T s = ( + p).
Following this strategy, in the last few years we have carried out extensive numerical computations in the SU(3) Yang-Mills theory also to prepare the ground for QCD. It is the aim of this letter to present the final results of this study.

Preliminaries
We regularize the SU(3) Yang-Mills theory on a finite fourdimensional lattice of spatial volume V = L 3 , temporal direction L 0 , and spacing a. The gauge field satisfies periodic boundary conditions in the spatial directions and shifted boundary conditions in the temporal direction where U µ (x 0 , x) ∈ SU(3) are the link variables, and the spatial vector ξ characterizes the moving frame in the Euclidean spacetime [8,9,10]. The action is discretized through the standard Wilson plaquette where the trace is over the color index, and g 0 is the bare coupling constant. The plaquette, defined as a function of the gauge links, is where µ, ν = 0, . . . , 3,μ is the unit vector along the direction µ, and x is the space-time coordinate. We are interested in the off-diagonal components of the energy-momentum tensor where the gluon field strength tensor is defined as [17] F a and , the minus sign standing for negative orientation. The T a are the Hermitian generators, see Ref. [16] for additional details.

The renormalization constant Z T
The field T 0k is multiplicatively renormalized by the finite renormalization constant Z T (g 2 0 ) which we have computed nonperturbatively in Ref. [16]. To have additional control on its discretization effects, we have also calculated the very same quantity at one loop in perturbation theory at finite lattice spacing [18]. By re-analyzing the data in Ref. [16] with the help of the one-loop improved definition, the final results for Z T (g 2 0 ) are well represented by the expression in the full range 0 ≤ g 2 0 ≤ 1. The error to be attached to this function, computed as in Ref. [16], is 0.4% up to g 2 0 ≤ 0.85 while it grows linearly from 0.4% to 0.7% in the range 0.85 ≤ g 2 0 ≤ 1. This new determination of Z T (g 2 0 ) is within one standard deviation from the one in Ref. [16], a fact which confirms that the systematics due to discretization effects is well within the quoted error.

The entropy density from the moving frame
When ξ is non zero, the entropy density at the temperature T = 1/(L 0 1 + ξ 2 ) can be written as [10] s(T ) Once Z T (g 2 0 ) has been determined, Eq. (7) provides a simple way for measuring the entropy density. The expectation value T 0k ξ can be measured by a single Monte Carlo simulation at each temperature and lattice spacing. No power divergences have to be subtracted, and the continuum limit is simply attained by increasing L 0 /a and by tuning the bare coupling constant g 0 so that the temperature stays unchanged in physical units. Since the observable is ultralocal, numerical simulations of lattices with large spatial sizes are not more expensive. The additional cost of updating a large lattice is compensated by the reduced statistical error due to volume averaging.
In our computation we opted for ξ = (1, 0, 0) for most of the temperatures, since very small lattice artifacts have been previously observed for the tree-level improved definition of T 0k ξ at this value of the shift [10,12,16]. We have used 2 ξ = (1, 1, 0)  for T/T c = 0.980, and ξ = (1, 1, 1) for T/T c = 0.904, 1.00, 1.061 and 2.30. In order to perform the continuum limit extrapolation of s(T )/T 3 with confidence, at each given temperature we have carried out numerical simulations with L 0 /a = 5, 6, 7, 8 and, sometimes, also 3, 4, 9 and 10. The inverse coupling constant 6/g 2 0 has been fixed as in Ref. [12]: from the Sommer scale r 0 /a [19] up to temperatures of 2 T c , while for higher temperatures from a quadratic interpolation in ln (L/a) of the data listed in Tables A.1 and A.4 of Ref. [20] corresponding to fixed values of the Schrödinger functional coupling constant g 2 (L). The critical temperature in units of the Sommer scale is r 0 T c = 0.750(4) [1,21]. For temperatures above 2T c , the accuracy on the temperature is about 2% [20] but its effect is negligible due to the very weak dependence of the thermodynamic quantities on the temperature. Error bars of the Monte Carlo data have been estimated using both the jackknife method and the binnining technique finding consistent results. In order to keep finite volume effects below the statistical errors, the lattice size in the spatial directions has been chosen to be L/a = 128 for L 0 /a up to 6 and L/a = 256 for larger values. Thus, LT ranges from 10 to 26, values where finite size effects are negligible within our statistical errors [10,12].
In Figure 1 we show the extrapolation to the continuum limit of the tree-level improved definition of s/T 3 at temperatures T/T c = 1.2, 1.5 and 2.0. As expected, lattice artifacts turn out to be very small. Results at other temperatures are qualitatively similar. Data generated with ξ = (1, 0, 0) are also compatible with a constant behaviour, i.e. no lattice spacing effects are observed within errors. Most of the uncertainty on s/T 3 comes from Z T , since in our simulations T 0k ξ is measured always with a statistical accuracy of permille or better. We extrapolated to the continuum limit s/T 3 linearly in (a/L 0 ) 2 . When discretization effects are not visible, we also fit the data to a constant taking into account the correlation among the values of the renormalization constant at different couplings. The fit   [22], each one including up to the order indicated in the legend. The continuous (red) curve is the O(g 6 ) perturbative prediction, but with a non-null unknown term fixed so to be able to reproduce the data in the temperature range considered.
are always similar to those showed in Figure 1. All results are reported in Table 1, and shown in Figures 2 and 3. Mainly for the computation of the pressure, see below, we have also carried out a rough measurement of the entropy density at T c in the deconfined phase. The continuum extrapolation gives s/T 3 c = 1.70 (24), where an estimate of the systematics is included in the error. The analogous rough measurement in the confined phase gives s/T 3 c = 0.37 (15). Above a few T c and within the statistical errors, s/T 3 exhibits a striking linear behaviour in ln(T/T c ) −1 over almost 2 orders of magnitude, see Figure 2. For T → ∞ data point straight to the Stefan-Boltzmann value, within errors, but with a slope almost 5 times smaller in magnitude than the leading-order perturbative prediction.

The equation of state
For T ≥ T c the entropy density is well represented by the Padé interpolating formula where w = ln(T/T c ), and the coefficients are: s 1 = 1.7015, s 2 = 77.757, s 3 = 232.33, s 4 = 19.033, s 5 = 32.200, s 6 = 6.9829, s 7 = −1.0348. The error attached to this curve at a given point can be safely estimated by interpolating linearly those of the two closest data points. This results in an uncertainty of 0.5% for T/T c > 1.5, which increases up to 1% going backward to 1.1. For T → ∞, the formula reproduces the Stefan-Boltzmann value within errors. Barring weird functional forms above 231.3 T c which cannot be logically excluded, Eq. (8) can be taken as a parametrization of the entropy for all temperatures above T c , with an error of 0.5% for T/T c ≥ 1.5. Within the rather large errors that we have for T ≤ T c , ln(s/T 3 ) can be fitted linearly in T c /T over the four data points in this range. The quality of the fit is very good.
Once the entropy density is known, the pressure is computed as For T ≤ T c , p is computed by integrating the curve resulting from the fit of ln(s/T 3 ) described above. For T ≥ T c , the integral is carried out as a sum of integrals between each couple of consecutive points at which s has been measured. Between any two such points, s/T 3 is interpolated by a quadratic curve in ln(T/T c ), with the three coefficients fixed by fitting the four data points closest to the integration region. The quality of the fits is always excellent, and the distance of the curve from the data is always a small fraction of the standard deviation. The statistical errors on p are computed by propagating linearly those on s, which in turn are dominated by the ones on the renormalization constant. The results for the pressure are reported in Table 1 and shown in the left plot of Figure 3. Being the entropy a very smooth function, the systematics due to the choice of the interpolating function is negligible with respect to the statistical error. As a further check, we have grouped the data in three samples of consecutive points and fitted them with independent Padé interpolants. By integrating the three of them in T , the results for the pressure are in perfect agreement with those in Table 1. Deviations are always a small fraction of the statistical error.
Analogously to the entropy density, for T ≥ T c the values for the pressure in Table 1 are well represented by the Padé interpolating formula where again w = ln(T/T c ), and the coefficients are: p 1 = 0.022288, p 2 = 2.0194, p 3 = 10.030, p 4 = 2.0941, p 5 = 5.6006, p 6 = 1.7469, p 7 = −0.79281 and p 8 = −0.30484. The error can be computed as for the curve of the entropy. It corresponds to 0.5% for T/T c ≥ 3, and it increases up to 5% going backward to 1.1. For T → ∞, the formula reproduces the Stefan-Boltzmann value within errors and can be taken as a parametrization of the pressure for all temperatures above T c , with an error of 0.5% for T/T c ≥ 3.
Once the entropy and the pressure are known, the energy density is computed by using the relation The results are again reported in Table 1 and shown in the left plot of Figure 3. The central values of can be computed by inserting Eqs. (8) and (10) into Eq. (11). The error is 0.5% for T/T c ≥ 1.3 while it grows up to 0.7% going backward to 1.1. The anomaly, computed as (T s − 4p)/T 4 , is shown in the right plot of Figure 3. It should be stressed once more that all results reported in this section refer to the continuum theory.

Discussion and conclusions
At our largest temperature, approximatively 230 T c , the entropy density still differs from the Stefan-Boltzmann value by roughly 3%. It raises linearly with ln(T/T c ) −1 following Eq. (8) which also connects well the data to the Stefan-Boltzmann value at T → ∞.
The perturbative expression is known to have a poor convergence rate [22]. If we use Λ MS r 0 = 0.586(48) [20,19], this can be seen in Figure 2 where the perturbative predictions at the various orders are shown. It must be said that if one fixes the O(g 6 ) undetermined coefficient so to be able to reproduce our results at large T , a reasonable description of the data above 5T c or so is achieved. This comes at the price of having an O(g 6 ) contribution at T = 231.3T c ∼ 68 GeV which is roughly 50% of the total correction to the entropy density given by the other terms. The perturbative formula is clearly of little help in determining the EoS of the theory in this interesting temperature range. The (essentially) linear functional form followed by the data in Figure 2, with a slope which is about 5 times smaller in magnitude than the leading pertubative prediction, may require sophisticated theoretical tools to be understood analytically, e.g. Refs. [23,24,25,26] and references therein.
For temperatures smaller than T c our results for the entropy density agree with those in Refs. [27,3,12]. Between T c and 3 T c they compare well with those in Refs. [1], that have significantly larger errors, and [12] while we observe a significant disagreement with the more precise ones in Ref. [3] as already reported in Ref. [12]. In particular we find a discrepancy of 4 to 6 standard deviations between T c and 1.5 T c , corresponding to a 2 to 4 percent effect, which becomes less than 2 standard deviations above 3 T c , see Ref. [28] for a detailed discussion. For temperatures larger than 3T c , our results are in good agreement with Refs. [1,3,12]. Similar discrepancies are propagated to pressure and energy density.
To understand the origin of the disagreement, it is instructive to compare directly the anomaly. This is the primary quantity computed in Ref. [3] from which all other potentials are derived. Results for the anomaly in Refs. [1] and [3] disagree significantly, see Figure 2 and 4 in [3] and [28] respectively. Between 1.1T c and 1.4T c , our results for (T s − 4p)/T 4 confirm those in [1] while they disagree with [3] by 2 to 4 standard deviations. Our continuum values, however, compare well with the finer lattice spacing results in Ref. [3], i.e. before the continuum limit extrapolation is carried out. We remark that similar disagreements with the data of Ref. [3] close to the peak of the trace anomaly have been reported also in Refs. [13] and [15]. Above 1.4T c no significant discrepancy is observed for the anomaly.
The study reported here demonstrates that lattice gauge theory with shifted boundary conditions offers a theoretically sound, simple and extremely powerful tool for an accurate determination of the EoS over several orders of magnitude in the temperature. Our results call for a non-perturbative determination of the EoS in QCD over an analogous range of temperatures where perturbation theory cannot help.

Acknowledgments
Simulations have been performed on the BG/Q Fermi and on the PC-clusters Galileo and Marconi at CINECA (CINECA-INFN and CINECA-Bicocca agreements), and on the PCcluster Wilson at Milano-Bicocca. We thankfully acknowledge the computer resources and technical support provided by these institutions.