New Way to Resum the Lattice QCD Taylor Series Equation of State at Finite Chemical Potential

Taylor expansion of the thermodynamic potential in powers of the (baryo)chemical potential $\mu_B$ is a well-known method to bypass the Sign Problem of Lattice QCD. Due to the difficulty in calculating the higher order Taylor coefficients, various alternative expansion schemes as well as resummation techniques have been suggested to extend the Taylor series to larger values of $\mu_B$. Recently, a way to resum the contribution of the first $N$ charge density correlation functions $D_1,\dots,D_N$ to the Taylor series to all orders in $\mu_B$ was proposed in Phys. Rev. Lett. 128, 2, 022001 (2022). The resummation takes the form of an exponential factor. Since the correlation functions are calculated stochastically, the exponential factor contains a bias which can be significant for large $N$ and $\mu_B$. In this paper, we present a new method to calculate the QCD equation of state based on the well-known cumulant expansion from statistics. By truncating the expansion at a maximum order $M$, we end up with only finite products of the correlation functions which can be evaluated in an unbiased manner. Although our formalism is also applicable for $\mu_B\ne0$, here we present it for the simpler case of a finite isospin chemical potential $\mu_I$ for which there is no Sign Problem. We present and compare results for the pressure and the isospin density obtained using Taylor expansion, exponential resummation and cumulant expansion, and provide evidence that the absence of bias in the latter actually improves the convergence.

Taylor expansion of the thermodynamic potential in powers of the (baryo)chemical potential µB is a well-known method to bypass the Sign Problem of Lattice QCD. Due to the difficulty in calculating the higher order Taylor coefficients, various alternative expansion schemes as well as resummation techniques have been suggested to extend the Taylor series to larger values of µB. Recently, a way to resum the contribution of the first N charge density correlation functions D1, . . . , DN to the Taylor series to all orders in µB was proposed in Phys. Rev. Lett. 128, 2, 022001 (2022). The resummation takes the form of an exponential factor. Since the correlation functions are calculated stochastically, the exponential factor contains a bias which can be significant for large N and µB. In this paper, we present a new method to calculate the QCD equation of state based on the well-known cumulant expansion from statistics. By truncating the expansion at a maximum order M , we end up with only finite products of the correlation functions which can be evaluated in an unbiased manner. Although our formalism is also applicable for µB = 0, here we present it for the simpler case of a finite isospin chemical potential µI for which there is no Sign Problem. We present and compare results for the pressure and the isospin density obtained using Taylor expansion, exponential resummation and cumulant expansion, and provide evidence that the absence of bias in the latter actually improves the convergence.

I. INTRODUCTION
The Equation of State (EoS) of strongly-interacting matter is an important input in the hydrodynamical modeling of heavy-ion collisions [1][2][3][4]. Unfortunately Lattice QCD, which is the preferred method of calculating observables in the non-perturbative regime of QCD, breaks down when the baryon chemical potential µ B is non-zero. This is the well-known Sign Problem of Lattice QCD [5]; despite recent progress [6][7][8][9][10], the current state-of-the-art results for the QCD EoS have been obtained by using either either analytical continuation from imaginary to real µ B [11,12], or by expanding the EoS in a Taylor series in the chemical potential µ B and calculating the first N coefficients [13,14]. In the latter case, a knowledge of the first several coefficients is necessary, not only to obtain the EoS for a fairly wide range of chemical potentials but also to determine the radius of convergence of the Taylor series beyond which the Taylor expansion must break down [15][16][17][18]. Unfortunately, the calculation of the higher order Taylor coefficients is computationally very challenging and it is natural to ask whether something can be learned about them from a knowledge of the first few Taylor coefficients. It turns out that this is indeed possible because the first N derivatives D 1 , . . . , D N of ln det M (µ B ), where M (µ B ) is the fermion matrix, also contribute to the higher order Taylor coefficients through products such as D 2 N , D N D 1 , etc. In fact, the contribution of the nth derivative D n to all higher orders can be shown to take the form of an exponential exp(D n µ n B /n!) [19]. Thus, if D n is known exactly, then its contribution to the Taylor series can be resummed to all orders through exponentiation. Exponential resummation, as we will refer to it from here on, can be shown to have several advantages compared to the original Taylor series: First, the resummed EoS converges faster than the Taylor series. Moreover, since the odd derivatives D 1 , D 3 , . . . are purely imaginary, the resummed expression directly gives us a phase factor whose expectation value approaches zero as µ B is increased, leading to a breakdown of the calculation. This breakdown is physical and related to the presence of poles or branch cut singularities of the QCD partition function in the complex µ B plane. The resummed expression for the partition function also makes it possible to calculate these singularities directly. Some of these advantages have been recently demonstrated through analytical calculations in a lowenergy model of QCD [20].
Despite its advantages, a technical drawback of exponential resummation is that the derivatives D 1 , . . . , D N are not known exactly in an actual Lattice calculation. As is easily seen from the identity ln det M = tr ln M , the D n can be expressed in terms of traces of various operators, all of which involve the inverse of the fermion matrix M . Since M is typically of size 10 8 or greater, calculating its exact inverse is prohibitively expensive. Instead the various traces, and hence the derivatives D n , are estimated stochastically using O(10 2 -10 3 ) random volume sources per gauge configuration. Now, the products of such stochastically estimated quantities e.g. D 2 N , need to be evaluated in an unbiased manner i.e. estimates coming from the same random vector must not be multiplied together. If D By contrast, the naïve Biased Estimate (BE) is given by (2) Eqs. (1) and (2) can both be readily generalized to any finite power or to the product of a finite number of traces.
In the Appendix, we present formulas for evaluating the unbiased estimate of such finite products in an efficient manner. However we do not know of any corresponding formula to calculate the unbiased estimate of an infinite series such as an exponential.
In this paper, we will present a new way of calculating the QCD EoS based on the well-known cumulant expansion from statistics. The cumulant expansion method is intermediate between a strict Taylor series expansion and exponential resummation in the sense that the contribution of Although the cumulant expansion method also works for µ B = 0, in this paper we will present the formalism for the simpler case of finite isospin chemical potential µ I instead. For µ I = 0, the fermion determinant is real and one has no Sign Problem. Thus one only works with real quantities which in turn simplifies the presentation. Moreover, the absence of the Sign Problem allows us to calculate observables for much larger values of µ I than would be possible for the µ B case, and it is precisely for these large values that bias can become significant. Lastly, the QCD phase diagram in the T -µ I plane is known from several studies to be interesting in its own right [22][23][24][25][26], and we hope to be able to apply this formalism to its study in the future.

II. FORMALISM
We consider Lattice QCD with 2+1 flavors of rooted staggered quarks. The partition function at non-zero isospin chemical potential µ I is given by where det M (T, µ I ) is shorthand for with m u = m d , µ u = −µ d = µ I and µ s = 0. The excess pressure ∆P (T, µ I ) ≡ P (T, µ I ) − P (T, 0) is given by where V is the spatial volume and T is the temperature. By employing the same arguments as in Ref. [19], we can write where the expectation value · is taken over a gauge field ensemble generated at µ u = µ d = µ s = 0, and The presence of only even powers is because the odd µ I derivatives vanish identically. Since even derivatives of the quark determinant are purely real, we see that det M (µ I ) is purely real and hence there is no Sign Problem. Note that this is true even when µ I is purely imaginary.
The D I n can be expressed as traces of various operators [27,28]. In Lattice calculations, the first N derivatives D I 1 , . . . , D I N are calculated stochastically using N rv random vectors, where N rv is typically of order 10 2 -10 3 . Then ∆P (T, µ I )/T 4 is approximately equal to (8) Here N σ and N τ are the number of lattice sites in the space and time directions respectively, whileD I 2n is the average of the N rv stochastic estimates of D I 2n . Eq. (8) is the N th order exponential resummation formula for ∆P (T, µ I )/T 4 . In the limit N rv → ∞, it accurately resums the contribution of the first N derivatives D I 1 , . . . , D I N to all orders in µ I [19]. For N rv < ∞ however, the formula contains bias. This is easily seen if one writes the exponential as an infinite series. The series expansion leads to terms such as (D I 2m ) p (D I 2n ) q · · · , and we have already seen that such products are biased due to multiplication of estimates coming from the same random vector.
The well-known cumulant expansion formula from statistics states that The coefficients C k (X) are known as the cumulants of X [29][30][31]. The first four cumulants are given by In our case t = 1, which we assume lies within the radius of convergence of the cumulant expansion, and Truncating Eq. (9) at k = M ≥ N/2 [32] gives us yet another way to estimate ∆P/T 4 , namely Eq. (12) may be compared to the familiar Taylor series expansion of ∆P/T 4 , which in our case is given by The restriction M ≥ N/2 in Eq. (12)  for evaluating these products efficiently in an unbiased manner. Thus, the cumulant expansion is free of the bias that can affect exponential resummation.
Finally, we will also present results for the net isospin density N (T, µ I ) which is given by The Taylor series expression N E N (T, µ I ) for the same is straightforward. The resummed and cumulant expansion expressions N R N (T, µ I ) and N C N,M (T, µ I ) can be obtained by differentiating Eqs. (8) and (12) respectively. We do not write down the explicit expressions here. Note however that the resummed formula, unlike the cumulant expansion expression, involves a ratio of expectation values.

III. RESULTS
To verify our formalism, we made use of the data generated by the HotQCD collaboration for their calculations of the finite-density EoS, finite-density chiral crossover temperature and conserved charge cumulants using Taylor series expansions [13,17,33]. The data consists of 2+1-flavor gauge configurations with N τ = 8, 12 or 16 and N σ = 4N τ in the temperature range 125 MeV T 178 MeV. The configurations were generated using a Symanzik-improved gauge action and the Highly Improved Staggered Quark (HISQ) action [34][35][36] for fermions. The lattice spacing was determined using both the Sommer parameter r 1 as well as the decay constant f K . The temperature values quoted in this paper were obtained using the f K scale. For each lattice spacing, the light and strange quark bare masses were tuned so that the pseudo-Goldstone meson masses reproduced the physical pion and kaon masses. A description of the gauge ensembles and scale setting can be found in Ref. [14].
The results presented here were obtained with around 20,000 N τ = 8 configurations for T = 135 MeV. On each gauge configuration, the first eight derivatives D f 1 , . . . , D f 8 for each quark flavor were estimated stochastically using around 2000 Gaussian random volume sources for D f 1 and around 500 sources for the rest.
We used the exponential-µ formalism [37] to calculate the first four derivatives, while the linear-µ formalism [38,39] was used in calculating all higher derivatives. We calculated the excess pressure ∆P (T, µ I ) and net isospin density N (T, µ I ) using (i) Taylor series expansion (Eq. (13)) for N = 2, 4 and 6; (ii) exponential resummation (Eq. (8)) for N = 2 and 4; and (iii) a fourth order cumulant expansion (Eq. (12) with M = 4) for N = 2 and 4. We calculated these observables for both real as well as imaginary µ I , in the range 0 |µ I /T | 2.
We present our results for ∆P (T, µ I ) and N (T, µ I ) in Figs. 1 and 2  Focussing first on the upper plots, we see that the second and fourth order Taylor expansion results (green and blue bands) start to differ significantly around |µ I /T | = 1. For real µ I , this difference is seemingly captured by the resummed results (red band) which agree with the fourth order Taylor results for both observables over nearly the entire range of µ I /T . For imaginary µ I however, the resummed results lie below the fourth order Taylor results. By contrast, the cumulant expansion results (blue inverted triangles) are in good agreement with the fourth order Taylor results for imaginary µ I , while they are only slightly less than the fourth order Taylor results for real µ I .
One explanation for the difference between the resummed and cumulant results is the higher order contributions that are present in the former but not in the latter. Another possibility is the bias that is present in the resummed but not in the cumulant results. To distinguish between the two possibilities, we recalculated the cumulant results using the biased formulas for the trace products rather than the unbiased ones (Eq. (A.2) rather than Eq. (A.3)). We found that the biased results (upright red triangles) agree well with the resummed results, thus suggesting that bias, rather than the contribution from higher orders, is responsible for the difference between the resummed and the cumulant expansion results.
To further confirm that this is the case, we also recalculated the resummed results using only 250 random vectors instead of 500. Since bias decreases as N rv is increased, conversely we should expect the bias to increase when we use fewer random vectors. From Figs. 1 and 2, we see that the N rv = 250 results (yellow band) lie further from the Taylor and unbiased cumulant results than the N rv = 500 results for |µ I /T | 1. Thus we see that the resummed results are indeed affected by bias for large values of the chemical potential.
The presence of bias must especially be accounted for when comparing higher order results (lower plots in Figs. 1 and 2). We see that the sixth order Taylor correction (blue band) to the fourth order results (green band) is small over the entire range of µ I considered here. By contrast the resummed results, although not containing the contributions of the operator D I 6 , nonetheless suggest that the higher-order contributions of D I 2 and D I 4 are large for both real as well as imaginary µ I . However, both the biased cumulant expansion results as well as the N rv = 250 resummed results once again suggest that at least some of this difference is due to bias. On the other hand, the unbiased cumulant results agree well with the Taylor series results for both real and imaginary µ I . We note that while the cumulant expansion too does not include the contribution of the operator D I 6 , it does contain higher order corrections that appear at O(µ 6 I , . . . , µ 16 I ). In fact, the (4, 4) cumulant expansion is exactly equal to a fourth order Taylor expansion, plus the contributions of the operators D I 2 and D I 4 to the Taylor coefficients χ I 6 , . . . , χ I 16 . The agreement between the unbiased cumulant and Taylor series results thus suggests that the contribution of D I 2 and D I 4 at higher orders is in fact small. All this goes to show that bias needs to be properly accounted for before one can identify the genuine higher-order corrections.

IV. CONCLUSIONS
In this paper, we presented a new way of resumming the QCD Taylor series EoS based on the well-known cumulant expansion from statistics. Our approach is a finite order truncation of the all orders resummation of the first N charge density correlation functions that was presented in Ref. [19]. The resummation presented there is susceptible to bias when the correlation functions are calculated stochastically. By contrast, the cumulant expansion contains only finite products of traces that can be evaluated efficiently in an unbiased manner. Moreover, while not an all orders resummation, the M th order cumulant expansion still captures the contributions of the lower order derivatives D 1 , . . . , D N to the higher order Taylor coefficients up to a maximum order χ M N .
Although our formalism is also applicable to µ B = 0, in this paper we presented results for finite isospin chemical potential µ I instead. Our reason for this was that the absence of a Sign Problem in the latter case meant that all quantities were real and this simplified the presentation. The µ I = 0 case is also of interest in its own right. We presented results for the excess pressure and net isospin density using Taylor series, resummation and cumulant expansion. We found evidence for bias in the resummed results at large (µ I /T ) 2 . We showed this by (a) calculating the cumulant expansion using biased rather than unbiased products, and (b) recalculating the resummed results using fewer random vectors. The cumulant expansion is a truncation of the resummed formula and when the terms of the expansion were calculated using biased products of traces, they were in good agreement with the resummed result, while they were closer to the Taylor series results when they were calculated using unbiased products.
There have been several proposals recently to extend the QCD EoS to larger values of the chemical potential [16,17,[40][41][42]. Exponential resummation is one such approach, which can be connected to reweightingbased approaches. By generalizing it to resummation in T and µ B , it can be also connected with the alternative expansion scheme proposed in Ref. [40]. The cumulant expansion approach that we have outlined here provides yet another way to extend the QCD EoS. The precise relation between this approach and the other proposals remains to be studied.
All data presented in the figures of this paper can be found in Ref. [45].

ACKNOWLEDGMENTS
We thank all the members of the HotQCD collaboration for their inputs and for valuable discussions, as well as for allowing us to use their data from the Taylor expansion calculations. CS was supported by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agree- In practice, the limit N rv → ∞ is usually dropped and the trace is estimated by the arithmetic mean over a finite number of estimates O i , which bears potential danger for the estimation of powers of the trace (tr O) n , with n ≥ 1.
The naïve estimate is a biased estimate since it receives contributions from products of the same estimates O k i , k ≤ n, which lead to a finite bias in the estimator (A.2). The bias can be quite significant in the regime of large n, even when we take N rv ∼ O(1000). An unbiased estimator for (tr O) n can be obtained by rejecting all diagonal contributions and considering only mutually distinct estimates in each product. We define In this case the sum is taken over all partitions of n, which makes the connection to the complete exponential Bell polynomials B n and is made explicit with the second equality. Below we give the final formulae in the form of Eq. (A.5), which are required up to n = 4, , , . (A.6) Although there are several (power) sums to be evaluated, these are not nested sums as in Eq. (A.3). Therefore, the unbiased product can be calculated in only O(N rv ) time.

Unbiased Estimators for Combinations of Traces of Multiple Operators
Quite often we encounter a situation where we need to estimate expressions involving combinations of traces of several operators O (1) , O (2) , . . . , O (m) , which are estimated on the same set of random vectors and are thus correlated. The evaluation of different operators on the same set of random vectors might seem to be avoidable at the first glance but could -as in the case of the derivative operators of the pressure discussed above -gain computational advantages, e.g., due to a recursive definition of the operators [27,44]. A biased estimator of a general expression of this kind is given by In order to construct an unbiased estimator for this more general case, we extend the framework of elementary symmetric polynomials and power sums presented above as follows: We introduce the multi-index notation α, β, γ ∈ N m with non-negative integer coefficients and define the metric |α| ≡ i α i . We further define the two types of symmetric polynomials as i β 1 +1 · · · O where the addition of multi-indices α + β is understood by components, the summation is over the product of all partitions of the components of γ and |α| α denote the multinomial coefficients. We can now define an unbiased version of (A.7) based on the symmetric polynomial e γ as UE tr O (1) γ1 · · · tr O (m) γm = |γ|! · e γ N rv · · · (N rv − |γ| + 1) .
(A.10) For the practical calculation of e γ we use the recursive definition (A.9). To reduce the computational efforts we manually cache previously computed values of p α and e β . We stress again that the main computational effort goes into the evaluation of the power sums. Each power sum has complexity O(N rv ). The number of distinct power sums given by i (γ i + 1) − 1. Even though the power sum number increases also drastically with the order |γ|, the values we encounter in this calculation are still of the order O(100) and thus an order of magnitude lower than the costs of a single power sum evaluation. For convenience we give in Tab. I some examples for cases with |γ| = 8, which we identified as the cases with the maximum number of power sums, for a given number of operators m. Note that for three-quark flavors (u, d, s) and |γ| ≤ 8, the maximum number of distinct operators we encounter in and unbiased estimator is m = 5.  Table I. The maximum number of distinct power sums that need to be evaluated for an unbiased estimator of order |γ| = 8, as a function of the number of distinct operators m (dimension of γ).