Bounding the QCD Equation of State with the Lattice

The equation of state of QCD matter at high densities is relevant for neutron star structure and for neutron star mergers and has been a focus of recent work. We show how lattice QCD simulations, free of sign problems, can provide an upper bound on the pressure as a function of quark chemical potentials. We show that at large chemical potentials this bound should become quite sharp; the difference between the upper bound on the pressure P-phase-quenched and the true pressure P is of order alpha^3 P. The corrections arise from a single Feynman diagram; its calculation would render remaining corrections of order alpha^4 P.


Introduction
The equation of state (EOS) of strongly interacting matter dictates the thermodynamics of any system ultimately composed of quarks and gluons.At high temperatures and low net baryon densities, the EOS can be computed directly from the partition function of Quantum Chromodynamics (QCD) using Monte-Carlo lattice techniques [1,2] and compared to experimental determinations of thermodynamic properties [3].However, at low temperatures and high net baryon densities, such techniques fail due to the well-known sign problem [4,5], and alternative methods are required to determine the EOS.
In this paper we argue that one more constraint (or family of constraints) can be added to this list of bounds on the dense QCD EOS.While lattice QCD cannot compute the EOS for the physical combination of chemical potentials, phase-quenched lattice QCD is free of sign problems.The pressure as a function of chemical potentials, calculated using phase quenching, P PQ (µ q ), is a strict upper bound on the true pressure at the same values of quark chemical potentials: P PQ (µ q ) ≥ P (µ q ).This fact was established by Cohen in Ref. [67], who used it to show how the case of an isospin chemical potential -opposite up and down quark chemical potentials -can be used to bound the more physically interesting case of equal up and down quark chemical potentials.Unfortunately, the combination of chemical potentials relevant for neutron stars is rather different: the down and strange chemical potentials should be equal µ d = µ s , while the up-quark chemical potential is somewhat smaller to accommodate charge neutrality, µ u < µ d .Bounding the equation of state for this combination of chemical potentials starting from lattice results for isospin chemical potentials requires the use of additional inequalities [68], as recently explored by Fujimoto and Reddy [69].
We argue here that the most effective way to use the lattice to bound the neutron-star equation of state is to perform new phase-quenched lattice simulations using the physical combination of up, down, and strange quark chemical potentials.Unlike an isospin chemical potential, the phase-quenched version of such a combination represents a completely unphysical system.But we expect that it will present the tightest bounds on the physical equation of state.In particular, we show here that, while at low densities the constraint is likely to be very loose, at high densities it should become ever sharper.In fact, we will demonstrate that, in the perturbative region, the relative difference is of order α 3 s .Furthermore, at this order in the coupling expansion, the difference arises from a single Feynman diagram, which we identify. 1n outline of this paper is as follows.In the next section we review the path integral for QCD at finite chemical potential, and how it is related to the phase-quenched version.The section reviews the (quite simple) proof that P PQ (µ q ) ≥ P (µ q ), and discusses a little more how one should interpret the phase-quenched calculation.In Section 3 we show how to represent the (unphysical) phase-quenched theory in Feynman diagrams, and we identify the unique O(α 3 s ) diagram which differs between P PQ (µ q ) and P (µ q ).Section 4 estimates the size of lattice artifacts in evaluating the pressure on the lattice, in order to give guidance for how small the lattice spacing must be at a given, large µ value.We end with a discussion, which lists the most important directions for further work.

Path integral and phase quenching
Here we review the proof that the pressure from the phase-quenched theory is a strict upper bound on the original theory.The proof as we formulate it is due to Cohen [67], though the underlying ideas are older [70].The partition function of QCD at a small finite temperature T = 1/β in a box with a large volume V and with quark chemical potentials Because / D + m q and µ q γ 0 have different γ 5 -Hermiticity [4], each determinant is complex.Phase quenching is the replacement of Z with That is, one uses the absolute values of the determinants, rather than the determinants themselves.Since the Euclidean gluonic action L E (G) is strictly real,3 the integrand is now real and positive and equals the absolute value of the integrand for Z.Since the integral of the absolute value of a function over a positive measure is greater than or equal to the integral of the original complex function, we have where P = ln(Z)/(βV ) is the pressure. 4et us investigate the interpretation of the phase-quenched theory.Applying γ 5 -Hermiticity, one finds that [4,74] Det / D + m q + µ q γ 0 * = Det / D + m q − µ q γ 0 (2.4) and therefore Therefore, the phase-quenched theory is equivalent to a theory with twice as many fermionic species, but where each is represented by the square root of a determinant -think of each as a half-species, which appear in pairs with equal mass but opposite chemical potential.This theory is clearly unphysical, since a half-species of fermion does not make sense as an external state.A well known exception is if two quark masses and chemical potentials are equal [70,[75][76][77].In particular, for m u = m d and µ u = µ d with µ s = 0, representing "standard" baryonic chemical potential, the phase-quenched version is equivalent, after re-labeling two of the half-species, to the case µ d = −µ u , corresponding to an isospin chemical potential.This among other things has motivated rather detailed lattice investigations of the case of finite isospin chemical potential, see for instance Refs.[78][79][80][81][82][83][84][85][86][87][88][89][90][91] or the comprehensive review by Mannarelli [92].In this case it is particularly clear that, for small chemical potentials, the two pressures are quite different.Tom Cohen has analyzed the case for an isospin chemical potential [74] and shown that nonzero µ q in Eq. (2.5) does not change the value of the determinant until µ q = m π /2, when eigenvalues of / D + m q + µ q γ 0 start to cross the imaginary axis.He refers to the constancy of the determinant up to this point as "the silver blaze problem."Applying his analysis for the general case of distinct µ u , µ d , µ s , one finds that P PQ first deviates from its vacuum value when either where m ss is the mass of the unphysical strange-antistrange pseudoscalar obtained by ignoring disconnected contributions to the correlation function.Therefore, at small chemical potentials m π /2 < µ q < m p /3, P PQ differs significantly from the physical pressure.
In the next section we will see that, for large chemical potentials, the two pressures become much more similar.

Perturbation theory
To construct a perturbation theory for Eq.(2.6), as usual for fermions, we rewrite The process of deriving Feynman rules from 1 2 Tr ln( / D + m ± µγ 0 ) is the same as without the 1  2 factor, except that each fermionic loop receives a factor of 1 2 .Therefore, in each place where a fermionic loop can appear, instead of performing a sum over n f fermions, each with mass m q and chemical potential µ q , we sum over 2n f terms; twice for each m q , once with µ q and once with −µ q but each with an overall factor of 1  2 .This is the same as averaging the loop over whether µ is positive or negative.
To see the impact of µ → −µ we remind the reader about Furry's theorem [93].Consider a fermion loop with n external gluons connected by n propagators, with incoming momenta q 1 , . . ., q n with q n = − n−1 j=1 q j by momentum conservation.The fermions are in a representation R (typically the fundamental representation) with Hermitian generators T A .Its group-theory factor and the numerator of its Feynman rule are where p 1 is the loop momentum and p i = p i−1 + q i .Inserting the charge conjugation matrix times its inverse, CC −1 , between each neighboring term and using Reversing the order of the symbols to get rid of the transposes on the gamma matrices, and using T ⊤ = T * for the Hermitian group generators, we find which is the same diagram,5 traversed in the opposite sense, with µ → −µ and with T A → (−T * A ) which is the generator of the conjugate representation -so if the T A generate the fundamental representation, the −T * A generate the antifundamental representation.Another way to state Furry's theorem is then that reversing the sign of µ is equivalent to considering fermions with the original µ value but in the conjugate representation (quarks with −µ are the same as antiquarks with +µ).
For the case of two external gluons we have and the fundamental and antifundamental representations give the same answer.The representations first differ when there are three gluons attached: where d ABC is a totally symmetric symbol which arises in groups of rank 2 or more.

Orientation A Orientation B
Figure 1.The lowest-order diagram which distinguishes between fundamental and antifundamental representations, with two relative orientations for the fermionic loops.
Because both f ABC d ABD = 0 and δ AB d ABC = 0, the first group-theoretical structure where d ABC can give a nonzero contribution, for the group SU(N c ), is This contraction can only appear in a diagram with at least two fermion loops with at least three gluon attachments each.The only such diagram at the α 3 s level is shown in Figure 1.Fermion number can traverse the two loops with two relative orientations, shown on the left and the right in the figure.We will write the two versions of the diagram, stripping off all group-theory factors, for two specific species (i, j), as A(µ i , µ j ) and B(µ i , µ j ).That is, A(µ i , µ j ) and B(µ i , µ j ) represent the value of each fermion orientation in the abelian version of the diagram.In this notation, the contribution of this diagram to the true pressure is In the abelian theory the group-theoretical coefficients on the first and second terms of the second line would be 1 and 0 respectively.Furry's Theorem means that B(µ i , µ j ) = −A(µ i , −µ j ) = −A(−µ i , µ j ).This allows us to rewrite Eq. (3.8) in terms of A only: (3.9) Averaging over µ j ↔ −µ j and/or µ i ↔ −µ i , as we are instructed to do in evaluating P PQ , eliminates the first expression and leaves only the second.Therefore, the difference between the phase-quenched and true pressure, at order α 3 s , is We can show that this combination is automatically positive; the proof, along with an explicit evaluation, will appear in a forthcoming publication.Evaluating this single combination of diagrams would determine the perturbative correction between the phase-quenched and true pressures, through to order α 3 s .Any quark with µ i = 0 does not contribute to Eq. (3.10).The contribution of soft gluon momenta to Eq. (3.10) is suppressed, since the three-gluon hard loop contains only the f ABC group-theory factor [94]. 6 Therefore Eq. (3.10) does not contain any logarithmically enhanced ∼ α 3 s ln(α s ) contributions.
As an aside, we comment on the behavior of these terms as a function of N c .In the large N c or t'Hooft limit, the group-theory factor on Orientation B in Eq. (3.8) is ∝ N 3 c while that for Orientation A is ∝ N c and is suppressed.In contrast, for the group SU(2), the group-theoretical factor on A+B in Eq. (3.8) vanishes and there is no distinction between the fundamental and antifundamental representations.In fact, because every SU(2) representation is equivalent to its conjugate representation, it is easy to show that any combination of chemical potentials is free of sign problems in 2-color QCD, which has motivated investigation of this theory at finite density [95][96][97][98][99][100][101].

Lattice considerations
The previous sections have made clear that a lattice calculation of P PQ would be valuable.Here we want to take one small step towards estimating how difficult such a lattice calculation would be.The larger µa is, the more rapidly a lattice calculation develops statistical power; since P ∝ µ 4 , we expect that the signal-to-noise should scale roughly as (µa) 4 .Therefore it is advantageous to work on lattices with the largest (µa) we can get away with.
But increasing µa increases systematic lattice-spacing effects.In general, for a lattice calculation to be accurate we need the lattice spacing to obey aΛ QCD ≪ 1.But µ introduces an additional scale, and for the regime µ ≫ Λ QCD , we expect to need the stronger condition µa ≪ 1.But what does µa ≪ 1 really mean -that is, how small does µa really need to be?To estimate this, we compute how different the vacuum-subtracted7 pressure P (µ) − P (0) is on the lattice from its continuum value, at lowest (zero) order in the strong coupling.Here we will only consider staggered quarks [102][103][104], since we expect this quark formulation to be used in practical calculations.
In the continuum, one free Dirac fermion with chemical potential µ provides a pressure of where the second term subtracts the µ = 0 "vacuum" pressure.Evaluating the trace, one finds Naively it appears that µ contributes to the pressure at any ⃗ p value including when ⃗ p 2 + m 2 > |µ|, but this is illusory.Separating the numerator and denominator of the log, one may deform the p 0 integration contour for the former, shifting it by −iµ.If p 2 + m 2 > |µ| then this contour deformation encounters no singularities, and the p 0 integral at nonzero and at zero µ are identical.For |µ| > ⃗ p 2 + m 2 there is a cut running from p 0 = 0 to p 0 = −i(µ − ⃗ p 2 + m 2 ) with discontinuity 2πi which the deformed contour must enclose, leading to as expected.The expressions for the energy density ε and for µN the chemical potential times the number density are the same but with the first factor, |µ| − ⃗ p 2 + m 2 , replaced by ⃗ p 2 + m 2 and by |µ| respectively, recovering the thermodynamical relation µN = ε + P .
Naive staggered fermions [104] represent four species of physical fermions with the limitation that each p µ component runs over [−π/(2a), π/(2a)] with a the lattice spacing 8 and with the substitution Here ĵ and 0 represent the unit vector in the j spatial direction and in the time direction respectively.Note that the effects of µγ 0 must be incorporated along with the time derivative because, in the staggered formulation, an insertion of γ 0 requires that the comparison be made between ψ and ψ which differ by an odd number of steps in the time direction.This is in any case the preferred way of introducing a chemical potential because it avoids quadratic-in-µ lattice artifacts, see [106].The fact that a staggered fermion represents four physical fermions is handled through the fourth-root trick and leads to a pressure for one physical fermion of There are two effects.First, the lattice, rather than physical, dispersion determines the energy E p .Second, the p 0 range is finite and periodic and p 0 appears inside a trigonometric function inside the log.Nevertheless, a contour deformation, p 0 → p 0 − iµ, is still possible, and it encounters a similar discontinuity in the log: We have only been able to compute the resulting ⃗ p integral numerically, but as expected, the leading small-µ corrections are of order µ 2 a 2 .Specifically, the flatter dispersion relation lowers E p and leads to a larger phase space region where fermionic states are occupied.
We show the results of a simple numerical evaluation of Eq. (4.6) in Figure 2. We have also checked that carrying out the integral shown in Eq. (4.5) leads to the same result.Unfortunately, if we set as a criterion for good convergence that the lattice and continuum pressures should differ by at most 10%, we are restricted to aµ < 0.42.To do better, we could use an improved fermionic action.Naik has advocated [107] that one replace Eq. (4.4) with As shown in Figure 2, this dramatically improves the match between the lattice and continuum pressure, such that aµ = 1 still has < 10% corrections.However, beyond aµ = 1.147, cosh(3aµ) > 9 cosh(aµ) and the "improvement" term in Eq. (4.7) starts to dominate over the nearest-neighbor term.This leads to additional cuts in the modified version of Eq. (4.5), and the performance of this implementation rather abruptly breaks down.With these results in mind, we feel that nearest-neighbor staggered quarks can only treat large µ accurately out to disappointingly small µa ∼ 0.4 or 0.5 -possibly a little higher with the help of extrapolation over a few lattice spacings.Improved quarks, such as ASQTAD [108,109] or HISQ [110] quarks, which use the Naik term, should be able to do better but it is very dangerous to venture beyond aµ = 1.1.
Another limitation of a lattice treatment is that one generically must compute at finite temperature.Therefore we should also estimate the size of thermal effects.We will assume that the thermal corrections are similar to those in the continuum.At the free theory level, the pressure of a single Dirac fermion with chemical potential µ at temperature T actually has a closed form: This implies µ > 14T to keep thermal effects at the 10% level.One can attempt to remove both lattice-spacing and temperature effects through extrapolation over multiple lattice spacings and box sizes.In the case of lattice spacing effects, an effective field theory analysis at the scale µ tells us that the lattice-spacing effects should scale as (µa) 2 up to anomalous dimension corrections.The coefficient will not necessarily equal the free-theory one, but the effect should definitely be a power law with power close to 2. More care is needed in extrapolating away temperature effects.While one might expect that the leading thermal effects are of order T 2 as in Eq. (4.8), this is really an assumption which can go wrong if, for instance, the theory develops a mass gap due to interactions between excitations near the Fermi surface.Therefore we expect that more care must be taken when extrapolating to small temperature.

Discussion
This paper has made three points.
• Lattice QCD can study arbitrary combinations of µ u , µ d , µ s using phase-quenching, providing P PQ (µ q ).Though this does not return the true QCD pressure, it returns a strict upper bound on the true pressure: P PQ (µ q ) ≥ P (µ q ), which could still be useful in constraining the equation of state for neutron star matter.
• The difference P PQ − P is small at weak coupling, in the sense that P PQ −P P ∝ α 3 s .Furthermore, the α 3 s contribution arises from a single diagram.If we could compute this diagram, we could use P PQ to determine P up to α 4 s corrections (up to logs and nonperturbative effects such as pairing gaps).
• Lattice calculations of P (µ) encounter lattice artifacts.We estimate the size of these artifacts and advocate that nearest-neighbor fermion formulations use aµ ≤ 0.4, while improved-dispersion fermions may be reliable at chemical potentials more than a factor of 2 larger.
The first step in utilizing this approach is to perform a lattice study at a series of chemical potentials.For neutron star physics one should choose µ d = µ s and µ u somewhat smaller to describe a charge-neutral system including equilibrated leptons.One might study a range of chemical potentials from µ s = 100 MeV to 1 GeV.Since a lattice approach will likely involve determining N at each µ value and integrating it to find the pressure using N q = dP/dµ q , it appears necessary to consider a tightly spaced series of µ values.One also needs a few lattice spacings in order to perform a continuum extrapolation.
The next step would be to use these P PQ (µ) values as constraints, when considering the high-density QCD equation of state.A proposed QCD equation of state is usually expressed as a curve in either the (P, N ), (N, ε), or (P, ε) plane.Standard thermodynamical relations can convert this into a curve in the P, µ plane; if the proposed EOS exceeds P PQ (µ) for any µ where we have data, it is excluded.
Next, the possibility of evaluating P PQ (µ) directly on the lattice at large µ can be used as a check on the performance of the perturbative expansion.Specifically, one can consider the perturbative expansion for P PQ (µ), which as we have argued is the same as the expansion for P (µ) at the currently known order of α 3 s ln(α s ) [46].By comparing this result to the latticedetermined P PQ (µ), one can assess how accurate the perturbative series is as a function of the scale µ.
Finally, one should perform an accurate evaluation of the single diagram which introduces O(α 3 s ) differences between P PQ and P .In any µ-region where the result is actually a small correction, this can be used together with a lattice-determined P PQ (µ) to provide an improved estimate of P (µ).The resulting estimate would also be perturbatively complete at O(α 3 s ), with only α 4 s ln(α s ), higher-order, and nonperturbative corrections remaining.

Figure 2 . 8 e 1 8 e
Figure2.Lattice-to-continuum pressure ratio as a function of the chemical potential in lattice units.The black curve is for nearest-neighbor staggered fermions, the blue curve is for 3'rd-neighbor-improved fermions.At the order of interest, the use of "fattened links" and other modifications to the fermionic action are not relevant.