The analytic structure of thermodynamic systems with repulsive interactions

Thermodynamic properties of systems with repulsive interactions, are considered in the grand canonical ensemble. The analytic structure of the excluded-volume model in the complex plane of the system chemical potential (fugacity) is elaborated, based on the fact that the pressure function can be given in terms of the Lambert W-function. Even though the excluded volume model has no phase transitions at real values of the chemical potential, it does exhibit a branch cut singularity in the complex plane, thus limiting the convergence range of the Taylor expansion in the chemical potential. Close similarities to analytic properties of the other models with repulsive interactions, such as a cluster expansion model, the mean-field model, and the ideal Fermi gas model, are pointed out. As an example, repulsive baryonic interactions in a hadron gas, with a focus on the fugacity/virial and Taylor expansion methods used in lattice QCD, are presented. The asymptotic behavior of the Fourier expansion coefficients in these various models suggests that the singular part of net baryonic density can to leading order be universally expressed in terms of polylogarithms.


I. INTRODUCTION
Properties of strongly interacting matter and determination of its different phases are the key questions which drive the heavy-ion collision experiments as well as finitetemperature lattice QCD simulations. The first-principle lattice methods are restricted to simulations at zero or imaginary chemical potentials due to the sign problem. The indirect lattice methods to probe finite baryon densities are based on extrapolations such as analytic continuation from imaginary chemical potential [1][2][3] or the Taylor expansion method [4][5][6][7]. Both methods are sensitive to the analytic properties of the pressure function in the complex chemical potential plane, in particular its singularities. These restrict the scope of the analytic continuation as well as the convergence radius of Taylor expansion. Knowledge of the possible singularities is thus useful to control the validity and accuracy of both these methods. Often the singularities of the pressure function are associated with phase transitions and critical phenomena. Important examples include the chiral phase transition in the chiral limit of QCD [8,9] and a suspected critical point at finite baryon density [9]. As we show below, however, the pressure function singularities are not necessarily connected to physical phase transitions.
Useful guidance is provided by phenomenological models, which incorporate various symmetries and physical mechanisms expected in a given region of the phase diagram. Here we employ hadron resonance gas (HRG) models which are used to provide a reasonable descrip-tion of the hadronic part of the QCD-matter phase diagram. A wide range of HRG applications includes the description of hadron yields in heavy-ion collisions [see, e.g., Refs. [10,11] for a review] and lattice QCD data at moderate temperatures [12,13]. Common extensions of the ideal HRG model include the incorporation of the repulsive interactions [14][15][16][17]. The relevance of repulsive baryonic interactions in the HRG equation of state has recently been established through an analysis of the lattice gauge theory data on baryon number susceptibilities at zero chemical potentials [18,19] and on Fourier coefficients of net baryon density at imaginary baryonic chemical potential [20] as well baryon number fluctuations in heavy-ion collisions [21,22].
In this paper we study the analytic properties of these models in the complex chemical potential plane. Certain related features, such as the distribution of the Lee-Yang zeros, have been studied within the excluded volume (EV) and mean field (MF) models long time ago [23].
Here we present a more complete picture, using the fact that the grand canonical thermodynamic functions in these two models can be expressed in terms of the Lambert W-function. We consider also the cluster expansion model from Ref. [24] and the ideal Fermi gas, in addition to the EV and MF models. All these distinct models are found to exhibit a quite similar analytic structure of their grand-canonical thermodynamic potentials. The results are discussed in light of possible applications to a reasonable analysis of lattice QCD data.
The paper is organized as follows. Section II goes in detail through the analytic solution of a single-component Maxwell-Boltzmann gas with EV interactions. Section III explores the analytic properties of the HRG with EV interactions. Section IV compares a number of distinct HRG models with repulsive interactions. Summary in Sec. V closes the article.

II. SINGLE COMPONENT EXCLUDED VOLUME MODEL
The pressure of a single-component Maxwell-Boltzmann gas with a van der Waals-type EV correction is given by where T and n are the system's temperature and particle number density, respectively, and b is the excluded volume parameter. The pressure (1) in the grand canonical ensemble (GCE) is presented in terms of the following transcendental equation [15]: Here λ = exp(µ/T ) is the fugacity, µ is the chemical potential, and Here d and m are particle's degeneracy factor and mass, respectively, and K 2 is the modified Bessel function. The pressure plays the role of the thermodynamical potential in the GCE, T and µ are the corresponding independent intensive variables. All thermodynamical functions can be calculated in terms of p(T, µ) and its partial derivatives. The solution of Eq. (2) can be written explicitly as [25]: in terms of the Lambert W-function [26] defined by the equation for any complex number z. Therefore, the representation (4) provides the analytic continuation of the EV model pressure function into the complex fugacity plane. W (z) is, in general, a multi-valued function. On the principal branch, W (z) is real for real values of the argument z, and W (z) ∼ z for small values of z. The principal branch therefore determines the physical behavior of the EV pressure at real, positive values of the fugacity λ. In the following we consider the principal branch of W (z) only. This principal branch has the following Taylor series representation: which follows from the Lagrange inversion theorem, applied to Eq. (5). The fugacity expansion around λ = 0 of the pressure in the EV model therefore reads The pressure of the ideal Boltzmann gas (b = 0), corresponds to the first term of the fugacity expansion (7).
It is instructive to consider the ratio R of the EV pressure to the ideal gas pressure: This ratio quantifies the deviations from the ideal gas case. It depends on the dimensionless variable z = b φ(T ) λ only, i.e. R ≡ R(z) = W (z)/z. The R(z) dependence is shown for real values of z > 0 in Fig. 1 (left panel). The EV effects are moderate (within 10%) for z 10 −1 . At z 1 the EV effects are quite strong, and the ideal gas picture simply breaks down. These observations are useful as a rule of thumb, to estimate the importance of the EV corrections in various settings, e.g. as in the application of the EV model to the HRG phenomenology.
The contour plot of |R(z)| in the complex z-plane is shown in Fig. 1 (right panel): R(z) exhibits a branch point at z = z br ≡ −e −1 , with a branch cut along the interval (−∞, −e −1 ), which follows from the analytic properties of the Lambert W-function. This branch cut is depicted by the black line. Note that the Lee-Yang zeroes of the EV model are distributed along this branch cut [23]. |R(z)| is a continuous function of the complexvalued z, but the imaginary part of R(z) flips its sign when crossing the branch cut. |R(z)| → e as z → z br .

A. Excluded volume HRG model
The EV approach is often applied to include repulsive interactions between hadrons in the HRG model. The HRG model with EV interactions between (anti)baryons (the EV-HRG model) was developed in Refs. [18,20,27]. This model treats the interactions between pairs of baryons and between pairs of anti-baryons, but not between any other pairs of hadrons, by excluded volume (EV) correctioná la van der Waals. These interactions are quantified by vdW-type eigenvolume parameter b. The pressure in the EV-HRG model reads where λ B = exp(µ B /T ) and µ B is the baryonic chemical potential. Here where ρ i (m) in Eq. The previous section has shown that the explicit form of the pressure in the EV-HRG model is given in terms of the Lambert W-function:

B. Taylor expansion properties
The branch points of the pressure function (15), are located exclusively at the negative real axis. Here The analytic structure of the EV-HRG model pressure function is depicted in the complex fugacity plane for (a) |λ br1 The blue and red lines with the points depict the branch cuts, the blue one corresponds to the branch cut in the 2nd term of Eq. (15) and the red one to the branch cut in the 3rd term of Eq. (15). The dashed curves correspond to purely imaginary values of the baryochemical potential in the range 0 < Im [µB/T ] < π, the integration contour in Eq. (25).
The locations of the distinct branch points are given (for k ∈ Z) in terms of the baryochemical potential: The pressure function (15) can now be written as a Taylor series expansion around µ B /T = 0: Here the coefficients of the expansion are the baryon number susceptibilities χ 2k (T ) The presentation (18) is quite general and is applied here for the QCD equation of state. The leading susceptibilities have been computed in lattice QCD simulations. Current data are available for susceptibilities up to χ B 8 [7,28]. Due to the CP-symmetry of QCD, all odd order susceptibilities vanish at µ B = 0.
The radius of convergence of the Taylor expansion (18) is determined by the singularity of the pressure function in the complex µ B /T plane, which is located in the closest to the expansion point, µ B /T = 0.
Thermodynamic singularities are often associated with phase transitions. For example, the critical endpoint of a first-order phase transition manifests itself as a singularity at real finite µ crit B , which limits the convergence of the Taylor expansion around µ B = 0 [29]. This fact has been used in various attempts to constrain the location of the critical point of QCD by numerical evaluation of a few leading coefficients with lattice QCD [4,30,31] or in effective models [9,32,33].
The EV-HRG model (17) These two branch points are symmetric with respect to µ B /T = 0, which reflects the symmetry between baryons and antibaryons. The radius of convergence, r µ/T , of the Taylor expansion in the EV-HRG model is given by the distance of these symmetric branch points to µ B /T = 0: The Taylor expansion (18) does converge only in the region |µ B |/T < r µ/T . For illustration, the behavior of the Taylor expansion (18) is studied in the EV-HRG model, where as an example T = 155 MeV and the model parameters from Ref. [20] are used: This yields the branch points (19): while the radius of convergence (20) becomes equal i.e. r µ 635 MeV. However, the behavior of the pressure function cannot be reliably described beyond the convergence radius by a truncated Taylor expansion, no matter how high is its order. Moreover, the agreement of the partial sums in Eq. (18) with the exact result becomes, with an increasing number of their terms, better at µ B /T < r µ/T , but worse at µ B /T > r µ/T outside the convergence radius, as can be seen in Fig. 3. Note that the divergence of the Taylor expansion at large real µ B /T > r µ/T does not at all indicate an emergence of physical effects. That observation does simply reflect the existence of complex chemical potential singularities, which limit the convergence range of a Taylor series.
The present results illustrate that the application of the Taylor expansion method in lattice QCD shall respect these findings and must be done carefully. The convergence ranges of the Taylor expansion method are often restricted just by pressure function singularities, which are not at all related to physical phase transitions.

C. Fourier coefficients
The QCD net baryon density n B can be written as a series in hyperbolic sines, This general representation is a consequence of the CPand Roberge-Weiss [34] symmetries of QCD. For purely imaginary chemical potentials µ B , this expansion becomes trigonometric Fourier series. Here the coefficients b k (T ) are the Fourier coefficients, which can be evaluated in the standard way: These Fourier coefficients have attracted considerable attention recently [24,[35][36][37], in particular in the context of lattice QCD simulations at imaginary µ B [20,38,39].
The four leading coefficients were analyzed in Ref. [24] within the EV-HRG model, in the context of lattice data. The analytic expression (15) determines the exact expressions for b k to arbitrary order in the EV-HRG model. The Taylor expansion of W (z) [Eq. (6)] yields The scaled net baryon density n B /T 3 = ∂(p/T 4 )/∂(µ B /T ) reads with the Fourier coefficients The four leading Fourier coefficients read They agree with those obtained in Ref. [20]. The closedform expression (28) suggests that the alternating sign structure of the Fourier coefficients in the EV-HRG model persists to asymptotically large k. The Stirling approximation k! ≈ √ 2πk (k/e) k yields the following large k asymptotics: The Fourier coefficients are exponentially damped, at large k, as long as the following condition is fulfilled: The corresponding analytic structure of the thermodynamic potential in this case is depicted in Fig. 2(a). In contrast, Eq. (33) implies an exponential growth of the coefficients at large k for |λ br1 B | < 1. Such a behavior contradicts the Riemann-Lebesgue lemma [40], which stipulates that Fourier coefficients of any function which is integrable on the imaginary µ B /T interval [0, π] vanish for large k, b k k→∞ → 0. This contradiction appears to be related to the divergence of the series in Eq. The Fourier coefficients can be evaluated numerically through Eq. (25) to cross-check these results with Eq. (28). Both results agree for |λ br1 B | ≥ 1 only, but they disagree for |λ br1 B | < 1. For the latter case, the numerically calculated b k , Eq. (25), show an asymptotic behavior b k ∼ (−1) k−1 /k.

IV. COMPARISON TO OTHER APPROACHES
The EV model is only one particular framework to treat repulsive interactions between particles. A comparison with the other approaches is instructive as it permits to establish the analytic properties of the generic features of all distinct repulsive interaction models presented here.

A. The mean-field approach
In the simplest version of a MF approach the interactions between particles are modeled through a common shift of the single-particle energies which is proportional to the number density by U = K n [14]. The relations K > 0 and K < 0 correspond to repulsive and attractive interactions, respectively. Such an approach has recently been used to model repulsive baryonic interactions in the HRG in the context of the lattice data on baryon number susceptibilities [19]. Similar results were achieved by the EV-HRG model [20]. In case of the Maxwell-Boltzmann statistics, the particle number density n of a single-component system in the GCE is given by the following transcendental equation: The similarity of Eq. (35) to the transcendental equation for the pressure (2) in the EV model is evident. The solution of (35) is given in terms of the Lambert Wfunction: The analytic properties of the MF model are determined by the analytic properties of the Lambert Wfunction, in close analogy to the EV model. The branch point of the MF-model thermodynamic potential is located at This singularity is located on the negative real axis for K > 0 (repulsive mean field) and on the positive real axis for K < 0 (attractive mean field). This result suggests that strong attractive interactions can lead to experimentally observable physical singularities.
The MF model can be used to model repulsive interactions between pairs of baryons and between pairs of antibaryons in the same fashion as was done in Sec. III for the EV model (see [19] for details). The resulting net baryon density reads (K > 0): Similar to Eq. (16) for EV interactions, the MF model (38) used here possesses two branch points located at the negative real axis. Here λ br1 B corresponds to baryons and λ br2 B to antibaryons, as in the EV-HRG model before. The Fourier coefficients of the net baryon density can be evaluated in the MF model using the Taylor series representation (6) The asymptotic behavior is the following: This asymptotic behavior is similar to the EV model. However, the MF model has a different power-law factor, namely k −3/2 , instead of k −1/2 which appears in the EV model. As in the EV model, Eqs. (40) and (41) are valid here for |λ br1 B | > 1.

B. The cluster expansion model
The cluster expansion model (CEM) for the equation of state of QCD-matter at finite baryon density has been introduced recently in Refs. [24,41]. Repulsive baryonic interactions are taken into account as well as the Stefan-Boltzmann limit of massless quarks at high temperatures. This provides a state-of-the-art description of the available lattice data on Fourier coefficients and baryon number susceptibilities. The CEM net baryon density reads are the Fourier coefficients as evaluated in the Stefan-Boltzmann limit of massless quarks. The analytic properties of the CEM are determined by the analytic properties of the polylogarithm. The branch points of the thermodynamic potential are located at The singularities are located on the negative real axis, ifb 1 /b 2 > 0. Lattice data suggestsb 1 > 0 andb 2 > 0 for T > 135 MeV [20] (lattice data are presently not yet available for T < 135 MeV). The Fourier coefficients in the CEM read (see [24]) with the following asymptotic behavior: This asymptotic behavior is similar to the EV and MF models discussed above, but has a power-law factor of k −1 , instead of the k −1/2 factor shown for the EV model or the k −3/2 factor in the MF model.

C. The ideal Fermi gas
The Fermi-Dirac and Bose-Einstein quantum statistical effects can be associated with "effective" repulsive (fermions) or attractive (bosons) interactions [29]. We analyze the analytic properties of the thermodynamic potential of both the ideal Fermi gases of the baryons and of the antibaryons. The GCE expression for the net baryonic number density of a relativistic ideal Fermi gas of degeneracy d and mass m is presented as [29] The series representation in the last line of Eq. (47) is valid for m > 0. The net baryonic density of the ideal Fermi gas has two singularities at located on the real negative fugacity axis, where the magnitude is determined by the mass of the particles, as follows from the integral representation in Eq. (47). Note that the relativistic ideal Bose gas, as for example an ideal gas of positively and negatively charged pions, exhibits singularities in the fugacity λ Q connected to the conserved electric charge. These are located at the positive axis Therefore, the ideal Bose gas does exhibit real physical singularities, which are connected to the Bose-Einstein condensation.
The behavior of ideal Fermi (Bose) gases is quite similar to the corresponding MF model with a repulsive (attractive) mean field [see Eq. (37)]. As stated, the singularities on the real axis evidently do correspond to the point of the onset of the Bose-Einstein condensation. For fermions, this issue is more subtle, as the singularities found do correspond to complex values of the chemical potential.
The expansion (47) allows to evaluate the Fourier coefficients of the net baryon density in an ideal gas of baryons and antibaryons: The asymptotic behavior of these Fourier coefficients is given by the following expression: This asymptotic behavior is exactly the same as the one  [20], the mean field HRG model (dashed black line) [19], and the cluster expansion model (blue symbols with error bars) [24].
found in the mean-field model (41).

D. Some remarks on the radius of convergence
All models with repulsive interactions considered show very similar analytic structure of the thermodynamic potential. In all cases, the branch points are located at the negative real fugacity axis. The radius of convergence of the Taylor expansion around µ B /T = 0 equals in all these models. Note that here (ln |λ 1 br |) 2 = (ln |λ 2 br |) 2 , i.e. both branch points lie at the same distance from µ B /T = 0. It is instructive to consider the behavior of r µ/T in these various models.
The radius of convergence in the ideal HRG model with quantum statistics is shown in Fig. 4 by the red line as a function of temperature. r µ/T is defined there by the singularity in the Fermi-Dirac distribution function for nucleons and its value is determined by the vacuum mass of nucleons.
The r µ/T values for EV-HRG and the MF-HRG models are shown in Fig. 4 by the black solid and dashed lines, respectively. Here we use b = 1 fm 3 for the EV-HRG model [20] and K = 350 MeV fm 3 for the MF-HRG model [19], reasonable parameter values suggested by comparisons to the lattice QCD data. Both models predict similar values of r µ/T ∼ 3 − 5 at T > 140 MeV, reaching the minimum value of r min µ/T = π at T 190 − 200 MeV. We do note that applicability of these hadron-based models might be questionable at high temperatures and our results there serve mainly for illustration purposes. Similar values of r min µ/T are predicted also by the CEM (blue symbols in Fig. 4) [24], where the lattice data for the two leading Fourier coefficients [20] were used as model input at each temperature value. The radius of convergence in the CEM tends to π at high temperatures, which may be associated with a Roberge-Weiss like transition [34]. The results presented suggest that Taylor expansion is likely to be divergent at µ B /T > 3 − 5 and T > 140 MeV, regardless of existence of the hypothetical chiral critical point of QCD.

E. Modeling the singular part of net baryon density with polylogarithms
The asymptotic behavior of the Fourier coefficients in all examples considered has the form of an exponential decay times a power-law damping: as long as |λ br1 B | > 1. This asymptotic behavior is determined by a singularity of the net baryon density. The corresponding singular part of n B (T, λ B ) can then be approximated to the leading order: as follows from the definition of the Fourier expansion for n B [see Eq. (24)]. Recalling the definition of the polylogarithm we arrive at as the leading order approximation of the singular part of the net baryon density in terms of the polylogarithm. This approximation can be improved further on by considering the higher-order terms in the asymptotic expansion (53) for b k , resulting in additional terms with polylogarithms of higher orders.
We considered an approximation of the Lambert Wfunction (see the EV and MF models) in terms of the polylogarithms as described above as an example. Namely, from an analysis of large k terms in Eq. (6) it follows that It is observed that a single polylogarithm Li 3/2 can approximate W (z) for |z| < 2 to relative accuracy of better than 15%, while the second-order approximation using two polylogarithms, Li 3/2 and Li 5/2 , improves this accuracy to within 2%. The presented resummation of complex chemical potential plane singularities using polylogarithms is useful for phenomenological studies of thermodynamic singularities in QCD.

V. SUMMARY
The analytic properties are studied within distinct approaches to treat repulsive interactions for the grandcanonical Maxwell-Boltzmann gas. Main results are based on an observation that the EV model pressure can be expressed in terms of the Lambert W-function. A single-component Maxwell-Boltzmann gas with an EV correction yields deviations from the ideal gas behavior which depends universally on the dimensionless parameter z = b φ(T ) λ, where b is the excluded-volume parameter, φ(T ) is the ideal gas density at zero chemical potential, and λ ≡ exp(µ/T ) is the fugacity.
The analytic properties of the EV model are fully determined by the properties of the Lambert W-function. The pressure function of the EV model has a regular behavior at all physical values of the fugacity/chemical potential, but exhibits a branch cut singularity in the complex domain, namely at λ br = [−b φ(T ) e] −1 . Therefore, the HRG model with baryonic eigenvolumes has a finite radius of convergence of its Taylor expansion around µ B /T = 0. This convergence radius is estimated to be r µ/T 4.1 for a crossover transition temperature (T ∼ 155 MeV), if a reasonable value is used for the baryonic excluded volume parameter, b 1 fm 3 .
The Lambert W-function is used to determine the explicit form of the Fourier coefficients of the net baryon density [Eq. (28)], which shows an alternating sign behavior in all orders. A number of other theories with repulsive interactions, such as the repulsive mean-field approach, the cluster expansion model, and the ideal gas of fermions, show strong similarities of their analytic properties to the EV model. In particular, the branch points of the pressure functions of all these approaches are all located on the negative real fugacity axis. The asymptotic behavior of the Fourier coefficients does for all these models exhibit the form of an exponential decay times a power-law damping: The magnitude of the exponential suppression is in all cases universally determined by the location of the branch point of the pressure function which can be directly associated with the repulsive interactions, whereas the power-law exponent γ is specific to each model. The alternating signs of the Fourier coefficients look the same in all considered examples and persist to asymptotically large n. The universal asymptotic form (58) allows to approximate the singular part of the net baryon density function in terms of polylogarithms, which is useful for phenomenological studies of thermodynamic singularities in QCD.
The present results are important in particular for the studies of the QCD phase structure, this concerns both the lattice-based methods such as the Taylor expansion of the pressure in µ B /T , as well as the Fourier expansion of the net baryon density at imaginary chemical potential.
In fact, a pressure function singularity which can not be related to a phase transition or a critical point does strongly restrict the convergence radius of the Taylorand/or Fourier expansion methods.
The present work focuses on theories with repulsive interactions only: Hence there is no possibility of a physical phase transition and/or a critical point. It will be interesting to extend these studies within more elaborate phenomenological models of QCD to determine the analytical structure of the pressure function for real values of the baryonic chemical potential, e.g. for the hypothetical case where a phase transition occurs.