The analytic structure of thermodynamic systems with repulsive interactions

TL;DR

This paper analyzes the complex analytic structure of thermodynamic systems with repulsive interactions, revealing singularities and limitations of series expansions, with applications to hadron gases and lattice QCD.

Contribution

It provides a detailed analysis of the analytic structure of the excluded-volume model using Lambert W-function and compares it with other models, highlighting universal features.

Findings

The pressure function can be expressed via Lambert W-function.

The model exhibits a branch cut singularity in the complex plane.

Universal behavior of net baryonic density in terms of polylogarithms.

Abstract

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.

Figures

Equations

$$p\text{leavevmode}=\text{leavevmode}\frac{nT}{1-bn}\text{leavevmode},$$

$$p(T,\lambda )=T\varphi (T)\lambda \exp \left(-\frac{bp}{T}\right).$$

$$\varphi (T)$$

$$p(T,\lambda )=\frac{T}{b}W[b\varphi (T)\lambda ]\text{leavevmode},$$

$$z\text{leavevmode}=\text{leavevmode}W(z)\exp [W(z)]$$

$$W(z)=\sum _{k=1}^{\infty }\frac{(-k{)}^{k-1}}{k!}{z}^k\text{leavevmode},$$

$$p(T,\lambda )=\sum _{k=1}^{\infty }\frac{(-k{)}^{k-1}}{k!}T{b}^{k-1}[\varphi (T){]}^k{\lambda }^k\text{leavevmode}.$$

$$p^{\mathrm{id}}(T,\lambda )\text{leavevmode}=\text{leavevmode}T\varphi (T)\lambda \text{leavevmode}=\text{leavevmode}nT\text{leavevmode},$$

$$R\equiv \frac{p(T,\lambda )}{p^{\mathrm{id}}(T,\lambda )}=\frac{W[b\varphi (T)\lambda ]}{b\varphi (T)\lambda }\text{leavevmode}.$$

$$p(T,{\lambda }_B)=p_M(T)+p_B(T,{\lambda }_B)+p_{\bar{B}}(T,{\lambda }_B),$$

$$p_{{}_M}(T)=T{\varphi }_M(T)\text{leavevmode},$$

$$p_{{}_B}(T,{\lambda }_B)=T{\varphi }_B(T){\lambda }_B\exp \left(\frac{-b{p}_{{}_B}}{T}\right)\text{leavevmode},$$

$$p_{{}_{\bar{B}}}(T,{\lambda }_B)=T{\varphi }_B(T){\lambda }_B^{-1}\exp \left(\frac{-b{p}_{{}_{\bar{B}}}}{T}\right)\text{leavevmode},$$

$${\varphi }_{M(B)}(T)=\sum _{i\in M(B)}\int dm{\rho }_i(m)\frac{d_i{m}^{2}T}{2{\pi }^2}{K}_2\left(\frac{m}{T}\right),$$

$$p(T,{\lambda }_B)=T{\varphi }_M(T)$$

$$+\frac{T}{b}\left\{W[b{\varphi }_B(T){\lambda }_B]+W[b{\varphi }_B(T){\lambda }_B^{-1}]\right\}\text{leavevmode}.$$

$${\lambda }_B^{br1,2}(T)=[-b{\varphi }_B(T)e{]}^{\mp 1}\text{leavevmode},$$

$$\frac{{\mu }_B^{\mathrm{br}}(T)}{T}=\pm \{\ln [b{\varphi }_B(T)]+1\}\pm i\pi (2k+1)\text{leavevmode}.$$

$$p(T,{\mu }_B)=p(T,{\mu }_B=0)+\sum _{k=1}^{\infty }\frac{{\chi }_{2k}(T)}{(2k)!}{\left(\frac{{\mu }_B}{T}\right)}^{2k}\text{leavevmode}.$$

$$\frac{{\mu }_B^{\mathrm{br}}(T)}{T}=\pm \{1+\ln [b{\varphi }_B(T)]\}\pm i\pi \text{leavevmode}.$$

$$r_{\mu /T}=\sqrt{\{1+\ln [b{\varphi }_B(T)]{\}}^2+{\pi }^2}\text{leavevmode}.$$

$$b=1\text{leavevmode}{\text{fm}}^3,b{\varphi }_B(T=155\text{leavevmode}\text{MeV})\simeq 0.026\text{leavevmode}.$$

$$\frac{{\mu }_B^{\mathrm{br}}}{T}\simeq \pm 2.634\pm i\pi \text{leavevmode},$$

$$r_{\mu /T}\text{leavevmode}\simeq 4.1,$$

$$\frac{n_B(T,{\mu }_B)}{T^3}=\sum _{k=1}^{\infty }{b}_k(T)\sinh \left(\frac{k{\mu }_B}{T}\right)\text{leavevmode}.$$

$$b_k(T)=\frac{2}{\pi }{\int }_0^{\pi }\text{Im}\left[\frac{n_B(T,i{\theta }_{B}T)}{T^3}\right]\sin (k{\theta }_B)d{\theta }_B\text{leavevmode}.$$

$$p(T,{\lambda }_B)=p_M(T)+\sum _{k=1}^{\infty }\frac{(-k{)}^{k-1}T{b}^{k-1}[{\varphi }_B(T){]}^k}{k!}{\lambda }_B^n$$

$$+\sum _{k=1}^{\infty }\frac{(-k{)}^{k-1}T{b}^{k-1}[{\varphi }_B(T){]}^k}{k!}{\lambda }_B^{-k}$$

$$=p_M(T)+\sum _{k=1}^{\infty }\frac{(-k{)}^{k-1}2T{b}^{k-1}[{\varphi }_B(T){]}^k}{k!}\cosh \left(\frac{k{\mu }_B}{T}\right)\text{leavevmode}.$$

$$\frac{n_B}{T^3}$$

$$b_k^{\mathrm{ev}}(T)=(-1{)}^{k-1}\frac{2k^k}{k!}\frac{{\varphi }_B(T)}{T^3}[b{\varphi }_B(T){]}^{k-1}\text{leavevmode}.$$

$$b_1^{\mathrm{ev}}(T)$$

$$b_2^{\mathrm{ev}}(T)$$

$$b_3^{\mathrm{ev}}(T)$$

Full text

complex chemical potential singularities, excluded volume model, radius of convergence, Fourier coefficients

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 finite-temperature 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 de Forcrand and Philipsen (2002); D’Elia and Lombardo (2003); Guenther et al. (2017) or the Taylor expansion method Allton et al. (2002); Gavai and Gupta (2008); Kaczmarek et al. (2011); Bazavov et al. (2017). 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 Pisarski and Wilczek (1984); Stephanov (2006) and a suspected critical point at finite baryon density Stephanov (2006). 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 description 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. Becattini (2009); Braun-Munzinger et al. (2016) for a review] and lattice QCD data at moderate temperatures Borsanyi et al. (2014); Bazavov et al. (2014). Common extensions of the ideal HRG model include the incorporation of the repulsive interactions Olive (1981); Rischke et al. (1991); Venugopalan and Prakash (1992); Yen et al. (1997). 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 Vovchenko et al. (2017a); Huovinen and Petreczky (2018) and on Fourier coefficients of net baryon density at imaginary baryonic chemical potential Vovchenko et al. (2017b) as well baryon number fluctuations in heavy-ion collisions Fukushima (2015); Albright et al. (2015).

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 Hauge and Hemmer (1963). 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. Vovchenko et al. (2018) 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

[TABLE]

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 Rischke et al. (1991):

[TABLE]

Here $\lambda=\exp(\mu/T)$ is the fugacity, $\mu$ is the chemical potential, and

[TABLE]

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 $\mu$ are the corresponding independent intensive variables. All thermodynamical functions can be calculated in terms of $p(T,\mu)$ and its partial derivatives.

The solution of Eq. (2) can be written explicitly as Noronha-Hostler et al. (2012):

[TABLE]

in terms of the Lambert W-function Lam defined by the equation

[TABLE]

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)\sim 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 $\lambda$. In the following we consider the principal branch of $W(z)$ only. This principal branch has the following Taylor series representation:

[TABLE]

which follows from the Lagrange inversion theorem, applied to Eq. (5). The fugacity expansion around $\lambda=0$ of the pressure in the EV model therefore reads

[TABLE]

The pressure of the ideal Boltzmann gas ($b=0$),

[TABLE]

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:

[TABLE]

This ratio quantifies the deviations from the ideal gas case. It depends on the dimensionless variable $z=b\,\phi(T)\,\lambda$ only, i.e. $R\equiv 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\lesssim 10^{-1}$. At $z\gtrsim 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_{\rm br}\equiv-e^{-1}$, with a branch cut along the interval $(-\infty,-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 Hauge and Hemmer (1963). $|R(z)|$ is a continuous function of the complex-valued $z$, but the imaginary part of $R(z)$ flips its sign when crossing the branch cut. $|R(z)|\to e$ as $z\to z_{\rm br}$.

III Hadron resonance gas with repulsive baryonic interactions

III.1 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. Vovchenko et al. (2017a); Satarov et al. (2017); Vovchenko et al. (2017b). 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

[TABLE]

where $\lambda_{B}=\exp(\mu_{B}/T)$ and $\mu_{B}$ is the baryonic chemical potential. Here

[TABLE]

where $\rho_{i}(m)$ in Eq. (14) takes into account the finite widths of the resonances while the sum runs over all mesonic (M) or baryonic (B) species.

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:

[TABLE]

III.2 Taylor expansion properties

The branch points of the pressure function (III.1),

[TABLE]

are located exclusively at the negative real axis. Here $\lambda_{B}^{\rm br1}$ corresponds to the branch point associated with the subsystem of baryons [the second term in Eq. (III.1)], while $\lambda_{B}^{\rm br2}$ corresponds to the subsystem of antibaryons [the third term in Eq. (III.1)]. The two singularities with positions related as $\lambda_{B}^{\rm br1}=\frac{1}{\lambda_{B}^{\rm br2}}$ emerge due to the presence of both, baryons and antibaryons, which leads to two different branch cuts, both located at the negative real fugacity axis. These branch cuts are depicted in Fig. 2 for three different cases:

1. (a)

$|\lambda_{B}^{\rm br1}|>1$: the branch cuts do not overlap; 2. (b)

$|\lambda_{B}^{\rm br1}|=1$: the two branch points coincide, $\lambda_{B}^{\rm br1}=\lambda_{B}^{\rm br2}=-1$; 3. (c)

$|\lambda_{B}^{\rm br1}|<1$: the branch cuts have a non-zero overlap.

The locations of the distinct branch points are given (for $k\in\mathcal{Z}$) in terms of the baryochemical potential:

[TABLE]

The pressure function (III.1) can now be written as a Taylor series expansion around $\mu_{B}/T=0$:

[TABLE]

Here the coefficients of the expansion are the baryon number susceptibilities $\chi_{2k}(T)=\partial^{2k}(p/T^{4})/\partial(\mu_{B}/T)^{2k}|_{\mu_{B}=0}$, evaluated at $\mu_{B}=0$.

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 $\chi_{8}^{B}$ Bazavov et al. (2017); Borsanyi et al. (2018). Due to the CP-symmetry of QCD, all odd order susceptibilities vanish at $\mu_{B}=0$.

The radius of convergence of the Taylor expansion (18) is determined by the singularity of the pressure function in the complex $\mu_{B}/T$ plane, which is located in the closest to the expansion point, $\mu_{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 $\mu_{B}^{\rm crit}$, which limits the convergence of the Taylor expansion around $\mu_{B}=0$ Greiner et al. (2012). 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 Allton et al. (2002); Gavai and Gupta (2005); Allton et al. (2005) or in effective models Stephanov (2006); Karsch et al. (2011); Skokov et al. (2011).

The EV-HRG model (17) exhibits no physical phase transition. Thus, it does not have singularities at real values of the baryochemical potential. Nevertheless, the model does contain branch point singularities in the complex plane, their locations are given by Eq. (17). The closest branch points to $\mu_{B}/T=0$ result by setting $k=0$ in Eq. (17):

[TABLE]

These two branch points are symmetric with respect to $\mu_{B}/T=0$, which reflects the symmetry between baryons and antibaryons. The radius of convergence, $r_{\mu/T}$, of the Taylor expansion in the EV-HRG model is given by the distance of these symmetric branch points to $\mu_{B}/T=0$:

[TABLE]

The Taylor expansion (18) does converge only in the region $|\mu_{B}|/T<r_{\mu/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. Vovchenko et al. (2017b) are used:

[TABLE]

This yields the branch points (19):

[TABLE]

while the radius of convergence (20) becomes equal

[TABLE]

i.e. $r_{\mu}\simeq 635$ MeV.

Figure 3 depicts the $\mu_{B}/T$ dependence of the subtracted scaled pressure $[p(T,\mu_{B})-p(T,0)]/T^{4}$, evaluated within the EV-HRG model at $T=155$ MeV. The full analytic solution [Eq. (10)] (solid black line) is shown, as well as the Taylor expansion (18), truncated at $\chi_{2}^{B}$ (dashed red line), $\chi_{4}^{B}$ (dot-dashed green line), $\chi_{6}^{B}$ (double-dot-dashed blue line), and $\chi_{20}^{B}$ (dotted grey line). The full analytic result is described fairly well by the Taylor expansion, if it is truncated at the $\mathcal{O}(\mu_{B}^{4})$ order or higher, for $\mu_{B}/T\lesssim r_{\mu/T}\simeq 4.1$. 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 $\mu_{B}/T<r_{\mu/T}$, but worse at $\mu_{B}/T>r_{\mu/T}$ outside the convergence radius, as can be seen in Fig. 3. Note that the divergence of the Taylor expansion at large real $\mu_{B}/T>r_{\mu/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.

III.3 Fourier coefficients

The QCD net baryon density $n_{B}$ can be written as a series in hyperbolic sines,

[TABLE]

This general representation is a consequence of the CP- and Roberge-Weiss Roberge and Weiss (1986) symmetries of QCD. For purely imaginary chemical potentials $\mu_{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:

[TABLE]

These Fourier coefficients have attracted considerable attention recently Vovchenko et al. (2018); Almasi et al. (2018); Bzdak and Koch (2019); Almási et al. (2019), in particular in the context of lattice QCD simulations at imaginary $\mu_{B}$ Lombardo (2006); Bornyakov et al. (2017); Vovchenko et al. (2017b).

The four leading coefficients were analyzed in Ref. Vovchenko et al. (2018) within the EV-HRG model, in the context of lattice data. The analytic expression (III.1) determines the exact expressions for $b_{k}$ to arbitrary order in the EV-HRG model. The Taylor expansion of $W(z)$ [Eq. (6)] yields

[TABLE]

The scaled net baryon density $n_{B}/T^{3}=\partial(p/T^{4})/\partial(\mu_{B}/T)$ reads

[TABLE]

with the Fourier coefficients

[TABLE]

The four leading Fourier coefficients read

[TABLE]

They agree with those obtained in Ref. Vovchenko et al. (2017b). The closed-form 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!\approx\sqrt{2\pi k}\,(k/e)^{k}$ yields the following large $k$ asymptotics:

[TABLE]

The Fourier coefficients are exponentially damped, at large $k$, as long as the following condition is fulfilled:

[TABLE]

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 $|\lambda_{B}^{\rm br1}|<1$. Such a behavior contradicts the Riemann-Lebesgue lemma Ryzhik and Gradshteĭn (1965), which stipulates that Fourier coefficients of any function which is integrable on the imaginary $\mu_{B}/T$ interval $[0,\pi]$ vanish for large $k$, $b_{k}\stackrel{{\scriptstyle k\to\infty}}{{\to}}0$. This contradiction appears to be related to the divergence of the series in Eq. (III.3) for purely imaginary values of the baryochemical potential, $|\lambda_{B}|=1$, used to evaluate $b_{k}$. In fact the integration endpoint $\theta_{B}=\pi$ in Eq. (25) for $|\lambda_{B}^{\rm br1}|<1$ lies on the branch cuts of both $W$-functions which enter Eq. (III.3) [see Fig. 2(c)]. Therefore, Eq. (28) is expected to coincide with the Fourier coefficients evaluated through (25) only when the condition (34) is fulfilled simultaneously.

The Fourier coefficients can be evaluated numerically through Eq. (25) to cross-check these results with Eq. (28). Both results agree for $|\lambda_{B}^{\rm br1}|\geq 1$ only, but they disagree for $|\lambda_{B}^{\rm br1}|<1$. For the latter case, the numerically calculated $b_{k}$, Eq. (25), show an asymptotic behavior $b_{k}\sim(-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.

IV.1 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$ Olive (1981). 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 Huovinen and Petreczky (2018). Similar results were achieved by the EV-HRG model Vovchenko et al. (2017b). 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:

[TABLE]

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 W-function:

[TABLE]

The analytic properties of the MF model are determined by the analytic properties of the Lambert W-function, in close analogy to the EV model. The branch point of the MF-model thermodynamic potential is located at

[TABLE]

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 Huovinen and Petreczky (2018) for details). The resulting net baryon density reads ($K>0$):

[TABLE]

Similar to Eq. (16) for EV interactions, the MF model (38) used here possesses two branch points

[TABLE]

located at the negative real axis. Here $\lambda_{B}^{\rm br1}$ corresponds to baryons and $\lambda_{B}^{\rm br2}$ 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) of $W$:

[TABLE]

The asymptotic behavior is the following:

[TABLE]

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 $|\lambda_{B}^{\rm br1}|>1$.

IV.2 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. Vovchenko et al. (2018, 2019). 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

[TABLE]

Here $\hat{b}_{1,2}=\displaystyle\frac{b_{1,2}(T)}{b_{1,2}^{\rm SB}}$, $x_{\pm}=\displaystyle-\frac{\hat{b}_{2}}{\hat{b}_{1}}\,\lambda_{B}^{\pm 1}$, $\textrm{Li}_{s}(z)=\displaystyle\sum_{k=1}^{\infty}\frac{z^{k}}{k^{s}}$ is the polylogarithm, and the

[TABLE]

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

[TABLE]

The singularities are located on the negative real axis, if $\hat{b}_{1}/\hat{b}_{2}>0$. Lattice data suggests $\hat{b}_{1}>0$ and $\hat{b}_{2}>0$ for $T>135$ MeV Vovchenko et al. (2017b) (lattice data are presently not yet available for $T<135$ MeV).

The Fourier coefficients in the CEM read (see Vovchenko et al. (2018))

[TABLE]

with the following asymptotic behavior:

[TABLE]

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.

IV.3 The ideal Fermi gas

The Fermi-Dirac and Bose-Einstein quantum statistical effects can be associated with “effective” repulsive (fermions) or attractive (bosons) interactions Greiner et al. (2012). 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 Greiner et al. (2012)

[TABLE]

The series representation in the last line of Eq. (IV.3) is valid for $m>0$.

The net baryonic density of the ideal Fermi gas has two singularities at

[TABLE]

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. (IV.3).

Note that the relativistic ideal Bose gas, as for example an ideal gas of positively and negatively charged pions, exhibits singularities in the fugacity $\lambda_{Q}$ connected to the conserved electric charge. These are located at the positive axis

[TABLE]

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 (IV.3) allows to evaluate the Fourier coefficients of the net baryon density in an ideal gas of baryons and antibaryons:

[TABLE]

The asymptotic behavior of these Fourier coefficients is given by the following expression:

[TABLE]

This asymptotic behavior is exactly the same as the one found in the mean-field model (41).

IV.4 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 $\mu_{B}/T=0$ equals

[TABLE]

in all these models. Note that here $(\ln|\lambda_{\rm br}^{1}|)^{2}=(\ln|\lambda_{\rm br}^{2}|)^{2}$, i.e. both branch points lie at the same distance from $\mu_{B}/T=0$. It is instructive to consider the behavior of $r_{\mu/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_{\mu/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_{\mu/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$ fm3 for the EV-HRG model Vovchenko et al. (2017b) and $K=350$ MeV fm3 for the MF-HRG model Huovinen and Petreczky (2018), reasonable parameter values suggested by comparisons to the lattice QCD data. Both models predict similar values of $r_{\mu/T}\sim 3-5$ at $T>140$ MeV, reaching the minimum value of $r_{\mu/T}^{\rm min}=\pi$ at $T\simeq 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_{\mu/T}^{\rm min}$ are predicted also by the CEM (blue symbols in Fig. 4) Vovchenko et al. (2018), where the lattice data for the two leading Fourier coefficients Vovchenko et al. (2017b) were used as model input at each temperature value. The radius of convergence in the CEM tends to $\pi$ at high temperatures, which may be associated with a Roberge-Weiss like transition Roberge and Weiss (1986). The results presented suggest that Taylor expansion is likely to be divergent at $\mu_{B}/T>3-5$ and $T>140$ MeV, regardless of existence of the hypothetical chiral critical point of QCD.

IV.5 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:

[TABLE]

as long as $|\lambda_{B}^{\rm br1}|>1$. This asymptotic behavior is determined by a singularity of the net baryon density. The corresponding singular part of $n_{B}(T,\lambda_{B})$ can then be approximated to the leading order:

[TABLE]

as follows from the definition of the Fourier expansion for $n_{B}$ [see Eq. (24)]. Recalling the definition of the polylogarithm

[TABLE]

we arrive at

[TABLE]

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 W-function (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

[TABLE]

It is observed that a single polylogarithm $\operatorname{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, $\operatorname{Li}_{3/2}$ and $\operatorname{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 grand-canonical 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\,\phi(T)\,\lambda$, where $b$ is the excluded-volume parameter, $\phi(T)$ is the ideal gas density at zero chemical potential, and $\lambda\equiv\exp(\mu/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 $\lambda^{\rm br}=[-b\,\phi(T)\,e]^{-1}$. Therefore, the HRG model with baryonic eigenvolumes has a finite radius of convergence of its Taylor expansion around $\mu_{B}/T=0$. This convergence radius is estimated to be $r_{\mu/T}\simeq 4.1$ for a crossover transition temperature ($T\sim 155$ MeV), if a reasonable value is used for the baryonic excluded volume parameter, $b\simeq 1$ fm3.

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:

[TABLE]

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 $\gamma$ 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 $\mu_{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 Taylor- and/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.

Acknowledgments

The work of M.I.G. is supported by the Alexander von Humboldt Foundation and by the Program of Fundamental Research of the Department of Physics and Astronomy of the National Academy of Sciences of Ukraine H.St. acknowledges the support through the Judah M. Eisenberg Laureatus Chair by Goethe University and the Walter Greiner Gesellschaft, Frankfurt.