Nt = 6 equation of state for two flavor QCD

We improve the calculation of the equation of state for two flavor QCD by simulating on Nt = 6 lattices at appropriate values of the couplings for the deconfinement/chiral symmetry restoration crossover. For amq = 0.0125 the energy density rises rapidly to approximately 1 GeV/fm3 just after the crossover mπ/mp ≈ 0.4 at this point). Comparing with our previous result for Nt = 4 [1], we find large finite Nt corrections as expected from free field theory on finite lattices. We also provide formulae for extracting the speed of sound from the measured quantities.

cc c University of Utah Institutional Repository Author Manuscript 1 The Nt = 6 equation of state for two flavor QeD C. Bernard a, T. Blum b, C.E. DeTar c, Steven Gottlieb d, U.M. Heller e, J.E. Hetrick b, L. Karkkainen f, K. Rummukainen d, R.L. Sugar g, D. Toussaint b , and M. Wingate h a Department of Physics, Washington University, St. Louis, MO 63130, USA bDepartment of Physics, University of Arizona, Tucson, AZ 85721, USA CPhysics Department, University of Utah, Salt Lake City, UT 84112, USA dDepartment of Physics, Indiana University, Bloomington, IN 47405, USA eSCRI, Florida State University, Tallahassee, FL 32306, USA fNordita, Blegdamsvej 17, DK-2100 Copenhagen 0, Denmark gDepartment of Physics, University of California, Santa Barbara, CA 93106, USA hphysics Department, University of Colorado, Boulder, CO 80309, USA We improve the calculation of the equation of state for two flavor QCD by simulating on Nt = 6 lattices at appropriate values of the couplings for the deconfinement/chiral symmetry restoration crossover. For amq = 0.0125 the energy density rises rapidly to approximately 1 GeV /fm3 just after the crossover(m7r/mp ~ 0.4 at this point). Comparing with our previous result for Nt = 4 [1] , we find large finite Nt corrections as expected from free field theory on finite lattices. We also provide formulae for extracting the speed of sound from the measured quantities. 1. INTRODUCTION In order to show the existence of the Quark­Gluon Plasma (QGP) in the aftermath of up com­ing heavy-ion collision experiments at RHIC and CERN or to understand the dynamics of the QGP in the early universe, one needs as input, among other things, the equation of state for QCD, i.e. the energy density and the pressure as a function of temperature and quark mass. We are continuing our program of computing the equation of state for two flavor QCD. Last year we reported results for Nt = 4 [1], and we are now working at Nt = 6. A similar program for the pure gauge theory is being pursued by the Biele­feld group[2]. The Nt = 6 simulations represent a significant increase in computational cost due to the increased lattice size, smaller quark masses, and corresponding smaller simulation step sizes. Step size (6.t) errors induced by the approximate integration of the gauge field equations of mo­tion are much more significant at Nt = 6 than at Nt = 4. They are handled by extrapolation of observables to 6.t = 0, and by running at small step sizes, with a corresponding large increase in the cost of simulation. We have surveyed the gauge coupling and quark mass plane for two flavor QCD in a region relevant to the nonzero temperature crossover in order to measure the nonperturbative pressure by integration. The interaction measure is also cal­culated and together with the pressure it yields the energy density. 2. THEORY A Euclidean N: x Nt lattice with periodic boundary conditions has a temperature T and volume V V = N:a3 , liT = Nta, (1) where a is the lattice spacing. Thermodynamic variables are derivatives of the partition function cc cc University of Utah Institutional Repository Author Manuscript 2 Z. In particular, the pressure p and energy den­sity e are given by p 8logZ T oV and 8logZ eV = - o(l/T) . (2) (3) The methods for computing the pressure and integration measure are discussed in Ref. [1]. The pressure is found by integrating either the plaque­tte or 1/J7jJ. The interaction measure is given by IV T 1 8log Z V8log Z T o(l/T) - 3 av o 0 8logZ = (-at- - as-)logZ = --- oat oas 8log a (4) (5) (6) where as and at are the spatial and temporallat­tice spacings. The scale dependence in the case of QCD with dynamical quarks leads to (7) The derivatives are the usual fJ function and the anomalous dimension of the quark mass. In the above the subscript "sym" refers to symmetric lattices with Nt = Ns which are used for vacuum subtraction. The knowledge of the non-perturbative pres­sure and interaction measure allows us to com­pute other bulk quantities. The energy and en­tropy s become e = I + 3p ,sT = I +4p . (8) 3. SOUND SPEED Acoustic perturbations travel III the system with a speed Cs : 1 de c; dp (9) One has to take the derivative keeping the phys­ical quark mass fixed, i. e. on the line of con­stant physics. Unfortunately, we know best only the variations of the energy density and pressure along lines of constant bare parameters. In or­der to measure the correct sound speed one has to use the fJ function to map the changes in bare parameters to physical changes of temperature. For QCD with zero chemical potential the only varying quantity is the temperature T . In the real world the masses of the quarks do not change and derivatives are taken at fixed quark mass. The volume is infinite and it is divided out of the equations that consider densities. Hence de: I d!T V,m7r /mp dp . dT lv,m7r/mp (10) Henceforth we drop the Im 7r/m p reminder. With only T varying, the fundamental relation of ther­modynamics becomes J(T) = e(T) - Ts(T) = -p (11) or e(T) = Ts(T) - p(T). (12) Then, :~ = Ts'(T) + s(T) - p'(T) = Ts'(T), (13) where we have utilized a Maxwell relation for the entropy s: s= ooTP lv=p'(T) . (14) In the Maxwell relation we interchanged pressure and free energy according to pV = logZ = -JV. T T Hence, Ts'(T) p'(T) dlog(s) dlogT' (15) (16) cc cc University of Utah Institutional Repository Author Manuscript where we have used Eq. (14) again. However, lattice simulations are always done at finite volume, and the volume varies with the lat­tice spacing a. To get the correct value one must take the derivative with respect to temperature with the volume constant. The derivative with respect to temperature is problematic. It requires asymmetric lattice spac­ing and leads to expressions with asymmetry co­efficients. These in turn, are poorly known in the regime of bare couplings, where simulations are currently feasible. It would be advantageous to find an expression that does not involve the asymmetry coefficients. Using Eq. (8), one can give the sound speed (16) with the interaction measure: 1 1 dI -c;= -s- d+T 3. (17) On the other hand, Eq. (6) allows us to write 1 1 o[I.. 8log z] _ = 3 _ V 8log a . c; s aT (18) Now, the T derivative is taken at constant V. Therefore, ~ = 3 _ ~ [~ Olog Z + T o[~]] . (19) c; s V Olog a V aT The last derivative can be given as o[8log z] 8loga aT lo[ei] T Ologa' (20) which can be seen by expanding the derivatives with respect to a and T as derivatives with re­spect to at and as: a a a Ologa ata+aS&' at as a 1 a aT --at-· (21) T oat We also used the definition of energy density, Eq. (3), in the form EV Olog Z - =-at---· T oat (22) 3 Then Eq. (19) becomes 1 1 [ T o[ei ]] c; = 3 + sT I - V Olog a . (23) Using sT = E + p and I = E - 3p this can be expressed in a rather compact form 1 1 [ T o[ei ]] c; = E + P 4E - V Olog a . (24) If the system is conformally invariant the energy density does not depend on the scale and the derivative term in (24) disappears: 1 4E - - - - 3 (25) C; - E+P- , the relativistic result for a free gas. For QeD with dynamical quarks an explicit form is The derivatives have to be taken on a line of con­stant physics. 4. SIMULATIONS We have measured (0) and \ ij;7jJ) on asymmet­ric lattices with Nt = 6 and Ns = 12, and on symmetric lattices with Nt = Ns = 12. The couplings in the simulations are appropri­ate for the Nt = 6 crossover [3]. At amq = 0.0125 and 0.025 we varied 6/ g2 between 5.37 and 5.53 and 5.39 and 5.53, respectively. These runs will allow us to extrapolate results to zero quark mass. At 6/g2 = 5.45 and 5.53, we have varied amq be­tween 0.01 and 0.1 and 0.0125 and 0.2, respec­tively. These runs also allow us to extrapolate to zero quark mass, but in addition they provide a cross check on the integrations over 6/ g2 and give information on the equation of state over a wide range of quark masses. Since we are simulating with two flavors of Kogut-Susskind fermions, the simulations are performed with the refreshed molecular dynam­ics R algorithm [4] . For the Nt = 6 lattices we cc cc University of Utah Institutional Repository Author Manuscript 4 ran at least 1800 trajectories after 200 trajecto­ries for thermalization. On the symmetric lattices we performed at least 800 trajectories after 200 trajectories for thermalization. Each trajectory had unit length in simulation time. The R algorithm induces an error in the ob­servables with a leading term proportional to the square of the step size, in simulation time, used to integrate the equations of motion of the gauge fields. In practice, the error is small for each ob­servable. However, the errors are different on the symmetric and asymmetric lattices, so they do not cancel from the vacuum subtractions. More­over, in many instances the errors are the same or­der of magnitude as the vacuum subtracted quan­tities. Therefore, they must be eliminated. In most cases, the step size errors are elimi­nated by extrapolation. For each gauge coupling and quark mass several simulations are run with two or more step sizes. If the step sizes are small enough, observables depend linearly on the step size squared. Once this linear dependence is ob­served, the quantities are extrapolated to zero step size. As a rule of thumb, the step size in the R algorithm should be less than or approx­imately equal to amq . Roughly speaking, this is because in the integration of the gauge mo­mentum, the fermion "force" to lowest order is proportional to 1/ amq . Thus the step taken in simulation time to update the momentum should be ~ amq (or smaller) to keep the change in the momentum less than 0(1). Thus we have used 0.007 ~ /::::,.t ~ 0.015 and 0.015 ~ /::::,.t ~ 0.03 for the runs with amq = 0.0125 and 0.025, re­spectively. For the larger quark mass runs at 6/ g2 = 5.45 and 5.53, we used 0.02 ~ /::::,.t ~ 0.03. Finally, for the smallest quark mass simulation, amq = 0.01 (6/ g2 = 5.45) , we use /::::,.t = 0.005. 5. EXTRAPOLATIONS In Fig. 1 we show the plaquette dependence on step size squared for amq = 0.0125. Similar results hold for amq = 0.025. Generally, the ef­fects are larger on the symmetric lattices and at smaller 6/g2 . For example, at 6/g2 = 5.39, lin­ear behavior sets in for much smaller step size on the cold lattice than on the hot lattice. In fact , 1\ D V ,ii-" --------405.53 1.65 )k--t: -------EO)--------1E') 5.47 i~~; ==~;;~O~5.~. ! 4~6===IlEll_ __ - _=€I=OO 5.45 :~======:E-l :~~ i!3 =1l5. 4 3 ~ m~ ~~~5.41 ,: ~5.40 ~ 1.60 5.39 1. 5 5 L...L.--'-----'---L---'----L-....l----"-----'-----"----'------"-----'-----'---..J 0.0000 0.0001 0.0002 bot 2 Figure 1. The plaquette as a function of /::::"t2 for amq = 0.0125. Lines are to guide the eye and are not fits. Labels for the couplings refer to the hot lattices (octagons). for /::::"t2 ~ 0.0001, one would conclude that the system is in the confined phase; only at smaller step size is a clear separation visible. From Fig. 1, the following is a reasonable extrapolation proce­dure. First, on the symmetric lattices for 5.39 ~ 6/ g2 ~ 5.43, use only the smallest two step sizes to extrapolate to /::::"t2 = O. At 6/g2 = 5.37 we effectively assume the system is in the cold phase and take the values at /::::,.t = 0.007 as the zero step size extrapolations. For 6/ g2 > 5.43 we use all /::::,.t values to do the extrapolations. On the hot lattices we do the following. For each 6/ g2 where there are multiple step sizes, we use all of them to extrapolate to /::::"t2 = O. For 6/g2 = 5.40 and 5.42 we interpolate the slope from the neighbor­ing points and use that along with the point at /::::,.t = 0.01 to extrapolate to zero step size. Note, for 6/ g2 ~ 5.43, the slopes are essentially zero. Therefore, for 6/ g2 > 5.43 where only one mea- cc cc University of Utah Institutional Repository Author Manuscript surement is available, we take that value as the zero step size value. For (1/J1jJ) the situation is similar to the pla­quette with the following exceptions. First, after vacuum subtraction, the relative step size errors are not as significant. Second, on the hot lattices for 6/g2 > 5.39, we find the slopes with respect to 6.t2 are zero within errors, so we take the zero step size (1/J1jJ) to be the value measured at the smallest 6.t for each of these couplings. For amq = 0.025, the situation is similar, but the smallest step size runs are still incomplete as of this writing. Because the cold lattices vary smoothly with 6/ g2, we made cold runs at values of 6/ g2 sepa­rated by 0.02, while the hot runs were separated by 6.6/ g2 = 0.01 near the transition. Cold ob­servables at the other couplings were obtained by interpolation. For 6/g2 ::; 5.47, quadratic fits to the zero step size plaquette and 1/J1jJ had X2 = 2.69 and 3.18 with three degrees of freedom respectively. These fits were used for the interpo­lated values. To get the symmetric observables at 6/ g2 = 5.53, one can either extrapolate in 6/ g2 for fixed amq , or extrapolate in amq for fixed 6/ g2. Both give the same result within errors, and we simply take the value from extrapolation in amq as it has a much smaller error. 6. PRESSURE The results for (1/J1jJ) and (0) from the previous section can now be integrated to yield the pres­sure. To begin consider (1/J1jJ) as a function of amq . Using Eq. (5), we find the pressure as a func­tion of amq (at 6/g2 = 5.45 and 5.53). The result is shown in Fig. 2. We also want the pressure at amq = 0 which is found by setting 1/J1jJ(0) = 0 and continuing the integration to amq = O. At 6/ g2 = 5.45, a linear fit to the data for amq ::; 0.025 that is constrained to go through the origin has X2 = 4.9 for three degrees of freedom. For 6/ g2 = 5.53, we only have mea­surements at amq = 0.0125, 0.025, and 0.05 for which the linear fit does not work well. However, a quadratic fit constrained to go through the ori­gin has X2 = 0.49 for one degree of freedom. The 5 pressure at amq = 0 calculated in this manner is also shown in Fig. 2. 4 ,- 3 ! "<j< :;, I t:2 o --- 6/g2= 5.53 p... o --- 6/g2= 5.45 \ I 1 Cj) tlI II IDII) II) 0 0.0 0.1 0.2 amq Figure 2. The pressure from integration of < 1/J1jJ > with respect to amq . The bursts are ex­trapolations to zero quark mass. Next we integrate the plaquette with respect to 6/ g2 to obtain the pressure as a function of 6/ g2 at fixed amq . The result for amq = 0.0125 is shown in Fig. 3. The pressure rises smoothly through the crossover region as it must. The re­sults from the quark mass integrations are also shown in Fig. 3, and are in agreement with the 6/ g2 integration. This is a good check on our analysis, since for the most part the two integra­tions are independent. In particular, the integra­tion with respect to 6/ g2 is sensitive to the step size extrapolations while the amq integration is not. 7. INTERACTION MEASURE The interaction measure is given by Eq. (7). For the f3 function we use our previous result cal­culated from the 7r and p masses at various values of 6/ g2 and amq [1] . The result is shown in Fig. 4. It rises sharply from zero in the cold phase to C C H ~ :> c rt :::J 0 H ~ ~ ~ C (fJ n .H.. ... ~ rt cc University of Utah Institutional Repository Author Manuscript 6 3 ! 2 II "<j< E-< ~ P-. !d 1 ~~ ~ iii 0 ID 5.35 5.40 5.45 5.50 6/g 2 Figure 3. The pressure from integration of (0) with respect to 6/ g2 (amq = 0.0125) . The crosses are from the (1{;'ljJ) integration. some maximum and then begins to drop off. For a noninteracting plasma, I is zero since E = 3p which is certainly not the case at the largest value of 6/ g2 in our simulations. From asymptotic free­dom we expect the system to approach a nonin­teracting plasma at very high temperature. The interaction measure at amq = 0 depends only on the plaquette since the anomalous dimen­sion of the quark mass is zero at amq = O. To extrapolate the plaquette to zero quark mass we use the fact that its slope is just its correlation with 1{;'ljJ, 8~~~q) = (01{;'ljJ) - (0) (1{;'ljJ). (27) Since 1{;'ljJ is discontinuous at the origin in the bro­ken phase and continuous in the chirally symmet­ric phase, we expect a cusp at the origin for the cold lattices and zero slope for the hot lattices [l]. At 6/ g2 = 5.45 a linear fit on the cold lattices for amq ::; 0.05 gives X2 = 1.6 with two degrees of freedom. On the hot lattice a quadratic fit con­strained to zero slope at the origin gives X2 = 0.74 with two degrees of freedom. At 6/ g2 = 5.53 the situation is less satisfactory. For the cold lattices, I T 1 10 I-- - "<j< E-< ~ ~ P-. C'J I 5 I-- W '-"' -, p O ~~~I~~~I~~~~I~~ 5.35 5.40 5.45 5.50 6/g2 Figure 4. The interaction measure for amq 0.0125. a quadratic fit to the data with amq ::; 0.1 has X2 = 0.09. Note the smallest quark mass here is 0.025, and we also use this fit for the cold ob­servables at amq = 0.0125. On the hot lattices the data seem to reach a maximum at nonzero quark mass, so we take the measurement at the smallest quark mass to be the extrapolated value. Similar behavior at 6/ g2 = 5.53 was observed in our earlier Nt = 4 study, which may indicate a systematic error. Naively one expects finite vol­ume effects to order the lattice which is opposite to the observed behavior. The data and the fit results are shown in Fig. 5. 8. EQUATION OF STATE In Fig. 6 we show the Nt = 6 equation of state as a function of 6/g2 for amq = 0.0125. The zero quark mass extrapolations are also shown. At the largest value of 6/ g2, E is still much larger than 3p; for amq = 0, 3p is approximately 75% of E. On the other hand, after a rapid rise, E/T4 more or less levels off for 6/ g2 ~ 5.43. This corre­sponds to an energy density of about 1 GeV /fm3 at 6/ g2 = 5.43 (we use the p mass to convert 6/ g2 to temperature). m7r / mp at this point is slightly C C H ~ :> C rt :::J 0 H ~ ~ ~ C (fJ n .H.. ... ~ rt cc /\ D V University of Utah Institutional Repository Author Manuscript I 1 'f- (;) (;) ~ ;;: ,j, CD [i] [i] CD CD 1.65 I-- [i] @ - JE em IIlCD CD IIlffin ffi CD 1.60 I-- § - I!I I I 0.0 0.1 0.2 amq Figure 5. The plaquette as a function of the quark mass. The lower(upper) curves are for 6/g2 = 5.45(5.53). The bursts are results from fits to the data except for the hot curve (oc­tagons) at 6/g2 = 5.53 where the measured value at amq = 0.0125 is taken as the zero quark mass result. less than 0.4. In Fig. 7 we compare the Nt = 4 and 6 equa­tions of state. There is a large finite size ef­fect as expected from the free field results on finite lattices (also shown) . For example, p/T4 for amq = 0 differs by about 15%. There ap­pears to be a sizable quark mass effect as well (the Nt = 4 results are for amq = 0.1 and 0.025). For Nt = 6 the approach to the Stefan-Boltzmann law is unclear since we do not have data at very high temperatures. The prominent peak in E /T4 for Nt = 4 just after the crossover is much smaller at Nt = 6. We have plotted the results versus temperature by using the p mass to set the scale. The crossover temperature is around 150 MeV and is rather insensitive to Nt and amq . In a re­cent preprint, Asakawa and Hatsuda have pointed out that many features of this equation of state are constrained by fundamental thermodynamic relations[5]. 7 15 IIIrI f ~ , ..,. 10 E--< ~ HI! · 0.. ~ (Y') ..,. E--< ~ Q) 5 ~ ~ T T T iji Cj) ilJ o ill ~~~~~~~~~~~~~~ 5.35 5.40 5.45 5.50 Figure 6. The Nt = 6 equation of state for amq = 0.0125 and extrapolations to amq O(bursts). The upper curve is E /T4 . 9. SOUND SPEED The sound speed is a quantity that depends on the second derivative of the partition function. Therefore it is more difficult to get its value than the values for the thermodynamic variables dis­cussed so far. Even if our formulation can avoid the asymmetry coefficients there is an awkward derivative of the energy density in the formula. To measure the change in the mass as accu­rately as possible we performed a set of simula­tions at amq = 0.09 with Nt = 4 lattices in addi­tion to our old data at amq = 0.1. The result is shown in Fig. 8. Close to the transition the error in the derivative of the energy density overwhelms our data and we have large error bars. For large temperatures, we see that the speed approaches the ideal gas value. 10. CONCLUSION We have calculated the equation of state for QCD with two flavors of quarks on Nt = 6 lat­tices. The algorithm used to generate gauge con­figurations introduces a step size error in observ­abIes that must be removed by extrapolation. cc cc University of Utah Institutional Repository Author Manuscript 8 20 o ~~~~~~--~~~~~~--~ 0. 1 0.2 T(GeV) 0.3 Figure 7. Comparison of the equation of state for Nt = 4 (solid lines) and 6 (dashed lines) . The results shown are for amq = 0.0125 (diamonds), 0.025 (octagons), and 0.1 (squares). Bursts are extrapolations to amq = o. The horizontal lines give the Stefan-Boltzmann law for Nt = 4, 6, and the continuum (lowest line) . The added computational cost is significant. For amq = 0.0125, we find the energy density just af­ter the crossover to be roughly 1 GeV /fm3 , and for T almost twice the critical value, the high­est that we simulated, three times the pressure is only 60% of the energy density. We find large ef­fects of nonzero lattice spacing from comparison with results at Nt = 4, as expected from free field theory. After the completion of runs at amq = 0.025, we will complete our extrapolation to zero quark mass. It still remains to eliminate remaining lat­tice size effects, include the effects of the strange quark, eliminate any effects of finite volume, and include the effect of nonzero net quark density. This work was supported by the US DOE and NSF. Computations were done at the San Diego Supercomputer Center, the Cornell Theory Cen­ter, and Indiana University. 40 C\l I rtl U 20 o 0.0 0.1 0.2 0.3 0.4 T (GeV) Figure 8. The inverse sound speed squared for the Nt = 4 system at mqa = 0.1. REFERENCES 1. T. Blum, S. Gottlieb, L. Karkhiinen, and D. Toussaint, Nucl. Phys. B (Proc. Suppl.) 42, 460, 1995; T . Blum, S. Gottlieb, L. Karkhiinen, and D. Toussaint, Phys. Rev. D. 51 (1995) 5153. 2. E. Laermann et al., Nucl. Phys. B (Proc. Suppl.) 42, (1995) 120; F. Karsch, hep­lat/ 9503010, review talk from "Quark Matter 95"; G. Boyd et al., hep-Iat/9506025. 3. C. Bernard, et al., Phys. Rev. D 45 (1992) 3854. 4. S. Gottlieb, W. Liu, R. L. Renken, R. L. Sugar and D. Toussaint, Phys. Rev. D 35 (1987) 2531. 5. M. Asakawa and T. Hatsuda, hep-ph/9508360.