Traces of the nuclear liquid-gas phase transition in the analytic properties of hot QCD

The nuclear liquid-gas transition at normal nuclear densities, $n \sim n_0 = 0.16$ fm$^{-3}$, and small temperatures, $T \sim 20$ MeV, has a large influence on analytic properties of the QCD grand-canonical thermodynamic potential. A classical van der Waals equation is used to determine these unexpected features due to dense cold matter qualitatively. The existence of the nuclear matter critical point results in thermodynamic branch points, which are located at complex chemical potential values, for $T>T_c \simeq 20$ MeV, and exhibit a moderate model dependence up to rather large temperatures $T \lesssim 100$ MeV. The behavior at higher temperatures is studied using the van der Waals hadron resonance gas (vdW-HRG) model. The baryon-baryon interactions have a decisive influence on the QCD thermodynamics close to $\mu_B = 0$. In particular, nuclear matter singularities limit the radius of convergence $r_{\mu_B/T}$ of the Taylor expansion in $\mu_B/T$, with $r_{\mu_B/T} \sim 2-3$ values at $T \sim 140-170$ MeV obtained in the vdW-HRG model.


I. INTRODUCTION
The thermodynamic properties of QCD at finite temperatures and densities are important issues of modern high-energy nuclear physics. Of particular interest are the phase structure of QCD matter and the nature of the hadron-parton transition. At zero baryon density, i.e. at µ B = 0, this transition is a crossover, according to lattice QCD simulations [1]. The nature of this transition at finite densities is not established yet. The experimental search for the hypothetical QCD chiral critical point (CP) [2] is performed at non-zero intermediate baryon densities using measurements of fluctuations in heavy-ion collisions [3][4][5][6] as well as indirect lattice gauge theory methods, such as a Taylor expansion around µ B = 0 [7,8] or analytic continuation from imaginary µ B [9,10]. Current high quality lattice QCD data at physical quark masses show no evidence or signatures of a chiral CP and disfavor the existence of a phase transition of first or second order at moderate baryon densities µ B /T π [11][12][13][14]. The location or even the existence of that CP is not settled to date. Possibilities for a phase transition at large baryon densities can be explored in heavy-ion collisions at moderate collision energies, such as the CBM experiment at FAIR [15], or through precision neutron star merger observations and gravitational wave astronomy [16,17].
In contrast to the chiral QCD CP, an existence of the nuclear liquid-gas phase transition with an associated CP at T c 20 MeV and µ c B 900 MeV is better established both theoretically [18][19][20][21][22] and experimentally [23][24][25] (see Ref. [26] for the current empirical estimates of the CP location). This transition is also accessible on the lattice through an effective theory [27]. Recently it has been pointed out that nuclear matter criticality has a sizeable influence on conserved charges susceptibilities in hot QCD matter, both in the vicinity of the crossover temperature region at µ B = 0 [28], and along the phenomenological freeze-out curve in heavy-ion collisions [29][30][31]. In spite of the fact that the nuclear phase transition ends at T c ∼ 20 MeV, its remnants appear to survive in certain observables to much higher temperatures.
A presence of a phase transition and a CP is imprinted in analytic properties of a thermodynamic potential. The pressure function, in particular, becomes a multi-valued function of the chemical potential, and exhibits branch cut singularities [32]. At subcritical temperatures these singularities correspond to spinodal instabilities, at T = T c the singularities merge at the CP, and at T > T c the singularities lie at complex values of the chemical potential [33]. Phase transitions are smoothed out in a finite volume, their remnants are characterized there by the Lee-Yang zeroes of the grand partition function [34,35].
The thermodynamic branch points associated with the nuclear liquid-gas transition are studied in detail in the present work. First, analytic results on the basis of the classical van der Waals (vdW) equation are presented in Sec. II. These are compared at intermediate temperatures (T 100 MeV) with numerical results obtained using quantum vdW, Walecka, and Skyrme models of nu-clear matter (Sec. III). An extrapolation to higher temperatures is achieved in the framework of the vdW-HRG model, with a focus on the influence of the nuclear matter LGPT singularities on convergence properties of the Taylor expansion in µ B /T around µ B = 0 (Sec. IV). Summary in Sec. V closes the article.

II. THERMODYNAMIC BRANCH POINTS OF A LIQUID-GAS PHASE TRANSITION
Let us first consider the system of interacting nucleons as a classical real gas described by the vdW equation. The pressure reads [36] where a > 0 and b > 0 correspond, respectively, to attractive and repulsive interactions. In the grand-canonical ensemble (GCE) the particle number density n(T, µ) is defined by a transcendental equation [37]: Here where d is the degeneracy factor and m is particle's mass 1 . Substituting n(T, µ) into Eq. (1) then allows one to reconstruct the GCE pressure function, i.e. a full thermodynamic potential in the GCE. At given values of T and (complex) µ, Eq. (2) may have more than a single solution, meaning that n(T, µ) is a multi-valued function. This multivalueness entails an existence of branch points. Early studies of the branch points for the classical vdW equation can be found in Ref. [38]. Here we present a systematic analysis of the behavior of branch points related to the nuclear liquidgas transition and their relevance for the QCD phase diagram.
Applied to Eq. (2) this yields Equation (5) is a cubic equation for n br defining the branch points. µ br is recovered by substituting n br into Eq. (2). The cubic equation (5) has three roots which are explicitly obtained using Cardano's formulas: where The third root, n br3 , given by Eq. (7), is real at all values of T and is larger than the limiting density of vdW excluded volume: n br3 > 1/b. Therefore, n br3 is not accessible in the region of physical solutions at any temperature and does not appear to be connected to the existence of the first-order phase transition in the vdW equation. n br3 will thus be omitted from consideration in the following.
The behavior of the two relevant roots (6) depends qualitatively on the value of the temperature. The two roots are real at subcritical temperatures, ∆ < 0 ⇔ T < T c = 8a/(27b). They correspond to the spinodal points of the subcritical isotherms, i.e. (∂p/∂n) T = 0.
At the critical temperature, ∆ = 0, T = T c , the two roots become degenerate. They coincide with the CP location, n br1 = n br2 = n c = 1/(3b).
At supercritical temperatures, ∆ > 0 ⇔ T > T c , the two roots correspond to a pair of complex conjugate numbers, i.e. the singularities lie in the complex plane. This is a manifestation of the so-called crossover transition [32].
In the present work we do only consider the analytic properties of pure phases and will not consider Maxwell's mixed phase construction.

III. BRANCH POINTS OF THE NUCLEAR LIQUID-GAS TRANSITION
In this section temperature dependence of the location of branch points associated with the nuclear liquid-gas transition is evaluated using different models of nuclear matter. The list of models considered is given below.

A. van der Waals
As the simplest model for the nuclear liquid-gas transition we take the classical van der Waals equation (1) for nucleons. We take the vdW parameter values a = 329 MeV fm 3 and b = 3.42 fm 3 from Ref. [39]. These parameter values yield the binding energy of 16 MeV in the nuclear ground state at n = n 0 = 0.16 fm −3 in the vdW model extended to include the Fermi statistics (see below). The locations of branch points are evaluated using Eqs. (2) and (5). The classical vdW equation predicts a nuclear liquid-gas transition with a CP at the following location: The model captures the qualitative features associated with a first-order phase transition but it is not accurate at small temperatures, where Fermi-Dirac statistics cannot be neglected. For the same reasons, the classical vdW equation overestimates the value of the critical temperature by about 10 MeV [39,40].
The quantum statistical effects are taken into account in the quantum van der Waals model (QvdW) [39]. The QvdW model is defined by the following equations: µ * = µ − b p id (T, µ * ) + 2 a n.
Here n id , p id are, respectively, the density and the pressure of the ideal Fermi gas. In the Boltzmann approximation Eqs. (10)- (12) reduce to the classical vdW equations (1) and (2).
Thermodynamic functions at fixed T and µ are usually determined by solving Eq. (12) numerically with respect to (w.r.t.) µ * , which then allows calculating all other quantities. The vdW parameter values are the same as for the classical vdW model above. The QvdW model predicts a CP at T c 19.7 MeV and n c 0.07 fm −3 (µ c 908 MeV).
One needs to evaluate the derivative (∂µ/∂n) T in order to determine the thermodynamic branch points. To do that we apply the derivative w.r.t. n at fixed T to Eqs. (11) and (12), which allows to determine (∂µ/∂n) T explicitly. The resulting equation (∂µ/∂n) T = 0 for the branch points reads 2 a n br T (1 − bn br ) 2 ω id (T, µ * br ) = 1.
Here ω id (T, µ * ) is the scaled variance of particle number fluctuations of an ideal quantum gas in the GCE: with η = +1 for fermions. In the Maxwell-Boltzmann approximation (η = 0) one has ω id = 1 and Eq. (13) reduces to Eq. (5) of the classical vdW model. Here we solve Eq. (13) numerically to determine µ * br 2 . At T = T c the solution of Eq. (13) corresponds to the CP. We use the CP as a starting point of the numerical procedure and move in small steps in temperature independently for T > T c (crossover) and T < T c (first-order phase transition), using the solution at the previous step as an initial guess for the next one. Figure 1 depicts the resulting chemical potential values corresponding to the branch points. At T < T c there are two real solutions which correspond to the spinodals of the first-order phase transition, as discussed in Sec. II for the classical vdW equation. These are depicted in Fig. 1 by two solid lines. At T = T c the two roots become degenerate at the CP. At T > T c , µ br have nonzero imaginary part, the branch points correspond to two complex conjugate roots. The behavior of the real part Fig. 1 by the dashed line.
A comparison between the classical and quantum vdW models can clarify the role of the Fermi statistics. The comparison is exhibited in Fig. 2, where the temperature dependencies of real and imaginary parts of µ br are depicted. Both models exhibit qualitatively similar behavior. The classical vdW model does not yield an accurate description of µ br at small temperatures. This is an expected artefact of neglecting the quantum statistics. Classical vdW model results approach the QvdW model at large temperatures, where effects of quantum statistics become negligible.

B. Skyrme model
To cross-check the robustness of results obtained in the framework of (quantum) vdW model we consider thermodynamic branch points in two alternative models of nuclear matter. In the Skyrme model of nuclear matter the attractive and repulsive interactions are modeled through a mean-field [41,42], which shifts the single-particle energy levels. Here the first term corresponds to intermediate-range attractive interactions and the second term to short-range repulsive interactions. The nucleon number density is given by a self-consistent equation, Here we use the following parameter values: n 0 = 0. 16 In practice, Eq. (17) is solved numerically for a quantity µ * br ≡ µ − u sk (n br ), as Eq. (16) gives n br as an explicit function of µ * br . As for the QvdW model, the CP location is used as a starting point of the numerical procedure to determine the temperature dependence of the branch points.

C. Walecka model
The last nuclear matter model under consideration is the Walecka model [20,43], which is one of the simplest examples of a relativistic mean field theory. The attractive and repulsive interactions are modeled through exchange of scalar σ and vector ω mesons, respectively. The mesonic fields are treated in a mean-field approximation. The interactions lead to an effective shift of the chemical potential µ → µ * and mass m → m * of nucleons, leading to the following form of the grand-canonical thermodynamic potential (pressure) 3 Here c 2 s > 0 and c 2 v > 0 are the coupling parameters corresponding to attractive and repulsive interactions, respectively. The effective chemical potential µ * and effective mass m * are determined from gap equations: Here n s id is the scalar density of an ideal Fermi gas of nucleons. The particle number density is The values of coupling parameters are determined from the nuclear ground state properties (see Ref. [44] for details): c 2 s = 14.6 fm 2 and c 2 v = 11.0 fm 2 . The model predicts nuclear matter CP at T c = 18.9 MeV, n c = 0.07 fm −3 (µ c 909 MeV).
The branch points are determined through Eq. (4). In order to evaluate (∂µ/∂n) T we first note that µ = µ * +   19) and (21). Therefore, (∂µ * /∂n) T is determined by applying the (∂/∂n) T derivative to the gap equations (19) and (20), and solving the resulting system of linear equations for (∂µ * /∂n) T and (∂m * /∂n) T : Here n * ≡ n id (T, µ * ; m * ) and n * s ≡ n s id (T, µ * ; m * ). The branch points equation (4) reads This equation is solved numerically to determine µ * br (the gap equation (20) is used to relate m * and µ * ). At temperatures T 100 MeV, which are probed by relativistic heavy-ion collisions and studied in finitetemperature lattice gauge theory, excitations of hadronic degrees of freedom other than nucleons cannot be neglected. The hot hadronic phase is typically modeled in the framework of the hadron resonance gas (HRG) model. The standard HRG model does not usually incorporate nuclear matter properties and the associated liquid-gas criticality.

D. Comparison between models
Here we employ a vdW-HRG model, which has been introduced in Ref. [28] as a "minimum" extension of the HRG model to incorporate the nuclear liquid-gas phase transition into a HRG picture. In the present paper we follow this "minimum" extension. The vdW-HRG model incorporates vdW interactions for all baryonbaryon (and, by symmetry, all antibaryon-antibaryon) pairs. The parameters a and b are taken to be the same for all baryon pairs. Their values are equal to the vdW parameters of nucleons (Sec. III). This ensures that the vdW-HRG model reduces to the vdW model of nuclear matter in Sec. III when the contributions of baryonic resonances become negligible, as is the case for the low T , large µ B nuclear matter region of the phase diagram.
We do not include vdW terms for baryon-antibaryon pairs as baryon-antibaryon interactions at short range are dominated by annihilations rather than by a repulsive core as in baryon-baryon interactions. Finally, most of the known meson-meson and meson-baryon scatterings are dominated by resonance formation. Such interactions are already incorporated in a HRG picture by including resonances as separate particles. Note also that our particle list has no resonances with |B| = 2, therefore, there is no double-counting of attractive interactions between baryon-baryon and antibaryon-antibaryon pairs.
The pressure in the vdW-HRG model reads: with where µ j = B j µ B + S j µ S + Q j µ Q is the chemical potential for baryon species j, with B j , S j , and Q j being its corresponding quantum numbers. n B and nB are total densities of baryons and antibaryons, respectively.
We neglect the quantum statistical effects for baryons and anti-baryons in the following. As was shown in Sec. III this is a good approximation for temperatures T 80 MeV (see Fig. 2). For µ Q = µ S = 0, the (anti)baryon densities n B(B) (T, µ) are defined by the transcendental equation: Here The sum in Eq. (30) runs over all baryons in the HRG. Densities n B(B) are multivalued functions of µ B . Both baryons and antibaryons lead to an appearance of branch points. Due to the charge conjugation parity symmetry, the corresponding branch points are related to each other through a transformation µ B → −µ B . The branch points of n B(B) are defined as The branch point coordinates are determined through the relations for the classical vdW equation (Sec. II) with a substitution φ(T ) → φ B(B) (T ). Figure 3 depicts the temperature dependence of the real and imaginary parts of the branch cut singularities associated with the nuclear liquid-gas transition, evaluated within the vdW-HRG model (solid black lines) and the vdW model with (anti)nucleons only (dashed red lines). Only limiting, i.e., closest to the µ B = 0 expansion point, singularities are presented. Two symmetric lines in (a) correspond to baryons and antibaryons. Imaginary parts of two complex-conjugated singularities, presented in (b), are equal for baryons and antibaryons. Circles represent the critical points of baryonic and antibaryonic matter. It is seen that the real part decreases with temperature and crosses zero at about T 180 MeV. This implies that vdW interactions become relevant even close to µ B = 0 at sufficiently large temperatures. This indeed was demonstrated for a number of thermodynamic quantities in Ref. [28]. The addition of the baryonic resonances leads to a faster decrease of µ R br /T towards zero. On the other hand, the resonances do not affect the behavior of the imaginary part µ I br /T , at least not within the vdW-HRG model used.

B. Taylor expansion
The presence of thermodynamic branch points leads to a number of consequences regarding the analytic properties of QCD. Of particular interest is the Taylor expansion of the QCD pressure: Here are the baryon number susceptibilities evaluated at µ B = 0. Expansion includes only even orders of chemical potential as follows from the charge conjugation parity symmetry of QCD.
The series (32) converges inside a circle in the complex µ B /T plane. The convergence is limited by a singularity closest to the expansion point, which lies on the border of the circle. A CP is an example of such singularity. A particular feature of a CP is that the singularity lies on the real axis, implying that Taylor expansion coefficients are asymptotically positive at the critical temperature. This fact is used in various attempts to constrain the location of the QCD CP using lattice QCD, by evaluating a number of leading order Taylor expansion coefficients at µ B = 0, verifying that all available coefficients are positive, and using various radius of convergence estimators [12,45,46]. Note that a divergent Taylor expansion can appear even without the presence of physical phase transitions, e.g. in systems with repulsive interactions only [47].
The thermodynamic singularities associated with the nuclear liquid-gas transition do limit the convergence range of Taylor expansion in the vdW-HRG model. The expected radius of convergence in the vdW-HRG model is given by Here µ br B is the location of the limiting singularity. At T > T c this corresponds to the crossover singularities [Eq. (6)], which both lie at the same distance from µ B = 0. At T = T c this is the nuclear matter CP. At T < T c the limiting singularity is the spinodal point which separates the gaseous and mechanically unstable nuclear phases (the right solid curve in Fig. 1). 4 Figure 4 depicts the temperature dependence of the radius of convergence in µ B and in dimensionless µ B /T variables. For crossover temperatures, T ∼ 140 − 170 MeV, the radius of convergence becomes as small as r µ/T ∼ 2 − 3. This indicates a real possibility that convergence properties of Taylor expansion in full QCD at crossover temperatures might be determined by the remnants of the nuclear liquid-gas transition, which manifest themselves in a form of singularities in the complex plane. 4 The branch point at the boundary of the liquid and the mechanically unstable phase (the left solid curve in Fig. 1) does not limit the radius of convergence of Taylor expansion, despite its smaller µ br B value. The reason is that this branch point lies on a Riemann surface different from the one where the expansion point µ B = 0 is. The radius of convergence is unaffected by the third root [Eq. (7)] for the same reason.  We cross-check our vdW-HRG model results by analyzing the convergence properties of Taylor expansion in this model directly. First, we analyze the convergence radius from the behavior of net baryon susceptibilities χ B k at µ B = 0. We calculate χ B k up to the order χ B 120 numerically, using an efficient algorithm described in the Appendix. The computed χ B k values are then used to determine r µ/T through various estimators.
The so-called ratio estimator, r RE n = |c n /c n+1 | 1/2 with c n ≡ χ 2n /(2n)!, fails to provide a useful estimate of r µ/T . This is a consequence of the fact that limiting singularity lies in the complex plane, with a non-zero imaginary part µ I br /T = 0. In such a case the ratio estimator does not converge [13,48].
An accurate estimate for the radius of convergence is obtained using the Mercer-Roberts estimator [49], and the so-called Domb-Sykes presentation [50,51] (see details in Ref. [13]). The resulting values of r µ and r µ/T using 120 Taylor expansion coefficients are depicted by blue symbols in Fig. 4. These values agree with the prior analytic expectations shown by the solid lines, indicating that the Mercer-Roberts estimator converges to the correct value. We also analyze how many coefficients are needed to obtain a meaningful estimate of r µ/T . The calculations suggest that r µ/T at T = 100 − 200 MeV can be estimated with a 10% accuracy using 5 − 10 nonzero Taylor expansion coefficients (see Tab This is shown in the 2nd row of Tab. I. Next, we study the convergence properties of the Taylor expansion by comparing the pressure isotherms evaluated using a truncated Taylor expansion around µ B = 0 and through a full numerical calculation. Figure 5 shows this comparison for the T = 150 MeV isotherm, where a subtracted scaled pressure [p(T, µ B ) − p(T, 0)]/T 4 as a function of µ B /T is analyzed. The truncated Taylor expansion describes well the full result for µ B /T < r µ/T as well as in a small region beyond r µ/T . We verified that this small region shrinks towards zero as more and more expansion terms are included, and that divergence of the series at µ B /T > r µ/T becomes more and more evident. Thus, for µ B /T > r µ/T the Taylor expansion can at best be only viewed as an asymptotic series.
The present calculation incorporates only hadronic degrees of freedom. Of course, mechanisms other than the nuclear liquid-gas transition are present in full QCD, which affect the analytic properties of the thermodynamic potential and which are not covered within the vdW-HRG model. This includes, for instance, the QCD transition to quark-gluon degrees of freedom and the associated chiral criticality. Another known QCD transition is the Roberge-Weiss transition at imaginary chemical potential [52], occurring at temperatures T > T RW , where T RW ∼ 208 MeV [53]. These mechanisms will yield additional restrictions on the radius of convergence of Taylor expansion, in addition to those due to nuclear liquid-gas transition alone that we study here. The work reported here, on the other hand, merely demonstrates how the presence of the nuclear liquid-gas transition alone can affect the convergence properties of the Taylor expansion. For similar reasons we do not perform here a comparison of the higher-order cumulants χ B 2k with the lattice data, but merely use these quantities to analyze the behavior of the Taylor-expanded pressure.
It should also be noted that the liquid-gas transition has been treated in this work on a mean-field level, and the associated critical behavior corresponds to the meanfield universality class. Going beyond the mean-field approximation can modify the universality class and the nature of the singularity associated with a critical point of a phase transition [32], which, in turn, will modify the behavior of the higher-order susceptibilities χ B 2k . Studying the thermodynamic singularities beyond the mean field level can, therefore, be an interesting future endeavor, which can be achieved e.g. using renormalization group methods.

V. SUMMARY
The presence of the nuclear liquid-gas transition at temperatures T 20 MeV in QCD leads to the emergence of thermodynamic branch points in the QCD grand potential. These branch points correspond to the spinodals of the first-order phase transition at T < T c , to the critical point at T = T c , and to crossover singularities in the complex µ B plane at T > T c . This qualitative result, obtained analytically within the classical vdW equation, is generic for any arbitrary mean-field description of nuclear matter, as follows from the universality argument for critical behavior. The van der Waals hadron resonance gas model analysis implies that signals from the nuclear liquid-gas transition are clearly visible in analytic properties of QCD even at crossover temperatures and moderate baryochemical potentials. In particular, the radius of convergence of a Taylor expansion reaches in the vdW-HRG model values as small as r µ/T ∼ 2 − 3 for temperatures T ∼ 140 − 170 MeV. Such high temperatures are typically assumed to only be associated with the chiral crossover transition at µ B = 0. However, the present results show that the radius of convergence of the Taylor expansion in QCD at these temperature exhibits clearly the remnants of the nuclear liquid-gas transition, at a region where we expected the signals of chiral criticality. If the hypothetical chiral critical point is located deeply in baryon-rich matter, as indicated by a recent analysis of QCD thermodynamics within the chiral mean-field approach [54], the attempts to locate the QCD CP by using the Taylor expansion method must take great care to distinguish the supposed signals of the conjectured chiral CP from the well established nuclear matter liquid-vapor CP. An even order baryon number susceptibility, χ 2m , is expressed at µ B = 0 through the (2m − 1) order derivative of the baryonic density, n B , w.r.t. µ B /T : Here we used Eq. (25), the fact that ∂p B(B) /∂µ B = ±n B and ∂ k nB/∂ To find an arbitrary order derivative ∂ k n B /∂(µ B /T ) k | µ B =0 we rewrite Eq. (29) for baryon density in the following form: where in the vdW-HRG model which links the n-th order derivative of n B with all of its lower order derivatives. For n = 1 we obtain the following from Eq. (39): By substituting (41) in Eq. (40) one can calculate the second order derivative, ∂ 2 n B /∂(µ B /T ) 2 . The procedure can then be applied iteratively to evaluate all derivatives of n B w.r.t µ B /T up to a desired order. The baryon number susceptibilities χ 2m are evaluated from Eq. (36).