Lattice QCD Equation of State for Nonvanishing Chemical Potential by Resumming Taylor Expansion

Taylor expansion in powers of baryon chemical potential ($\mu_B$) is an oft-used method in lattice QCD to compute QCD thermodynamics for $\mu_B>0$. Based only upon the few known lowest order Taylor coefficients, it is difficult to discern the range of $\mu_B$ where such an expansion around $\mu_B=0$ can be trusted. We introduce a resummation scheme for the Taylor expansion of the QCD equation of state in $\mu_B$ that is based on the $n$-point correlation functions of the conserved current ($D_n$). The method resums the contributions of the first $N$ correlation function $D_1,\dots,D_N$ to the Taylor expansion of the QCD partition function to all orders in $\mu_B$. We show that the resummed partition function is an approximation to the reweighted partition function at $\mu_B\ne0$. We apply the proposed approach to high-statistics lattice QCD calculations using 2+1 flavors of Highly Improved Staggered Quarks with physical quark masses on $32^3\times8$ lattices and for temperatures $T\approx145$-176 MeV. We demonstrate that, as opposed to the Taylor expansion, the resummed version not only leads to improved convergence but also reflects the zeros of the resummed partition function and severity of the sign problem, leading to its eventual breakdown. We also provide a generalization of our scheme to include resummation of powers of temperature and quark masses in addition to $\mu_B$, and show that the alternative expansion scheme of [S. Bors\'anyi et al., Phys. Rev. Lett. 126, 232001 (2021).] is a special case of this generalized resummation.

Introduction.-Lattice Quantum Chromodynamics (QCD) results for the QCD equation of state (EoS) plays a critical role in the dynamical modeling of heavy-ion collisions [1][2][3][4] and, thereby, in the experimental explorations of the QCD phase diagram in the T -µ B plane. Due to the fermion sign problem it is difficult to carry out lattice QCD computations directly at µ B = 0. Despite some recent progress [5][6][7][8][9][10], direct lattice computations of the QCD EoS µ B = 0 with physical quark masses, fine lattice spacings and large lattice volumes have remained elusive. Instead, the present state-of-the-art lattice QCD EoS at µ B > 0 has been obtained using the Taylor expansion [11,12] and the analytic continuation [13,14] methods. In the Taylor expansion method one expands the pressure in powers of µ B around µ B = 0 and directly computes the Taylor coefficients at µ B = 0. For the analytic continuation, one avoids the fermion sign problem using simulations at purely imaginary values µ B , fits these results with a power series in µ B to determine the Taylor coefficients at µ B = 0 and then provides the EoS at real µ B > 0 based on these Taylor coefficients. On the other hand, it is well-known that the applicability of the Taylor expansion as well as the analytic continuation should be limited by the zeros, nearest to µ B = 0, of the partition function in the entire complex-µ B plane [15][16][17]. In principle, it is possible to gain some knowledge about the locations of the zeros of the partition function by re-expressing the power series in real or imaginary µ B in terms of Padé approximants [12] or in a power series of the fugacity [18][19][20]. Armed, in reality, with only the few lowest order Taylor coefficients, this becomes a very difficult task and, in practice, one just restricts the EoS to {T, µ B } that avoids any pathological nonmonotonicity in the truncated Taylor series [11,14]. Furthermore, these methods provide very little guidance on the severity of the fermion problem, i.e. how rapidly the phase of the partition function fluctuates as µ B is increased. It is possible to determine the zeros of the partition function as well as its average phase by reweighting the fermion determinant to µ B = 0 [21][22][23][24][25]. However, due to the computational cost associated with exact evaluation of the fermion determinant, at present this method is restrained within coarse lattice spacings and small lattice volumes.
In this work, we introduce a method for the calculation of the lattice QCD EoS that genuinely resums the truncated Taylor series to all orders in µ B and whose breakdown encodes the severity of the sign problem and zeros of the resummed partition function.
The resummation method.-The Taylor expansion to O µ N B of the excess pressure, ∆P (T, µ B ) ≡ P (T, µ B ) − P (T, 0), is given by where the Taylor coefficients are defined as Here, the QCD partition function is denoted as Z = e −S det[M ]DU , V is the spatial volume, U is the SU ( and the · denotes average over gauge field ensembles at µ B = 0, i.e. O = Oe −S det[M (T, 0)]DU/Z. The physical interpretation of D n is simple for the continuum theory: D n = dx 1 · · · dx n J 0 (x 1 ) · · · J 0 (x n ) is the integrated n-point correlation function of the 0 th component of the conserved (baryon) current J 0 (x) at a space-time point x. Note that, due to CP symmetry of QCD all D n for n = odd(even) are purely imaginary(real) and only the n = even terms contribute to Eq. 1. In practice, lattice QCD computations of the χ B N involve computations of all D n for n ≤ N as intermediate steps, and χ B N are obtained from combinations of D n and their powers.
Contributions of various combinations of D n to the few lowest order Taylor coefficients are sketched in FIG. 1. If one considers the factorials and the powers of µ B /T associated with each D n in the sum of Eq. 1, it is not difficult to realize that all contributions of each D n to ∆P E can be resummed into exponential forms. For example, contributions of D n 1 from all χ B n in Eq. 1 can be resummed as exp D 1 (µ B /T ) . Similarly, contributions of all D n 2 can be resummed as exp D 2 (µ B /T ) 2 , and so on. Also it is easy to see that the contributions of the mixed terms like D n 1 D m 2 arise from exp D 1 (µ B /T ) × exp D 2 (µ B /T ) 2 . Thus, it is possible to write down a resummed version of Eq. 1, viz.
providing the EoS up to infinite orders in µ B . The ∆P R N can be considered as a µ B -dependent effective action obtained by resumming up to N -point correlation functions of the conserved current. Expansion of ∆P R N in powers of µ B /T yields an infinite series in µ B /T , in addition to the truncated Taylor series: where i, j = 0, . . . , N satisfying 1 · i + · · · + N · j = n. The Taylor expanded (N E N ) and the resummed (N R N ) net baryon-number densities can be straightforwardly obtained as a single µ Bderivative of ∆P E and ∆P R in Eq. 1 and Eq. 4, respectively.
The resummed version in Eq. 4 also highlights the connection between the Taylor expansion and the reweighting method. In the reweighting method Since CP symmetry dictates that the even(odd) D n are purely real(imaginary) and the partition function must be real, a measure of the severity of the sign problem is given by the average phase factor for Z R N (with µ B real), An expansion of cos Θ R N in µ B /T leads to the Taylor expanded measure of the average phase of the partition function [23,26], which we will denote by Θ E N . As the sign problem becomes more severe the average phase cos Θ R N ≈ 0 and resummed results will also show signs of breakdown. Furthermore, although ∆P E N can be evaluated for any complex value of µ B , ∆P R N becomes undefined when Re[Z R N ] ≤ 0 for a given N and statistics, leading to a natural breakdown of the resummed results. The location of the zeros of Z R N in the complex-µ B plane will indicate the µ B region where such resummation can be applicable. Obviously, for any given N the region of applicability of ∆P E N cannot exceed the same for ∆P R N . Lattice QCD computations.-For this work, we used the data for χ B n and D n generated by the HotQCD collaboration for calculations of the QCD EoS [11] and the chiral crossover temperature [28] at µ B > 0 using the Taylor expansion method. The HotQCD ensembles were generated with 2+1-flavors of Highly Improved Staggered Quarks Comparisons between the Taylor expanded and resummed results for different orders for the excess pressure (top) and net baryon-number density (bottom) at T = 157 MeV. Results for real and imaginary µB/T are plotted on the positive and negative x-axis respectively. and the tree-level improved Symanzik gauge action [29][30][31]. Bare quark masses were chosen to reproduce, within a few percent, the physical value of the kaon mass and a pseudo-Goldstone pion mass of 138 MeV in the continuum limit at T = µ B = 0 and the lattice spacing were calibrated against the physical value of the kaon decay constant [32]. We present lattice QCD results from a single lattice size 32 3 × 8 and for 6 temperatures T = 145, 151, 157, 166, 171, 176 MeV. About 475K, 520K, 716K, 522K, 232K and 152K gauge field configurations were used to measure D n at these temperatures respectively. The gauge field configurations were separated by 10 Rational Hybrid Monte Carlo trajectories of unit length. The D n were calculated within the formalism adopted in Refs. [11,28], i.e. using the exponential-µ formalism [33] for n ≤ 4 and the linear-µ formalism [34,35] for n > 4. The expressions for D n in terms of the traces involving the inverse of the staggered fermion matrix and its µ B -derivatives are well-known [26,36]. Each trace was calculated stochastically for each configuration by employing 2000 random Gaussian volume sources for the trace D 1 and 500 random sources for the rest [36].
Results.-To demonstrate the superiority of the resummation method over the Taylor expansion, we chose the temperature where we had the largest statistics, i.e. T = 157 MeV, which is also closest to the QCD crossover temperature [28]. In Fig. 2  isons are shown both for real as well as imaginary values of µ B , corresponding to positive and negative values of (µ B /T ) 2 , respectively. The ∆P R N and N R N show very good convergence between different orders N = 2, 4, 6, 8. The Taylor-expanded results seem to approach their respective resummed results as contributions from higher orders in µ B are included; however the convergence of the Taylor-expanded results is slow due the alternating signs of the higher order χ B n near the QCD crossover [11][12][13]. The resummation method overcomes this problem by including contributions from all orders in µ B and shows markedly improved convergence. In contrast to the Taylor expansion, the resummed results break down for |µ B /T 1.5|. For |µ B /T 1.5|, Re[Z 8 ] ≤ 0 and N R 8 becomes divergingly large. We checked that such a breakdown is not a mere statistical issue by repeating the calculations using only parts of the gauge configurations available at this temperature. Similar breakdown for µ B /T 1.5 was also observed in Refs. [12,37,38] when the EoS was reconstructed from the Padé approximants of the Taylor series in µ B . While Padé-based continuations of the QCD crossover temperature from imaginary values µ B did not encounter such breakdowns [39,40], the same in the case of the EoS seemed to break down due to singularities in the complex-µ B plane [41].
To investigate the origin of this breakdown, we computed the average phase as a function of real µ B , c.f. Eq. 6. The results are shown in FIG. 3 (top). Also, cos Θ R N ≈ 0 for µ B /T 1.5, which shows that the sign problem is uncontrollably severe where the EoS calculations broke down. The resummation method thus faithfully captures the severity of the sign problem, as opposed to the Taylor expansion. The phase factor cannot be calculated exactly within the Taylor series approach. Its Taylor series expansion too converges very slowly, as the bands plotted in FIG. 3 (top) show. Further, we searched for the zeros of resummed partition function, c.f. Eq. 5, in the complex-µ B plane. We solved for Z R N = 0 using the Newton-Raphson algorithm with initial guesses chosen from a uniform distribution over a grid 0 ≤ {Re(µ B /T ), Im(µ B /T )} ≤ 2.5. The results are shown in FIG. 3 (bottom). The zeros of Z R 6 and Z R 8 are more or less consistent with each other and appears only for |µ B /T | 1.5. The exact nature of the singularity responsible for breakdown of the resummation method is certainly of great interest, i.e. whether it is associated with the Yang-Lee edge singularity of the QCD chiral transition [15,17] or the QCD critical point and approaches the real axis [12,21,22,37,38] etc. This will need detailed quantitative studies involving careful finite-volume scaling analyses using more sophisticated techniques [18,20,42] and is beyond the scope of the present work. But our results demonstrate that the breakdown of the resummation method reflects the associated singularities of the partition function, at least qualitatively.
Finally and N R 6 with the corresponding ∆P E 6 and N E 6 . As in the case of T = 157 MeV, ∆P R and N R show improved convergence over ∆P E and N E at all temperatures. Again, in contrast to the Taylor expansion the resummation method shows signs of breakdown for µ B 200 − 250 MeV, depending on the temperature. As before, we checked that in all cases, these breakdowns reflect the severity of the sign problem and the singularities of the partition function in the complex-µ B plane.
Generalization to multi-parameter and joint expansion in T, µ B .-Akin to multi-parameter reweighting [21][22][23][24] in bare gauge coupling, ∆β = β − β 0 , and quark mass, ∆m = m − m 0 , our resummation scheme also can be extended to obtain Z R N (T, µ B ) starting from a different temperature T 0 (β 0 ) and bare quark mass m 0 , where the expectation value is taken over gauge fields associated with {β 0 , m 0 , 0}. Here, S G is the pure gauge action and Note,Ḡ i0 =D i (Eq. 3),Ḡ 0j are the chiral condensate and higher order chiral susceptibilities, and gen-eralḠ ij are µ B -derivatives of various chiral observables [26,36,43,44]. This generalization can possibly mitigate the overlap problem that one might encounter while resumming only in µ B . Further, a systematic expansion of the logarithm of Eq. 7 in powers of ∆β, ∆m and µ B yields the expansion of the pressure difference, P (T, µ B )/T 4 − P (T 0 , µ B )/T 4 0 , in powers of ∆T = T − T 0 and µ B ; particular choice of T 0 (µ B ) defined by a line of constant physics in the T -µ B -plane reproduces the expansion scheme used in Ref. [45] by resumming up to N -point baryon-current correlations to all orders in µ B and ∆T [46]. Thus, our method also generalizes the alternative expansion scheme of Ref. [45].
Conclusions.-We have introduced a new method to compute lattice QCD EoS by resumming contributions of up to N -point baryon-current correlations to all orders in µ B . When expanded in powers of µ B this resummed partition function exactly reproduces the Taylor expansion up to O µ N B , plus an infinite series in µ B capturing all possible contributions involving only the n ≤ N -point baryon-current correlations. This resummation method also amounts to an approximate reweighting method, thereby bridging two traditional lattice QCD techniques for µ B = 0. With illustrative high-statistics lattice QCD computations we have demonstrated that the resummation method show improved convergence over the Taylor expansion method. The method also faithfully captures the severity of the sign problem as well as reflects the singularities in the complex-µ B plane that are responsible for its eventual breakdown. Thus the resummation method not only provides a more convergent lattice QCD EoS but also a more reliable one by enabling us to judge its validity with increasing µ B . Although the resummation method is more general and powerful than the Taylor expansion, computationally it is somewhat simpler. The resummation method relies on the computations of D n which come as an intermediate step in the computations of the Taylor coefficients. Comparison with the resummed results and the direct lattice QCD simulations for purely imaginary µ B will help us to decide up to what values Im(µ B ) an analytic continuation using only the power series of µ B is justified and whether Padé-type analytic continuations [39][40][41] are necessary to avoid singularities in the complex-µ B plane. We have also introduced a generalized multi-parameter version of the resummation, Eq. 7, and shown that the method of Ref. [45] is a special case of this-Taylor expansion of Eq. 7 in T and µ B along a specific line in the T -µ B -plane.
Acknowledgments.-We are indebted to members of the HotQCD collaboration for letting us reuse the data they had generated for the Taylor expansion calculations as well as for several valuable discussions. Our starting point for the multi-parameter expansion is Eq. 7. There, the expectation value is taken over gauge fields associated with {β 0 , m 0 , µ B = 0}, where β is the QCD gauge coupling, S G is the pure gauge action and theḠ ij are given by Eq. 8.
For brevity, in this section, we will use notationŝ µ B ≡ µ B /T ,m = m/T and ∆m ≡m −m 0 , and provide explicit demonstration of the generalized resummation by keeping only leading order terms, i.e. O(∆β), O(∆m) and O μ 2 B , in all expansions. Extensions to higher or-ders are straightforward. We consider the case where the temporal extent (N τ ) of the lattice is kept fixed. In this case, T is changed by varying the bare gauge coupling β, and the bare quark mass m must be tuned with β to keep vacuum hadron masses constant. Thus, T (β, m) and T 0 (β 0 , m 0 ). Applying chain rule for derivatives as well as expanding T (β, m) around (β 0 , m 0 ) one gets to the leading order The goal is to obtain Z R N (T, µ B ) by expanding around Z(T 0 , 0), while resumming contributions of up to N -point baryon-current correlations to all orders in µ B and ∆T . Following multi-parameter reweighting technique where the expectation value is with respect to a gauge field ensemble generated for {β 0 , m 0 , µ B = 0}. Simultaneously expanding in µ B and ∆m the determinant ratio in Eq. 12 can be written as whereḠ ij are defined through Eq. 8. By plugging Eq. 13 back into Eq. 12 and truncating the sum at i + j = N , we obtain Eq. 7. Next, by Taylor expanding Eq. 7 in powers of ∆β, ∆m andμ B we get The pressure difference is given by where N s is the spatial extent of the lattice. Using Eq. 14, expanding the logarithm in powers of ∆β, ∆m,μ B and keeping only the real part one obtains N 3 s N 3 τ ∆ P T 4 = Ḡ 01 ∆m − S G ∆β + D 2 +D 2 1 /2 μ 2 B − S G (D 2 +D 2 1 /2) − S G (D 2 +D 2 1 /2) μ 2 B ∆β + Ḡ 01 (D 2 +D 1 /2) − Ḡ 01 (D 2 +D 2 1 /2) + Ḡ 21 +Ḡ 11D1 μ 2 B ∆m + . . . .
Noting that and using Eq. 9 of the main paper, it is easy to identify that Eq. 15 is nothing but a joint Taylor expansion of P (T, µ B ) in T and µ B around (T 0 , 0), Thus, the generalized version given by Eq. 7 genuinely resums contributions of up to N -point baryon current in the Taylor expansion of EoS to all orders in T, µ B . Following Ref. [45], the generalized resummation of Eq. 7 can be made even more powerful by choosing the expansion point T 0 along some physically motivated line in the T -µ B -plane, i.e. by choosing some physically motivated β 0 (µ B ) and m 0 (µ B ). The one-to-one correspondence between the Taylor expansion of Eq. 7 and alternative expansion scheme presented in Ref. [45] can be readily observed. By including the O(μ 4 B ) in Eq. 17 and taking aμ B -derivative we get If one chooses in Eq. 18 then one arrives at the starting point of Ref. [45], namely χ B 1 (T 0 ,μ B ) =μ B χ B 2 (T 0 , 0). Hence, the method used in Ref. [45] is a special Taylor-expanded case of the generalized resummation Eq. 7.