Cluster Expansion Model for QCD Baryon Number Fluctuations: No Phase Transition at $\mu_B / T<\pi$

A Cluster Expansion Model (CEM), representing a relativistic extension of Mayer's cluster expansion, is constructed to study baryon number fluctuations in QCD. The temperature dependent first two coefficients, corresponding to the partial pressures in the baryon number $B = 1$ and $B = 2$ sectors, are the only model input, which we fix by recent lattice data at imaginary baryochemical potential. All other coefficients are constructed in terms of the first two and required to match the Stefan-Boltzmann limit at $T \to \infty$. The CEM allows calculations of the baryon number susceptibilities $\chi_k^B$ to arbitrary order. We obtain excellent agreement with available lattice data for the baryon fluctuation measures $\chi_2^B$, $\chi_4^B$, $\chi_6^B$ and predict higher order susceptibilities, that are not yet available from Lattice QCD. The calculated susceptibilities are then used to extract the radius of convergence of the Taylor expansion of the pressure. The commonly used ratio test fails due to the non-trivial asymptotic behavior of the Taylor coefficients. At the same time, a more elaborate estimator provides finite convergence radii at all temperatures and in agreement with the singularities of Pad\'e approximants. The associated singularities lie in the complex $\mu_B/T$-plane and appear smoothly connected to the Roberge-Weiss transition at high temperatures and imaginary chemical potential. No evidence for a phase transition at $\mu_B / T \lesssim \pi$ is found.

The grand canonical thermodynamic potential of QCD is an even function of the real baryon chemical potential µ B because of the CP-symmetry. Taking the quantization of the baryon charge into account, the QCD pressure has a fugacity expansion This relation can formally be viewed as a relativistic extension of Mayer's cluster expansion in fugacities [1]. A similar formalism is also employed in the canonical approach [2][3][4][5], where the fugacity expansion is applied to the partition function Z itself (see Refs. [6,7] for recent developments). The net baryon density reads where b k (T ) ≡ k p k (T ). Analytic continuation to an imaginary chemical potential yields a purely imaginary ρ B /T 3 , with b k (T ) becoming its Fourier expansion coefficients. The analytic continuation to/from imaginary µ is one of the methods used to study QCD at finite net baryon density [8][9][10][11][12][13][14][15][16]. The leading four Fourier coefficients b k were studied in Refs. [17,18] using lattice simulations with physical quark masses at imaginary µ B , as well as phenomenological models. At low temperatures, below the crossover transition, a behavior similar to an uncorrelated gas of hadrons, i.e. |b k (T )/b 1 (T )| 1 for k ≥ 2, is seen in lattice QCD (LQCD) [17,18]. Lattice simulations at T > 135 MeV yield negative b 2 (T ) values, which could be interpreted in terms of repulsive baryonic interactions. The first four LQCD Fourier coefficients, up to T 185 MeV, can indeed be described rather well by a hadron resonance gas (HRG) model with excluded-volume (EV) interactions between the (anti)baryons, with one canonical "eigenvolume" parameter b 1 fm 3 for all (anti)baryon species [18].
QCD in the high-temperature limit resembles an ideal gas of massless quarks and gluons. In this Stefan-Boltzmann (SB) limit, the quarks carry a fractional baryon charge of 1/3, and the coefficients b k read [18] b The fugacity series (2)  treats the non-interacting hadron limit at low temperatures and the non-interacting quark limit at high temperatures in the same framework. In this work a model is constructed which allows to calculate the coefficients b k at all intermediate temperatures between these two limiting cases. This Cluster Expansion Model (CEM) is based on the following assumptions: • The first coefficient b 1 (T ) -the QCD partial pressure in the |B| = 1 sector -is taken as input. It is interpreted as a temperature dependent density of "free" excitations with B = ±1.
• The second coefficient, b 2 (T ), is also taken as input. In the spirit of a cluster expansion it parametrizes the baryon-baryon interactions. In the HRG-EV model b 2 is rewritten as where b(T ) is a temperature dependent "coupling" parameter.
• Mayer's cluster expansion assumes two-baryon interactions only, expected to be a good approximation at sufficiently low density or high temperature, i.e. moderate µ B /T . The higher-order coefficients b k (T ) are then expressed in terms of the first two, motivated by a HRG-EV-type system with two-particle hard core interactions [18]: the α k are temperature independent parameters.
• The model is constrained by the SB limit (3) of massless quarks and gluons at high tempera- , this condition fixes the coefficients α k : Eqs.  In what we term CEM-LQCD, b 1 (T ) and b 2 (T ) are fixed by recent (2+1)-flavor, N τ = 12 lattice QCD simulations at imaginary µ B of the Wuppertal-Budapest collaboration [18]. In an alternative CEM-HRG, b 1 (T ) and b 2 (T ) are taken from the HRG-EV model with a constant b(T ) = 1 fm 3 value [18][19][20].
Note that, for a calculation of the pressure using the CEM, also the partial pressure p 0 (T ) in the |B| = 0 sector is required as input. Here we only study baryon number fluctuations for which this is not needed.
Temperature dependences of the first four coefficients b k (T ), as calculated in lattice QCD simulations [18], the CEM-LQCD model, and the CEM-HRG model, are shown in Fig. 1 by the circles, the stars, and the dashed lines, respectively. The CEM-LQCD parametrization reproduces the lattice data for Calculations within CEM-LQCD and CEM-HRG (b = 1 fm 3 ) models are depicted by the red stars and dashed black lines, respectively. Lattice QCD results of the Wuppertal-Budapest [22,23] and HotQCD [24,25] collaborations are shown, respectively, by the blue and green bands where available.
in Eq. (2) within CEM is given by the ratio test: This condition is violated at sufficiently high µ B /T values. At µ B = 0, the condition is satisfied at all temperatures considered for the CEM-LQCD. In the CEM-HRG at µ B = 0, it is satisfied only up to T 195 MeV, which correlates with the breakdown of CEM-HRG shown in Fig. 1.
In the following, the net baryon number susceptibilities χ B 2n (T ) are calculated at zero baryochemical potential. They are defined as follows: The χ B 2n are proportional to the coefficients of the Taylor expansion of the QCD pressure with respect to µ B /T . The odd-order susceptibilities, The temperature dependence of the net baryon number susceptibilities, up to χ B 12 , is shown in Fig. 2 by red stars (CEM-LQCD) and dashed black lines (CEM-HRG), respectively. The error bars shown for the CEM-LQCD calculations were obtained by standard statistical propagation of the uncertainties in the input data, i.e. by calculating the derivatives of the observables with respect to b 1 (T ) and b 2 (T ). We cross checked the validity by varying b 1 (T ) and b 2 (T ) within their uncertainties and observing consistent variations in the observables.
In our calculations, the sum (8) is truncated at a sufficiently high value k = k max , such that the terms with k > k max have negligible contributions to χ B 2n at a given temperature T . Fig. 2(a) shows that the first two terms (k max = 2, i.e. the lattice input), reproduce the full result for χ B 2 only at moderate temperatures, T 160 MeV. Hence, a higher number of terms is required to correctly calculate the susceptibilities at higher temperatures: For k max = 20, the full result for χ B 2 is reproduced up to T 200 MeV, while k max 80 is required to calculate χ B 2 at T = 230 MeV. In order to cope with the large round-off errors, which arise in the numerical calculations, the arbitrary precision arithmetic provided by the Mathematica package is employed.
Lattice QCD results for χ B 2 , χ B 4 /χ B 2 , and χ B 6 /χ B 2 of the Wuppertal-Budapest [22,23] and HotQCD [24,25] collaborations are also shown in Fig. 2, where available. Comparing the CEM results with the lattice data elucidates the excellent predictive power of the present CEM-LQCD approach. A precise calculation of the χ B 2n values requires summation over tens and,  in some cases, hundreds of b k (T ) coefficients. All of them, except the first two, are predicted by CEM. The CEM-LQCD predictions for χ B 8 , χ B 10 , and χ B 12 are shown in Fig. 2(d)-(f). The comparison with future lattice data will be able to confirm (or refute) the validity of the CEM approach presented here.
The CEM-HRG model results, as shown by the dashed lines in Fig. 2, agree very well with CEM-LQCD calculations, up to T 180 MeV, for all considered observables. Hence, the drastic temperature dependence of the baryon number fluctuations in this temperature range, as well as the particularly strong deviations from the ideal HRG baseline -the Skellam distribution -are convincingly interpreted in terms of repulsive baryonic interactions (see also [18,19,21]).
The ability of the CEM-formalism to calculate baryon number susceptibilities to very high order provides a unique opportunity to analyze the radius of convergence of the Taylor expansion of the QCD pressure, The radius of convergence, r µ/T , of this series at a given temperature corresponds to the distance to the nearest singularity in the complex µ B /T plane and this has been used in various attempts to constrain the location of the critical point of QCD by numerical evaluation of a few leading coefficients in lattice QCD [26][27][28] or in effective models [29][30][31]. Deriva-tives of the pressure series expansion may be used equally well. In the present work, estimates based on the Taylor series of p/T 4 , χ B 2 , and χ B 4 are analyzed. First the ratio estimator, r n = |c n /c n+1 | 1/2 , is used. The square root in this estimator [as well as the extra square root in Eq. (10)] appears due to the fact that the Taylor expansion (9) is actually in (µ B /T ) 2 rather than just in µ B /T . Here c n = χ 2n /(2n)! for the p/T 4 expansion, c n = χ 2n /(2n − 2)! for the χ B 2 expansion, and c n = χ 2n /(2n − 4)! for the χ B 4 expansion. The n → ∞ limit of r n , if it exists, is the same for all three expansions and corresponds to the true radius of convergence. This limit can be determined with the Domb-Sykes presentation [32], by plotting 1/r 2 n−1 versus 1/n for a finite number of terms, and then extrapolating the result linearly to 1/n = 0. To illustrate the behavior of the r µ/T estimators we show T = 160 MeV as an example, the behavior at all other temperatures investigated is similar. The Domb-Sykes plot for the Taylor series of p/T 4 , as obtained within the CEM-LQCD model at T = 160 MeV by using the first 200 terms of the Taylor expansion, is depicted in Fig. 3a by the open symbols. (The plots for χ B 2 and χ B 4 are similar and not shown). Note how the different orders jump between several branches of 1/r 2 n−1 as 1/n approaches zero, with no unique limiting value in sight. This behavior is caused by the irregular asymptotic structure of the Taylor coefficients. Convergence of a Domb-Sykes plot for the ratio test requires the coefficients to asymptotically be of the same sign or to alternate in sign. Neither of the two scenarios is realized in the CEM-LQCD: even at very high order, at least two positive-and at least two negative coefficients appear regularly in a row. Therefore, the ratio estimator does not give a correct estimate of r µ/T since the limit lim n→∞ r n simply does not exist. (Note that the ratio estimator is commonly used in the lattice QCD studies of the Taylor expansion [24,33,34]).
More elaborate estimators do exist which deal with the irregular asymptotic structure of the Taylor coefficients. Consider the Mercer-Roberts estimator [35], The corresponding 1/r 2 n−1 vs 1/n plot is shown by the full symbols in Fig. 3b. For all three Taylor expansions, the Mercer-Roberts estimators appear to converge to the same point as 1/n → 0. Linear extrapolations to 1/n → 0 give a value for the radius of convergence r µ/T . The behavior of both estimators shown in Fig. 3 is similar at all considered temperatures. Various QCD critical point estimates [44][45][46][47][48] are shown by the black symbols.
The temperature dependence of the radius of convergence, r µ/T , as calculated within the CEM-LQCD (red stars) and the CEM-HRG model (dashed line) models using the Mercer-Roberts procedure, is presented in Fig. 4. r µ/T is a smooth function of T and it is finite, at all temperatures considered. The corresponding limiting singularities lie at complex µ B /T values, as follows from the absence of a regular asymptotic behavior of the Taylor expansion coefficients. r µ/T decreases with temperature and it approaches the asymptotic value of r µ/T = π at higher temperatures, T > 190 MeV. This value can be identified with the Roberge-Weiss (R-W) transition [36], which was predicted to appear at sufficiently high temperatures at imaginary chemical potential values of Im [µ B /T] c = π(2k + 1), and studied quite extensively in LQCD simulations [37][38][39][40][41][42][43]. This transition is a consequence of the R-W periodicity of the QCD partition function, Z(µ B ) = Z(µ B + i2πT ), due to the center symmetry [36], which is fully respected by the CEM.
We have cross-checked our results for r µ/T by constructing Padé approximants [49,50]  It is interesting that numerical lattice studies at purely imaginary µ indicate T RW = 208 ± 5 MeV for the endpoint temperature of the R-W transition [43], a temperature value where r µ/T is already almost indistinguishable from π in CEM-LQCD. We conclude that the radius of convergence of the Taylor series at T > 135 MeV is only determined by the singularities in the complex plane which appear to be smoothly connected to the R-W transition at high temperatures, a scenario suggested in Refs. [8,10]. The CEM-LQCD "knows" about the the spontaneous breaking of the center symmetry at the high temperature R-W transition indirectly, being matched to baryonic excitations at low temperatures and to quark degrees of freedom at high temperatures. The exact nature and relation to the R-W transition of the singularities at intermediate temperatures still need to be clarified. We note that CEM also inherits aspects of the chiral symmetry restoration, in the form of the input coefficients b 1 (T ) and b 2 (T ) taken from the lattice.
In any case, our analysis within CEM-LQCD and CEM-HRG shows no evidence for the existence of a phase transition or a critical point at real µ B /T < r µ/T , with r µ/T ≥ π at all temperatures considered. This is consistent with all available lattice results at zero and imaginary chemical potential, but in contrast to various other QCD critical point estimates available in the literature: these are based on lattice reweighting techniques [44], experimental finitesize scaling analyses [45], the Dyson-Schwinger [46] or holographic [47,48] approaches, which are also shown in Fig. 4. We note that CEM is not full QCD, therefore we do not rule out conclusively these other estimates. Note also that our results at T < 135 MeV are based on the HRG extrapolation of the lattice data, and therefore should be treated with care.
The particular CEM formulation presented here is simple and powerful, but it has limitations. The relation (5) expressing the higher-order Fourier coefficients through the first two is likely to get modified whenever effects of genuine many-body interactions become important. We therefore expect the model to break down at large µ B /T values, e.g., in the dense nuclear matter region. Note that the formalism itself can accommodate any pressure function periodic under the µ B → µ B + i 2πT transformation, as required by the Z(3) symmetry of QCD. The CEM model can thus be extended once new and possibly contradicting lattice data become available. However, given that CEM is consistent with all presently available lattice data we conclude that its range of applicability is at least as large as that of current lattice methods.
To summarize, a novel Cluster Expansion Model for the QCD equation of state has been developed and applied to calculate the baryon number susceptibilities at µ = 0, to very high order. The only model inputs are the partial pressures in the |B| = 1 and |B| = 2 sectors, taken from the lattice simulations at imaginary µ B . The model yields excellent agreement with the available lattice data for χ B 2 , χ B 4 /χ B 2 , and χ B 6 /χ B 2 . The extended model predictions for χ B 8 , χ B 10 , and χ B 12 shall be verified by future lattice data. The commonly used ratio estimator is unable to determine the radius of convergence of the Taylor series of the pressure in µ B /T , due to a non-trivial asymptotic behavior of the Taylor coefficients. The radius of convergence is instead determined with the more elaborate Mercer-Roberts estimator, which provides finite values of the convergence radii at all temperature values considered, 135 < T < 230 MeV, in full agreement with the singularities of Padé approximants. These singularities lie in the complex plane and appear to be smoothly connected to the R-W transition at high temperatures and imaginary (baryo)chemical potential. The analysis within CEM shows no evidence for the existence of a phase transition or a critical point at real values of the baryochemical potential at µ B /T π for temperatures above 135 MeV.
The CEM model can be straightforwardly extended to calculate the equation of state of QCD at finite µ B /T , by supplying the B = 0 partial pressure p 0 (T ) as additional model input. Furthermore, the CEM formalism is rather flexible, and the model assumptions and input can be modified if new and contradicting lattice data becomes available. 2