The $N_t=6$ equation of state for two flavor QCD

TL;DR

This paper improves the calculation of the equation of state for two-flavor QCD by using $N_t=6$ lattices, revealing significant finite lattice spacing effects and providing formulas for the speed of sound.

Contribution

It presents the first $N_t=6$ lattice results for the QCD equation of state, highlighting finite lattice corrections and offering methods to determine the speed of sound.

Findings

Energy density rises rapidly after crossover

Large finite $N_t$ corrections observed

Formulas for speed of sound provided

Abstract

We improve the calculation of the equation of state for two flavor QCD by simulating on $N_t=6$ lattices at appropriate values of the couplings for the deconfinement/chiral symmetry restoration crossover. For $am_q=0.0125$ the energy density rises rapidly to approximately 1 ${\rm GeV/fm^3}$ just after the crossover($m_\pi/m_\rho\approx 0.4$ at this point). Comparing with our previous result for $N_t=4$~\cite{eos}, we find large finite $N_t$ corrections as expected from free field theory on finite lattices. We also provide formulae for extracting the speed of sound from the measured quantities.

Equations

$$V=N_s^3{a}^3,$$

$$1/T=N_{t}a,$$

$$\frac{p}{T}=\frac{\partial \log Z}{\partial V}$$

$$\epsilon V=-\frac{\partial \log Z}{\partial (1/T)}.$$

$$p{a}^4=2{\int }_{\mathrm{cold}}^{6/g^2}[⟨\square ⟩-{⟨\square ⟩}_{\mathrm{sym}}]d(6/g^{\prime 2}).$$

$$p{a}^4={\int }_{\mathrm{cold}}^{m_{q}a}[⟨\bar{\psi }\psi ⟩-{⟨\bar{\psi }\psi ⟩}_{\mathrm{sym}}]d(m_q^{\prime }a).$$

$$\frac{IV}{T}=-\frac{1}{T}\frac{\partial \log Z}{\partial (1/T)}-3V\frac{\partial \log Z}{\partial V}$$

$$=(-a_t\frac{\partial }{\partial a_t}-a_s\frac{\partial }{\partial a_s})\log Z=-\frac{\partial \log Z}{\partial \log a}$$

$$I{a}^4$$

$$-$$

$$\epsilon =I+3p$$

$$,sT=I+4p.$$

$$\frac{1}{c_s^2}=\frac{d\epsilon }{dp}.$$

$$\frac{1}{c_s^2}=\frac{\frac{d\epsilon }{dT}{|}_{V,m_{\pi }/m_{\rho }}}{\frac{dp}{dT}{|}_{V,m_{\pi }/m_{\rho }}}.$$

$$f(T)=\epsilon (T)-Ts(T)=-p$$

$$\epsilon (T)=Ts(T)-p(T).$$

$$\frac{d\epsilon }{dT}=T{s}^{\prime }(T)+s(T)-p^{\prime }(T)=T{s}^{\prime }(T),$$

$$s=\frac{\partial p}{\partial T}{|}_V=p^{\prime }(T).$$

$$\frac{pV}{T}=\log Z=\frac{-fV}{T}.$$

$$\frac{1}{c_s^2}=\frac{T{s}^{\prime }(T)}{p^{\prime }(T)}=\frac{d\log (s)}{d\log T},$$

$$\frac{1}{c_s^2}=\frac{1}{s}\frac{dI}{dT}+3.$$

$$\frac{1}{c_s^2}=3-\frac{1}{s}\frac{\partial [\frac{T}{V}\frac{\partial \log Z}{\partial \log a}]}{\partial T}.$$

$$\frac{1}{c_s^2}=3-\frac{1}{s}\left[\frac{1}{V}\frac{\partial \log Z}{\partial \log a}+\frac{T}{V}\frac{\partial [\frac{\partial \log Z}{\partial \log a}]}{\partial T}\right].$$

$$\frac{\partial [\frac{\partial \log Z}{\partial \log a}]}{\partial T}=\frac{1}{T}\frac{\partial [\frac{\epsilon V}{T}]}{\partial \log a},$$

$$\frac{\partial }{\partial \log a}$$

$$\frac{\partial }{\partial T}$$

$$\frac{\epsilon V}{T}=-a_t\frac{\partial \log Z}{\partial a_t}.$$

$$\frac{1}{c_s^2}=3+\frac{1}{sT}\left[I-\frac{T}{V}\frac{\partial [\frac{\epsilon V}{T}]}{\partial \log a}\right].$$

$$\frac{1}{c_s^2}=\frac{1}{\epsilon +p}\left[4\epsilon -\frac{T}{V}\frac{\partial [\frac{\epsilon V}{T}]}{\partial \log a}\right].$$

$$\frac{1}{c_s^2}=\frac{4\epsilon }{\epsilon +p}=3,$$

$$\frac{1}{c_s^2}=\frac{1}{\epsilon a^4+p{a}^4}[4\epsilon a^4-(\frac{\partial (6/g^2)}{\partial \log a})\frac{\partial \epsilon a^4}{\partial (6/g^2)}-$$

$$(\frac{\partial m_{q}a}{\partial \log a})\frac{\partial \epsilon a^4}{\partial (m_{q}a)}].$$

$$\frac{\partial ⟨\square ⟩}{\partial (a{m}_q)}=⟨\square \bar{\psi }\psi ⟩-⟨\square ⟩⟨\bar{\psi }\psi ⟩.$$

Full text

1 INTRODUCTION

In order to show the existence of the Quark-Gluon Plasma (QGP) in the aftermath of upcoming 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 $N_{t}=4$ [1], and we are now working at $N_{t}=6$. A similar program for the pure gauge theory is being pursued by the Bielefeld group[2]. The $N_{t}=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 ($\Delta t$) errors induced by the approximate integration of the gauge field equations of motion are much more significant at $N_{t}=6$ than at $N_{t}=4$. They are handled by extrapolation of observables to $\Delta 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 calculated and together with the pressure it yields the energy density.

2 THEORY

A Euclidean $N_{s}^{3}\times N_{t}$ lattice with periodic boundary conditions has a temperature $T$ and volume $V$

[TABLE]

where $a$ is the lattice spacing. Thermodynamic variables are derivatives of the partition function $Z$. In particular, the pressure $p$ and energy density $\varepsilon$ are given by

[TABLE]

and

[TABLE]

The methods for computing the pressure and integration measure are discussed in Ref. [1]. The pressure is found by integrating either the plaquette or $\bar{\psi}\psi$.

[TABLE]

[TABLE]

The interaction measure is given by

[TABLE]

where $a_{s}$ and $a_{t}$ are the spatial and temporal lattice spacings. The scale dependence in the case of QCD with dynamical quarks leads to

[TABLE]

The derivatives are the usual $\beta$ function and the anomalous dimension of the quark mass. In the above the subscript “sym” refers to symmetric lattices with $N_{t}=N_{s}$ which are used for vacuum subtraction.

The knowledge of the non-perturbative pressure and interaction measure allows us to compute other bulk quantities. The energy and entropy $s$ become

[TABLE]

3 SOUND SPEED

Acoustic perturbations travel in the system with a speed $c_{s}$:

[TABLE]

One has to take the derivative keeping the physical quark mass fixed, i.e. on the line of constant physics. Unfortunately, we know best only the variations of the energy density and pressure along lines of constant bare parameters. In order to measure the correct sound speed one has to use the $\beta$ 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

[TABLE]

Henceforth we drop the $|_{m_{\pi}/m_{\rho}}$ reminder. With only $T$ varying, the fundamental relation of thermodynamics becomes

[TABLE]

or

[TABLE]

Then,

[TABLE]

where we have utilized a Maxwell relation for the entropy $s$:

[TABLE]

In the Maxwell relation we interchanged pressure and free energy according to

[TABLE]

Hence,

[TABLE]

where we have used Eq. (14) again.

However, lattice simulations are always done at finite volume, and the volume varies with the lattice 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 spacing and leads to expressions with asymmetry coefficients. 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:

[TABLE]

On the other hand, Eq. (6) allows us to write

[TABLE]

Now, the $T$ derivative is taken at constant $V$. Therefore,

[TABLE]

The last derivative can be given as

[TABLE]

which can be seen by expanding the derivatives with respect to $a$ and $T$ as derivatives with respect to $a_{t}$ and $a_{s}$:

[TABLE]

We also used the definition of energy density, Eq. (3), in the form

[TABLE]

Then Eq. (19) becomes

[TABLE]

Using $sT=\varepsilon+p$ and $I=\varepsilon-3p$ this can be expressed in a rather compact form

[TABLE]

If the system is conformally invariant the energy density does not depend on the scale and the derivative term in (24) disappears:

[TABLE]

the relativistic result for a free gas.

For QCD with dynamical quarks an explicit form is

[TABLE]

The derivatives have to be taken on a line of constant physics.

4 SIMULATIONS

We have measured $\left\langle\Box\right\rangle$ and $\left\langle\bar{\psi}\psi\right\rangle$ on asymmetric lattices with $N_{t}=6$ and $N_{s}=12$, and on symmetric lattices with $N_{t}=N_{s}=12$.

The couplings in the simulations are appropriate for the $N_{t}=6$ crossover [3]. At $am_{q}=0.0125$ and 0.025 we varied $6/g^{2}$ 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/g^{2}=5.45$ and 5.53, we have varied $am_{q}$ between 0.01 and 0.1 and 0.0125 and 0.2, respectively. These runs also allow us to extrapolate to zero quark mass, but in addition they provide a cross check on the integrations over $6/g^{2}$ 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 dynamics R algorithm [4]. For the $N_{t}=6$ lattices we ran at least 1800 trajectories after 200 trajectories 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 observables 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 observable. However, the errors are different on the symmetric and asymmetric lattices, so they do not cancel from the vacuum subtractions. Moreover, in many instances the errors are the same order of magnitude as the vacuum subtracted quantities. Therefore, they must be eliminated.

In most cases, the step size errors are eliminated 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 observed, 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 approximately equal to $am_{q}$. Roughly speaking, this is because in the integration of the gauge momentum, the fermion “force” to lowest order is proportional to $1/am_{q}$. Thus the step taken in simulation time to update the momentum should be $\simeq am_{q}$ (or smaller) to keep the change in the momentum less than ${\cal O}(1)$. Thus we have used $0.007\leq\Delta t\leq 0.015$ and $0.015\leq\Delta t\leq 0.03$ for the runs with $am_{q}=0.0125$ and 0.025, respectively. For the larger quark mass runs at $6/g^{2}=5.45$ and 5.53, we used $0.02\leq\Delta t\leq 0.03$. Finally, for the smallest quark mass simulation, $am_{q}=0.01$ $(6/g^{2}=5.45)$, we use $\Delta t=0.005$.

5 EXTRAPOLATIONS

In Fig. 1 we show the plaquette dependence on step size squared for $am_{q}=0.0125$. Similar results hold for $am_{q}=0.025$. Generally, the effects are larger on the symmetric lattices and at smaller $6/g^{2}$. For example, at $6/g^{2}=5.39$, linear behavior sets in for much smaller step size on the cold lattice than on the hot lattice. In fact, for $\Delta t^{2}\geq 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 procedure. First, on the symmetric lattices for $5.39\leq 6/g^{2}\leq 5.43$, use only the smallest two step sizes to extrapolate to $\Delta t^{2}=0$. At $6/g^{2}=5.37$ we effectively assume the system is in the cold phase and take the values at $\Delta t=0.007$ as the zero step size extrapolations. For $6/g^{2}>5.43$ we use all $\Delta t$ values to do the extrapolations. On the hot lattices we do the following. For each $6/g^{2}$ where there are multiple step sizes, we use all of them to extrapolate to $\Delta t^{2}=0$. For $6/g^{2}=5.40$ and 5.42 we interpolate the slope from the neighboring points and use that along with the point at $\Delta t=0.01$ to extrapolate to zero step size. Note, for $6/g^{2}\geq 5.43$, the slopes are essentially zero. Therefore, for $6/g^{2}>5.43$ where only one measurement is available, we take that value as the zero step size value.

For $\left\langle\bar{\psi}\psi\right\rangle$ the situation is similar to the plaquette with the following exceptions. First, after vacuum subtraction, the relative step size errors are not as significant. Second, on the hot lattices for $6/g^{2}>5.39$, we find the slopes with respect to $\Delta t^{2}$ are zero within errors, so we take the zero step size $\left\langle\bar{\psi}\psi\right\rangle$ to be the value measured at the smallest $\Delta t$ for each of these couplings.

For $am_{q}=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/g^{2}$, we made cold runs at values of $6/g^{2}$ separated by $0.02$, while the hot runs were separated by $\Delta 6/g^{2}=0.01$ near the transition. Cold observables at the other couplings were obtained by interpolation. For $6/g^{2}\leq 5.47$, quadratic fits to the zero step size plaquette and $\bar{\psi}\psi$ had $\chi^{2}=2.69$ and 3.18 with three degrees of freedom respectively. These fits were used for the interpolated values. To get the symmetric observables at $6/g^{2}=5.53$, one can either extrapolate in $6/g^{2}$ for fixed $am_{q}$, or extrapolate in $am_{q}$ for fixed $6/g^{2}$. Both give the same result within errors, and we simply take the value from extrapolation in $am_{q}$ as it has a much smaller error.

6 PRESSURE

The results for $\langle\bar{\psi}\psi\rangle$ and $\langle\Box\rangle$ from the previous section can now be integrated to yield the pressure.

To begin consider $\langle\bar{\psi}\psi\rangle$ as a function of $am_{q}$. Using Eq. (5), we find the pressure as a function of $am_{q}$ (at $6/g^{2}=5.45$ and 5.53). The result is shown in Fig. 2. We also want the pressure at $am_{q}=0$ which is found by setting $\bar{\psi}\psi(0)=0$ and continuing the integration to $am_{q}=0$. At $6/g^{2}=5.45$, a linear fit to the data for $am_{q}\leq 0.025$ that is constrained to go through the origin has $\chi^{2}=4.9$ for three degrees of freedom. For $6/g^{2}=5.53$, we only have measurements at $am_{q}=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 origin has $\chi^{2}=0.49$ for one degree of freedom. The pressure at $am_{q}=0$ calculated in this manner is also shown in Fig. 2.

Next we integrate the plaquette with respect to $6/g^{2}$ to obtain the pressure as a function of $6/g^{2}$ at fixed $am_{q}$. The result for $am_{q}=0.0125$ is shown in Fig. 3. The pressure rises smoothly through the crossover region as it must. The results from the quark mass integrations are also shown in Fig. 3, and are in agreement with the $6/g^{2}$ integration. This is a good check on our analysis, since for the most part the two integrations are independent. In particular, the integration with respect to $6/g^{2}$ is sensitive to the step size extrapolations while the $am_{q}$ integration is not.

7 INTERACTION MEASURE

The interaction measure is given by Eq. (2). For the $\beta$ function we use our previous result calculated from the $\pi$ and $\rho$ masses at various values of $6/g^{2}$ and $am_{q}$ [1]. The result is shown in Fig. 4. It rises sharply from zero in the cold phase to some maximum and then begins to drop off. For a noninteracting plasma, $I$ is zero since $\varepsilon=3p$ which is certainly not the case at the largest value of $6/g^{2}$ in our simulations. From asymptotic freedom we expect the system to approach a noninteracting plasma at very high temperature.

The interaction measure at $am_{q}=0$ depends only on the plaquette since the anomalous dimension of the quark mass is zero at $am_{q}=0$. To extrapolate the plaquette to zero quark mass we use the fact that its slope is just its correlation with $\bar{\psi}\psi$,

[TABLE]

Since $\bar{\psi}\psi$ is discontinuous at the origin in the broken phase and continuous in the chirally symmetric phase, we expect a cusp at the origin for the cold lattices and zero slope for the hot lattices[1]. At $6/g^{2}=5.45$ a linear fit on the cold lattices for $am_{q}\leq 0.05$ gives $\chi^{2}=1.6$ with two degrees of freedom. On the hot lattice a quadratic fit constrained to zero slope at the origin gives $\chi^{2}=0.74$ with two degrees of freedom. At $6/g^{2}=5.53$ the situation is less satisfactory. For the cold lattices, a quadratic fit to the data with $am_{q}\leq 0.1$ has $\chi^{2}=0.09$. Note the smallest quark mass here is 0.025, and we also use this fit for the cold observables at $am_{q}=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/g^{2}=5.53$ was observed in our earlier $N_{t}=4$ study, which may indicate a systematic error. Naively one expects finite volume 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 $N_{t}=6$ equation of state as a function of $6/g^{2}$ for $am_{q}=0.0125$. The zero quark mass extrapolations are also shown. At the largest value of $6/g^{2}$, $\varepsilon$ is still much larger than $3p$; for $am_{q}=0$, $3p$ is approximately $75\%$ of $\varepsilon$. On the other hand, after a rapid rise, $\varepsilon/T^{4}$ more or less levels off for $6/g^{2}\geq 5.43$. This corresponds to an energy density of about 1 ${\rm GeV/fm^{3}}$ at $6/g^{2}=5.43$ (we use the $\rho$ mass to convert $6/g^{2}$ to temperature). $m_{\pi}/m_{\rho}$ at this point is slightly less than 0.4.

In Fig. 7 we compare the $N_{t}=4$ and 6 equations of state. There is a large finite size effect as expected from the free field results on finite lattices (also shown). For example, $p/T^{4}$ for $am_{q}=0$ differs by about $15\%$. There appears to be a sizable quark mass effect as well (the $N_{t}=4$ results are for $am_{q}=0.1$ and $0.025$). For $N_{t}=6$ the approach to the Stefan-Boltzmann law is unclear since we do not have data at very high temperatures. The prominent peak in $\varepsilon/T^{4}$ for $N_{t}=4$ just after the crossover is much smaller at $N_{t}=6$. We have plotted the results versus temperature by using the $\rho$ mass to set the scale. The crossover temperature is around 150 MeV and is rather insensitive to $N_{t}$ and $am_{q}$. In a recent preprint, Asakawa and Hatsuda have pointed out that many features of this equation of state are constrained by fundamental thermodynamic relations[5].

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 discussed 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 accurately as possible we performed a set of simulations at $am_{q}=0.09$ with $N_{t}=4$ lattices in addition to our old data at $am_{q}=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 $N_{t}=6$ lattices. The algorithm used to generate gauge configurations introduces a step size error in observables that must be removed by extrapolation. The added computational cost is significant. For $am_{q}=0.0125$, we find the energy density just after the crossover to be roughly 1 ${\rm GeV/fm^{3}}$, and for $T$ almost twice the critical value, the highest that we simulated, three times the pressure is only $60\%$ of the energy density. We find large effects of nonzero lattice spacing from comparison with results at $N_{t}=4$, as expected from free field theory.

After the completion of runs at $am_{q}=0.025$, we will complete our extrapolation to zero quark mass. It still remains to eliminate remaining lattice 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 Center, and Indiana University.