Full result for the QCD equation of state with 2+1 flavors

We present a full result for the 2+1 flavor QCD equation of state. All the systematics are controlled, the quark masses are set to their physical values, and the continuum extrapolation is carried out. This extends our previous studies [JHEP 0601:089 (2006); 1011:077 (2010)] to even finer lattices and now includes ensembles with Nt = 6,8,10,12 up to Nt = 16. We use a Symanzik improved gauge and a stout-link improved staggered fermion action. Our findings confirm our earlier results. In order to facilitate the direct use of our equation of state we make our tabulated results available for download.


Introduction
The early universe went through a rapid transition from a phase dominated by colored degrees of freedom to a phase dominated by color neutral degrees of freedom (hadrons). The same transition is now routinely reproduced in heavy ion collisions at the Large Hadron Collider (LHC, CERN) and at the Relativistic Heavy Ion Collider (RHIC, Brookhaven National Lab.). The only systematic theoretical approach to determine many of the features of this transition is lattice QCD (for recent reviews see, e.g., [1,2,3]).
Our most important qualitative knowledge about the transition is its nature. We know from lattice calculations that the transition is analytic. Therefore, without singular behavior, the system evolves smoothly from one phase to the other. (Since there is no real phase transition, the word "phase" merely indicates the dominant degrees of freedom.) This result [4] of the Wuppertal-Budapest Collaboration was obtained through a finite size scaling analysis of the continuum extrapolated observables, computed using physical quark masses. Since there is a theoretically debated technical procedure 1 required in any N f = 2 + 1 flavor staggered calculation, as a future project it is desirable to reproduce the existing result using another -preferably chiral -lattice formalism.
One of the most important quantitative parameters of the transition is its absolute scale, the transition temperature T c . In the absence of a real phase transition, there is, however, no 2 Action, its physical motivation and simulation setup We use a tree-level Symanzik improved gauge action with 2-step stout-link improved staggered fermions. The precise definition of the action can be found in ref. [18]. Though the naive power counting tells us that the gauge part is more improved (a 4 ) than the fermionic part (a 2 ), this is no expensive overkill: the gauge part is relatively cheap and one wants to avoid a situation in which both sectors have the same order cutoff effects, but the computationally cheaper gauge part has accidentally a much larger prefactor. Indeed, our experiences show that using stout-link smearing the scaling features are very good [7].
The main advantages of this action are threefold.
a.) It is computationally fast (even faster than the completely unimproved action). This feature allows one to go to finer lattices (and thus closer to the continuum limit) at given computational costs than with essentially any other action in the literature.
b.) The action has a very well behaving continuum extrapolation. Our action approaches the continuum value of the Stefan-Boltzmann limit in the T → ∞ limit slower than actions with p4 or Naik terms (the latter is an additional fermionic term in the asqtad and hisq actions). Nevertheless, our action is monotonous and reaches the asymptotic a 2 behavior quite "early". Extrapolations from moderate temporal extents, e.g., using N t ≥ 8, provide an accuracy on the percent level, which is the typical accuracy one aims to reach. Furthermore, stout link smeared actions are ultralocal, both in fermion space as also in the gauge background, with a small exponential locality range in the gauge background only [15]. According to experience, these features allow one to carry out a smooth continuum extrapolation. Furthermore, applying simple tree-level improvement factors for the bulk thermodynamic observables leads to results, which are already close to the continuum ones. Since the action is cheep to simulate, one can have several lattice spacings, which enables a controlled continuum extrapolation using many points. Other improved actions, with larger locality extent (p4 or Naik-type asqtad/HISQ) can have non-monotonic behavior (consider, e.g., the N t dependence of the free energy density for the Naik term) [17]. Furthermore, for T = 0 simulations, required to compute the LCP and also to renormalize the T > 0 data points, these actions have O(a 2 ) cutoff effects. Therefore, their improvement which is motivated by T → ∞ studies does not remove all O(a 2 ) lattice artifacts at T = 0.
c.) Staggered fermions are cheap to simulate, however, the correct spectrum is recovered only in the continuum limit. In the 2+1 flavor framework there are 3/16 pseudo-Goldstone bosons instead of 3 (these are the pions in continuum QCD). To compensate for the deficit there is a tower of much heavier non-Goldstone pseudoscalars. As we approach the continuum limit, these heavier states merge with the 3/16 pseudo-Goldstone bosons and finally form the 3 pions that we need. This so-called "taste-violation" at non-vanishing lattice spacing can be characterized by the splitting between the pseudo-Goldstone and the lowest level mass states in the pseudoscalar tower and also by the splittings within the tower. The smaller the splitting becomes (i.e. the closer the non-Goldstone get to the pseudo-Goldstones bosons) the faster one reaches the continuum limit. This effect turned out to be more important for staggered QCD thermodynamics than improving the action in the T → ∞ limit and motivated our choice of the two-stout smeared action for our large scale staggered QCD thermodynamics projects (see Figure 1 of ref. [18] or Figure 2 of ref. [8]). As of today, the new HISQ action possesses an even smaller taste violation (see, e.g., Figure 4 of ref. [9]), though at higher computational costs. Taste violation and/or heavier than physical quark masses increase the height of the trace anomlay's peak. This fact is illustrated in Fig. 16 of Ref. [14].
In the following four paragraphs we summarize the improvements over our previous EoS paper [14].
a.) In this paper the tree-level Symanzik improved gauge action with 2-step stout-link improved staggered fermions is used on T>0 lattices with five lattice spacings corresponding to N t =6, 8, 10, 12, and 16 temporal extensions. The finest lattice N t =16 is used to verify the earlier finding about the peak's height of the trace anomaly and to determine the additive renormalization (see below). This huge data set allows a fully controlled continuum extrapolation. In contrast, in ref. [14] we used N t =6, 8, and 10 up to T =350 MeV above which only N t =6 and 8 were used. At three characteristic temperatures also N t =12 was included and the continuum extrapolation was carried out (with relatively large errors). b.) The T→0 limit is particularly difficult to reach, since for a given N t lower and lower temperatures correspond to larger and larger lattice spacings, thus larger and larger taste violating effects. Nevertheless, this limit is important, since according to the standard choice, the renormalization is done at zero temperature: p(T =0)=0. A mismatch at T=0 leads to a shift in the whole EoS. In ref. [14] we calculated the difference in the pressure between the physical theory and its counterpart with 720 MeV heavy pions at a selected temperature (100 MeV) on N t = 6, 8 and 10 lattices. At this low temperature the latter theory has practically zero pressure, thus the difference gives the pressure of the physical theory p(T = 100 MeV) with the desired β fit to f K scale fit to w 0 scale scale from w 0 based step scaling scale from w 0 along LCP Figure 1: The lattice spacing a for the β range used in this paper. Down to a = 0.047 fm two quantities were used to determine the scale, f K and w 0 . For higher couplings a step-scaling approach was applied.
normalization. Using this technique in ref. [14] we obtained about half (p/T 4 =0. 16) the value of the prediction of the hadron resonance gas model (p/T 4 =0.27). This difference was included in the systematic error. Though this mismatch is really tiny compared to the Stefan-Boltzmann value of 5.209, it was obviously a suboptimal solution. Now we use all five lattice spacings (including N t =16) to fix the additive term in the pressure, arriving at a complete agreement with the hadron resonance gas model at low temperatures. c.) We determined the scale and the line of constant physics (LCP) with higher accuracy than in ref. [14]. For the scale in [14] a step-scaling technique was used on lattice spacings smaller than 0.073 fm. Here we have a direct determination of the lattice spacing based on the w 0 scale [19] down to 0.047 fm. Below that we again used a variant of the step-scaling method [20] to determine the β(a) function. To do this one needs a renormalized quantity which has a significant volume dependence even for very small physical volumes. We chose to use the derivative of the Yang-Mills flow of E(t) = F µν F µν [21] given by t · dt 2 E(t)/dt. The flow time was coupled to the box size as t = 0.01L 2 , so our renormalized observable was O = t · dt 2 E(t)/dt t=0.01L 2 . In the first step we used a physical volume corresponding to L = 24 · a min where a min is the smallest lattice spacing we could reach with our direct approach. Using this lattice and two coarser ones, 16 4 and 20 4 , both corresponding to the same physical volume, we extrapolated O to a finer lattice spacing a 1 = 24/32 · a min . Then, using simulations on 32 4 lattices we searched for the β 1 coupling which provides this value of O. The pair (β 1 , a 1 ) is the first new point of our scale function. In every further step this procedure was repeated: in the i-th step the physical box size was chosen as L = 24·a i−1 , the observable O was extrapolated using 16 4 , 20 4 and 24 4 lattices to a i = 24/32a i−1 and the 32 4 lattice was used to find β i which provides the extrapolated value of O and thus corresponds to a i . The difference between the f K and w 0 scale settings is included into our systematic error estimate. We show the scale in Figure 1. The LCP is defined by fixing the kaon decay constant to pion mass ratio f K /M π and the light to strange quark mass ratio m s /m ud to their physical values. (For the latter we used, similarly to our previous studies, a value of 28.15 as determined in ref. [7]. Most recent studies (e.g., refs. [22,23]) lead to a 2% lower value. Note that this is in the same ballpark as the accuracy of the LCP and, furthermore, this ratio has far less than this 2% influence on the EoS. We studied this dependency of the EoS on the mass parameters in [14].) d.) In ref. [14] we explicitly pointed out that "for a rigorous continuum extrapolation one would need N t =12 for the entire temperature region". With the present paper we fulfill this condition. Furthermore, we extend our analysis procedure to control the different sources of systematic errors, using our histogram method [15]. We considered various fit methods (each of which is in principle a completely valid approach), calculated the goodness of fit Q and weights based on the Akaike information criterion AICc [24] of that fit and looked at the unweighted or weighted (based on Q or AICc) distribution of the results. The median gives our central value, whereas the central region containing 68% of all the possible methods gives an estimate on the systematic uncertainties. This procedure provides very conservative errors. In the present case we had four basic types of continuum extrapolation methods (with or without tree level improvement for the pressure and with a 2 alone or a 2 and a 4 discretization effects) and two continuum extrapolation ranges (including or excluding the coarsest lattice N t =6 in the analysis). We used seven ways to determine the subtraction term at T=0 (subtracting directly at the same gauge coupling β or interpolating between the β values with various orders of interpolation functions). We applied two scale setting procedures; one based on the kaon decay constant and one, as described above, using the w 0 scale. Finally, we had eight options to determine the final trace anomaly by choosing among various spline functions. This gives altogether 4·2·7·2·8=896 methods. Note that using either an AICc or Q based distribution changed the result only by a tiny fraction of the systematic uncertainty. Furthermore, the unweighted distribution always delivered consistent results within systematical errors. The systematic error procedure clearly demonstrates the robustness of our final result. Even in the case of applying or not applying tree level improvement, where the data points at finite lattice spacing change considerably, the agreement between the continuum extrapolated results, and hence the contribution to the systematic error, is on the few percent level. Figure 2 shows the lattice extents and the collected statistics for our T = 0 runs. Asymmetric ones were used for the scale setting, symmetric lattices for renormalization (and, additionally, for w 0 scale setting). We had the highest number of trajectories (67k) for the 48 4 lattice which was used to renormalize the N t = 16 data at T = 214 MeV. At finite temperature essentially the same ensembles were used as in ref. [25], with additional simulations on 32 3 × 6, 32 3 × 8 lattices and ∼ 13×10 3 or 50×10 3 trajectories, respectively. To reduce the potential finite-volume effects we also added six ensembles of 48 3 × 12 lattices in the range T = 220 . . . 335 MeV with ∼ 3 · 10 4 trajectories, each.

Results
At high temperatures we face two technical challenges. a.) If the lattice geometry is kept constant the physical volume will drop and relevant scales might be absent from the lattice. To prevent this, we increased the volumes once more and generated ensembles with lattice sizes 32 3 × 6, 48 3 × 8, 64 3 × 10, and 64 3 × 12, with 5, 40, 10 and 12 ×10 3 trajectories, respectively.
b.) The renormalization runs at high temperatures require extremely fine lattices (below 0.05 fm lattice spacing), where T = 0 simulations are no longer feasible due to algorithmic limitations. Simulations are, however, still possible in the deconfined phase. Starting from a temperature of 335 MeV, we follow the strategy described in refs. [26,27]. Firstly, we subtract the value of the trace anomaly at the same coupling but doubled time extent (and thus a . Right: comparison of the result with the parallel effort using the HISQ action by the hotQCD collaboration (as it was presented at the Lattice 2012 conference [28], with f K scale setting) and the related parametrization 's95p-v1' of [29]. A comparison to the Hadron Resonance Gas model's prediction and our result [14] from 2010 ("WB 2010") is also shown.
Adding to this result the value of the trace anomaly at T /2 and the same N t , we get the total trace anomaly. For the half-temperature subtractions we generated ensembles on 48 3 × 16, 64 3 × 20 and 64 3 × 24 lattices with matching parameters and statistics to their finite temperature counterparts. The continuum extrapolated trace anomaly is shown in Figure 3 (left). In parallel to our investigations the hotQCD group is pursuing a similar strategy to calculate the continuum extrapolated equation of state using the HISQ action. As of the lattice conference 2012 the results appear to be inconsistent [28]. The situation might improve, however, when the HISQ analysis becomes complete with physical quark masses, a continuum extrapolation and a systematic error estimate.
The apparent discrepancy in Figure 3 is strongest in the peak region between 200 and 230 MeV, thus we have selected a temperature (T ≈ 214 MeV) where we use an N t = 16 data point to demonstrate the continuum scaling for our action. On the r.h.s. of Figure 4 we give the trace anomaly both with and without tree level improvement. The extrapolated values are consistent with each other, however, tree level improvement results in smaller uncertainties.
An ongoing project of the Wuppertal-Budapest collaboration is the determination of the equation of state with N f = 2 + 1 + 1 flavors, i.e. including the contribution of the charm quark. This is done using a different action with 4 steps of stout smearing with taste breaking artifacts roughly equal to the HISQ action. Our 4-stout action has been tuned independently from our previous efforts by bracketing the physical point to ±2% in the quark masses, in boxes with Lm π > 4. The scale was set using the pion decay constant f π = 130.41 MeV. The tuning strategy was pursued to sub-percent accuracy down to a lattice spacing of a ≈ 0.077 fm, which corresponds to T ≈ 214 MeV temperature at N t = 12. Since the contribution of the charm quark is expected to be far below the present accuracy at T = 214 MeV [30], a continuum limit result for the trace anomaly at this temperature, computed using this new action, provides an independent cross-check of our 2-stout results. Therefore, in order to verify our 2+1 flavor continuum extrapolation against a different action and scale setting technique we made dedicated simulations using this 4-stout action at T = 214 MeV. To also check for the volume dependence here we use the aspect ratio LT = 4. The result for the trace anomaly is shown in on the r.h.s of Figure 4. Lattices up to N t = 12 are used. Both with and without tree level improvement, the continuum extrapolated trace anomaly is in perfect agreement with that of the 2-stout action. This cross-check with an independent regularization provides confidence on the reliability of our continuum extrapolation. The pressure is obtained via integration from the trace anomaly. As discussed before, a crucial step is to have the correct additive normalization to satisfy the p(T = 0) = 0 condition. We used our N t = 16 dataset to improve on the normalization presented in ref. [14] in the following way.
As it has also been discussed in ref. [14] the normalized pressure is directly related to the partition function through p/T 4 = N 3 t /N 3 s log Z in the thermodynamic limit. Since log Z itself is not an observable one calculates derivatives instead, such as the trace anomaly which can then be integrated to give the pressure. There are many other possible choices for derivatives, e.g., the bare mass parameter, which is equal to the chiral condensate. In order to get the correct additive renormalization for the pressure one has to integrate from a starting point where the pressure is equal to the T = 0 pressure. We chose to integrate in the quark masses along fixed gauge couplings, selected, for each N t , such that they correspond to T * = 214 MeV at the physical point. These same gauge couplings give a temperature that is deep in the confined phase for inifinitely heavy quark masses. Therefore integrating down from sufficiently heavy quark masses along these fixed couplings one gets the correctly normalized pressure at T * = 214 MeV. The pressure at all other temperatures is normalized to this point. In Figure 5 we show the derivatives with respect to the light and strange quark masses as we vary the quark masses on a logarithmic scale. Each magenta data point represent the chiral condensate coming from two simulations, a finite temperature run and a dedicated renormalization run. Here sea and valence quark masses are kept equal. For the same pair of runs we show the strange contribution in turquoise.
We used this continuum result p(T * ) in Figure 5 as a starting point of the trace anomaly integration: The pressure is plotted in Figure 6 (left) together with the predictions of the hadron resonance gas (HRG) model at low temperatures. There is a perfect agreement with HRG in the hadronic phase. The energy and entropy densities as well as the speed of sound are shown in the right panel of Figure 6.
There is good agreement with our earlier result [14]. The only obvious difference is that in   h 0 h 1 h 2 f 0 f 1 f 2 g 1 g 2 this work 0.1396 -0.1800 0.0350 1.05 6.39 -4.72 -0.92 0.57 2010 [14] 0.1396 -0.1800 0.0350 2.76 6.79 -5.29 -0.47 1.04 Table 1: Constants for our parametrization of the trace anomaly in Eq. (2). the high temperature region (T>350 MeV), where we used only two lattice spacings previously (corresponding to N t =6 and 8) the N t =10 and 12 data changed the EoS by a few percent (the difference is on the 1.2 sigma level). This difference means that we should provide a new parametrization for the trace anomaly. The analytic function we suggest to use is of the same form: with the parameters of the fit are slightly changed. Table 1 contain the parametrization of ref. [14] and the parametrization of our present result. The maximal difference between our old parametrization [14] and the new one is only 2.8% of the Stefan-Boltzmann value for the energy density. Note that though the two results differ only on the percent level, the parameters in the new parametrization changed more (these changes merely reflect some flat directions in the parameter space).

Conclusions
We have presented a full result for the N f = 2 + 1 QCD equation of state. Our contiuum extrapolated results are completely consistent with our previous continuum estimate based on coarser lattices. The main advancement of the present work is the complete control over all systematic uncertainties. We presented a parametrization of our result which makes it easy to use in other calculations and provide our tabulated results for download.