Net-Baryon Probability Distributions from Lattice Simulations

Abstract

We use the results of lattice simulations of the net-baryon number density at imaginary baryon chemical potential in $N_f=2$ QCD to construct the equation of state of dense and hot strong-interacting matter both above the Roberge–Weiss temperature $T>T_{RW}$ and below the critical temperature $T<T_c$. For these cases, we also evaluate probability distributions of the net-baryon number, as well as the respective cumulants and moments. The consequences of the asymptotic behavior of these probability distributions for the problem of reconstruction of the net-baryon probability distributions from cumulants are discussed.

1. Introduction

The properties of nuclear matter at low temperature and high density and the location of the phase transition to chral-symmetric and/or deconfined quark matter are subjects of both experimental and theoretical studies. With large baryon densities, the only approach based on first principles useful in nonperturbative computations—that is, lattice QCD—is plagued by the so-called sign problem.

For this reason, the generally accepted knowledge in the physics of dense, strong-interacting matter largely comes from low-energy effective models, which are consistent with the first-order chiral phase transition at low temperatures and high baryon densities, as is shown in Figure 1; some experimental evidence for this transition can be found in [1].

There are several methods to obviate the sign problem in lattice studies. In this work, we discuss Taylor expansion at zero baryon chemical potential ${\mu }_B=0$ and the analytical continuation of the quantities computed at imaginary values of ${\mu }_B$, where there is no sign problem.

Our attention is focused on the net-baryon number density, which can be evaluated on a lattice at imaginary ${\mu }_B$. Given the net-baryon number density, one can find the respective pressure and the dependence of the grand canonical partition function on the baryon chemical potential, as well as, therefore, the probability that the system has a definite net-baryon number.

The probability mass function describing the distribution of fireballs produced in heavy-ion collisions in the net-baryon number can be measured experimentally, modulo some plausible assumptions about thermodynamical equilibrium in fireballs. The cumulants associated with this distribution have received considerable attention in the literature [2]. In this regard, considerable theoretical and experimental work has been undertaken on higher-order cumulants corresponding to the distribution of the fireballs produced in heavy-ion collisions in the net-baryon charge.

We discuss the problem of reconstruction of the probability mass function from the cumulants in connection with the phase transition.

On the basis of our lattice results obtained in a simplified version of QCD, we analyze the asymptotic behavior of the distribution ${\mathcal{P}}_n$ in the net-baryon number and study its dependence on the temperature. We find that the decrease in ${\mathcal{P}}_n$ with n as $n\to \infty$ becomes progressively slower as the temperature decreases. We analyze the consequences of this observation for the relation between the cumulants and the probability mass function at low temperatures, where the phase transition is expected.

The results of our lattice simulations agree well with the hadron resonance gas (HRG) model. We try to find signatures of the phase transition that should appear in the behavior of the baryon density at imaginary chemical potentials, where it can be evaluated on a lattice. We discuss in detail possible corrections to the the grand canonical partition function of the HRG model, which are small at the imaginary baryon chemical potential; therefore, however, they are consistent with the first-order phase transition at a real value of ${\mu }_B$. We also consider the equation of state (EoS)— that is, the relation between the pressure and the net-baryon number density— that follows from the results of our lattice simulations. We discuss its changes with temperature and outline a scenario of phase transition in terms of the EoS.

Our study relies on the results of simulations described in [3], which use lattices of size $N_s^3\times 4$ with $N_s=16,20,40$ at several temperature values over the range $0.84T_c<T<1.35T_c$ with $m_{\pi }/m_{\rho }=0.8$.

2. Definitions and Notation

The thermodynamical state of strong-interacting matter is described by the grand canonical partition function

$$Z_{GC}(\theta ,T,V)=\mathrm{T}\mathrm{r}exp\left(-\frac{\widehat{H}-{\mu }_B\widehat{B}}{T}\right),$$

where $\widehat{H}$ is the Hamiltonian, $\widehat{B}$ is the operator of the baryon charge, ${\mu }_B$ is the baryon chemical potential, and ${\displaystyle \theta =\frac{{\mu }_B}{T}}$. The grand canonical partition function can be expressed in terms of the canonical partition functions,

$$Z_C(n,T,V)=\sum _j\langle j|exp\left(-\frac{\widehat{H}}{T}\right)|j\rangle ,$$

(1)

where $|j\rangle$ runs all states with the net-baryon number $B=n$. This expression is as follows:

$$Z_{GC}(\theta ,T,V)=\sum _{n=-\infty }^{\infty }{Z}_C(n,T,V){\xi }_B^n,{\xi }_B=e^{{\mu }_B/T},{\mu }_B=N_c{\mu }_q,\theta ={\theta }_R+ı{\theta }_I,$$

(2)

It is referred to as the fugacity expansion. It should be emphasized that $Z_C(n,T,V)\ge 0$ and, since the baryon charge takes integer values, the grand canonical partition function is periodic in ${\theta }_I$ [4]: $Z_{GC}(\theta ,T,V)=Z_{GC}(\theta +2ı\pi ,T,V)$. The inverse of the fugacity expansion is provided by the Fourier transform

$$Z_C(n,T,V)={\int }_{-\pi }^{\pi }\frac{d{\theta }_I}{2\pi }{e}^{ın{\theta }_I}{Z}_{GC}(\theta =ı{\theta }_I,T,V).$$

(3)

As was shown in [5], the canonical partition functions are of phenomenological significance because they are related to the probabilities ${\mathcal{P}}_n$ that the net-baryon number at given values of ${\mu }_B$ and T equals n:

$${\mathcal{P}}_n\left(\theta \right)=\frac{Z_C(n,T,V){\xi }_B^n}{Z_{GC}(\theta ,T,V)},$$

(4)

We introduce special notation for these probabilities at ${\mu }_B=0$: ${\mathcal{P}}_n\left(0\right)\equiv {\mathbf{P}}_n$, which are of particular interest. The probabilities ${\mathcal{P}}_n\left(\theta \right)$ can be expressed in terms of ${\mathbf{P}}_n$ by the formula

$${\mathcal{P}}_n\left(\theta \right)=\frac{{\mathbf{P}}_n{e}^{n\theta }{Z}_{GC}\left(0\right)}{Z_{GC}\left(\theta \right)},$$

(5)

where

$$Z_{GC}\left(\theta \right)=Z_{GC}\left(0\right)\sum _{n=-\infty }^{\infty }{\mathbf{P}}_n{e}^{n\theta }$$

(6)

In addition to these probabilities, we also consider the equation of state; that is, the relation between the pressure

$$p(\theta ,T)=\frac{T}{V}ln{Z}_{GC}(\theta ,T,V)$$

(7)

and the net-baryon number density

$$\rho (\theta ,T)=\frac{1}{V}\frac{\partial ln{Z}_{GC}(\theta ,T,V)}{\partial \theta };$$

(8)

We also use the respective dimensionless quantities

$$\widehat{\rho }=\frac{\rho }{T^3},\widehat{p}=\frac{p}{T^4},\mathrm{and}\nu =V{T}^3.$$

(9)

Considering ${\mathcal{P}}_n\left(0\right)$ defined in Equation (4) as a discrete probability mass function, one can conclude that the partition function represents the moment-generating function (up to a constant factor); namely

$$\mathfrak{M}\left(t\right)=\sum _n{\mathcal{P}}_n\left(\theta \right)e^{nt}=\sum _{k=0}^{\infty }\frac{t^k}{k!}\sum _n{\mathcal{P}}_n\left(\theta \right)n^k=\sum _{k=0}^{\infty }\frac{{\mu }_k\left(\theta \right)t^k}{k!},$$

(10)

where ${\displaystyle {\mu }_k\left(\theta \right)=\sum _{n=-\infty }^{\infty }{\mathcal{P}}_n{n}^k}$ are the moments of the random variable n. Now, we use Equations (2) and (4) to obtain the relation

$$\mathfrak{M}\left(t\right)=\frac{Z_{GC}(t+\theta ,T,V)}{Z_{GC}(\theta ,T,V)}.$$

(11)

The cumulants of this distribution are given by the coefficients of the Taylor expansion of the cumulant-generating function

$${\mathcal{K}}_{\theta }\left(t\right)=ln{Z}_{GC}(\theta +t,T,V)-ln{Z}_{GC}(\theta ,T,V)=\frac{\left(p\right(\theta +t,T)-p(\theta ,T\left)\right)V}{T}=\frac{V}{T}\sum _{j=1}^{\infty }\frac{{\kappa }_j\left(\theta \right)t^j}{j!},$$

(12)

We also consider ${\varkappa }_j\equiv {\kappa }_j\left(0\right)$; from parity conservation, it follows that ${\mu }_{2n-1}\left(0\right)=0$ and ${\varkappa }_{2n-1}=0$. Of particular interest are a few of the lowest-order cumulants

$$\begin{array}{ccc}& {\kappa }_2={\mu }_2-{\mu }_1^2& {\kappa }_4={\mu }_4-4{\mu }_3{\mu }_1-3{\mu }_2^2+12{\mu }_2{\mu }_1^2-6{\mu }_1^4\hfill \\ \hfill & {\kappa }_3={\mu }_3-3{\mu }_2{\mu }_1+2{\mu }_1^3& {\kappa }_5={\mu }_5-5{\mu }_4{\mu }_1-10{\mu }_3{\mu }_2+20{\mu }_3{\mu }_1^2+30{\mu }_2^2{\mu }_1-60{\mu }_2{\mu }_1^3+24{\mu }_1^5.\hfill \end{array}$$

(13)

3. Lattice Results

First, we write down the lattice results we need for the net-baryon number density $\widehat{\rho }$ at imaginary values of $\theta$ from [3]. It should be noticed that we presented the net-baryon number density in contrast to the net-quark number density in [3].

In this work, we consider the cases $T>T_{RW}$ and $T<T_{RW}$, where $T_{RW}$ is the Roberge–Weiss temperature [4] (see Figure 2 illustrating the physical sense of $T_{RW}$) and $T_c$ is, strictly speaking, the pseudocritical temperature of the chiral crossover at ${\mu }_B=0$. The range $T_c<T<T_{RW}$ is omitted for simplicity.

3.1. $T>T_{RW}$

It was shown in [3] that one can fit the lattice data for the baryon number density ${\rho }_I$ on the segment $-\pi <{\theta }_I<\pi$ with a low-degree polynomial. We used the polynomial

$$-ı{\widehat{\rho }}_I=a_1{\theta }_I-a_3{\theta }_I^3$$

(14)

in the case where $T=1.35T_c$ and

$$-ı{\widehat{\rho }}_I=a_1{\theta }_I-a_3{\theta }_I^3-a_7{\theta }_I^7,$$

(15)

in the case where $T=1.2T_c$. Formula (14) is the lowest-degree polynomial that fits the data at $T=1.35T_c$ well. It is interesting to note that it turns into the net-baryon density formula for free massless quarks when ${\displaystyle a_1=\frac{N_f}{3N_c}}$ and ${\displaystyle a_3=\frac{N_f}{3{\pi }^2{N}_c^3}}$. At $T=1.2T_c$, two-parameter polynomial-fit formulas do not provide $p$-values greater than 0.05, whereas the three-parameter fit formulas (Formula (15) and $a_1{\theta }_I-a_3{\theta }_I^3+a_5{\theta }_I^5$) give rather good fit quality ($p>0.5$). Therefore, the latter formula leads to unphysical behavior at ${\theta }_I=0,{\theta }_R\to \infty$ ($\widehat{\rho }<0$ and $\widehat{p}<0$), and thus we choose the fit Formula (15). We obtain the following values for the fit parameters:

$$\begin{array}{ccc}\hfill T=1.20T_c:& & a_1=0.4989\left(6\right),a_3=0.0127\left(4\right),a_7=1.43\left(13\right)·{10}^{-5}\hfill \\ \hfill T=1.35T_c:& & a_1=0.5189\left(3\right),a_3=0.0122\left(4\right),\hfill \end{array}$$

(16)

3.2. $T<T_c$

At $T<T_c$ and $\theta =ı{\theta }_I$, the net-baryon number density can be fitted well by the trigonometric series

$$-ı\widehat{\rho }(ı{\theta }_I)\simeq f_{1}sin{\theta }_I+f_{2}sin\left(2{\theta }_I\right)+f_{3}sin\left(3{\theta }_I\right)+\dots$$

(17)

in which only a few terms survive. The coefficients $f_n$ can be evaluated in some models. It has been shown [3,6,7] that, over the range $T_c<T<T_{RW}$, the cluster expansion model (CEM) agrees with the lattice data well, and one of the parameters of this model, corresponding to the ratio $f_{n+1}/f_n$ at $n\to \infty$, tends to zero as the temperature becomes lower than $T_c$. We consider temperatures $T=0.99T_c$ and $T=0.93T_c$ in the confinement domain in order to study the temperature dependence of the probabilities ${\mathbf{P}}_n$ where the CEM reaches its applicability limit. The values of $f_n$ from [3] are as follows:

$$\begin{array}{ccc}\hfill T=0.93T_c:& f_1=0.0869\left(3\right)& f_2=0.0\widehat{p}\left(0\right)=0.569;\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill T=0.99T_c:& f_1=0.2541\left(8\right)& f_2=-0.0053\left(7\right)\widehat{p}\left(0\right)=0.692.\hfill \end{array}$$

(19)

Here, $\widehat{p}(0,T)$ are the values of the pressure at ${\mu }_B=0$ evaluated with the formulas from [8]; these values are needed for the EoS. Then, the pressure $\widehat{p}(\theta ,T)$ can be readily determined from $\widehat{\rho }$ by integration with respect to ${\theta }_I$.

4. Equation of State

To derive the EoS in the case $T>T_{RW}$, we consider the fit function shown in Formula (14) at both $T=1.35T_c$ and $T=1.2T_c$. Thus, we obtain the EoS at $T>T_{RW}$ in the parametric form,

$$\begin{array}{ccc}\hfill \widehat{\rho }& =& a_1\theta +a_3{\theta }^3\hfill \\ \hfill \widehat{p}& =& \frac{a_1}{2}{\theta }^2+\frac{a_3}{4}{\theta }^4+{\widehat{p}}_0,\hfill \end{array}$$

(20)

where ${\widehat{p}}_0=\widehat{p}(0,T)$.

In the case $T<T_c$, we consider the situation when $f_1\ne 0$ and employ analytical continuation. Thus, we arrive at

$$\widehat{\rho }\left(\theta \right)\simeq f_{1}sinh\theta$$

(21)

with real values for the baryon chemical potential. In this case, the dependence of the pressure $p(\theta ,T)$ on $\theta$ is given by

$$\widehat{p}(\theta ,T)=\widehat{p}(0,T)+f_1\left(T\right){\int }_0^{\theta }sinh\left(x\right)dx=\widehat{p}(0,T)+f_1\left(T\right)(cosh\left(\theta \right)-1)=\frac{ln{Z}_{GC}(\theta ,T,V)}{V{T}^3}.$$

(22)

From Formulas (17) and (22), we readily obtain the relation

$$\widehat{p}(\rho ,T)=\widehat{p}(0,T)+\left(\sqrt{{\widehat{\rho }}^2+f_1^2\left(T\right)}-f_1\left(T\right)\right),$$

(23)

where $\widehat{p}(0,T)$ is the dimensionless pressure corresponding to $\rho =0$. This dependence of the pressure of the density at a given temperature furnishes the equation of state (EoS) for interacting quark–gluon matter obtained in lattice QCD for temperatures well below $T_c$, where only one term in the Fourier expansion of the baryon density survives.

In the case of two-parameter fit, which works at $T=0.99T_c$, we obtain the respective EoS in the parametric form (here, $s=sinh\left(\theta \right)$) as follows:

$$\begin{array}{ccc}\hfill \widehat{\rho }& =& f_{1}s+2f_{2}s\sqrt{s^2+1};\hfill \\ \hfill \widehat{p}& =& f_1(\sqrt{s^2+1}-1)+f_2{s}^2+{\widehat{p}}_0.\hfill \end{array}$$

(24)

We plot these dependencies in Figure 3.

In the left panel, we see that, at the densities of the order of the saturation density, the dependence of the pressure on the net-baryon number is qualitatively the same over the range $0.93T_c<T<1.2T_c$, and the values of the pressure at $\theta =0$ are taken from [8]; the pressure gradually decreases with the temperature. The difference between the EoSs corresponding to the two-parameter fit (Formula (24)) and the one-parameter fit (Formula (23)) at $T=0.99T_c$ is shown in the right panel. The dark blue solid curve corresponds to the case where $f_2$ in Formula (24) is set equal to zero, and the cyan dashed curve corresponds to the negative value of $f_2$, as was obtained in lattice simulations. It deflects upwards from the $f_2=0$ case, which is a consequence of the EoS (Formula (24)) that follows from the expansion shown in Formula (17). The physical reason for such deflection can be explained as follows. In a wide class of statistical models, such as the excluded-volume hadron resonance gas (EV-HRG) model or the CEM, the sign of the coefficients $f_n$ alternates, which is a consequence of the excluded volume, as was shown in [9]. The excluded-volume effect can be formulated in terms of the EoS as follows: the pressure of the gas of finite-volume particles increases with an increase in the density more rapidly than in the case of point-like particles. This is precisely what one sees in the case of a van der Waals gas, as well as in Figure 3.

The (hypothetical) case where $f_2>0$ is shown for $T=0.93T_c$ by a dash-dotted curve; the isotherm deflects downwards as compared with the case $f_2=0$ (corresponding to our data and shown by the green solid curve). It should be noticed that the Fourier expansions of the density with all positive coefficients were obtained in [10] in the Nambu–Jona-Lasinio model, where some evidence for the phase transition was found. The (hypothetical) isotherm of the van der Waals type that would correspond to this phase transition is shown by the dark red dash-dotted line; $AB$ is the physical segment of this isotherm associated with a two-phase mixed state. Thus, in the right panel, we have outlined a plausible scenario for the phase transition in terms of the isotherms. It is important in what follows that positive higher-order coefficients in the Fourier expansion in Formula (17) are associated with the isotherms that can be naturally deformed to the isotherms of the van der Waals type.

5. Asymptotic Behavior of the Net-Baryon Number Distributions

Given the pressure $p({\theta }_I,T)$, we can find the grand canonical partition function,

$$Z_{GC}=exp\left(\frac{p(\theta ,T)V}{T}\right),$$

(25)

with imaginary values for the baryon chemical potential and then compute the canonical partition functions using Formula (3). We present the results of the outlined procedure at high and low temperatures separately.

5.1. $T>T_{RW}$

First, we consider in detail the case where only $a_1$ and $a_3$ are nonzero. A numerical estimate of $Z_C(n,T,V)$ based on high-precision computations leads to negative values of $Z_C(n,T,V)$ at $n>1.6V{T}^3$, which is unphysical. Moreover, the series in Formula (2) with the coefficients calculated with Formula (3) using $Z_{GC}(\theta ,T,V)$ obtained by integration of Equation (8) with the baryon density (Formula (14)) diverges everywhere in the fugacity plane except the unit circle. This divergence stems from the discontinuities of the fit function shown in Formula (14) at ${\theta }_I=(2n+1)\pi ,n\in \mathbb{Z}$ and the general properties of Fourier series.

Yet another way to estimate the canonical partition functions associated with the baryon density $-ı{\widehat{\rho }}_I=a_1{\theta }_I-a_3{\theta }_I^3$ is to perform an asymptotic expansion of the integral in Formula (3) in the limit where

$$n\to \infty andV\to \infty ,\mathrm{whereas}T\mathrm{and}\varrho =\frac{n}{\nu }\mathrm{are}\mathrm{kept}\mathrm{fixed}.$$

(26)

Such an estimate was performed in [11], where the leading and next-to-leading terms of the asymptotic expansion of the canonical partition function in $\nu$ were taken into account, and it was found that

$$Z_C(n,T,V)\simeq exp\left(\nu \left[\frac{a_1}{2}{\theta }_s^2+\frac{a_3}{4}{\theta }_s^4-\varrho {\theta }_s\right]-\frac{1}{2}ln\left(4\pi \nu (a_1+3a_3{\theta }_s^2)\right)\right),$$

(27)

where

$${\theta }_s=\sqrt{\frac{a_1}{3a_3}}\left[{(\sqrt{x^2+1}+x)}^{1/3}-{(\sqrt{x^2+1}-x)}^{1/3}\right]$$

(28)

and ${\displaystyle x={\varrho }_n\frac{\sqrt{27a_3}}{2a_1^{3/2}}}$. At high $\varrho$, we arrive at ${\displaystyle {\theta }_s\sim \sqrt[3]{\frac{\varrho }{a_3}}}$ and

$$Z_C(n,T,V)\simeq exp\left(-\frac{3}{4}\sqrt[3]{\frac{n^4}{\nu a_3}}\right).$$

(29)

It should be noticed that the canonical partition functions defined by this formula give the series shown in Formula (2), which converges everywhere in the fugacity plane and is periodic in ${\theta }_I$ with the period $2\pi$. Therefore, the precision of the approximation of the grand canonical partition function

$$\frac{Z_{GC}(\theta ,T,V)}{Z_{GC}(0,T,V)}=exp\left(\frac{\nu a_1}{2}{\theta }^2+\frac{\nu a_3}{4}{\theta }^4\right)$$

(30)

using the formula associated with the coefficients (Formula (27)) should be investigated in great detail; however, we leave this problem for other research. Here, we only check whether the asymptotic behavior of $Z_{GC}(\theta ,T,V)$ at ${\theta }_I=0,{\theta }_R\to \infty$ defined by Formula (30) corresponds to the asymptotic behavior (Formula (30)) of $Z_C(n,T,V)$ at $n\to \infty$.

The estimate of the grand canonical partition function at large $\theta >0$ associated with the asymptotic behavior of the canonical partition functions (Formula (27)) can be obtained as follows:

$$Z_{GC}(\theta ,T,V)\simeq \sum _{n=-\infty }^{\infty }exp\left(-\frac{3}{4}\sqrt[3]{\frac{n^4}{\nu a_3}}+\theta n\right)\simeq {\int }_0^{\infty }exp\left(-\frac{3}{4{\left(\nu a_3\right)}^{1/3}}{n}^{4/3}+\theta n\right)dn,$$

(31)

where the integral with respect to n is evaluated using the saddle-point approximation. The result is

$$Z_{GC}(\theta ,T,V)\simeq exp\left(\frac{\nu a_3{\theta }^4}{4}\right),$$

(32)

which is in agreement with Formula (30) (the factor of ${\displaystyle exp\left(\frac{\nu a_1{\theta }^2}{2}\right)}$ is associated with the subleading contribution to the pressure). Such agreement gives some evidence for the validity of our conclusion that the asymptotic estimate (Formula (27)) of the canonical partition functions in the limit shown in Formula (26) makes physical sense.

Therefore, of particular importance is the third method of estimating $Z_C(n,T,V)$ or, more precisely, ${\mathcal{P}}_n$ given that the pressure takes the form

$$\widehat{p}=\frac{a_1}{2}{\theta }^2+\frac{a_3}{4}{\theta }^4.$$

(33)

In terms of probability theory, $a_1={\varkappa }_2$ and $a_3={\varkappa }_4$ furnish the second and the fourth cumulants corresponding to the probability mass function ${\mathbf{P}}_n$ to be determined.

To do this, one should begin with finding the moments of the random variable n (net-baryon number), given that only the two cumulants ${\varkappa }_2$ and ${\varkappa }_4$ differ from zero. In view of the general relation ${\displaystyle \mathfrak{C}\left(t\right)=ln\mathfrak{M}\left(t\right)}$ between the cumulant-generating function ${\displaystyle \mathfrak{C}\left(t\right)=\sum _{j=1}^{\infty }\frac{{\kappa }_j{t}^j}{j!}}$ and the moment-generating function (notice that ${\mu }_0\left(\theta \right)=1$) ${\displaystyle \mathfrak{M}\left(t\right)=\sum _{k=0}^{\infty }\frac{{\mu }_k{t}^k}{k!}}$, we obtain

$$\begin{array}{ccc}\hfill \mathfrak{C}\left(t\right)& =& exp\left(\frac{{\varkappa }_2{t}^2}{2}\right)exp\left(\frac{{\varkappa }_4{t}^4}{4}\right)\hfill \\ \hfill & =& \sum _{j=0}^{\infty }\frac{t^{2j}}{2^j}\sum _{m=0}^{\left[\frac{j}{2}\right]}\frac{{\varkappa }_4^m{\varkappa }_2^{j-2m}}{m!(j-2m)!}=\sum _{j=0}^{\infty }\frac{t^{2j}}{2^j}\frac{{ı}^j\sqrt{{\varkappa }_4^j}}{j!}{H}_j\left(-\frac{ı{\varkappa }_2}{2\sqrt{{\varkappa }_4}}\right),\hfill \end{array}$$

(34)

where $H_n\left(x\right)$ is the Hermite polynomial. Using the asymptotic formula for the Hermite polynomials,

$$H_n{\left(x\right)|}_{n\to \infty }\simeq \sqrt{2}exp\left(\frac{x^2}{2}\right){\left(\frac{2n}{e}\right)}^{\frac{n}{2}}cos\left(x\sqrt{2n}-\frac{\pi n}{2}\right){\left(1-\frac{x^2}{2n+1}\right)}^{-\frac{1}{4}},$$

(35)

at large n, we arrive at the asymptotic behavior of the moments

$${\mu }_{2n}\simeq {\left(\frac{8{\varkappa }_4{n}^3}{e}\right)}^{\frac{n}{2}}exp\left(-{\left(\sqrt{n}-\frac{{\varkappa }_2}{\sqrt{8{\varkappa }_4}}\right)}^2\right).$$

(36)

According to the Carleman criterion [12], the probability mass function ${\mathcal{P}}_n$ exists and is unique provided that the series ${\sum }_{n=1}^{\infty }{\left({\mu }_{2n}\right)}^{-\frac{1}{2n}}$ diverges. In the case under consideration, ${\displaystyle {\left({\mu }_{2n}\right)}^{-\frac{1}{2n}}\sim n^{-\frac{3}{4}}}$, and thus the Carleman criterion is met. Thus, we can use the ansatz

$${\mathrm{P}}_j\sim exp(-x{j}^y)$$

(37)

in order to determine the asymptotic behavior of the probability distribution at $n\to \infty$ related to the asymptotic behavior of moments (Formula (36)). We find x and y as follows: we determine the asymptotic behavior of the moments that follows from the ansatz (Formula (37)) using Laplace’s method,

$${\mu }_{2n}={\mathbf{P}}_0+2\sum _{j=1}^{\infty }{j}^{2n}{\mathbf{P}}_j\simeq \int exp\left(-x{j}^y+2nln\left(j\right)\right)dj,$$

(38)

and then compare the obtained formula with Formula (36). The result ${\displaystyle x=\frac{3}{4\sqrt[3]{\nu {\varkappa }_4}}}$, ${\displaystyle y=\frac{4}{3}}$ coincides with Formula (29).

We consider the asymptotic of the moments of the ${\mathbf{P}}_n$ distribution when the term $a_7{\theta }_I^7$ in Formula (15) does not vanish. From the above reasoning, we arrive at the asymptotic behavior as follows: ${\displaystyle {\mathbf{P}}_n\sim exp\left(-\frac{7}{8}\sqrt[7]{\frac{n^8}{\nu a_7}}\right)}$. In the general case, when the quark density is fitted by a polynomial of degree q, we obtain

$${\mathbf{P}}_n\sim exp\left(-\frac{q}{q+1}\sqrt[q]{\frac{n^{q+1}}{\nu a_q}}\right).$$

(39)

As the temperature $T>T_{RW}$ decreases, the degree of the polynomial needed to fit the quark density at imaginary $\theta$ increases. Therefore, the respective probability mass function ${\mathbf{P}}_n$ decreases with n progressively more slowly as the temperature decreases over the range $T>T_{RW}$.

5.2. $T<T_c$

With $T<T_c$ and $\theta =\rho {\theta }_I$, the net-baryon number density can be fitted by Formula (17), in which only two (with $T=0.99T_c$) or one (with $T=0.93T_c$) terms survive. First, we consider the case where $f_1\ne 0$ and employ analytical continuation. Thus, we arrive at $\widehat{\rho }\left(\theta \right)\simeq f_{1}sinh\theta$ with real values for the baryon chemical potential, derive the expression for the pressure (22), and, thus, obtain the grand canonical partition function as well. Then, Formula (3) gives

$${\mathbf{P}}_n=e^{n\theta -\nu f_{1}cosh\theta }{I}_n\left(\nu f_1\right),$$

(40)

where $I_n\left(z\right)$ is the Infeld function of z, which coincides with the Skellam distribution

$$𝒮(n;b,\overline{b})=e^{-(b+\overline{b})}{\left(\frac{b}{\overline{b}}\right)}^{\frac{n}{2}}{I}_n\left(2\sqrt{b\overline{b}}\right)$$

(41)

with the parameters ${\displaystyle b=\frac{\nu f_1}{2}{e}^{\theta }}$, ${\displaystyle \overline{b}=\frac{\nu f_1}{2}{e}^{-\theta }}$.

The HRG model also predicts the Skellam distribution of the net-baryon number [13,14], and our results agree with this prediction provided that b and $\overline{b}$ are the average baryon and antibaryon multiplicities. In this case, our fit parameter is interpreted as follows:

$$f_1\nu =2\sqrt{b\overline{b}},$$

(42)

It should be emphasized that the product $b\overline{b}$ is independent of ${\mu }_B$. This excellent agreement between lattice and HRG predictions gives some evidence that a “naive” analytic continuation of the grand canonical partition function and other thermodynamical quantities to real values of chemical potential with Formula (22) is valid when T is well below $T_c$ so that the baryon density can be fitted by trigonometric series with a few terms only. Therefore, the dependence of the pressure and baryon density on $\theta$ is rather smooth, implying that there is no transition in the ${\mu }_B-T$ plane at $T<T_s$ associated with singular behavior for the baryon density.

The agreement of the lattice data for several low-order cumulants with the HRG model was found in [6]; however, our approach demonstrates the agreement of the net-baryon number probability distribution, which is more informative.

The asymptotic behavior of the probabilities ${\mathbf{P}}_n$ at $n\to \infty$ in the case where exactly J coefficients in Formula (17) do not vanish was obtained in [3],

$${\mathbf{P}}_n\sim \frac{{\left(\nu f_J\right)}^{n/J}}{{\displaystyle \Gamma \left(\frac{n}{J}+1\right)}},$$

(43)

However, this estimate was proven for positive values of J only. In the case of alternating coefficients $f_n$, the Formula (43) is not justified; however, in the case of alternating coefficients $f_n$, the probability mass function ${\mathbf{P}}_n$ declines even more rapidly than in the case where $f_n>0\forall n$.

To justify this statement, one can consider the CEM, the EV-HRG model, etc., where the net-baryon number density at imaginary chemical potentials is described by Formula (17) with an alternating sign for $f_n$ and $|f_n|$ decreasing as a geometric progression or faster. The probabilities ${\mathbf{P}}_{\mathbf{n}}$ corresponding to such quark densities can be evaluated numerically using the algebraic procedures described in [7], both in the case $f_n=\left|f_n\right|$ and in the case of an alternating sign for $f_n$. Such estimates show that the probabilities ${\mathbf{P}}_{\mathbf{n}}$ in the case of an alternating sign for $f_n$ decrease with n more rapidly than those in the case of a constant sign for $f_n$.

It should also be noticed that the alternating sign for $f_n$ was obtained in lattice simulations in QC₂D at $T_c<T<T_{RW}$. In [7], it was shown that fitting the CEM-parameterized fit-function to the lattice data on the imaginary part of the net-baryon number density over the segment $-\pi <{\theta }_I<\pi$ gives alternating coefficients $f_n$.

Using expression (43) and assuming the positivity of the coefficients of the Fourier expansion, we find that the second term of the expansion shown in Formula (17) can be neglected as compared with the first term provided that

$$n\lesssim \frac{e{f}_1^2\nu }{2f_2},$$

(44)

where $\nu =V{T}^3$. Then, assuming that, for the fireball, $\nu \gtrsim 10$ and $f_n$ are given by Formula (18) (for $T=0.99T_c$), we find that deviations from the HRG model should be dramatic at $n\sim 160$, whereas experimental studies deal with $n\sim 40÷60$ [15].

6. Discussion

From Formula (43), it follows that, when we deal with infinite series of positive coefficients $f_n$, the sequence ${\mathbf{P}}_n$ decreases more slowly than any function of the type ${\displaystyle exp\left(-\frac{n}{J}lnn\right)}$, $J\in \mathbb{N}$. The positivity of the coefficients of the Fourier expansion of the density obtained in [10] at low temperatures suggests that ${\mathbf{P}}_n$ shows a very slow decrease with n, though the corresponding fugacity expansion converges. It may occur that the cumulants obtained from such slowly decreasing probability mass functions correspond to a large variety of functions; that is, given the complete set of the cumulants, one cannot unambiguously reconstruct the probability mass function. To be more precise, we should indicate the conditions for a unique reconstruction of ${\mathbf{P}}_n$ from its moments ${\mu }_k$, since the moments are readily determined by the cumulants. The Krein criterion [12] states that the problem of moments becomes indeterminate when

$$\int dx\frac{ln\phi \left(x\right)}{(1+x^2)}>-\infty ,$$

(45)

where $\phi \left(x\right)$ is the probability density function. The rate of decrease in ${\mathbf{P}}_n$ at low temperatures is very close to the line of demarcation between probability mass functions generating determinate and indeterminate moment problems. Therefore, at low temperatures, the problem of reconstruction of the probability distribution of the net-baryon number by the measured cumulants is highly probably an ill-posed problem. In any case, a comprehensive study of possible unambiguity in the reconstruction of the net-baryon number probability distribution is needed.

Therefore, such a distribution should also be measured experimentally. This point of view was advocated in [5], where, on the basis of Formula (4) and the C-parity conservation represented by the relation $Z_C(n,T,V)=Z_C(-n,T,V)$, a procedure for the measurement of the baryon chemical potential was suggested. Since ${\mathcal{P}}_n$ are experimentally measurable quantities and ${\displaystyle \frac{{\mathcal{P}}_n}{{\mathcal{P}}_{-n}}={\xi }^{2n}}$, one may conclude that theoretical expectations of the quantities

$${\mu }_B^{\left(n\right)}=\frac{T}{2n}ln\left(\frac{{\mathcal{P}}_n}{{\mathcal{P}}_{-n}}\right)$$

(46)

coincide with each other; that is, ${\mu }_B^{\left(n\right)}={\mu }_B$, provided that thermodynamical equilibrium in the fireball is approached.

Thus, it should be checked experimentally whether the data on ${\mu }_B^n$ considered as a function of n can be fitted by a constant function. If they can, this constant should be interpreted as ${\mu }_B$. Given ${\displaystyle \theta =\frac{{\mu }_B}{T}}$, one can determine the empirical cumulative distribution function (ECDF) of the net-baryon number at $\theta =0$ using the formula

$${\mathbf{P}}_n=\frac{{\mathcal{P}}_n{e}^{-n\theta }}{{\sum }_{k=-\infty }^{\infty }{\mathcal{P}}_k{e}^{-k\theta }}.$$

(47)

Thus, one can combine data obtained at different values of ${\mu }_B$ (however, at the same temperature) to increase statistics for the determination of the probability mass function ${\mathbf{P}}_n$.

7. Conclusions

The results of the analytical continuation of the net-baryon number density evaluated in $N_f=2$ lattice QCD at imaginary baryon chemical potentials were used to estimate both the equation of state and the asymptotic behavior of the net-baryon number probability distribution ${\mathbf{P}}_n$ when $n\to \infty$.

It was found that ${\mathbf{P}}_n$ decreases sufficiently fast for the fugacity expansion to converge; however, the rate of decrease becomes gradually slower with the decrease in the temperature.

At low temperatures, the rate of decrease becomes so slow that the respective moment problem may become indeterminate—the used lattice data do not imply a definite conclusion. That is, different probability distributions of the net-baryon number may correspond to the same equation of state provided that at least one of them decreases sufficiently slowly such that the respective moment problem becomes indeterminate.

In any case, the problem of reconstruction of ${\mathbf{P}}_n$ from the respective cumulants is either ill posed, even from the mathematical point of view, or highly sensitive to small variations in higher-order cumulants. For this reason, the analysis of experimental data based on the net-baryon number distribution involves additional information compared to that extracted from the set of cumulants.

The dependence of the equation of state on the temperature and fit parameters has also been used to formulate a possible scenario of emergence of the van der Waals isotherms corresponding to the first-order chiral phase transition. However, a check of this hypothesis, as well as the relation of the ambiguity of the reconstruction of the net-baryon number probability distribution from the cumulants to the first-order phase transition, should be the subject of further studies.