Critical point signatures in the cluster expansion in fugacities

The QCD baryon number density can formally be expanded into a Laurent series in fugacity, which is a relativistic generalization of Mayer's cluster expansion. We determine properties of the cluster expansion in a model with a phase transition and a critical point at finite baryon density, in which the Fourier coefficients of the expansion can be determined explicitly and to arbitrary order. The asymptotic behavior of Fourier coefficients changes qualitatively as one traverses the critical temperature and it is connected to the branch points of a thermodynamic potential associated with the phase transition. The results are discussed in the context of lattice QCD simulations at imaginary chemical potential. We argue that the location of a branch point closest to the imaginary chemical potential axis can be extracted through an analysis of an exponential suppression of Fourier coefficients. This is illustrated using the four leading coefficients both in a toy model as well as by using recent lattice QCD data.


I. INTRODUCTION
Identification of the phases and structure of strongly interacting matter at finite baryon densities is one of the outstanding issues in modern nuclear physics. It has been established in the framework of lattice QCD that the quark-hadron transition at vanishing baryon density is a smooth crossover [1]. On the other hand, it is expected (although not proven) that a first-order quark-hadron transition takes place at sufficiently high baryon density, with the associated QCD critical point (CP) [2]. The search for a critical behavior at finite µ B is performed using measurements of fluctuations in heavy-ion collisions [3][4][5][6] and also using indirect lattice methods such as Taylor expansion around µ B = 0 [7][8][9][10] or analytic continuation from imaginary µ B [11][12][13]. Currently, the available lattice QCD results show little (if any) hints for a CP [10,14].
In the present work we will consider the above questions in the framework of an expansion which represents the QCD grand canonical potential as a Laurent series in baryon number fugacity λ B ≡ e µ B /T . Formally, it can be viewed as a relativistic extension of Mayer's cluster expansion in fugacities [15]. 1 The net baryon density, ρ B = (∂p/∂µ B ) T , 1 This expansion is often called the "relativistic virial expan- The cluster expansion (2) of net baryon density is particularly interesting in the context of lattice QCD simulations at imaginary µ B . Indeed, ρ B attains a form of a trigonometric Fourier series for a purely imaginary µ B = i θ B T [12,17,18]. The cluster expansion coefficients b k become Fourier coefficients: The four leading Fourier coefficients have been computed at the physical point in Ref. [19]. Recent applications of the Fourier expansion method include a construction of a crossover equation-of-state for finite baryon densities [20][21][22], a determination of the net baryon number distribution in heavy-ion collisions at the LHC [23], and analysis of the scaling properties of b k related to the chiral phase transition [24,25] or repulsive interactions [26].

arXiv:1909.02276v2 [hep-ph] 23 Dec 2019
In the present work we explore how a a critical endpoint of a first-order phase transition at finite baryon density affects the properties of the cluster expansion, in particular the asymptotic behavior of the Fourier coefficients. To this end, we develop a toy model containing a phase transition and a critical point for which one can evaluate all coefficients b k of the cluster expansion explicitly. The behavior of Fourier coefficients associated with the phase transition criticality is elaborated. Based on universality argument, the results obtained are expected to be quite generic for any first-order phase transition with a critical endpoint at finite density. This is additionally demonstrated in Appendix B for a Nambu-Jona-Lasinio description. We then explore a possibility of extracting the location of thermodynamic singularities from a number of leading Fourier coefficients and show that such a procedure is feasible under certain circumstances.

II. TRIVIRIAL MODEL
For simplicity, we consider first a single-component Maxwell-Boltzmann gas of interacting particles. In the context of QCD, the particles can be regarded as abstract baryonic degrees of freedom. The scaled particle number density, n/T 3 , has the following cluster expansion form n(T, µ) The ideal gas limit corresponds to truncating the series at the first term. This is the case for the partial pressure of baryons in the ideal hadron resonance gas model. The CP-symmetry of QCD can be recovered by adding an antisymmetric contribution of antibaryons to Eq. (4). In such a case b k corresponds to the Fourier coefficients defined in Eq. (3). Here we would like to determine how a phase transition at finite density influences the behavior of b k . To achieve this goal, we are looking for a theory containing a phase transition where one can evaluate b k explicitly.
Before proceeding to a model calculation, it is worthwhile to point out the expected large-k behavior of b k based on generic features of power series expansions. The series (4) converges for all complex values of λ inside a circle around the origin which has a radius of |λ r |. |λ r | is the radius of convergence of the power series (4), which corresponds to the distance from the origin to the nearest point λ r in the complex λ plane where n/T 3 cannot be defined as a holomorphic function of λ. We will refer here to such a point as a singularity. This singularity is located on the above-mentioned circle of radius |λ r |. The radius of convergence is encoded in the asymptotic behavior of the expansion coefficients. The general definition of |λ r | is A simple possibility which satisfies (5) is an exponential dependence of b k on λ r in the large-k limit: Below we demonstrate the validity of Eq. (6) in an explicit model calculation.

A. Model definition
Perhaps the simplest theory with a critical point of a first-order phase transition is the van der Waals (vdW) equation, which is given in terms of the pressure as a function of the temperature and particle number density: The grand canonical formulation of the vdW equation was considered in Refs. [27][28][29] to study particle number fluctuations associated with criticality. On the other hand, an explicit vdW model determination of the coefficients b k (T ) to arbitrary order k does not appear to be straightforward. For this reason we consider here a slightly different model, which is obtained by expanding the vdW equation (7) in power series in n and truncating the series at the 3rd order: In this model the pressure is represented as a thirdorder polynomial in the particle number density. For this reason we will call this equation of state the trivirial model (TVM). The qualitative behavior of the TVM isotherms coincides with the one in the standard vdW model: above a certain critical temperature T c the pressure isotherms are monotonically decreasing functions of the specific volume v = n −1 while at T < T c they contain non-monotonic wiggles (see Fig. 1). This implies an existence of a first-order liquid-gas transition with a critical point (CP) in the TVM. The CP location is determined from equations (∂p/∂n) T = 0 and (∂ 2 p/∂n 2 ) T = 0: where we picked only the solution with T c > 0.

B. Grand canonical ensemble
Equation (8) defines the model pressure in the canonical ensemble. As T and n are not the natural variables of the pressure function, Eq. (8) so far does not define the thermodynamic potential. This is achieved through a transformation to the grand canonical ensemble (GCE). This will allow the analysis of the cluster expansion coefficients b k . In order to achieve that, we follow the procedure done in Ref. [27]. First, the free energy F (T, V, N ) is determined from the thermodynamic relation p = −(∂F/∂V ) T,N . Integrating the pressure function (8) and requiring that the free energy reduces to that of an ideal gas in the limit N/V → 0 one obtains with Here d and m are particle's degeneracy factor and mass, respectively, and K 2 is the modified Bessel function. The chemical potential, µ ≡ (∂F/∂N ) T,V , reads The fugacity reads Equation (13) [or, equivalently, (12)] defines the particle number density n(T, µ) in the GCE, i.e. density as a function of T and µ. Substituting n(T, µ) into Eq. (8) then allows one to reconstruct the GCE pressure. Relation (13) is a transcendental equation for the GCE density.
At given values of T and (complex) µ, Eq. (13) may have more than a single solution, meaning that n(T, µ) is a multivalued function. This multi-valuedness translates to the analytic properties of the GCE thermodynamic potential p(T, µ) and is expected to determine the convergence properties of the cluster expansion (4).

C. Branch points
The particle number density n, as defined by Eq. (13), is a multi-valued function of the (complex) fugacity λ at a fixed temperature T . Therefore, n(λ; T ) has branch points, which correspond to the zeroes of the derivative of the inverse function, i.e. (∂λ/∂n) T = 0. Thus, λ br = λ(n br ) are the branch points where n br satisfies with the solution Here ν ≡ a/(bT ) and ν c = a/(bT c ) = 2/( √ 3 − 1). For the subcritical temperatures, T < T c , the branch points n br1,2 are real, whereas for the supercritical temperatures, T > T c , they correspond to a pair of complex conjugate numbers. Expressing ν = ν c + ∆ν where ∆ν < 0 for T > T c and ∆ν > 0 for T < T c , one has A more detailed look shows that the two real roots at T < T c correspond to the spinodal points on the isotherms, see the blue line in Fig. 1. In fact, one can see that Eq. (14) is equivalent to the equation ∂ p(T, n)/∂n = 0 defining the spinodal points. These points correspond to the boundaries separating the mechanically unstable part of the isotherm (∂p/∂n < 0) from the phases of metastable gas (n br1 ) and metastable liquid (n br2 > n br1 ). At T = T c the two roots become degenerate, this corresponds to the critical point. At T > T c the two complex conjugate roots are located in the complex plane away from the real axis. This corresponds to the so-called crossover singularities [30]. A similar behavior of thermodynamic singularities associated with a phase transition has earlier been reported for the vdW equation of state [31]. The locations of the branch points in the fugacity plane are the following where we used Eq. (14). For the complex chemical potential plane one has: Here the second term appears due to the periodicity of the chemical potential in the imaginary axis direction.

D. Cluster expansion coefficients
Let us define the principal branch of the GCE density n(T, λ) by the condition that n is real for real values of λ and that it reduces to the ideal gas density in the dilute limit, i.e. n → φ(T ) λ as λ → 0. This principal branch can then be expressed in a Taylor series of the form (4) around λ = 0. Here we determine the coefficients b k (T ) of the expansion.
The function n(T, λ) is defined by Eq. (13) implicitly, namely as the inverse of the function λ(T, n). Given that λ(T, n) is analytic at the point n = 0 and λ n (T, n = 0) = 0, one can apply the Lagrange inversion theorem [32] to evaluate the series coefficients b k . One obtains: Let us make a variable substitution ω = x 2/(3k): One can recognize the generating function of Hermite polynomials in the r.h.s. of Eq. (20), with t = −[1 − a/(bT )] 2k/3. The higher-order derivatives evaluated at x = 0 in the r.h.s. of Eq. (20) correspond to the Taylor coefficients of the generating function of Hermite polynomials. One obtains The four leading coefficients read The first three coefficients evaluated in the TVM coincide with the vdW model result (see Appendix in Ref. [19]). This is not surprising as the pressures in the two models coincide up to the 3rd power of the particle number density. Starting from the 4th coefficient (26), however, the two models differ.

E. Asymptotic behavior of b k
The asymptotic behavior of the cluster integrals b k determines the convergence properties of the cluster expansion (4). This behavior is determined by the asymptotic properties of Hermite polynomials, which are known. Extra care should be taken here, as both the index and the argument of the Hermite polynomials in (22) tend to large values as k → ∞. In such a case the asymptotic behavior depends on the relative increase rate of the Hermite polynomial index and its argument. These different behaviors were studied in Ref. [33]. There are three cases relevant for our analysis: 1. For x > √ 2n the Hermite polynomials H n (x) admit the asymptotic representation (Theorem 1 in [33]) In our case n = k − 1 and The condition x > √ 2n corresponds to the subcritical temperatures, T < T c . Recalling ν ≡ a/(bT ) and setting ν = ν c + ∆ν one obtains the following for the The asymptotic behavior of the cluster expansion coefficients has the form of an exponential damping superimposed on a power-law suppression. The exponential suppression is determined by the branch point λ br1 , located on the real fugacity axis. 2 2. For x ≈ √ 2n the Hermite polynomials H n (x) admit the asymptotic representation (Theorem 3 in [33]) The notation x ≈ √ 2n here means that x → √ 2n in the limit n → ∞, however x can be different from √ 2n for a finite value of n. The condition x ≈ √ 2n corresponds in the TVM to the critical temperature, T = T c . One obtains: The asymptotic behavior of b k at T = T c has the form of an exponential suppression superimposed on a power-law damping. The exponential part is determined by the critical fugacity value λ c which corresponds to the critical point location.
3. For |x| < √ 2n (i.e. T > T c ), the Hermite polynomials H n (x) have the asymptotic representation (Theorem 5 in [33]) where − π 2 < θ < π 2 . Here n = k − 1 and sin θ = First, one observes that the fugacity values λ br1,2 (17) at the branch points of the thermodynamic potential can be written with The asymptotic behavior of b k at T > T c corresponds to a damped oscillator superimposed on a powerlaw decay. The branch points that define this behavior correspond to the so-called crossover singularities. Denoting λ br = e µ br /T and µ br = µ R br ± i µ I br one can see that |λ br | = e µ R br /T and θ br = µ I br /T . Thus, the real part µ R br /T of the chemical potential at the branch point determines the exponential suppression of the magnitude of b k at T > T c whereas the imaginary part µ I br /T defines the period of oscillations.
The obtained TVM results can be summarized as follows: Here µ sp1 corresponds to the spinodal point of the first-order phase transition which delineates the metastable gaseous phase and the mechanically unstable phase at T < T c , µ c corresponds to the critical point at T = T c , and µ R crs and µ I crs are, respectively, the real and imaginary parts of the chemical potential corresponding to the crossover branch points at T > T c . Figure 2 depicts the k-dependence of b k in the TVM for five different temperatures: two temperatures above the critical one, T = 1.3 T c and T = 1.1T c , the critical temperature, T = T c , and two temperatures below the critical one, T = 0.7 T c and T = 0.9 T c . The coefficients here are divided by the expected power-law, exponential, and amplitude factors from Eqs. (32)- (34). We also set bφ(T ) = 1 in this calculation. The large k behavior of the computed coefficients is consistent with the expected asymptotics. For T < T c the asymptotic behavior is approached monotonically. The rate of approach depends on the distance of an isotherm to the critical one. For example, the reduced coefficients are within 10% of the asymptotic limit at k = 7 for T = 0.9T c and already at k = 2 for T = 0.7T c . For T = T c the large-k limit is reached considerably slower. For T > T c the coefficients exhibit an oscillatory behavior. The period of oscillations is large at temperatures slightly above the critical one and decreases with increasing temperatures. This behavior reflects the increase of the imaginary part µ I br /T of the crossover singularity chemical potential with the temperature at T > T c , in accordance with Eq. (34).
The asymptotic behavior (32)- (34) has been obtained here in the framework of the TVM. Nevertheless, due to a universality of the critical behavior this result is expected to be the same for any theory with a phase transition and a critical point at finite density which belongs to the mean-field universality class, with nonuniversal constants A − , A c , A + , and θ 0 . In particular, we checked numerically that this holds for the vdW model [Eq. (7)], through an evaluation of a large number of leading b k coefficients in that model. Additionally, in Appendix B we analyze the behavior of b k in a Nambu-Jona-Lasinio (NJL) model through numerical calculations at an imaginary chemical potential. These calculations confirm that the asymptotic behavior (32)-(34) also holds in NJL. It should be noted that a phase diagram of a statistical system may contain richer structures than those given solely by a critical point that are studied here. These could be, for example, a tricritical point and a line of second-order phase transitions, as expected for QCD in the chiral limit [34,35], or inhomogeneous phases [36][37][38]. How the phase structures of such systems are related to the b k asymptotics is not obvious. The TVM presented here is not suitable in such cases, as the model is only suited to determine features associated with a critical endpoint of a first-order phase transition. Therefore, an analysis of Fourier coefficients within other manageable models possessing these involved phase structures is an interesting future possibility. One possible choice is the Gross-Neveu model in large N f limit [39], which has been used in the past to test proposals to analyze the phase structure of QCD [40], in particular using imaginary chemical potentials [41].
A behavior of Fourier coefficients of net baryon density similar to Eqs. (32)-(34) have recently been obtained in Ref. [25] in the framework of Landau theory of phase transitions applied to the chiral phase transition in the limit of small quark masses, as well as using scaling relations. One notable difference to the present work is a pre-exponential factor of k −2 in the case of a chiral crossover in the mean-field approximation, which is different from the k −3/2 factor obtained here [see Eq. (34)]. The apparent reason for this difference is that Ref. [25] considers the chiral criticality at µ B = 0, where the perturbation in µ B is coupled to the temporal Ising variable t, but where a coupling of µ B to the magnetic field variable h is forbidden by symmetry. For a critical point at finite baryon density the situation is different: µ B is coupled to both the t and h Ising variables, and the variable h is expected to dominate the scaling near the critical point [30]. This behavior is reflected in the TVM, leading to the pre-exponential factor k −3/2 instead of k −2 .
Another important remark is related to the universality class of the critical behavior associated with a critical point of the phase transition. As mentioned above, in the TVM this is the mean-field universality class. The expected universality class for the QCD critical point is Z(2) (3D-Ising) [35,42], which is characterized by somewhat different critical exponents [15]. Therefore, for a universality class different from mean-field one expects a similar asymptotic behavior to the one given in Eqs. (32)-(34), but with corrections to the power-law exponents.
Based on the considerations above, we expect the following asymptotic behavior of b k in a general case: Here µ br = µ R br ±i µ I br is a singularity (a branch point) of the thermodynamic potential which determines the asymptotic behavior of the cluster expansion coefficients. Evidently, this has to be the singularity located the closest to the imaginary µ B axis, as contributions from all other singularities will have a stronger exponential suppression, rendering their contributions to b k subleading. This singularity may not necessarily be connected to a critical point of a phase transition at finite density studied here.
The exponent α depends on the nature of the singularity (universality class, critical point, spinodal or crossover, etc.). While the only singularities in the TVM are those related to the phase transition, in a more general case the form (35) can also accommodate singularities not related to physical phase transitions (see Ref. [26] for a number of examples).

III. EXTRACTING THERMODYNAMIC SINGULARITIES FROM FOURIER COEFFICIENTS
The TVM introduced above can be used to model a hypothetical phase transition and a critical point at finite baryon density. In such a case one can associate interacting particles in the TVM with abstract baryonic degrees of freedom. The net baryon density reads 3 where which implies The form (38) coincides with the relativistic cluster expansion [Eq. (2)] meaning that b k correspond to the Fourier coefficients of net baryon density at imaginary µ B . The large k behavior of Fourier coefficients associated with a phase transition at finite baryon density is given by Eqs. (32)- (34). Leading Fourier coefficients can in principle be calculated using lattice QCD simulations at imaginary µ B through the Fourier transform. In fact, this has already been done for the four leading coefficients in Ref. [19] on N τ = 12 lattices. The question that we want to address is the following: can one extract useful information about QCD thermodynamic singularities from a number of leading Fourier coefficients based on 3 Here we consider the simplest relativistic generalization of the TVM, where baryon-antibaryon interaction terms are neglected.
the known expected asymptotic behavior? We argue that the answer to this question is affirmative. Moreover, we show that some useful information can be extracted from the already available lattice data. Based on the general asymptotics (35) of Fourier coefficients one notices that by far the strongest effect on the overall magnitude of Fourier coefficients is exerted by the real part µ R br of the branch point closest to the real µ B axis. More specifically, one has where the strongest k-dependence is in the third term. It can be reasonable to expect the appearance of strong exponential suppression of Fourier coefficients already in the leading coefficients. If that is the case, µ R br can be extracted by fitting the absolute magnitudes of a number of the leading b k 's with an ansatz As a proof of concept, we take the TVM for baryons with model parameters fixed in such a way as to obtain a critical point at T c = 120 MeV and µ c = 528 MeV (a = 328 MeV fm 3 , b = 1 fm 3 , d = 10, and m = 938 MeV/c 2 ). We fit the four leading Fourier coefficients in the TVM with the ansatz (40) at three temperatures: T = 100 MeV, T ≈ T c = 120 MeV, and T = 150 MeV. Parameters A and µ R br are fitted while α is fixed to its expected value of 3/2.
The fit results are depicted in Fig. 3 by the blue lines for the k-dependence of ln |b k |. The fitted function provides a reasonable description of the higherorder b k 's that were not used in the fitting procedure. values reproduce the true ones to a fairly good precision (10% or better), even at T > T c where the exponential suppression of b k is superimposed on an oscillatory behavior. The procedure can be improved by including higher-order coefficients into the fit. Omitting a number of leading coefficients from the fit could be helpful as well, as the asymptotic form (39) is less justified for these coefficients than for the higher-order ones.
The exercise shows that even only four leading Fourier coefficients might be sufficient to extract the  real part µ R br of the limiting thermodynamic singularity under certain circumstances. Unfortunately, the method yields no conclusive answer with regards to the nature of the extracted singularity, in particular to the possible presence of an imaginary part. Extraction of µ I br requires an analysis of the possible oscillatory behavior of b k 's which might require knowledge of a considerably larger number of coefficients. Nevertheless, the extracted µ R br in all likelihood serves as a lower bound on the value of the critical chemical potential at a given temperature. Indeed, a presence of a singularity which is closer to the imaginary axis would imply a weaker exponential damping of Fourier coefficients.
We now repeat the fit procedure using the real lattice QCD data for b 1 , b 2 , b 3 , and b 4 [19] for the temperature range 135 < T < 230 MeV. Results are depicted in Fig. 4. In addition to α = 3/2, here we also consider two additional cases: α = 1 and α = 2. This allows to asses the sensitivity of the results to the choice of α, which can be different depending on the nature of singularity. While α ≈ 3/2 is expected for a (crossover) phase transition at finite baryon density, one can also have α = 1 in the case of a Roberge-Weiss transition [43] or α ≈ 2 for a chiral crossover singularity for small quark masses [25]. These considerations lead us to assume the interval 1 ≤ α ≤ 2 as a reasonable bound on α.
The left panel of Fig. 4 depicts the fit results for three different temperatures: T = 135, 170, and 215 MeV. One can see that the exponential suppression is stronger at lower temperatures. At T = 135 MeV fits are not sensitive to the chosen value of α. At T = 175 MeV the differences are seen more clearly: while α = 1 and α = 3/2 fits still describe the data reasonably well, α = 2 does a noticeably worse job. For T = 215 MeV, the α = 2 case appears to be ruled out. The temperature dependence of the extracted µ R br /T values is depicted in the right panel of Fig. 4. Fits at T 200 MeV are characterized by small (α = 1) or vanishing (α = 3/2, 2) values of µ R br /T , indicating a possibility of a power-law suppression of Fourier coefficients instead of an exponential one. Such a scenario would correspond to a singularity at a purely imaginary value of the chemical potential. This may be an indication of the Roberge-Weiss transition at imaginary chemical potential at T > T RW [44], where T RW 208 MeV according to lattice QCD estimates [45].
For comparison, we also depict in Fig. 4 the predictions of the cluster expansion model (CEM) of Ref. [20]. The CEM describes the available lattice data on b k within errors, the asymptotic behavior of the Fourier coefficients in this model matches the general form given in Eq. (35) with α = 1. This model predicts µ R br /T = ln b 1 (T )/b 2 (T ) withb 1,2 (T ) being the lattice data for the two leading Fourier coefficients scaled by the high-temperature Stefan-Boltzmann limiting values. The CEM predictions agree quite well with fit results performed for α = 1, especially at larger temperatures, as seen by comparing the open and full circles in Fig. 4.

IV. DISCUSSION AND CONCLUSIONS
We have determined properties of the cluster expansion in fugacities that are associated with a presence of a phase transition and a critical point at finite baryon density. This has been achieved through the trivirial model (TVM) -a model which does contain  The nontrivial behavior of b k associated with the phase transition present in the TVM is encoded in the properties of Hermite polynomials. The asymptotic behavior of b k changes qualitatively as one traverses the critical temperature: at T < T c one observes a monotonic behavior characterized by an exponential suppression of b k at large k [Eq. (32)], which is superimposed on a power-law damping. At T = T c the behavior is similar, with a modification to the powerlaw exponent [Eq. (33)]. The magnitude of the exponential suppression at T = T c is determined by the value of the critical chemical potential µ c which corresponds to the critical point, b k ∼ e −kµc/Tc . At T > T c the thermodynamic branch points move into the complex µ B plane, which leads to an emergence of an oscillatory behavior in addition to the exponential and power-law decays in magnitude [Eq. (34)]. An appearance of negative values of b k above a certain temperature may therefore signal reaching the crossover temperature region, T > T c .
In all cases, the asymptotic behavior of b k is determined in the TVM by the location of the phase transition branch points, and is given in most general case by Eq. (35). Given the universality of the critical behavior, our results obtained in the framework of the TVM are expected to be qualitatively generic for any critical endpoint of a phase transition. In Appendix B we supplement these results with similar findings obtained within a Nambu-Jona-Lasinio description in a mean-field approximation.
In QCD, the cluster expansion properties are particularly interesting in the context imaginary chemical potentials. There, the cluster expansion coefficients become Fourier expansion coefficients of net baryon density. First-principle lattice QCD simulations are free of sign problem at imaginary µ B , and a calculation of a number of the leading Fourier coefficients appears to be feasible. In fact, the four leading coefficients have already been computed on the lattice for temperatures 135 < T < 230 MeV, although the error bars at the smaller temperatures of that range are still quite large.
According to Eq. (35), the magnitude of b k drops exponentially with k, with the slope proportional to the real part µ R br /T of the chemical potential of the phase transition branch point, be it a spinodal point at T < T c , a critical point at T = T c , or a crossover branch point at T > T c . An analysis of this exponential suppression appears to be a fairly reliable way to extract the value of µ R br /T , even just four leading Fourier coefficients might be sufficient, as we show in Sec. III. Our analysis of the available lattice data suggests µ R br /T ≤ 2 − 3 at T > 135 MeV, with decreasing values at higher temperatures (see Fig. 4). These values can serve as a reliable lower bound on the critical point location at T > 135 MeV. Furthermore, given that the available lattice data contains negative Fourier coefficients b k < 0 at all temperatures where the data are available (135 < T < 230 MeV), this disfavors the existence of the critical point at these temperatures. At T 200 MeV the µ R br /T values are small, in some cases even vanishing. This might serve as an indication of the Roberge-Weiss transition at purely imaginary µ B .
In summary, we have presented an analysis of a hypothetical phase transition and a critical point at finite baryon density using the coefficients of the cluster expansion in fugacities. This provides a complementary approach in the hunt for the QCD critical point, in addition to the commonly used methods based on conserved charges susceptibilities, which are either calculated at µ B = 0 in lattice QCD or measured at non-zero µ B in heavy-ion collisions. The results obtained are useful for future lattice QCD simulations at imaginary µ B , which will hopefully yield more accurate values of Fourier coefficients at temperatures where the QCD critical point can be expected. Analysis of the structure of these Fourier coefficients will be able to yield new bounds on the possible location (or even existence) of the QCD critical point. The TVM introduced in this paper (Sec. II) can be viewed as a variant of a real gas equation of state, constructed for a system of particles with short-range repulsive (excluded volume) and intermediate range attractive (mean field) interactions. A generic framework of real gas models, including the effects of quantum statistics, was developed in Ref. [46] and applied to model the nuclear matter. The free energy of a real gas in this framework is the following: Here f (η) is an available volume fraction. It models the short-range repulsive interactions in a form of a generalized excluded volume procedure. η ≡ (bN )/4V and b is the excluded volume parameter. u(n) with n ≡ N/V is an attractive mean-field. Comparing Eq. (41) with the free energy expression (10) in the TVM allows to obtain the explicit TVM expressions for f (η) and u(n): u tvm (n) = −a n.
The real gas formulation of the TVM brings new possible applications. For example, the TVM can be used to model nucleon-nucleon interactions and the nuclear liquid-gas transition. Following the generic procedure described in Ref. [46] one can fix the parameters a and b of nucleons to reproduce the saturation density n 0 = 0. 16 which is in a reasonable agreement with empirical estimates [47]. The TVM can also be used to incorporate the baryon-baryon interactions and the associated nuclear liquid-gas criticality into a hadron resonance gas model, in a similar way as it was done in Ref. [48] for the vdW equation.

B. Fourier coefficients in an NJL model
This Appendix presents the behavior of Fourier coefficients in an NJL model [49,50]. The NJL model is a low-energy effective theory of QCD [51], which has been used in the past to study the phase structure of QCD, in particular that associated with the critical behavior [52,53]. The model exhibits a chiral critical point and a first-order phase transition at finite net quark number densities. Therefore, it can be interesting to consider the behavior of the Fourier coefficients in NJL, in particular to verify the asymptotic behavior (32)-(34) of b k associated with the critical point, which was obtained in the framework of the TVM and which we expect to be model-independent.
We take a mean-field variant of the NJL model for 2 flavors and 3 colors and neglect the vector repulsion. The quark chemical potential µ q plays the role of the chemical potential µ. The grand potential reads [51] Ω(T, µ; M ) = − 12 2π 2 Λ 0 k 2 dk k 2 + M 2 The model parameters are the momentum cut-off Λ, the bare quark mass m 0 , and the scalar coupling G S . The constituent quark mass M at given T and µ is determined by minimizing the grand potential. This is defined by the gap equation: Two solutions of the gap equation may merge at a branch point. The branch points are defined from ∂ 2 Ω ∂M 2 = 0, (48) and have the same physical meaning as the thermodynamic branch points introduced in Sec. II within the TVM. The branch point coordinates at a given temperature can be determined in the NJL model by solving numerically Eqs. (47) and (48). Two branch points merge at the critical point. This corresponds to Equations (47) The behavior of Fourier coefficients b k is studied by computing the Fourier integrals of net quark number density at imaginary µ through a numerical integration: b k (T ) = 2 π π 0 Im ρ(T, iθ T ) T 3 sin(k θ) dθ . (51) The evaluation of the net quark number density ρ(T, µ) at imaginary µ is done in two steps. First, the effective mass M at a given µ is computed by minimizing the grand potential Ω. This is achieved by solving the gap equation (47). Then, the density is computed as ρ(T, µ) = −(∂Ω/∂µ) T . Calculation results for b k are depicted in Fig. 5, for three different temperatures: a supercritical temperature of T = 360 MeV, the critical temperature, T = T c = 120 MeV, and a subcritical temperature of T = 100 MeV. The coefficients are scaled by the expected asymptotic power-law and exponential factors from Eqs. (32)- (34), for which we determine µ R crs , µ c , and µ sp1 numerically, by solving Eqs. (47) and (48): 368 MeV for T = 100 MeV. The scaled coefficients quickly flatten for T < T c and T = T c (the two lower panels in Fig. 5) whereas at T > T c they show an oscilla-tory behavior (dashed line in Fig. 5), as predicted by the TVM [Eq. (34)]. The numerical NJL model results thus confirm the analytic TVM predictions for the asymptotic behavior of b k .