The QCD Equation of State with almost Physical Quark Masses
M. Cheng, N. H. Christ, S. Datta, J. van der Heide, C. Jung, F. Karsch, O. Kaczmarek, E. Laermann, R. D. Mawhinney, C. Miao, P. Petreczky, K. Petrov, C. Schmidt, W. Soeldner, T. Umeda
aa r X i v : . [ h e p - l a t ] J a n BNL-NT-07/38, BI-TP 2007/20, CU-TP-1181
The QCD Equation of State with almost Physical Quark Masses
M. Cheng a , N. H. Christ a , S. Datta b , J. van der Heide c , C. Jung d , F. Karsch c , d , O. Kaczmarek c ,E. Laermann c , R. D. Mawhinney a , C. Miao c , P. Petreczky d , e , K. Petrov f , C. Schmidt d , W. Soeldner d and T. Umeda g a Physics Department,Columbia University,New York, NY 10027, USA b Department of Theoretical Physics,Tata Institute of Fundamental Research,Homi Bhabha Road, Mumbai 400005, India c Fakult¨at f¨ur Physik, Universit¨at Bielefeld,D-33615 Bielefeld, Germany d Physics Department, Brookhaven National Laboratory,Upton, NY 11973, USA e RIKEN-BNL Research Center,Brookhaven National Laboratory, Upton, NY 11973, USA f Niels Bohr Institute, University of Copenhagen,Blegdamsvej 17, DK-2100 Copenhagen, Denmark g Graduate School of Pure and Applied Sciences,University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan (Dated: October 29, 2018)We present results on the equation of state in QCD with two light quark flavors and a heavierstrange quark. Calculations with improved staggered fermions have been performed on latticeswith temporal extent N τ = 4 and 6 on a line of constant physics with almost physical quark massvalues; the pion mass is about 220 MeV, and the strange quark mass is adjusted to its physicalvalue. High statistics results on large lattices are obtained for bulk thermodynamic observables, i.e. pressure, energy and entropy density, at vanishing quark chemical potential for a wide rangeof temperatures, 140 MeV ≤ T ≤
800 MeV. We present a detailed discussion of finite cut-offeffects which become particularly significant for temperatures larger than about twice the transitiontemperature. At these high temperatures we also performed calculations of the trace anomaly onlattices with temporal extent N τ = 8. Furthermore, we have performed an extensive analysis of zerotemperature observables including the light and strange quark condensates and the static quarkpotential at zero temperature. These are used to set the temperature scale for thermodynamicobservables and to calculate renormalized observables that are sensitive to deconfinement and chiralsymmetry restoration and become order parameters in the infinite and zero quark mass limits,respectively. PACS numbers: 11.15.Ha, 11.10.Wx, 12.38Gc, 12.38.Mh
I. INTRODUCTION
Reaching a detailed understanding of bulk thermodynamics of QCD, e.g. the temperature dependence of pressureand energy density as well as the equation of state, p ( ǫ ) vs. ǫ , is one of the central goals of non-perturbative studiesof QCD on the lattice. The equation of state clearly is of central importance for the understanding of thermalproperties of any thermodynamic system. It provides direct insight into the relevant degrees of freedom and theircorrelation in different phases of strongly interacting matter. We have some understanding of the equation of state inlimiting cases of high and low temperatures from perturbation theory [1, 2, 3, 4, 5] and hadron gas phenomenology[6], respectively. In the transition region from the low temperature hadronic regime to the high temperature quarkgluon plasma, however, one has to rely on a genuine non-perturbative approach, lattice regularized QCD, to studythe non-perturbative properties of strongly interacting matter.Lattice studies of bulk thermodynamics are particularly demanding as the most interesting observables, pressureand energy density, are given in terms of differences of dimension 4 operators. These differences are particularlydifficult to evaluate because both terms being subtracted contain the pressure or energy density of the vacuum, anunphysical quantity that is approximately 1 / ( aT ) larger than the sought-after difference. Numerical signals thusrapidly decrease with the fourth power of the lattice spacing, a , when one tries to approach the continuum limit atfixed temperature ( T ). For this reason improved actions, which allow one to perform calculations on rather coarselattices with relatively small lattice discretization errors, are quite useful in thermodynamic calculations. Indeed,the early calculations of bulk thermodynamics with standard staggered [7] and Wilson [8] fermion discretizationschemes have shown that at high temperature bulk thermodynamic observables are particularly sensitive to latticediscretization errors. This closely follows observations made in studies of the thermodynamics of SU(3) gauge theories[9]. In order to minimize discretization errors at high temperature, improved staggered fermion actions - the p4-action[10] and the asqtad action [11] - have been used to study the QCD equation of state. Recent studies, performed withthe asqtad action with almost physical quark mass values on lattices with two different values of the lattice cut-off[11], indeed show much smaller discretization errors than similar studies performed with the 1-link, stout smearedstaggered fermion action [12]. Another source for cut-off errors arises, however, from the explicit breaking of flavorsymmetry in the staggered fermion formulation. While this is not of much concern in the chirally symmetric hightemperature phase of QCD, it leads to cut-off dependent modifications of the hadron spectrum and thus may influencethe calculation of thermodynamic observables in the low temperature hadronic phase of QCD. Techniques to reduceflavor symmetry breaking through the introduction of so-called ’fat links’ are thus generally exploited in numericalcalculations with staggered fermions [14].In this paper we report on a calculation of bulk thermodynamics in QCD with almost physical light quark massesand a physical value of the strange quark mass. Our calculations have been performed with a tree level Symanzik-improved gauge action and an improved staggered fermion action, the p4-action with 3-link smearing (p4fat3), whichremoves O ( a ) cut-off effects at tree-level and also leads to small cut-off effects in O ( g ) perturbation theory [13].At each temperature, we perform simulations with two degenerate light quark masses and a heavier strange quarkmass for two different values of the lattice cut-off, corresponding to lattices with temporal extent N τ = 4 and 6.In these calculations we explore a wide range of temperatures varying from about 140 MeV to about 800 MeV.This corresponds to the temperature interval relevant for current experimental studies of dense matter in heavy ioncollisions at RHIC as well as the forthcoming experiments at the LHC. Bare quark masses have been adjusted tokeep physical masses approximately constant when the lattice-cut off is varied. At high temperatures, T > ∼
350 MeV,we also performed calculations on lattices with temporal extent N τ = 8 to get control over cut-off effects in the hightemperature limit.We will start in the next section by reviewing basic thermodynamic relations in the continuum valid for thermody-namic calculations on such lines of constant physics (LCP). In section III we outline details of our calculational set-upwith improved staggered fermions. In section IV we present our zero temperature calculations needed to define theline of constant physics and the temperature scale deduced from properties of the static quark-antiquark potential.Section V is devoted to the presentation of our basic result, the difference between energy density and three times thepressure from which we obtain all other thermodynamic observables, e.g. the pressure, energy and entropy densitiesas well as the velocity of sound. Section VI is devoted to a discussion of the temperature dependence of Polyakov loopexpectation values and chiral condensates which provides a comparison between the deconfining and chiral symmetryrestoring features of the QCD transition. We finally present a discussion of our results and a comparison with otherimproved staggered fermion calculations of bulk thermodynamics in Section VII. II. THERMODYNAMICS ON LINES OF CONSTANT PHYSICS
To start our discussion of QCD thermodynamics on the lattice we recall some basic thermodynamic relations in thecontinuum. For large, homogeneous media the basic bulk thermodynamic observables we will consider here can beobtained from the grand canonical partition function with vanishing quark chemical potentials, Z ( T, V ). We introducethe grand canonical potential, Ω(
T, V ), normalized such that it vanishes at vanishing temperature,Ω(
T, V ) = T ln Z ( T, V ) − Ω , (1)with Ω = lim T → T ln Z ( T, V ). With this we obtain the thermal part of the pressure ( p ) and energy density ( ǫ ) p = 1 V Ω( T, V ) , ǫ = T V ∂ Ω( T, V ) /T∂T , (2)which by construction both vanish at vanishing temperature. Using these relations one can express the differencebetween ǫ and 3 p , i.e the thermal contribution to the trace of the energy-momentum tensor Θ µµ ( T ), in terms of aderivative of the pressure with respect to temperature, i.e. Θ µµ ( T ) T ≡ ǫ − pT = T ∂∂T ( p/T ) , (3)In fact, it is Θ µµ ( T ) which is the basic thermodynamic quantity conveniently calculated on the lattice. All otherbulk thermodynamic observables, e.g. p/T , ǫ/T as well as the entropy density, s/T ≡ ( ǫ + p ) /T , can be deducedfrom this using the above thermodynamic relations. In particular, we obtain the pressure from Θ µµ ( T ) throughintegration of Eq. 3, p ( T ) T − p ( T ) T = Z TT d T ′ T ′ Θ µµ ( T ′ ) . (4)Usually, the temperature for the lower integration limit, T , is chosen to be a temperature sufficiently deep in thehadronic phase of QCD where the pressure p ( T ), receives contributions only from massive hadronic states and isalready exponentially small. We will discuss this in more detail in Section V. Eq. 4 then directly gives the pressureat temperature T . Using p/T determined from Eq. 4 and combining it with Eq. 3, we obtain ǫ/T as well as s/T .This makes it evident, that there is indeed only one independent bulk thermodynamic observable calculated in thethermodynamic (large volume) limit. All other observables are derived through standard thermodynamic relations sothat thermodynamic consistency of all bulk thermodynamic observables is insured by construction!We stress that the normalization introduced here for the grand canonical potential, Eq. 1, forces the pressure andenergy density to vanish at T = 0. As a consequence of this normalization, any non-perturbative structure of theQCD vacuum, e.g. quark and gluon condensates, that contribute to the trace anomaly Θ µµ (0), and would lead toa non-vanishing vacuum pressure and/or energy density, eventually will show up as non-perturbative contributionsto the high temperature part of these thermodynamic observables. This is similar to the normalization used, e.g. inthe bag model and the hadron resonance gas, but differs from the normalization used e.g. in resummed perturbativecalculation at high temperature [15, 16] or phenomenological (quasi-particle) models for the high temperature phase ofQCD [17]. This should be kept in mind when comparing results for the EoS with perturbative and model calculations.We also note that ambiguities in normalizing pressure and energy density at zero temperature drop out in a calculationof the entropy density which thus is the preferred observable for such comparisons. III. LATTICE FORMULATION
In a lattice calculation, temperature and volume are given in terms of the temporal ( N τ ) and spatial ( N σ ) latticeextent as well as the lattice spacing, a , which is controlled through the lattice gauge coupling β ≡ /g , T = 1 N τ a ( β ) , V = ( N σ a ( β )) . (5)As all observables that are calculated on the lattice, are functions of the coupling, β , we may rewrite Eq. 3 in termsof a derivative taken with respect to β rather than T . Furthermore we adopt the normalization of the pressure asintroduced in Eq. 1. This insures a proper renormalization of thermodynamic quantities and, as a consequence, forcesthe pressure to vanish in the vacuum, i.e. at T = 0.Let us write the QCD partition function on a lattice of size N σ N τ as Z LCP ( β, N σ , N τ ) = Z Y x,µ d U x,µ e − S ( U ) , (6)where U x,µ ∈ SU (3) denotes the gauge link variables and S ( U ) = βS G ( U ) − S F ( U, β ) is the Euclidean action, which iscomposed out of a purely gluonic contribution, S G ( U ), and the fermionic part, S F ( U, β ), which arises after integrationover the fermion fields. We will specify this action in more detail in the next section but note here that we will usetree level improved gauge and fermion actions. Although it would be straightforward to introduce one-loop or tadpoleimprovement factors in the action the setup used here greatly simplifies the analysis of thermodynamic observablesand in some cases also gives a more direct relation to corresponding observables in the continuum.When using only tree level improvement the gluonic action does not depend on the gauge coupling, β , and thefermion action depends on β only through the bare light ( ˆ m l ) and strange ( ˆ m s ) quark masses. The subscript LCP in Eq. 6 indicates that we have defined the partition function Z LCP on a line of constant physics (LCP), i.e. whenapproaching the continuum limit by increasing the gauge coupling ( β → ∞ ) the bare quark masses ( ˆ m l ( β ) , ˆ m s ( β ))in the QCD Lagrangian are tuned towards zero such that the vacuum properties of QCD remain unchanged. Thequark masses thus are not independent parameters but are functions of β which are determined through constraintsimposed on zero temperature observables; e.g. one demands that a set of hadron masses remains unchanged whenthe continuum limit is approached on a LCP.We now may rewrite Eq. 3 in terms of observables calculable in lattice calculations at zero and non-zero temperature,Θ µµ ( T ) T = − R β ( β ) N τ (cid:18) N σ N τ (cid:28) d S d β (cid:29) τ − N σ N (cid:28) d S d β (cid:29) (cid:19) . (7)Here h ... i x , with x = τ, N σ N τ , with N τ ≪ N σ , and zero temperature lattices, i.e. on lattices with large temporal extent, N σ N τ with N τ ≡ N > ∼ N σ ,respectively. Furthermore, R β denotes the lattice version of the QCD β -function which arises as a multiplicativefactor in the definition of Θ µµ ( T ) because derivatives with respect to T have been converted to derivatives withrespect to the lattice spacing a on lattices with fixed temporal extent N τ , R β ( β ) ≡ T d β d T = − a d β d a . (8)We note that in the weak coupling, large β limit, R β approaches the universal form of the 2-loop β -function of3-flavor QCD, R β ( β ) = 12 b + 72 b /β + O ( β − ) , (9)with b = 9 / π and b = 1 / π .We analyze the thermodynamics of QCD with two degenerate light quarks ( ˆ m l ≡ ˆ m u = ˆ m d ) and a heavier strangequark ( ˆ m s ) described by the QCD partition function given in Eq. 6. For our studies of bulk thermodynamics we usethe same discretization scheme which has been used recently by us in the study of the QCD transition temperature[22], i.e. we use an O ( a ) tree level improved gauge action constructed from a 4-link plaquette term and a planar6-link Wilson loop as well as a staggered fermion action that contains a smeared 1-link term and bent 3-link terms.We call this action the p4fat3-action; further details are given in Ref. [10] where the p4fat3 action was first used instudies of the QCD equation of state on lattices with temporal extent N τ = 4 and larger quark masses. With thisaction, bulk thermodynamic quantities like pressure and energy density are O ( a ) tree level improved; corrections tothe high temperature ideal gas limit only start at O (1 /N τ ) and are significantly smaller than for the Naik action or thestandard staggered action which suffers from large O (1 /N τ ) cut-off effects at high temperature. An analysis of cut-offeffects in the ideal gas limit and in O ( g ) lattice perturbation theory [13] shows that deviations from perturbativeresults are already only a few percent for lattices with temporal extent as small as N τ = 6.Following the notation used in Ref. [22] the Euclidean action is given as S ( U ) = βS G ( U ) − S F ( U, β ) , (10)with a gluonic contribution, S G ( U ), and a fermionic part, S F ( U, β ). The latter can be expressed in terms of thestaggered fermion matrices, D ˆ m l ( ˆ m s ) , for two light ( ˆ m l ) and a heavier strange quark ( ˆ m s ), S F ( U, β ) = 12 Tr ln D ( ˆ m l ( β )) + 14 Tr ln D ( ˆ m s ( β )) . (11)Here we took the fourth root of the staggered fermion determinant to represent the contribution of a single fermionflavor to the QCD partition function .We also introduce the light and strange quark condensates calculated at finite ( x = τ ) and zero temperature ( x = 0),respectively, h ¯ ψψ i q,x ≡
14 1 N σ N x (cid:10) Tr D − ( ˆ m q ) (cid:11) x , q = l, s , x = 0 , τ , (12)as well as expectation values of the gluonic action density, h s G i x ≡ N σ N x h S G i x . (13)All numerical calculations have been performed using the Rational Hybrid Monte Carlo (RHMC) algorithm [21] withparameters that have been optimized [22] to reach acceptance rates of about 70%. Some details on our tuning ofparameters of the RHMC algorithm have been given in [23].For the discussion of the thermodynamics on a line of constant physics (LCP) it sometimes is convenient toparametrize the quark mass dependence of S F in terms of the light quark mass ˆ m l and the ratio h ≡ ˆ m s / ˆ m l ratherthan ˆ m l and ˆ m s separately. We thus write the β -dependence of the strange quark mass as, ˆ m s ( β ) = ˆ m l ( β ) h ( β ). In There is a controversy regarding the validity of the rooting approach in numerical calculations with staggered fermions. For furtherdetails we refer to recent reviews presented at Lattice conference [18, 19, 20] and references therein. the evaluation of ( ǫ − p ) /T we then will need to know the derivatives of these parametrizations with respect to β .We define R m ( β ) = 1ˆ m l ( β ) d ˆ m l ( β )d β , R h ( β ) = 1 h ( β ) d h ( β )d β . (14)With these definitions we may rewrite Eq. 7 as ǫ − pT = T dd T (cid:16) pT (cid:17) = R β ( β ) ∂p/T ∂β = Θ µµG ( T ) T + Θ µµF ( T ) T + Θ µµh ( T ) T , (15)with Θ µµG ( T ) T = R β [ h s G i − h s G i τ ] N τ , (16)Θ µµF ( T ) T = − R β R m (cid:2) m l (cid:0) h ¯ ψψ i l, − h ¯ ψψ i l,τ (cid:1) + ˆ m s (cid:0) h ¯ ψψ i s, − h ¯ ψψ i s,τ (cid:1)(cid:3) N τ , (17)Θ µµh ( T ) T = − R β R h ˆ m s (cid:2) h ¯ ψψ i s, − h ¯ ψψ i s,τ (cid:3) N τ . (18)We will show in the next section that to a good approximation h ( β ) stays constant on a LCP. R h thus vanishes onthe LCP and consequently the last term in Eq. 15, Θ µµh , will not contribute to the thermal part of the trace anomaly,Θ µµ ( T ). The other two terms stay finite in the continuum limit and correspond to the contribution of the thermalparts of gluon and quark condensates to the trace anomaly. We note that the latter contribution vanishes in thechiral limit of three flavor QCD ( ˆ m l , ˆ m s → µµG ( T ) andthe observables entering the calculation of bulk thermodynamic quantities in the chiral limit of QCD would reduceto those needed also in a pure SU (3) gauge theory [9]. In fact, we also find that for physical values of the quarkmasses the trace anomaly is dominated by the gluonic contribution, Θ µµG ( T ). As will become clear in Section VΘ µµF ( T ) contributes less than 10% to the total trace anomaly for temperatures large than about twice the transitiontemperature.We also note that the prefactor in Eq. 17 will approach unity in the continuum limit as R m attains a universalform up to 2-loop level which is similar to that of R − β and is only modified through the anomalous dimension of thequark mass renormalization [24]. For the relevant combination of β -functions that enters the fermionic part of thetrace anomaly, one has − ( R β ( β ) R m ( β )) − loop = 1 + 16 b β . (19) IV. STATIC QUARK POTENTIAL AND THE LINE OF CONSTANT PHYSICSA. Construction of the line of constant physics
We will calculate thermodynamic observables on a line of constant physics (LCP) that is defined at T = 0 asa line in the space of light and strange bare quark masses parametrized by the gauge coupling β . Each point onthis line corresponds to identical physical conditions at different values of the lattice cut-off which is tuned towardsthe continuum limit by increasing the gauge coupling β . We define the line of constant physics by demanding (i)that the ratio of masses for the strange pseudo-scalar and the kaon mass, m ¯ ss /m K , stays constant and (ii) that m ¯ ss expressed in units of the scale parameter r stays constant. The latter gives the distance at which the slope of the zerotemperature, static quark potential, V ¯ qq ( r ), attains a certain value. We also introduce the scale r , which frequentlyis used on finer lattices to convert lattice results expressed in units of the cut-off to physical scales, (cid:18) r d V ¯ qq ( r )d r (cid:19) r = r = 1 . , (cid:18) r d V ¯ qq ( r )d r (cid:19) r = r = 1 . . (20)We checked that (i) and (ii) also hold true, if we replace m ¯ ss by the mass of the light quark pseudo-scalar meson, m π . However, errors on m π r and ( m π /m K ) are generally larger which, in particular at large values of β makes theparametrization of the LCP less stringent.Leading order chiral perturbation theory suggests that the ratio ( m ¯ ss /m K ) is proportional to ˆ m s / ( ˆ m l + ˆ m s ). Onethus expects this ratio to stay constant for fixed h = ˆ m s / ˆ m l . This is, indeed fulfilled in the entire regime of couplings, β , explored in our calculations (see Table I). The first condition for fixing the LCP parameters thus, in practice, hasbeen replaced by choosing h = ˆ m l / ˆ m s to be constant. As a consequence we find R h ( β ) = 0, which simplifies thecalculation of thermodynamic quantities.In order to define a line of constant strange quark mass, as a second condition for the LCP we demand that theproduct m ¯ ss r stays constant. For our LCP we chose 1.59 as the value for the product. Here one should notethat m ¯ ss determined in our calculations only receives contributions from connected diagrams and does not includedisconnected loops. In order to compare our value (1.59, see discussion below) to a physical one, we therefore followthe argumentation of Ref. [25] and adopt m ¯ ss = p m K − m π = 686 MeV as the physical mass of our strangepseudoscalar. Together with the scale r = 0 . r withlevel splittings of the charmonium system [27], this yields m ¯ ss r ≃ .
63. Of course, there is some ambiguity in thischoice as current determinations of r differ by about 10% [26, 28]. This introduces some systematic error in thedefinition of the physical LCP . The main reason for deviation from the physical LCP in the present calculation,however, is due to the choice of the light quark masses which are about a factor two too large.Fixing the light and strange pseudo-scalar masses in units of r required some trial runs for several β values. Wethen used the leading order chiral perturbation theory ansatz m ss ∼ ˆ m s (or m π ∼ ˆ m l ) to choose ˆ m s and ˆ m l ≡ ˆ m s / β values at which high statistics simulations have been performed. It turned out thatthese values are best fitted by m ¯ ss r = 1 .
59. We thus use this value rather than the value 1 .
63 mentioned above, todefine our LCP. For all other simulations we then used the results of these zero temperature calculations to determinethe quark mass values that belong to a line of constant physics characterized by:LCP : ( i ) m ¯ ss r = 1 . , ( ii ) h ≡ ˆ m s / ˆ m l = 10In general our calculations are thus performed at parameter values close to the LCP which is defined by the abovecondition. The parameters of all our zero temperature calculations performed to determine the LCP, results formeson masses and parameters of the static quark potential are summarized in Table I. As can be seen, at ouractual simulation points the results for m ¯ ss r fluctuate around the mean value by a few percent. We also checkedthe sensitivity of the meson masses used to determine the LCP to finite volume effects. At β = 3 .
49 and 3 .
54 weperformed calculations on 32 lattices in addition to the 16 ·
32 lattices. As can bee seen from Table I results for m ¯ ss and m K agree within statistical errors and volume effects are at most on the level of 2% for the light pseudo-scalar.The LCP is furthermore characterized by m π /m K = 0 . m ¯ ss /m K = 1 . r = 0 . m π ≃ m ¯ ss ≃ m K ≃ B. The static quark potential and the scale r On the LCP we determine several parameters, e.g. the short distance scale r and the linear slope parameter, thestring tension σ , that characterize the shape of the static quark potential calculated at T = 0 in a fixed range ofphysical distances. The distance r , defined in Eq. 20, is used to define the temperature scale for the thermodynamicscalculations.The static quark potential, V ¯ qq ( r ), has been calculated from smeared Wilson loops as described in [22] for allparameter sets listed in Table I. We checked that the the smeared Wilson loops project well onto the ground stateat all values of the cut-off by verifying the independence of the extracted potential parameters on the number ofsmearing levels used in the analysis. The set of gauge couplings, β ∈ [3 . , . a ≃ . a ≃ .
05 fm. When analyzingthe static potential over such a wide range of cut-off values one should make sure that the potential is analyzed inapproximately the same range of physical distances. The fit interval [( r/a ) min , ( r/a ) max ] for fits with a Cornell typeansatz for the static potential thus has been adjusted for the different values of gauge couplings such that it coversapproximately the same range of physical distances, r / < ∼ r < ∼ r , or 0 . < ∼ r< ∼ r also independently by using spline interpolations which arenot biased by a particular ansatz for the form of the potential.The left hand part of Fig. 1 shows the static quark potential for several of our parameter sets. We have renormalized β
100 ˆ m l N σ · N τ m π a m ¯ ss a m K a r /a √ σa c ( g ) r ·
32 0.3410( 2) 1.0474( 1) 0.7854( 2) 1.467(72) 0.75(18) 0.97(12)3.210 1.000 16 ·
32 0.3262( 1) 0.9988( 1) 0.7496( 1) 1.583(36) 0.685(75) 1.118(68)3.240 0.900 16 ·
32 0.3099( 2) 0.9485( 2) 0.7125( 3) 1.669(31) 0.658(36) 1.243(67)3.277 0.765 16 ·
32 0.2881( 7) 0.8769( 5) 0.6599( 6) 1.797(19) 0.612(53) 1.362(53)3.290 0.650 16 ·
32 0.2667( 8) 0.8104( 7) 0.6101( 8) 1.823(16) 0.623(32) 1.362(29)3.335 0.620 16 ·
32 0.2594( 3) 0.7884( 2) 0.5941( 5) 1.995(11) 0.5668(73) 1.504(22)3.351 0.591 16 ·
32 0.2541( 7) 0.7692( 5) 0.5800( 7) 2.069(12) 0.551(11) 1.594(24)3.382 0.520 16 ·
32 0.2370( 6) 0.7194( 5) 0.5422( 5) 2.230(14) 0.5100(82) 1.718(57)3.410 0.412 16 ·
32 0.2098( 4) 0.6371( 6) 0.4796( 8) 2.503(18) 0.440(10) 2.073(49)3.420 0.390 24 ·
32 0.2029( 8) 0.6177( 5) 0.4675( 5) 2.577(11) 0.4313(56) 2.124(33)3.430 0.370 24 ·
32 0.1986( 6) 0.6000( 3) 0.4529( 5) 2.6467(81) 0.4225(53) 2.178(17)3.445 0.344 24 ·
32 0.1909( 7) 0.5749( 4) 0.4335( 5) 2.813(15) 0.3951(68) 2.388(35)3.455 0.329 24 ·
32 0.1833(10) 0.5580( 6) 0.4204( 8) 2.856(20) 0.3895(68) 2.375(42)3.460 0.313 16 ·
32 0.1808(16) 0.5443(11) 0.4102(11) 2.890(16) 0.3831(84) 2.391(55)3.470 0.295 24 ·
32 0.1686(19) 0.5233( 8) 0.3940(12) 3.065(18) 0.3592(75) 2.617(41)3.490 0.290 16 ·
32 0.1689(14) 0.5115(11) 0.3842(11) 3.223(31) 0.3423(66) 2.757(59)3.490 0.290 32 ·
32 0.1525(40) 0.4740(20) 0.3554(22) 3.423(61) 0.322(14) 2.934(92)3.540 0.240 16 ·
32 0.1495(24) 0.4458(20) 0.3358(19) 3.687(34) 0.3011(46) 3.128(51)3.540 0.240 32 ·
32 0.1347(53) 0.4053(18) 0.3028(23) 4.009(26) 0.2743(38) 3.414(47)3.630 0.170 24 ·
32 0.1126(20) 0.3386( 7) 0.2537( 8) 4.651(41) 0.2352(44) 3.939(59)3.690 0.150 24 ·
32 0.1020(90) 0.2960(20) 0.2230(30) 5.201(48) 0.2116(36) 4.320(63)3.760 0.130 24 · ·
48 0.0857(32) 0.2530(16) 0.1894(16) 6.050(61) 0.1810(29) 4.984(73)3.820 0.125 32 ·
32 0.0830(40) 0.2310(38) 0.1744(50) 6.835(44) 0.1701(21) 5.541(106)3.920 0.110 32 ·
32 0.0750(70) 0.2020(10) 0.1550(20) 7.814(83) 0.1423(24) 6.037(72)4.080 0.081 32 ·
32 0.0700(70) 0.1567(36) 0.1220(50) 10.39(23) 0.1060(35) 7.710(183)TABLE I: Light quark and strange pseudo-scalar meson masses and parameters of the static quark potential calculated on zerotemperature lattices of size N σ N τ . The last column gives the renormalization constants times r needed to renormalize theheavy quark potential to the string potential at distance r/r = 1 . these potentials by matching all potentials at a large distance, r/r = 1 .
5, to a common value that is taken to beidentical to the large distance string potential, V string ( r ) = − π/ r + σr . The result of this matching is shown in thelower part of Fig. 1(left) and the constant shifts needed to obtain these renormalized potentials are listed in Table I.The good matching of all the potential data obtained at different values of the cut-off already gives a good idea of thesmallness of finite cut-off effects in this observable. We note that this matching procedure provides renormalizationconstants for the static quark potential, which we also will use later to renormalize the Polyakov loop expectationvalue.To further analyze the shape of the static quark potential we determined the scale parameter r /a as well as thesquare root of the string tension in lattice units, √ σa . These parameters have been obtained from three and fourparameter fits. As described in [22] the latter fit ansatz has been used to estimate systematic errors in our analysisof the scale parameters.Results for r /a and √ σa are given in Table I. We note that the product r √ σ stays constant on the LCP andchanges by less than 2% in the entire range of couplings β in which the lattice cut-off changes by a factor 6. For a ≤ .
15 fm we used a quadratic fit ansatz, ( r √ σ ) a = r √ σ + c ( a/r ) , to fit 10 data points. The asymptotic valuefor r √ σ coincides within errors with a simple average over all values of ( r √ σ ) a in this interval. This confirms that O ( a ) corrections indeed are small for this product. Similarly we determined the scale parameter r frequently used Further details on the matching of the zero temperature heavy quark potentials, its application to the renormalization of the Polyakovloop and finite temperature free energies will be published elsewhere. -2 0 2 4 6 8 0 0.5 1 1.5 2 2.5 3 r/r r V qq (r) - string a/r a [fm] r /r r σ FIG. 1: The static quark potential in units of the scale r versus distance r/r (left) and dimensionless combinations of thepotential shape parameters r /r and r √ σ extracted from fits to these potentials (right). The left hand figure shows potentialsfor several values of β taken from our entire simulation interval, β ∈ [3 .
15 : 4 . r = 0 .
469 fm. to set the scale in calculations performed on finer lattices. Both fits for r √ σ and r /r yield χ /dof ≃ .
7. Fromthis analysis we obtain the parameters characterizing the shape of the heavy quark potential at masses in the vicinityof the LCP, r √ σ = 1 . ,r /r = 1 . . (21)We note that the result obtained here for r /r is in good agreement with the corresponding continuum extrapolatedvalue, r /r = 1 . a ≃ .
12 fm and a ≃ .
09 fm, respectively [36]. We show resultsfor r /r and r √ σ calculated at parameter sets close to the LCP in Fig. 1 (right).Despite the good scaling behavior of dimensionless combinations of scale parameters deduced from the staticpotential, one expects, of course, to still find substantial deviations from asymptotic scaling relations that are controlledby universal 2-loop β -functions. For the scale parameter r /a we parametrize deviations from asymptotic scaling usinga rational function ansatz, ˆ r ≡ r a = 1 + e r ˆ a ( β ) + f r ˆ a ( β ) a r R ( β ) (1 + b r ˆ a ( β ) + c r ˆ a ( β ) + d r ˆ a ( β )) , (22)where R ( β ) = exp (cid:18) − β b (cid:19) (cid:18) b β (cid:19) − b / (2 b ) (23)denotes the 2-loop β -function of QCD for three massless quark flavors and ˆ a ( β ) = R ( β ) /R (3 . β -function R β entering all basic thermodynamic observables, R β ( β ) = r a (cid:18) d r /a d β (cid:19) − . (24)Furthermore, we need a parametrization of the β -dependence of the bare quark masses to determine the second β -function entering the thermodynamic relations, i.e. R m ( β ) defined in Eq. 14. For this purpose we use a parametriza-tion of the product of the bare light quark mass, ˆ m l and ˆ r that takes into account the anomalous scaling dimensionof quark masses [24], ˆ m l ˆ r = a m (cid:18) b β (cid:19) / P ( β ) , (25) b m c m d m e m f m g m -2.149(121) 1.676(178) -0.365(144) -2.290(162) 1.829(425) -0.356(335) a r b r c r d r e r f r r in lattice units based on the ansatz given in Eq. 22 (lower half) andthe fit of the renormalization group invariant combination of light quark masses and r (Eq. 25) on the line of constant physics(upper half). The χ /d.o.f for these fits are 1 . .
84 for 18 degrees of freedom (upperhalf). r /a β RG fit 0.0100.0150.0200.0250.0300.0350.040 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 m l r ( β /12b ) β FIG. 2: The scale parameter ˆ r ≡ r /a versus β = 6 /g (left) and its product with the bare light quark mass on the LCP(right). The two curves shown in the left hand part of this figure correspond to two different fit ans¨atze. As explained in thetext in addition to the renormalization group motivated ansatz given in Eq. 22 the result from a 3-interval fit is shown. Thecurve in the right hand part of the figure shows a fit based on the ansatz given in Eqs. 25 and 26. with a m being related to the renormalization group invariant quark mass in units of r and P ( β ) being a sixth orderrational function that parametrizes deviations from the leading order scaling relation for the bare quark mass, P ( β ) = 1 + b m ˆ a ( β ) + c m ˆ a ( β ) + d m ˆ a ( β )1 + e m ˆ a ( β ) + f m ˆ a ( β ) + g m ˆ a ( β ) . (26)This ansatz insures that the parametrization for the two β -functions as well as the parametrization of their product, R β ( β ) R m ( β ), reproduces the universal 2-loop results given in Eqs. 9 and 19.In Fig. 2 we show our results for ˆ r = r /a and ˆ m l ˆ r together with the fits described above. The fit parametersdefining the quark masses on the LCP have been obtained from χ -fits in the interval β ∈ [3 . , . a m = 0 . . r /a differ from the actually calculated values given in Table II by less than one percent.Like in the pure gauge theory calculations of the equation of state, we also find for QCD with light dynamical quarksthat, in the parameter range of interest for finite temperature calculations, β -functions deviate significantly from theasymptotic scaling form. In particular, we find a dip in R β at β ≃ .
43. For small values of N τ , the interestingparameter range thus includes the crossover region from the strong to weak coupling regime.We use the interpolating fits for ˆ r and ˆ m l ˆ r , to determine the two β -functions R β and R m that enter the calculationsof thermodynamic quantities. As all basic thermodynamic observables are directly proportional to R β , we should checkthe sensitivity of R β on the particular interpolation form used. We thus have used a completely different interpolationthat restricts the renormalization group motivated ansatz to the small coupling regime, β ≥ .
52, and uses purelyrational functions to piecewise fit 2 intervals at smaller β . We find that results for R β are sensitive to the fit ansatzonly for small β -values, i.e β< ∼ .
25, where the dependence of ˆ r on β becomes weak. As discussed later the uncertaintyon R β at small values of the coupling only affects the three smallest temperatures used for the analysis of the equationof state on the N τ = 4 lattices.0 β R β β -R β R m FIG. 3: The β -function on the LCP (Eq. 24) (left) and the product R β R m (right). The horizontal lines show the weak couplingbehavior given in Eqs. 9 and 19. The two curves result from two different fits of ˆ r as discussed in the text. Using the parametrizations of ˆ r and ˆ m l ˆ r given in Eq. 22 and Eq. 25 as well as the above discussed piecewiseinterpolation of ˆ r we now can derive the two β -functions R β ( β ) and R m ( β ). In Fig. 3 we show R β as well as thecombination − R β R m which enter the calculation of the gluonic and fermionic contributions to ( ǫ − p ) /T . For ˆ r as well as for the two β -functions, we show result obtained with our two different fit ans¨atze. As can be seen, thedifferent fit forms lead to differences in the fit result at the edges of the parameter range analyzed. We take care of thisin our analysis of the equation of state by averaging over the results obtained with the two different fit ans¨atze andby including the difference of both fit results as a systematic error. We note that the β -function R β has a minimumat β ≃ .
43. This characterizes the transition from strong to weak coupling regions and is similar to what is knownfrom β -functions determined in pure gauge theory [9] as well as in QCD with heavier quark masses [10]. The detailsof this region will differ in different discretization schemes as the QCD β -functions are universal only up to 2-looporder in perturbation theory. In order to understand the origin of cut-off effects in thermodynamic observables it is,however, important to have good control over R β in this non-universal regime as well, as R β enters the calculation ofall relevant lattice observables as an overall multiplicative factor . V. BULK THERMODYNAMICSA. The trace anomaly: ( ǫ − p ) /T The basic lattice observables needed to determine the QCD equation of state with our tree level improved gauge andfermion actions are expectation values of the gauge action as well as the light and strange quark chiral condensatescalculated on the LCP on finite ( N τ ≪ N σ ) and zero ( N τ > ∼ N σ ) temperature lattices. We have performed finitetemperature calculations on lattices with temporal extent N τ = 4, 6 and 8. In all cases the spatial extent of thelattices ( N σ ) was at least four times larger than the temporal extent ( N τ ), i.e. most finite temperature calculationshave been performed on lattices of size 16
6, respectively. In particular at high temperature, we found itimportant to increase the spatial volume in our calculations on N τ = 6 lattices to check for possible finite volumeeffects and also to add a few calculations on N τ = 8 lattices to get control over the cut-off dependence seen in thetrace anomaly. In these cases, calculations on 32
32 and 24
32. In a few caseswe used lattices of size 24 · ·
48 as well as 32 . The length of individual calculations on the finite temperaturelattices varied between 6500 and 35000 trajectories on the N τ = 4 lattices and 5000 to 17600 iterations on the N τ = 6lattices, where Metropolis updates were done after hybrid Monte Carlo evolutions of trajctory length τ MD = 0 . We note that in order to insure thermodynamic consistency the β -function used in the definition of thermodynamic quantities has tobe determined from the cut-off dependence of the observable used to set the temperature scale, i.e. r /a in our study. T > ǫ − p ) /T , of below 20% at all temperatures. In fact, they are below 10% in thetemperature interval T ∈ [180MeV , T ∈ [195MeV , T [MeV] Tr ( ε -3p)/T p4: N τ =468 FIG. 4: The trace anomaly Θ µµ ( T ) ≡ ǫ − p in units of T versus temperature obtained from calculations on lattices withtemporal extent N τ = 4, 6, and 8. The temperature scale, T r (upper x-axis) has been obtained using the parametrizationgiven in Eq. 22, and T [MeV] (lower x-axis), has been extracted from this using r = 0 .
469 fm.
The basic zero and finite temperature observables needed to calculate the trace anomaly in units of the fourthpower of the temperature, Θ µµ ( T ) /T = ( ǫ − p ) /T , from Eq. 7 are summarized in Tables III, IV, V and VI. Toextract Θ µµ ( T ) /T one furthermore needs to know the derivatives of bare couplings and quark masses, R β and R m .Their calculation from zero temperature observables has been discussed in the previous section. With this input, weobtain the result for Θ µµ ( T ) /T , shown in Fig. 4 for the entire range of temperatures explored by us. Here, and inall subsequent figures, the temperature scale has been determined from our results for r /a , which characterizes theslope of the static quark potential and has been extracted from the zero temperature potential as discussed in theprevious section. On lattices with temporal extent N τ we then have T r ≡ ˆ r /N τ . Whenever we show in the followingtemperatures in units of MeV we use r = 0 .
469 fm [26] to convert
T r to a MeV-scale. We will, however, show in allfigures both scales which should allow us to compare the results presented here unambiguously with any other latticecalculation performed within a different regularization scheme.In QCD with light (u, d)-quarks and a heavier strange quark the trace anomaly receives, in addition to the gluoniccontribution to the trace over the energy-momentum tensor, also contributions from the light and heavy quark chiralcondensates (Eqs. 16, 17). In the chiral limit only the former contributes and all fermionic contributions enterindirectly through modifications of the gauge field background. It thus is interesting to check the relative importanceof direct contributions from the chiral condensates to ( ǫ − p ) /T . In Fig. 5 we show the fermion contribution Θ µµF /T to the total trace anomaly shown in Fig. 4. The right hand part of this figure shows the relative magnitude of thelight and strange quark contributions. As can be seen they are of similar size close to the transition temperature.With increasing temperature, however, the importance of the light quark contribution rapidly drops and becomessimilar to the ratio of light to strange quark masses at about twice the transition temperature. As can be seen inFig. 5(left) the total fermionic contribution shows a significant cut-off dependence. This partly arises from the largechange of the product of β -functions, R β R m that still deviates a lot from the asymptotic weak coupling value in therange of couplings relevant for the N τ = 4 and 6 calculations, respectively (see Fig. 3(right)). The influence of thiscut-off dependence on the calculation of the total trace anomaly, however, is strongly reduced as the contribution ofΘ µµF /T only amounts to about 20% in the transition region and already drops below 10% at about 1 . T c .As all other thermodynamic observables will eventually be deduced from ( ǫ − p ) /T using standard thermodynamicrelations, we should analyze its structure carefully. Bulk thermodynamics of QCD in different temperature intervalsaddresses quite different physics. This includes (i) the low temperature regime, which in the vicinity of the transitiontemperature often is compared with the physics of a resonance gas and which at lower temperatures is sensitive toproperties of the hadron spectrum controlled by chiral symmetry breaking; (ii) the genuine non-perturbative physicsin the transition region and at temperatures above but close to the crossover region which is probed experimentally atRHIC and presumably is a still strongly interacting medium with a complicated quasi-particle structure; and (iii) thehigh temperature regime, which eventually becomes accessible to resummed perturbative calculations. In numerical2 β
100 ˆ m l N σ · N τ h s G i h ¯ ψψ i l, h ¯ ψψ i s, ·
32 4544 4.82564(21) 0.28727(11) 0.392677(53)3.210 1.000 16 ·
32 5333 4.68944(27) 0.25284(14) 0.358813(80)3.240 0.900 16 ·
32 5110 4.61441(29) 0.23156(16) 0.333957(88)3.277 0.765 16 ·
32 3408 4.51660(41) 0.20232(17) 0.29834(12)3.290 0.650 16 ·
32 3067 4.47696(37) 0.18807(19) 0.27506(14)3.335 0.620 16 ·
32 3689 4.36044(25) 0.15429(17) 0.24425(10)3.351 0.591 16 ·
32 7005 4.31880(34) 0.14175(20) 0.23045(13)3.382 0.520 16 ·
32 5051 4.23499(26) 0.11515(14) 0.19922(11)3.410 0.412 16 ·
32 5824 4.15990(43) 0.09013(27) 0.16256(20)3.420 0.390 24 ·
32 2448 4.13616(20) 0.08303(17) 0.15304(12)3.430 0.370 24 ·
32 1849 4.11217(29) 0.07606(15) 0.14364(11)3.445 0.344 24 ·
32 1707 4.07770(23) 0.06650(10) 0.130718(86)3.455 0.329 24 ·
32 2453 4.05605(36) 0.06098(24) 0.12314(18)3.460 0.313 16 ·
32 2513 4.04471(35) 0.05733(25) 0.11734(17)3.470 0.295 24 ·
32 3079 4.02346(18) 0.05237(10) 0.109388(88)3.490 0.290 16 ·
32 4300 3.98456(31) 0.04424(22) 0.10072(15)3.510 0.259 16 ·
32 2279 3.94649(29) 0.03657(21) 0.08764(14)3.540 0.240 16 ·
32 4067 3.89302(37) 0.02816(22) 0.07513(17)3.570 0.212 24 ·
32 2400 3.84392(17) 0.021767(89) 0.062829(68)3.630 0.170 24 ·
32 3232 3.75291(10) 0.013176(93) 0.045175(67)3.690 0.150 24 ·
32 2284 3.669908(81) 0.008740(85) 0.035734(47)3.760 0.130 24 · ·
48 2538 3.580005(77) 0.005781(55) 0.027805(20)3.820 0.125 32 ·
32 2913 3.508124(74) 0.004467(68) 0.024666(37)3.920 0.110 32 ·
32 4677 3.396477(51) 0.002967(69) 0.019635(15)4.080 0.081 32 ·
32 5607 3.234961(31) 0.001546(43) 0.012779(16)TABLE III: Expectation values of the pure gauge action density, light and strange quark chiral condensates calculated on zerotemperature lattices of size N σ N τ . Also given is the number of trajectories generated at each value of the gauge coupling β with light quarks of mass ˆ m l and bare strange quark mass ˆ m s = 10 ˆ m l . calculations on a lattice these, three regimes also deserve a separate discussion as discretization effects influencelattice calculations in these regimes quite differently. Before proceeding to a calculation of other bulk thermodynamicobservables we therefore will discuss in the following three subsections properties of ( ǫ − p ) /T in three temperatureintervals: (i) T < ∼
200 MeV or
T < ∼ T c , (ii) 200 MeV < ∼ T < ∼
300 MeV or 1 . < ∼ T /T c < ∼ . T > ∼
300 MeV or
T > ∼ . T c .
1. Trace anomaly at low temperatures
In Fig. 6 we show the low temperature part of ( ǫ − p ) /T obtained from our calculations with the p4fat3 action onlattices with temporal extent N τ = 4 and 6 and spatial size N σ /N τ = 4. We compare these results with calculationsperformed with the asqtad action [11] for N τ = 6. These latter calculations have been performed on lattices withsmaller spatial extent, N σ /N τ = 2, and results are based on lower statistics. These calculations are, however, consistentwith our findings. We also note that results obtained for two different values of the lattice cut-off, N τ = 4 and 6, arecompatible with each other.In the transition region from high to low temperature it is generally expected that thermodynamic quantities canbe described quite well by a hadron resonance gas (HRG) [6]; the freeze-out of hadrons in heavy ion experimentstakes place in this region and observed particle abundances are, in fact, well described by a HRG model [29]. Also acomparison of lattice results for the EoS with heavier quarks with a resonance gas model ansatz was quite satisfactory[30] but required the use of a suitably adjusted hadron mass spectrum. As we now can perform lattice calculationswith almost physical quark mass values a more direct comparison using the HRG model with physical quark massvalues should be appropriate.We use a HRG model constructed from all resonances, with masses taken from the particle data book up to a3 β
100 ˆ m l N σ h s G i τ h ¯ ψψ i l,τ h ¯ ψψ i s,τ ( ǫ − p ) /T p/T N τ = 4. The last two columns give the trace anomaly, ǫ − p , and the pressure, p , in units of T .(*) Note that at β = 3 .
82 simulations on N τ = 4 and 6 lattices have been performed at slightly different quark masses. T [MeV] Tr Θ F µµ /T p4: N τ =468 T [MeV] Tr m l (< ψψ > l,0 - < ψψ > l,T )2 m s (< ψψ > s,0 - < ψψ > s,T ) p4: N τ =468 FIG. 5: The fermionic contribution to the trace anomaly (left) and the ratio of the light and strange quark contributions toΘ µµF /T (right). maximal value m max = 2 . (cid:18) ǫ − pT (cid:19) low − T = X m i ≤ m max d i π ∞ X k =1 ( − η i ) k +1 k (cid:16) m i T (cid:17) K ( km i /T ) . (27)Here different particle species of mass m i have degeneracy factors d i and η i = − ǫ − p ) /T quite well, although the lattice data seem to drop somewhat faster at low temperature. Whether this points towardsa failure of the HRG model at lower temperatures, or is due to difficulties in correctly resolving the low energy hadron4 β
100 ˆ m l N σ h s G i τ h ¯ ψψ i l,τ h ¯ ψψ i s,τ ( ǫ − p ) /T p/T N τ = 6. The last two columns give the trace anomaly, ǫ − p , and the pressure, p , in units of T . β
100 ˆ m l N σ h s G i τ h ¯ ψψ i l,τ h ¯ ψψ i s,τ ( ǫ − p ) /T N τ = 8. The last column gives the trace anomaly, ǫ − p , in units of T . spectrum in the current calculations on still rather coarse lattices, will require more detailed studies on finer ( N τ = 8)lattices in the future. We will return to this question in Section VII.We also note that the current lattice calculations are performed with light quark masses that are a factor two largerthan the physical ones. Reducing the light quark masses to their physical values will shift the lattice data to smallertemperatures and will thus improve the comparison with the HRG model. From the known systematics of the quarkmass dependence of other thermodynamic quantities, e.g. the transition temperature, chiral condensates, or Polyakovloop expectation values [22] one can estimate this shift to be less than 5 MeV. Moreover, we note that the scale r used to convert lattice results to physical units has an error of about 2%. This is indicated in Fig. 6 by a horizontalerror bar for the data. Within this error all data may be shifted coherently.The low temperature region of the QCD EoS clearly deserves more detailed study in the future.
2. Trace anomaly in the strongly non-perturbative regime
At temperatures just above the transition temperature, ( ǫ − p ) /T shows the largest deviations from the conformallimit, ǫ = 3 p . The peak in ( ǫ − p ) /T at a temperature T max that is only slightly larger than the transition temperature T c constitutes a prominent structure of the trace anomaly which is relatively easy to determine in a lattice calculation.It is closely related to the softest point in the QCD equation of state [31], i.e. the minimum of p/ǫ as function of theenergy density. T max thus plays an important role for the construction of model equations of state that are consistentwith lattice calculations and may be used in hydrodynamic models for the expansion of dense matter created in heavy5 T [MeV] Tr ( ε -3p)/T p4: N τ =46asqtad: N τ =6 FIG. 6: Comparison of the low temperature part of ( ǫ − p ) /T calculated on lattices with temporal extent N τ = 4 and 6 witha resonance gas model that includes all resonances up to mass 2 . N τ = 6 data obtained with the p4fat3 action. Data for calculations with the asqtad action are taken from [11]. ion collisions. As T max and, in particular ( ǫ − p ) /T max , are easily determined they may also serve as consistencychecks between different lattice calculations.In Fig. 7 we show results for ( ǫ − p ) /T in the intermediate temperature interval 180 MeV < T <
300 MeV. Alsoshown here are results from calculations performed with the asqtad action on lattice with temporal extent N τ = 6[11]. As can be seen these calculations are in quite good agreement with the results obtained with the p4fat3 actionon lattice with the same temporal extent but larger spatial volume. Estimates for T max and the peak height onlattices with temporal extent N τ = 4 and 6 are given in Table VII. Also given there are estimates for the transitiontemperature T c r obtained previously in a dedicated analysis of the transition temperature on N τ = 4 and 6 lattices[22]. The values quoted in the Table give the fit results obtained from a joint fit of transition temperatures on bothlattice sizes for different quark mass values evaluated at the pseudo-scalar mass values corresponding to our LCP. Wenote that the temperature T max is only about 3% larger than the transition temperatures determined from peaks inthe chiral susceptibility.On the coarse N τ = 4 lattices the analysis of ( ǫ − p ) /T in the transition region is still quite sensitive to thenon-perturbative structure of the β -functions, R β and R β R m shown in Fig. 3; this region is still close to the strongcoupling regime below and in the vicinity of the dip in R β shown in Fig. 3(left). This seems to be the main reason forthe large differences seen in the peak height for ( ǫ − p ) /T between the N τ = 4 and 6 lattices. In the latter case thetransition and peak region is already in the regime where the lattice β -functions smoothly approach the continuumresults. We thus expect that these results are much less affected by this source of lattice artifacts. Nonetheless, abetter control over the cut-off dependence in this region clearly is needed and does require calculations on a largerlattice in order to control the continuum extrapolations of T max r as well as ( ǫ − p ) /T max .
3. Trace anomaly at high temperatures
In Fig. 8 we show results for ( ǫ − p ) /T in the high temperature regime, T > ∼ . T c . A comparison with dataobtained with the asqtad action on lattices with temporal extent N τ = 6 shows that the results obtained here withthe p4fat3 action are compatible with the former for T < ∼
400 MeV ( ∼ T c ). The current analysis performed with thep4fat3 action, however, has been extended to much larger temperatures, T ∼ T c , i.e. into the temperature regimeaccessible to heavy ion experiments at the LHC.For temperatures larger than T max the trace anomaly rapidly drops. Eventually, when the high temperatureperturbative regime is reached, the temperature dependence is expected to be controlled by the logarithmic running6 T [MeV] Tr ( ε -3p)/T p4: N τ =46asqtad: N τ =6 FIG. 7: The trace anomaly in the vicinity of the transition temperature calculated on lattices with temporal extent N τ = 4and 6 on lattices with aspect ratio N σ /N τ = 4. Data for calculations with the asqtad action are taken from [11] which havebeen performed on finite temperature lattices with a smaller physical volume corresponding to an aspect ratio N σ /N τ = 2. N τ T c r T c [MeV] ( T r ) max T max [MeV] ( ǫ − p ) /T max ǫ − p ) /T and its value calculated on lattices with different values of the temporal extent N τ on a line of constant physics that corresponds to a pion mass of about 220 MeV. Errors on the peak positions have beenestimated from cubic fits in the peak region by varying the fit intervals. The second and third columns show the transitiontemperature determined on the LCP used for this study of the EoS. For N τ = 4 this had been determined in [22] and for N τ = 6 in this analysis (see discussion in Section VI). of the QCD coupling constant. To leading order in high temperature perturbation theory Θ µµ ( T ) for massless quarksis given by [1], ǫ − pT = 13 b (cid:18) n f (cid:19) g ( T ) + O ( g ) , (28)with n f = 3 for massless 3-flavor QCD, which corresponds to the high temperature limit for our (2+1)-flavor QCDcalculations performed on a LCP with fixed non-zero quark mass values.For temperatures larger than about 2 . T c results for Θ µµ ( T ) /T obviously are sensitive to lattice cut-off effects.The results on N τ = 6 lattices drop significantly slower with temperature than the N τ = 4 results. In order to makesure that this effect does not superimpose with possible finite volume effects, we increased in this temperature regionthe spatial lattice size from 24 to 32 . No statistically significant volume effects have been observed for Θ µµ ( T ) /T ,although we observe a sensitivity of the zero temperature light and strange quark chiral condensates on the volume;as the condensates contribute less than 10% to the trace anomaly at these high temperature values (see Fig. 5)modifications of the condensates by a few percent contribute insignificantly to finite volume effects in Θ µµ ( T ) /T .Moreover, as the entire fermionic contribution, Θ µµF ( T ) /T , to the total trace anomaly is small for T > ∼
400 MeV, it isobvious that the contribution of the fermion condensates is not the source for the cut-off effects at high temperature.The cut-off dependence seen in Fig. 8 arises from the gluonic sector of Θ µµ ( T ) /T , which of course also receivescontributions from virtual quark loops.In the high temperature region we also added calculations on lattices with temporal extent N τ = 8 at 3 differentvalues of the temperature. Results from these calculations are summarized in Table VI and are also shown in Fig. 8.As can be seen in this figure results obtained for the trace anomaly on the N τ = 8 lattice are in good agreement withthe N τ = 6 results suggesting that remaining cut-off effects in this temperature range are small for N τ ≥ µµ ( T ) /T at high temperature also lead to larger values for the pressure, which isobtained from an integral over the trace anomaly, and also results in larger values for the energy and entropy densities, i.e. these quantities approach the Stefan-Boltzmann limit more rapidly on the N τ = 6 lattices than they did on the7 T [MeV] Tr ( ε -3p)/T p4: N τ =468asqtad: N τ =6 FIG. 8: The high temperature part of ( ǫ − p ) /T calculated on lattices with temporal extent N τ = 4, 6 and 8. The curvesshow fits to the N τ = 4 and 6 data with the ansatz given in Eq. 29. N τ = 4 lattice. It thus is important to get good control over cut-off effects at high temperatures and obtain furtherconfirmation of the results obtained in our N τ = 6 lattices, and through further calculations, on the N τ = 8 latticesat higher temperatures.As discussed previously, Θ µµ ( T ) /T contains a contribution from the vacuum quark and gluon condensates thatgets suppressed by a factor T at high temperature. In the case of a pure gauge theory it has, however, been notedthat up to temperatures a few times the transition temperature the dominant power-like correction to the perturbativehigh temperature behavior is O ( T − ) rather than O ( T − ) [32, 33]. These qualitative features also show up in ourresults for Θ µµ ( T ) /T at temperatures T > ∼ . T c . In Fig. 8 we show a comparison of the lattice results with such aphenomenologically motivated polynomial fit ansatz, (cid:18) ǫ − pT (cid:19) high − T = 34 b g + bT + cT . (29)Here we used the parametric form of the leading order perturbative result given in Eq. 28 with a tempera-ture independent coupling g to characterize the high temperature behavior of ( ǫ − p ) /T in the fit interval T ∈ [300 MeV ,
800 MeV]. For N τ = 4 we only performed a 2-parameter fit as it turned out that the fit doesnot require the contribution from a constant term ( g ≡ T ≥
300 MeV are given in Table VIII. We note that the vacuum condensate contribution ( ∼ c/T ) is small com-pared to the genuine thermal part. The present analysis, however, does not yet allow us to disentangle logarithmicfrom power-like (quadratic) corrections.Of course, it is tempting to relate the coefficient of the quartic term to a bag constant or zero temperature quarkand gluon condensate contribution, c = 4 B . This yields quite a reasonable value, B / = 247(25) MeV. Nonetheless,it seems that a more detailed analysis of the scaling behavior in the high temperature region and better control overcut-off effects is needed before a proper running of the gauge coupling can be established that would unambiguouslyallow to single out power-like (quadratic) corrections in the high temperature regime which then would also allow oneto establish a connection to the perturbative regime for the trace anomaly.8 N τ g b [GeV ] c [GeV ]4 – 0.101(6) 0.024(1)6 2.3(7) 0.16(6) 0.013(6)TABLE VIII: Fit parameters for 2-parameter fits ( N τ = 4) and 3-parameter fits ( N τ = 6) to ( ǫ − p ) /T in the region T ≥
300 MeV using the ansatz given in Eq. 29.
B. Pressure, energy and entropy density
As indicated in Eq. 4 we obtain the pressure difference,∆ p ( T ) ≡ p ( T ) T − p ( T ) T , (30)by integrating over the trace anomaly weighted with an additional factor of T − in the interval [ T , T ]. We havestarted our integration at T = 100 MeV, or T r ≃ .
24, by setting the trace anomaly to zero at this temperature.As discussed in the previous section, this leaves us with an uncertainty for the value of the pressure at T , which weestimate to be of the order of the pressure in a hadron resonance gas, i.e. [ p ( T ) /T ] HG = 0 . p ( T ) from our lattice calculations for the pressure at higher temperatures thus yield p/T up to asystematic uncertainty on p ( T ) /T . We also note again that the normalization at T does not take care of the overallnormalization of the pressure at T = 0.To calculate ∆ p ( T ) by integrating the numerical results obtained for Θ µµ ( T ) /T from Eq. 4, we have used straightline interpolations of our results for Θ µµ /T at adjacent values of the temperature. We also used stepwise interpolationsobtained by fitting quadratic polynomials to the data in small intervals that are matched to fits in the previous interval.Results of the latter approach are then used to perform the integration in the various regions analytically. Differencesbetween this approach and the straight line interpolations are nowhere larger than 1.5%. We then used the smoothpolynomial interpolations to determine the pressure and combined this result with that for Θ µµ ( T ) to obtain theenergy density. Both are shown in the left hand part of Fig. 9. The uncertainty arising from the normalization ofthe pressure at T is indicated as a small vertical bar in the upper right part of this figure. We note that at T ∼ T c results for p/T and ǫ/T stay about 10% below the ideal gas value.In particular, for applications to heavy ion phenomenology and for the use of the QCD equation of state inhydrodynamic modeling of the expansion of matter formed in heavy ion collisions, it is of importance to eliminatethe temperature in favor of the energy density and thus obtain the pressure as function of energy density. The ratio p/ǫ is shown in the right hand part of Fig. 9. As can be seen at low temperature, in the vicinity of the minimum in p/ǫ , results are consistent with values extracted for this quantity from a hadron resonance gas model. We also notethat in the high temperature regime it has been found in [34] that the ratio p/ǫ shows little dependence the baryonnumber density when evaluated on lines of constant entropy per baryon number.The density dependence of p/ǫ is related to the square of the velocity of sound c s = d p d ǫ = ǫ d p/ǫ d ǫ + pǫ . (31)In the high temperature limit as well as in the transition region where the derivative d( p/ǫ ) / d ǫ vanish, c s is directlygiven by p/ǫ . We therefore find that the velocity of sound is close to the ideal gas value, c s = 1 /
3, for energydensities ǫ> ∼
100 GeV/fm and drops by a factor of 4 to a minimal value of about ( c s ) min ≃ .
09 that is reached at ǫ> ∼ (1 −
2) GeV/fm . The dependence of p/ǫ on the energy density can be parametrized in the high temperature regionwith a simple ansatz [34], pǫ = 13 (cid:18) C − A Bǫ fm / GeV (cid:19) , (32)which then also allows a simple calculation of the velocity of sound, using Eq. 31. We find that the above parametriza-tion yields a good fit of the N τ = 6 data in the interval 1 . ≤ ǫ / / (GeV / fm ) / ≤ χ /dof of 1.3. For thefit parameters we obtain, C = 0 . A = 1 . B = 0 . ǫ ≃ the lattice calculations indicate a rise of p/ǫ as expected in hadronresonance gas models. However, the current resolution and accuracy of lattice calculations in this regime clearly isnot yet sufficient to allow for a detailed comparison between both.9
0 2 4 6 8 10 12 14 16 100 200 300 400 500 600 700 0.4 0.6 0.8 1 1.2 1.4 1.6
T [MeV] Tr ε SB /T ε /T : N τ =463p/T : N τ =46 0.000.050.100.150.200.250.300.35 0 1 2 3 4 5 6 ε [(GeV/fm ) ] p/ ε SBN τ =46fit: p/ ε HRG: p/ ε c s2 FIG. 9: Energy density and three times the pressure as function of the temperature (left) and the ratio p/ǫ as function of thefourth root of the energy density (right) obtained from calculations on lattices with temporal extent N τ = 4 and 6. Temperatureand energy density scales have been obtained using the parametrization of r /a given in Eq. 22 and r = 0 .
469 fm. The smallvertical bar in the left hand figure at high temperatures shows the estimate of the systematic uncertainty on these numbersthat arises from the normalization of the pressure at T = 100 MeV. The dashed curve (HRG) in the right hand figure showsthe result for p/ǫ in a hadron resonance gas for temperatures T <
190 MeV.
0 5 10 15 20 100 200 300 400 500 600 700 0.4 0.6 0.8 1 1.2 1.4 1.6
T [MeV] s/T Tr s SB /T p4: N τ =46asqtad: N τ =6 FIG. 10: Entropy density as function of the temperature obtained from calculations on lattices with temporal extent N τ = 4and 6. Temperature and energy density scales have been obtained using the parametrization of r /a given in Eq. 22 and r = 0 .
469 fm. The small vertical bar in the left hand figure at high temperatures shows the estimate of the systematicuncertainty on these numbers that arises from the normalization of the pressure at T = 100 MeV. As pointed out in Section II the non-perturbative vacuum condensates of QCD show up at high temperature aspower-like corrections to temperature dependence of the trace anomaly and consequently also to pressure and energydensity. These vacuum condensate contributions drop out in the entropy density which is shown in Fig. 10. It thusis an observable most suitable for comparisons with (resummed) perturbative calculations [15]. Like energy densityand pressure, the entropy also deviates from the ideal gas value by about 10% at T ∼ T c .We note that for T < ∼ T c the results obtained with the asqtad action [11] for the entropy density are in goodagreement with the results obtained with the p4fat3 action, although at least in the high temperature limit the cut-offdependence of both actions is quite different. This suggests that at least up to temperature T ≃ T c non-perturbativecontributions dominate the properties of bulk thermodynamic observables like the entropy density. It also gives rise tothe expectation that additional cut-off effects are small. Nonetheless, the result presented in this section on propertiesof bulk thermodynamic observables clearly need to be confirmed by calculations on lattices with larger temporalextent.0 VI. RENORMALIZED POLYAKOV LOOP AND CHIRAL CONDENSATES
As part of our analysis of bulk thermodynamic observables we have gathered a lot of information on the static quarkpotential at zero temperature. This has been discussed in Section IV and results obtained for V ¯ qq ( r ) have been usedthere to determine a temperature scale for our thermodynamic calculations. Furthermore, we have obtained a lot ofinformation on the chiral condensates at zero temperature that entered our calculation of thermodynamic quantities.Together with corresponding results on heavy quark free energies and chiral condensates at finite temperature thisallows us to analyze the deconfining properties as well as the change of chiral properties of the finite temperaturetransition in terms of observables which are related to exact order parameters for deconfinement and chiral symmetryrestoration in the infinite quark mass and vanishing quark mass limits of QCD, respectively.As discussed in the previous sections, the deconfining aspect of the finite temperature transition, i.e. the suddenliberation of partonic degrees of freedom in QCD, is reflected in the rapid change of bulk thermodynamic observables.This is also reflected in the rapid change of the static quark free energy which characterizes the response of a thermalmedium to the addition of static quark sources. The static quark free energy, F q , is related to the Polyakov loopexpectation value, h L i ∼ exp( − F q ( T ) /T ), h L i = * N σ X ~x L ~x + with L ~x = 13 Tr N τ Y x =1 U ( x ,~x ) , ˆ0 . (33)It may more rigorously be defined through the asymptotic large distance behavior of static quark-antiquark correlationfunctions [35], h L i = lim | ~x − ~y |→∞ h L ~x L † ~y i . (34)The Polyakov loop needs to be renormalized in order to attain a physically meaningful value in the continuum limit.To construct the renormalized Polyakov loop from the bare Polyakov loop expectation values, h L i , calculated onlattices with temporal extent N τ at a temperature controlled by the gauge coupling β , L ren ( T ) = Z N τ ren ( β ) h L i , (35)we can make use of our extensive calculations of the static potential at zero temperature. As outlined in SectionIV we have extracted renormalization constants, ( c ( β ) a ), from the matching of the static potential to the stringpotential. These renormalization constants are given in Table I in terms of the product c ( β ) r . With this we obtainthe renormalization constants for the Polyakov loop as, Z ren ( β ) = exp( c ( β ) a/ L ren on lattices with temporal extent N τ = 4 and 6 is small, which is in agreement with results obtained in studies of L ren in pure SU (3) gauge theories [35]. A similar renormalization of the Polyakov loop obtained in calculations with the1-link, stout smeared staggered fermion action has been used in [37]. The large cut-off dependence of L ren observedin this case mainly seems to be due to the choice of observable ( f K ) that has been used to set the temperature scale.In fact, when using r instead of f K to determine the lattice spacing, and thus the temperature, most of the cut-offdependence of L ren is removed in the data shown in [37].Another important aspect of the QCD transition is, of course, the change of chiral properties with temperature.This is generally reflected in the temperature dependence of the chiral condensate or related susceptibilities. Alsothe chiral condensates need to be renormalized to obtain finite, well defined quantities in the continuum limit. Toeliminate the quadratic divergences in the linear quark mass dependent correction to the chiral condensates [24]we calculate a suitable combination of light and strange quark condensates at finite temperature. We furthermorenormalize this quantity by the corresponding combination of condensates calculated at zero temperature at the samevalue of the lattice cut-off, i.e. at the same value of the gauge coupling β ,∆ l,s ( T ) = h ¯ ψψ i l,τ − ˆ m l ˆ m s h ¯ ψψ i s,τ h ¯ ψψ i l, − ˆ m l ˆ m s h ¯ ψψ i s, . (36)This eliminates multiplicative renormalization factors.In Fig. 11(right) we show results obtained for ∆ l,s ( T ) on the LCP for N τ = 4 and 6. In this figure, as well as in thecorresponding figure for L ren ( T ) shown on the left hand side, we also give estimates for the pseudo-critical temperatureextracted from the position of the peak in the disconnected part of the light quark chiral susceptibility. For N τ = 4this value has been determined previously by us [22] as the quark mass parameters for the LCP used here are closeto those used in [22] to determine T c on the N τ = 4 lattices for ˆ m l / ˆ m s = 0 .
1. For N τ = 6 the choice of the strange1 T [MeV] Tr L ren N τ =468 T [MeV] Tr ∆ s,l N τ =46 FIG. 11: Renormalized Polyakov loop on lattices with temporal extent N τ = 4, 6 and 8 (left) and the normalized difference oflight and strange quark chiral condensates defined in Eq. 36. The vertical lines show the location of the transition temperaturedetermined in [22] on lattices with temporal extent N τ = 4 (right line) and in this analysis for N τ = 6 (left line). quark mass differs slightly from the one used in that earlier study. We therefore performed a new determination ofthe transition temperature for the N τ = 6 lattice and the parameters of the LCP used here. From the peak positionsof the disconnected parts of the light and strange quark susceptibilities we find β c ( N τ = 6) = 3 . r /a quoted for this value of the coupling in Table I we find T c r = 0 . T c = 196(3).We note that the region of most rapid change in the subtracted and normalized chiral condensate, ∆ l,s ( T ), is ingood agreement with the region where the Polyakov loop expectation value as well as bulk thermodynamic quantities,e.g. the energy and entropy densities change most rapidly. VII. DISCUSSION AND CONCLUSIONS
We have presented here a detailed analysis of the QCD equation of state with an almost physical quark massspectrum. The current calculations have been performed with a physical strange quark mass value and two degeneratelight quark masses that are about a factor two larger than the physical average light quark mass value. In a widetemperature range, results have been obtained on large spatial lattices close to the thermodynamic limit for twodifferent values of the lattice cut-off, corresponding to lattices of temporal extent N τ = 4 and 6. At high temperaturesadditional calculations on lattices with temporal extent N τ = 8 have been performed, which allow us to controlapparent cut-off effects in this temperature range. All finite temperature calculations have been supplemented withcorresponding zero temperature calculations to perform necessary vacuum subtractions and to accurately set thetemperature scale.At high temperature, T > ∼ T c , bulk thermodynamic observables such as pressure, energy and entropy density deviatefrom the continuum Stefan-Boltzmann values only by about 10% and show little cut-off dependence. This weak cut-off dependence could only be achieved through the use of O ( a ) improved gauge and fermion actions. On the otherhand, a closer look at the trace anomaly, ( ǫ − p ) /T , from which these quantities are derived, clearly unravels cut-offeffects when comparing results obtained for the N τ = 4 and 6 lattices; for temperatures T > ∼ . T c or equivalently T > ∼
500 MeV results for ( ǫ − p ) /T on the N τ = 4 lattices are systematically lower than for N τ = 6. Additionalcalculations performed on N τ = 8 lattices in this high temperature region are consistent with the results obtained on N τ = 6 lattices and thus suggest that cut-off effects are small on lattice with temporal extent N τ ≥
6. Of course, thisshould be confirmed through additional calculations on lattices with temporal extent N τ = 8 at larger temperatures.On these fine lattices it also will be interesting to analyze in more detail the contribution of charm quarks to theequation of state [38, 39]. Our earlier analysis for m l = 0 . m s on a 16 T c r = 0 . T c = 201(2). ǫ − p ) /T at high temperature clearly is important whenone wants to make contact between lattice calculations for e.g. the entropy density and high temperature perturbationtheory. Although our present high statistics analysis seems to have achieved good control over the cut-off dependenceof ( ǫ − p ) /T in this high temperature regime, a more extended analysis of the temperature dependence on N τ = 8lattices is still needed to make firm contact with perturbative or resummed perturbative calculations.At low temperatures, T < ∼
200 MeV, the influence of cut-off effects is less apparent. We observe that at a givenvalue of the temperature results for ( ǫ − p ) /T obtained on the N τ = 6 lattice are systematically larger than thoseobtained on the N τ = 4 lattice. This is, in fact, expected and is consistent with the cut-off dependence observed incalculations of the transition temperature on lattices with temporal extent N τ = 4 and 6 [22]. Also on lattices withtemporal extent N τ = 8 indications for this to happen have been found in preliminary studies of chiral and quarknumber susceptibilities as well as Polyakov loop expectation values [40]. We thus expect that with increasing N τ , i.e. closer to the continuum limit, the region where ( ǫ − p ) /T and all other thermodynamic observables will start to riserapidly, continues to shift towards smaller temperatures. A further, although smaller shift of the transition regiontowards smaller values of the temperature will arise from an extrapolation to physical quark masses. Judging fromthe known temperature dependence of the transition temperature [22] and other thermodynamic observables, like thePolyakov loop expectation value or quark number susceptibilities [40], this will amount to a shift of the scale by a fewMeV. In fact, extrapolations of the transition temperature in quark mass and lattice spacing to the physical pointhave been performed by several groups for staggered as well as Wilson fermions [22, 41, 42]. These extrapolationsconsistently show that the quark mass dependence of the transition temperature is weak. We take the quark massdependence of the transition temperature as indicator for the shift of the transition region one has to expect in futurecalculations on finer lattices with physical values of the quark masses. We also should stress that current estimatesfor the bare lattice parameters that correspond to physical values of the light quark mass, ˆ m q ≃ .
04 ˆ m s , of course,are based on studies of the spectrum on lattices with finite cut-off. Eliminating these systematic effects will requirefurther calculations on finer lattices. This may also improve the comparison with model calculations of the equationof state in the low temperature phase of QCD. Acknowledgments
This work has been supported in part by contracts DE-AC02-98CH10886 and DE-FG02-92ER40699 with the U.S.Department of Energy, the Helmholtz Gesellschaft under grant VI-VH-041, the Gesellschaft f¨ur Schwerionenforschungunder grant BILAER and the Deutsche Forschungsgemeinschaft under grant GRK 881. Numerical simulations havebeen performed on the QCDOC computer of the RIKEN-BNL research center, the DOE funded QCDOC at BNL,the apeNEXT at Bielefeld University and the new BlueGene/L at the New York Center for Computational Science(NYCCS). We thank the latter for generous support during the burn-in period of NYBlue. [1] J. I. Kapusta, Nucl. Phys. B , 461 (1979).[2] C. X. Zhai and B. M. Kastening, Phys. Rev. D , 7232 (1995).[3] P. Arnold and C. x. Zhai, Phys. Rev. D , 1906 (1995).[4] A. Vuorinen, Phys. Rev. D , 054017 (2003).[5] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, Phys. Rev. D , 105008 (2003);F. Di Renzo, M. Laine, V. Miccio, Y. Schroder and C. Torrero, JHEP , 026 (2006)[6] P. Braun-Munzinger, K. Redlich and J. Stachel, in Quark Gluon Plasma 3, eds. R.C. Hwa and Xin-Nian Wang, WorldScientific Publishing, Singapore 2004, nucl-th/0304013.[7] C. W. Bernard et al. [MILC Collaboration], Phys. Rev. D , 6861 (1997).[8] A. Ali Khan et al. [CP-PACS collaboration], Phys. Rev. D , 074510 (2001).[9] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier and B. Petersson, Nucl. Phys. B , 419 (1996).[10] F. Karsch, E. Laermann and A. Peikert, Phys. Lett. B , 447 (2000).[11] C. Bernard et al. , Phys. Rev. D , 094505 (2007).[12] Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP , 089 (2006).[13] U. M. Heller, F. Karsch and B. Sturm, Phys. Rev. D , 114502 (1999).[14] T. Blum et al. , Phys. Rev. D , 1133 (1997).[15] J. P. Blaizot, E. Iancu and A. Rebhan, Phys. Lett. B , 181 (1999) and Phys. Rev. D , 065003 (2001).[16] J. O. Andersen, E. Braaten, E. Petitgirard and M. Strickland, Phys. Rev. D , 085016 (2002).[17] M. Bluhm, B. Kampfer, R. Schulze, D. Seipt and U. Heinz, Phys. Rev. C , 034901 (2007)[18] S. R. Sharpe, PoS LAT2006 , 022 (2006).[19] M. Creutz, PoS
LAT2006 , 208 (2006) and PoS (Lattice 2007) 007. [20] A.Kronfeld, PoS (Lattice 2007) 016.[21] I. Horv´ath, A. D. Kennedy and S. Sint, Nucl. Phys. Proc. Suppl. , 834 (1999);M. A. Clark, A. D. Kennedy and Z. Sroczynski, Nucl. Phys. Proc. Suppl. , 835 (2005).[22] M. Cheng et al. , Phys. Rev. D , 054507 (2006).[23] M. Cheng et al. , Phys. Rev. D , 034506 (2007).[24] J. Gasser and H. Leutwyler, Phys. Rept. , 77 (1982).[25] C. Bernard et al., Phys. Rev. D , 054506 (2001).[26] A. Gray et al., Phys. Rev. D72 , 094507 (2005).[27] C.T.H. Davies et al., Phys. Rev. lett. , 022001 (2004).[28] T. Ishikawa et al, PoS LAT2006, 181 (2006).[29] A. Andronic, P. Braun-Munzinger and J. Stachel, Nucl. Phys. A , 167 (2006).[30] F. Karsch, K. Redlich and A. Tawfik, Eur. Phys. J. C , 549 (2003).[31] C.M. Hung and E. Shuryak, Phys. Rev. Lett. , 4003 (1995).[32] R. D. Pisarski, Prog. Theor. Phys. Suppl. , 276 (2007), hep-ph/0612191.[33] E. Megias, E. Ruiz Arriola and L. L. Salcedo, Phys. Rev. D , 105019 (2007).[34] S. Ejiri, F. Karsch, E. Laermann and C. Schmidt, Phys. Rev. D , 054506 (2006).[35] O. Kaczmarek, F. Karsch, P. Petreczky and F. Zantow, Phys. Lett. B , 41 (2002);O. Kaczmarek and F. Zantow, Phys. Rev. D , 114510 (2005);P. Petreczky and K. Petrov, Phys. Rev. D , 054503 (2004).[36] C. Aubin et al. , Phys. Rev. D , 094505 (2004).[37] Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B , 46 (2006).[38] M. Laine and Y. Schroder, Phys. Rev. D , 085009 (2006).[39] M. Cheng [RBC-Bielefeld Collaboration], PoS (Lattice 2007) 173.[40] F. Karsch, PoS (CPOD07), 026 (2007), arXiv:0711.0656 [hep-lat]C. DeTar and R. Gupta, PoS (Lattice 2007) 179.[41] C. Bernard et al. [MILC Collaboration], Phys. Rev. D71