BI-TP 2003/11

SWAT/03/374

The Equation of State for

Two Flavor QCD

at Non-zero Chemical Potential

C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. Schmidt

Department of Physics, University of Wales Swansea, Singleton Park,

Swansea SA2 8PP, U.K.

Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany.

ABSTRACT

We present results of a simulation of QCD on a lattice with 2 continuum flavors of p4-improved staggered fermion with mass . Derivatives of the thermodynamic grand potential with respect to quark chemical potential up to fourth order are calculated, enabling estimates of the pressure, quark number density and associated susceptibilities as functions of via Taylor series expansion. Discretisation effects associated with various staggered fermion formulations are discussed in some detail. In addition it is possible to estimate the radius of convergence of the expansion as a function of temperature. We also discuss the calculation of energy and entropy densities which are defined via mixed derivatives of with respect to the bare couplings and quark masses.

May 2003

## 1 Introduction

Non-perturbative studies of QCD thermodynamics with small but non-zero baryon charge density by numerical simulation of lattice gauge theory have recently made encouraging progress [1] - [4]. In particular, it has proved possible to trace out the pseudocritical line separating hadronic and quark-gluon plasma (QGP) phases in the plane out to MeV [2] - [4], where the quark chemical potential is the appropriate thermodynamic control variable in the description of systems with varying particle number using the grand canonical ensemble. In addition, the first estimate has been made of the location along this line of the critical endpoint expected for light quark flavors, where the crossover between hadron and QGP phases becomes a true first order phase transition [2]. As well as being of intrinsic theoretical interest, such studies are directly applicable to the regime under current experimental investigation at RHIC, where corrections to quantities evaluated at are both small and calculable. In this respect it is worth reminding the reader that in a relativistic heavy ion collision of duration s, thermal equilibration is possible only for processes mediated by the strong interaction, rather than the full electroweak equilibrium achievable, say, in the core of a neutron star. This means that each quark flavor is a conserved charge, and conditions at RHIC are thus approximately described by

(1) |

with MeV [5] when we relate the quark and baryon number chemical potentials via . In this paper we will present numerical results for the equation of state, i.e. pressure and quark number density , obtained from a lattice QCD simulation with , which should give a qualitatively correct description of RHIC physics, and provide a useful warm-up exercise for the physical case of 2+1 flavors with realistic light and strange quark masses.

Direct simulation using standard Monte Carlo importance sampling is hampered because the QCD path integral measure , where is the Euclidean space fermion kinetic operator, is complex once . In the studies which have appeared to date, two fundamentally distinct approaches to this problem have emerged. In the reweighting method results from simulations at are reweighted on a configuration-by-configuration basis with the correction factor yielding formally exact estimates for expectation values. Indeed, it is found that if reweighting is performed simultaneously in two or more parameters, convergence of this method on moderately-sized systems is considerably enhanced [6]. This method has been used on lattice sizes up to with to map out the pseudocritical line and estimate the location of the critical endpoint at MeV, MeV [2]. More recently the equation of state in the entire region to the left of the endpoint has been calculated this way [7]. However, it remains unclear whether the thermodynamic limit can be reached using this technique.

Analytic approaches use data from regions where direct simulation is possible, either by calculating derivatives with respect to (or more properly with respect to the dimensionless combination ) to construct a Taylor expansion for quantities of interest [1, 3, 8, 9], or more radically by analytically continuing results from simulations with imaginary (for which the integration measure remains real) to real . The second technique has been used to map for QCD with both [4] and [10], in the latter case finding evidence that the line is first order in nature. Fortunately, the pseudocritical line found in [4] is in reasonable agreement with that found by reweighting; moreover the radius of convergence within which analytic continuation from imaginary is valid corresponds to [11]. The leading non-trivial term of quadratic order in the Taylor expansion appears to provide a good approximation throughout this region. In general though, while analytic approaches have no problem approaching the thermodynamic limit, it is not yet clear if and how they can be extended into the region around the critical endpoint (but see [12]), and to observables that vary strongly with like eg. the pressure or energy density.

In our previous paper [3] we used a hybrid of the two techniques, by making a Taylor series estimate of the reweighting factor to . Since this is considerably cheaper than a calculation of the full determinant, we are able to explore a larger system, and also exploit an improved action in both gauge and fermion [14] sectors, thus dramatically reducing discretisation artifacts on what at MeV is still a coarse lattice. Our results yield a curvature of the phase transition line , consistent with the other approaches. Although our simulation employs quark masses , not yet realistically light, the results also suggest that any dependence on quark mass is weak. In the current paper we extend the Taylor series to the next order but this time remain entirely within the analytic framework, using derivatives calculated at to evaluate non-zero density corrections to the pressure and quark number susceptibility , as well as the quark number density itself. In fact, since the correction can be evaluated at fixed temperature, it turns out to be considerably easier to calculate than the equation of state at [15, 16]. Since we now have the first two non-trivial terms in the Taylor expansion, we are also able to estimate its radius of convergence as a function of , and confirm that close to the results of our previous study for the critical line curvature can be trusted out to , whereas at higher temperatures a considerably larger radius of convergence is likely to be found. Finally we consider mixed derivatives with respect to both and the other bare parameters and , which are required to estimate energy and entropy densities. Due to the presence of a critical singularity, these latter quantities appear considerably harder to calculate in the critical region using this approach.

Sec. 2 outlines the formalism used in the calculation and specifies which derivatives are required. In Sec. 3 we present a calculation of the cutoff dependence of the terms in the expansion for standard, and for Naik- and p4-improved staggered lattice fermions, showing that both improvements result in a dramatic reduction of discretisation effects. Our numerical results are presented in Sec. 4, and a brief discussion in Sec. 5. Two appendices contain further details on the calculation of the required derivatives and the cutoff dependence.

## 2 Formulation

In the grand canonical ensemble pressure is given in terms of the grand partition function by

(2) |

Note that we have been careful to express both sides of this relation in dimensionless quantities. Since the free energy and its derivatives can only be calculated using conventional Monte Carlo methods at , we proceed by making a Taylor expansion about this point in powers of the dimensionless quantity :

(3) | |||||

where derivatives are taken at . Note that calculating is considerably easier than itself, because whereas must be estimated by integrating along a trajectory in the bare parameter plane [15, 16], its derivatives can be related to observables which are directly simulable at fixed , where is the gauge coupling parameter and the bare quark mass. Only even powers appear in (3) because as shown in [3], odd derivatives of the free energy with respect to vanish at this point. Note also that we will work throughout with fixed bare mass, implying that our computation of is strictly valid along a line of fixed .

For QCD with staggered fermions the partition function may be written

(4) |

with denoting the gauge field variables, the link action and the kinetic operator describing a single staggered fermion, equivalent to continuum flavors. On a lattice of size with physical lattice spacing , so that , we define a dimensionless lattice chemical potential variable . Equation (3) then becomes

(5) |

The derivatives may be expressed as expectation values evaluated at :

(6) |

(7) | |||

All expectation values are calculated using the measure and in deriving (6) and (7) we used the fact that for odd. To evaluate these expressions we need the following explicit forms:

(8) | |||

(9) | |||

(10) | |||

(11) | |||

The traces can be estimated using the stochastic method reviewed in [3]. Since is a local operator, expressions containing powers of require operations of matrix inversion on a vector.

Next we discuss the quark number density and its fluctuations. Starting from the Maxwell relation

(12) |

where is the thermodynamic grand potential and the net quark number, we can write an equation for the quark number density analogous to (5):

(13) |

It is also possible to interpret derivatives of with respect to in terms of the various susceptibilities giving information on number density fluctuations [8, 17]. We define quark number , isospin and charge susceptibilities as follows:

(14) | |||||

(15) | |||||

(16) |

Quark and baryon number susceptibilities are related by . Any difference between and is due to correlated fluctuations in the individual densities of and quarks. With the choice , , which approximates the physical conditions at RHIC, can then be expanded in terms of quantities already used in the calculation of and :

(17) |

whereas the expansion of is determined by the following expectation values:

(18) | |||

(19) |

where we have explicitly set . Charge fluctuations are then given by the relation

(20) |

where the third term vanishes for , .

Finally we discuss the energy density , most conveniently extracted using the conformal anomaly relation

(21) |

where and are the bare coupling and quark mass respectively. In fact, for the derivation of this expression needs careful discussion. Start from the defining relation

(22) |

where is entropy. For a Euclidean action defined on an isotropic lattice of spacing we have the identity

(23) |

It follows that

(24) | |||||

(25) |

implying

(26) |

where we have allowed for the dependence of the lattice action on all bare parameters. Since however , and a parameter multiplying a conserved charge experiences no renormalisation, the third terms on each side cancel leaving the relation (21).

Taylor expansion of (21) about leads to the expression

(27) | |||||

The beta function may be estimated by measurements of observables at ; the factor is poorly constrained by current lattice data but vanishes in the chiral limit, so is frequently neglected. In order to assess the magnitude of the resulting error, it is nonetheless useful to calculate all the derivative terms. They may be estimated using the formulæ

(28) | |||||

(29) |

The derivative is, of course, simply the combination of plaquettes comprising the gauge action itself, and derivatives with respect to can be evaluated using

(30) |

The implementation of the second square bracket in (27) in terms of lattice operators is straightforward but unwieldy; for reference the non-vanishing terms are listed in Appendix A.

## 3 Analyzing the cut-off dependence

In this section we discuss the influence of a non-zero chemical potential on the cut-off effects present in calculations of bulk thermodynamic observables on a lattice with finite temporal extent . For this issue has been discussed extensively for both gluonic and fermionic sectors of QCD. In particular, it has been shown that the use of improved actions is mandatory if one wants to ensure that discretisation errors in the calculation of quantities like the pressure or energy density are below the 10% level on moderately sized lattices [14]. We now want to extend these considerations to the case , which affects the quark sector only. Following [14] we will concentrate on an evaluation of the pressure. As we will be evaluating thermodynamic quantities using a Taylor expansion in we want to understand the cut-off dependence of and its expansion coefficients in terms of .

In the limit of high temperature or density, due to asymptotic freedom thermodynamic observables like or are expected to approach their free gas, i.e. Stefan-Boltzmann (SB) values. In this limit cut-off effects become most significant as the relevant momenta of partons contributing to the thermodynamics are and thus of similar magnitude to the UV cut-off . Short distance properties thus dominate ideal gas behaviour and cut-off effects are controlled by the lattice spacing expressed in units of the temperature, .

In the continuum the pressure of an ideal gas of quarks and anti-quarks is given by

(31) |

with the fugacity and the relativistic single particle energies . For massless quarks one finds from an evaluation of the integral the pressure as a finite polynomial in :

(32) |

For non-zero the pressure is a series in the fugacity:

(33) |

where is a Bessel function. Of course, Eq. (33) can also be reorganised as a power series in .

It is well known that the straightforward lattice representation of the QCD partition function in terms of the standard Wilson gauge and staggered fermion actions leads to a systematic cut-off dependence of physical observables. In the infinite temperature limit this gives rise to deviations of the pressure from the SB value (32);

(34) |

Using improved discretisation schemes it is possible to ensure that cut-off effects only start to contribute at [13], or to considerably reduce the magnitude of the coefficient relative to the standard discretisation scheme for staggered fermions [14].

For the pressure of free staggered fermions on lattices with infinite spatial volume () but finite temporal extent is given by

(35) | |||||

In the first term the sum runs over all discrete Matsubara modes, i.e. , whereas in the second term we have an integral over which gives the vacuum contribution. For quarks of mass the function is given by . Here we have introduced functions to describe the momentum dependence of the propagator for the standard, Naik [13] and p4 staggered fermion actions [14]:

(36) | |||||

(37) | |||||

(38) |

The introduction of a non-zero chemical potential is easily achieved by substituting every temporal momentum by . The integrals in (35) have been evaluated numerically for different . Results for different values of and are shown for the different fermion actions in Fig. 1.

For the standard action cutoff effects remain % out to , whereas both improved actions are hard to distinguish from the continuum result by . We note that lines for different values but the same quark mass fall almost on top of each other. Cutoff effects are thus almost independent of . The effect of on the cutoff dependence of the pressure is even smaller than the effect of quark mass .

As can be seen from (32) for moderate values of the -dependence of the continuum ideal gas pressure is dominated by the leading contribution. In order to control the cut-off dependence of the various expansion terms we have expanded the integrand of (35) up to order . For the standard action the series starts with

(39) | |||||

Here we use the shorthand notation . The remaining orders as well as the series for Naik and p4 actions are given in Appendix B. A common feature of these expansions is that the odd terms are pure imaginary and the integral and sum over of those terms vanish due to a factor which always appears. To be more precise, this factor always forms the pattern which can be shown to vanish, either after summation over the discrete set of values, or integration from to , for . Performing the momentum integration and the summation over Matsubara modes we obtain the coefficients of the expansion of the pressure;

(40) |

We checked numerically that with increasing the coefficients , and do indeed converge to their corresponding SB values,

(41) |

Fig. 1b shows , and for the standard, Naik and p4 actions with massless quarks, normalized to the corresponding SB value.

We see here again that the cutoff dependence of the pressure at is qualitatively the same as at .

For massless quarks the coefficient should vanish with increasing , as checked in Fig. 2a. It is expected that this term will approach zero like in the large limit. In order to define the numerical factor, we plot over . A fit yields for the standard action. This is at least an order of magnitude larger than for the p4 improved action, for which the dominant cut-off dependence seems to be as for the Naik action.

In the case of massive quarks the expansion (40) no longer terminates at . After expanding (31) in terms of and performing a numerical integration we find for the expansion coefficients up to the values given in Table 1.

0.967 | 0.827 | |||

0.976 | 0.863 | |||

0.999 | 0.976 | |||

— | — |

The mass value is the value we use in our numerical calculations, corresponding to and . We note that the coefficient no longer vanishes. As shown in Fig. 2b, for finite there are large deviations from the continuum value. Even at , however, the absolute value of this coefficient is still a factor of about smaller than the leading term . The deviations thus do not show up in the calculation of the complete expression for the pressure shown in Figure 1a. These terms, however, become more important in higher derivatives of the pressure such as the quark number susceptibility . In summary, for a gas of free quarks we find that the expansion up to is sufficient for and . In the continuum the deviation from the full expression over this range is smaller than 0.01%. On the lattice, however, cut-off effects lead to deviations of approximately 10% on coarse () lattices.

## 4 Numerical Results

We applied the formalism outlined in Sec. 2 to numerical
simulations of QCD with quark flavors on a lattice,
using both Symanzik improved gauge and p4-improved
staggered fermion actions. The simulation method is exactly as presented in
[3] ^{a}^{a}aThe coefficient of the knight’s move hopping term
was incorrectly reported to be 1/96 in [3]; its correct value is 1/48..
The bare quark mass was for which the
pseudocritical point for zero chemical potential is estimated to be
. In order to cover a range of temperatures on
either side of the critical point we examined 16 values in the
range . The simulation employed a hybrid molecular dynamic
‘R’-algorithm with discrete timestep ,
and measurements were performed on equilibrated
configurations separated by . In general for each 500 to
800 configurations were analysed, with 1000 used in the critical region
. On each configuration
50 stochastic noise vectors were used to estimate the
required fermionic operators. For each noise vector, 7 matrix inversions
are required to estimate the required operators
(8-11)
and (49-52).

Following the procedure used for the equation of state at [15], we translate to physical units using the following scaling ansatz [18]:

(42) |

where , being the two-loop perturbative scaling function appropriate for two light flavors. Using string tension data at the best values for the fit parameters corresponding to a reference are , [15]. We find that our simulations span a temperature range , where is the critical temperature at .

In Fig. 3 we show the first two coefficients, (a) and (b) , of the Taylor expansion of introduced in (3) as functions of . Also shown are the corresponding SB limits: (a) and (b) , where the coefficients are defined in (40), with values relevant for both the lattice used () and the continuum limit () plotted. Both and vary sharply in the critical region, but except in the immediate vicinity of the transition the quadratic term dominates the quartic. This is consistent with the results of [7] where data at varying obtained by reweighting was found lie on an almost universal curve when plotted as a fraction of the SB prediction. The asymptotic value of appears to be approached from above.

A notable feature is that in the high- limit our data lies closer to the continuum SB prediction rather than their values corrected for lattice artifacts, assuming 80% of the continuum value for whereas is almost coincident with its continuum value. By contrast recent calculations with unimproved staggered fermions [7, 9] find that the high- limit of the data lie close to the lattice-corrected SB value. This situation can be modelled by making the coefficient of the correction appearing in (34) temperature dependent. In thermodynamic calculations performed with pure unimproved SU(3) lattice gauge theory [19], where extrapolations to the continuum limit are currently practicable, it is found that for , becoming even smaller closer to . The behaviour of and we have observed using p4 fermions is broadly consistent with this behaviour.

In Fig. 4 we plot defined in
(3) as a function of
for various temperatures.
In most cases
and the
relation thus almost linear. The strongest departures from linearity are for
, but even here the quadratic term is dominant for
, corresponding to MeV.
Given enough terms of the Taylor expansion in , one could determine
its radius of convergence via^{b}^{b}bThe argument of [4] that
due to the presence of a phase transition as imaginary
chemical potential is increased beyond this value [20] does not hold
for calculations with real; in this case the pressure
corresponding to the unit element of the Z(3) sector is always the maximum and
hence dominates the partition function in the thermodynamic limit – hence the
issue concerns the analytic properties
within this physical unit sector.

(43) |

Data from the pressure at [15] and the current study enable us to plot the first two estimates and on the plane along with the estimated pseudocritical line found in [3] in Fig. 5. Also shown are the corresponding values from the SB limit (32). For one finds that increases markedly as increases from 0 to 2; if the SB limit is a good predictor for the QGP phase we might expect to be very small, and the next estimate correspondingly very large in this regime. Close to the transition line, however, the thermodynamic singularities appear to restrict ; this in turn gives an approximate lower bound for the position of the critical endpoint. From the figure we deduce , not inconsistent with the result of [2]. The new results at are important because they justify in retrospect our neglect of fourth order reweighting factors in our earlier calculation of the critical line [3]. Indeed, simulations with imaginary suggest that neglect of these terms in the analytic continuation to physical is justified for MeV [4].

Fig. 7 plots the dimensionless correction to the equation of state as a function of both and . In the latter case the correction rises steeply across the transition and peaks for , before rapidly approaching a form characteristic of the SB limit, with the coefficient having 82% of the continuum SB value. Comparison with the equation of state results at from Ref. [15] suggest that the correction will give a significant correction to the pressure for , , but will decrease in importance as rises further. The curves of Fig. 7b are in good qualitative agreement with those of Refs. [7, 9], although we consider any quantitative agreement to be somewhat accidental as the numerical data obtained in [7, 9] with unimproved actions have large discretisation errors which have been corrected for by renormalising the raw data with the known discretisation errors in the infinite temperature limit. Experience gained in calculations of thermodynamic quantities in the pure SU(3) gauge theory suggests that in the temperature range of a few times this procedure overestimates the importance of cut-off effects by a factor two or so [19].

Fig. 7a shows the quark number density evaluated using Eqn. (13). As increases, rises steeply as the QGP phase is entered; for reference, if the quark number density in nuclear matter is denoted , then the ratio . Our results are numerically very similar to those obtained using exact reweighting in [7], where a mass for the light quark flavors was used. Note that a significant quark mass dependence for was observed in [3], and indeed is present even in the SB limit as described in Section 3; however analysis of the SB limit suggests that the difference between the chiral limit and is about 4%. In Fig. 7b we show the result of eliminating in favour of via