Probing QCD Phase Structure by Baryon Multiplicity Distribution

The canonical partition functions $Z_n$ and the number distributions $P_n$ which are obervable in experiments, are related by a single parameter, the fugacities $\xi=\exp(\mu/T)$. With the charge parity invariance, $Z_n$ and $\xi$ can be determined. Thermodynamic quantities such as the number density susceptibility and the kurtosis are then calculated from the grand canonical partition function $Z(\xi,T)=\sum Z_n(T) \xi^n$, ($n=-N_{\rm max},\cdots, N_{\rm max}$), for any chemical potential $\mu$, although the region over which the results are reliable for these quantities is constrained by $N_{\rm max}$. We then calculate the Lee-Yang zeros, which are the zeros of $Z(\xi)$ in the complex fugacity plane, as poles of $d\log Z(\xi)/d\xi$ by using the Cauchy integral theorem. With the help of a multiple precision library, this method provides any precision required without misidentification of the zeros. We analyse $Z_n$ from the net-proton number distributions recently measured at the Relativistic Heavy Ion Collider (RHIC) by assuming the net-proton number is approximately propotional to that of the baryon after the freeze-out, and calculate the moments. We also evaluate the Lee-Yang zero structures obtained from RHIC data and compare them with those obtained from lattice quantum chromodynamics (QCD) calculations. Possible regions of QCD phase transition lines are estimated from the thermodynamics quantities and the Lee-Yang zeros. We discuss how the limited $N_{\rm max}$ in both experimental and numerical studies affects the reliability of the thermodynamic results and Lee-Yang zeros.


Introduction
When temperature and density are varied, QCD is expected to have a rich phase structure [1]. One of the most important challenges in particle and nuclear physics is to experimentally discover the phases that only appear under extreme conditions and to theoretically understand their nature. Such achievements would not only deepen our understanding of QCD but also extend our knowledge of the early universe and compact stars.
The Relativistic Heavy Ion Collider (RHIC) was built to explore the properties of QCD matter [2]. Recently net-proton multiplicity measurements at RHIC are gaining attention [3,4] because they provide valuable information about the QCD phase diagram [5]. In these measurements, the colliding energy is varied, and trajectories of the produced hot matter in (T, µ B ) plane may pass near the critical region. Event-by-event fluctuations are expected to encompass the critical point, where the correlation length rapidly changes [6,7,8]. In particular, conserved quantities such as the charge or baryon number may reveal possible correlations that existed inside the system before hadronization. See Ref. [9] and references therein.
Usually, data obtained at a given colliding energy is assigned with a set of temperature T and chemical potential µ, which are referred to as chemical freeze-out point, due to the success of a thermal statistical model for hadrons in heavy ion collisions [10].
In this paper, we propose a method by which, we can obtain information of the QCD phase diagram not only at the experimental (T, µ/T ) point, but other values of µ/T . This may seem to be magical, but is possible because we take into account all the multiplicities. Basic idea is the use of CP symmetry to extract canonical partition functions from measured number distribution according to a fundamental equation in statistical mechanics. This idea also provides a model independent method to determine µ/T only using CP symmetry. As we will show later, µ/T obtained in the present method is consistent with those obtained from the chemical freeze-out data, for wide range of colliding energies.
In addition, even without direct lattice QCD calculations in physical chemical potential regions, we can calculate the canonical partition functions, which helps us to understand the QCD at the finite density. The method provides us an approach beyond the Taylor expansion method, namely it is possible to calculate large µ/T regions.
In section 2, we describe how to extract the canonical partition functions, Z n , from experimental data, and show the obtained results. In section 3, we construct the grand partition function Z(µ, T ) from these Z n , and calculate the moments as a function of µ/T . In section 4, we show that we can calculate the Lee-Yang zero structure from the canonical partition functions obtained by the lattice QCD and the RHIC experiment. To our knowledge, this is the first calculation of the Lee-Yang zero for the high energy nuclear collisions. Section 5 is devoted to the concluding remarks. In Appendix, we give detailed numerical data Z n and Moments λ 4 /λ 2 for RHIC data. Part of the results here were reported in several proceedings [11,12].

Constructing the canonical partition functions
The grand partition function Z and the canonical partition functions Z n are related as  Figure 1: Experimental multiplicity data P n , and Z n = P n /ξ n for √ s N N = 39 GeV. ξ is tuned so that Z n = Z −n , and ξ=1.88336.
where ξ = exp(µ/T ) is fugacity, and Here we assume that the number operatorN commutes with H, that is,N is a conserved quantity.N can be any conserved number operators, such as baryon, charge and strangeness. In the following we consider the baryon number as a concrete example. Because of the charge-parity symmetry, Z n defined by Eq.(2) satisfy The multiplicity distributions P n observed in experiments are related to Z n as P n (ξ) = Z n ξ n .
Using Eqs. (3) and (4), we can determine ξ and Z n . In Fig.1 and Table 1, we show as an example P n (ξ) and Z n ξ n for √ s N N = 39 GeV [4]. The data correspond to the 0-5% centrality in Au+Au collisions. Here ξ=1.88336 as given in Table 1.
Figure2 shows the obtained ξ together with that obtained by freeze-out analysis in Refs. [13] and [14]. Errors in ξ due to the breaking of the relation Eq.(3) are less than one percent for all colliding energies. Note that we use here only the multiplicity data and the charge-parity symmetry relation (3).
We assume that the net-proton multiplicity data are approximately proportional to those of the baryon 1 . This approximation is justified if (i) after the chemical freeze-out, the net-proton number is effectively constant, or (ii) a created fireball is approximately isoneutral. See also Sec.3 in Ref. [15].

Moments as a function of µ/T and the QCD transition boundary
From the grand partition function given by Eq.(1), the moments are evaluated using The quantity N max in Eq.(1) is finite because of the measurement statistics (simulation statistics and finite volume) in the experiments (lattice QCD). The finite nature of N max places an upper bound on the chemical potential for which the calculation is reliable. To estimate the effect of the finite N max , we test two cases: 1. the values of the final three Z n (n = N max − 2, N max − 1, N max ) are increased by 15% 2. the final two Z n (n = N max − 1, N max ) are set to zero.
As an example, we plot the number susceptibility λ 2 /N max in Fig.3 as a function of µ/T at √ s NN = 27 GeV for these two cases together with the one constructed from all the Z n . Let us suppose that, as µ increases, we encounter a phase transition or a cross-over, and cross it. In this case, we would expect the structure of the moments to be rough in this area. At the peak position of the lower curve, indicated by an arrow in Fig.3, the center line continues to increase. We write this value asμ(T ), and any transition may occur for µ ≥μ. In other words, this positionμ is a candidate for the lower bound of the real susceptibility peak. We then investigate the behavior of λ 3 /λ 2 and λ 4 /λ 2 . Although higher moments have large effects of the finite N max , they may be a good tool for detecting the transition region of the QCD phase diagram [6]. The former ratio does not indicates a significant structure for µ ≤μ. Around the freeze-out points, λ 4 /λ 2 ∼ 1 (Poisson) and it becomes negative as µ increases. See Fig.4. In Refs. [16,17,8], it is argued that a negative value of λ 4 /λ 2 indicates that the phase transition has been reached. We give λ 4 /λ 2 at √ s N N = 11.5, 19.6, 27, 39, 62.4 and 200 GeV in the appendix.
We take the points where the lower curve of λ 4 /λ 2 becomes negative (indicated by an arrow in Fig.4) as an another candidate for the lower boundμ(T ) and show them in Fig.5.

Lee-Yang zero structure
Next we extend the fugacity ξ to complex values. Lee-Yang zeros (LYZs) α l are zeros of the grand partition function Z(ξ) in the complex fugacity plane; that is,   The distribution of these zeros reflects the phase structure of the corresponding statistical system. Lee and Yang argued that for a finite system, no zeros appear on the real axis. In the thermodynamic limit, the number of zeros becomes infinite and the zeros coalesce onto one-dimensional curves [18]. For the firstorder (second-order) phase transition, the coalescing zeros cross (pinch) the real positive ξ axis [19], whereas for the crossover, they do not reach the real axis. The first pioneering work to calculate the LYZs of the lattice QCD was carried out by Barbour and Bell [20], who found that the zero distribution is considerably different in the confinement and deconfinement phases. In Ref. [21] the authors calculated the LYZs in the complex β = 6/g 2 plane to distinguish between a crossover and a second-order phase transition by investigating the volume dependence. Ejiri pointed out that this approach is very difficult with finite statistics since the phase of det ∆(µ) is mixed with the complex β [22].
Although there have been many phenomenological approaches to extracting information on the QCD phase [9,6,8], no study has yet been carried out to employ experimental data to investigate the QCD phase through the LYZs.
For this purpose it is important to reliably determine the LYZs; that is, all zeros must be found without ambiguity, and their positions in the complex fugacity plane must be determined with high accuracy. We obtain the LYZs as follows. We first map the problem to the calculation of the residue of Z ′ /Z: The left hand side of Eq. (7) is integrated along a contour, and the residues inside the contour are summed according to Cauchy's theorem. Because of the 1 1 Figure 6: Schematic of the cBK contours in the divide-and-conquer search for residues.
symmetry +µ ↔ −µ, if α l is a solution, then 1/α l is also a solution. Therefore, we only need to search for residues inside the unit circle. Figure 6 shows the cut Baum-Kuchen (cBK) shape contours used in the study. Starting from 0 ≤ θ < 2π, and 0 < r ≤ 1 in polar coordinates, the region is divided into four pieces, and the Cauchy integral is evaluated over each section. This divide-and-conquer process is iterated as many times as required (here 20 times). When no residue is found in a section, no further divisions are applied to that section. At each divide-and-conquer level, we check the conservation of the residue sum. The technical details of this approach will be presented elsewhere, including the parallelization.
All calculations were performed using the multiple-precision package, FMlib [23] and the number of significant digits was 50 -100. With this algorithm, we can construct LYZ diagrams from Z n .

Lattice calculations
First, we study the LYZs obtained by lattice QCD simulation. Here we do not distinguish between u and d quarks. Details of calculating Z n by lattice QCD simulation are given in [24], where the fugacity expansion formula [25] plays an essential role in obtaining Z n . We update 11000 trajectories including 3000 thermalization trajectories. The measurement is performed every 10 (20) trajectories for a 8 3 × 4(10 3 × 4) lattice size. A Monte Carlo update is performed with the fermion measure at µ = 0; thus, we avoid the sign problem due to the complex fermion determinant. However, an overlap problem may still exist.
The LYZ diagram of the lattice QCD above the phase-transition temperature is shown in Fig.7. For Z n , we use data with (signal/noise) ≥ 2. The condition Z n = 0 for n = 1, 2 (mod 3) is imposed on the lattice data for Z n [26,27] which guarantees 2π/3-translational invariance for Im(µ)/T . We evaluate the LYZs B . In Ref. [28] it is shown that two widely used methods, the multi-parameter reweighting and Taylor expansion, are consistent for the EoS and number density up to µ q /T ∼ 0.8 and for number susceptibility up to µ q /T ∼ 0.6. This implies that the current lattice QCD calculation should not be considered reliable in the higher density regions. In the right panel of Fig.7, the circle corresponding to |ξ| = exp(−3µ q /T ) with µ q /T = 0.6 is displayed. GeV. "Full" result is calculated by using all Z n . "-1" is calculated by removing a pair with Z n of n = ±N max . For "-2", we remove further pair of Z n . The inset is an enlarged view near the positive real axis, which gives the predicted region of the QCD phase transition in Fig.5.
Because the Z n obtained in the confinement region suffer from significant noise, the LYZ diagram should be considered qualitative. However, despite this, distinctive differences are observed above and below the confinement/deconfinement transition temperature. At T ∼ 1.2T c , a line of the zero accumulation appears at arg(ξ) = ±π/3, which is consistent with the Roberge-Weiss (RW) phase transition. Of course, in order to confirm this is a real RW phase transition, we must go to large volume (large N max ) and check that zeros are accumulated.
Roberge and Weiss discussed the regions of pure imaginary chemical potential, and found that at Im(µ)/T = π/3 + 2kπ/3 (k = 0, 1), a phase transition occurs for T ≥ T RW [29]. The transition is the first-order transition above T RW , and it is easy to detect at experiments. See Ref. [30] for detailed discussions on the order of the RW phase transition. See Refs. [31,32] for more detailed discussion how to occur the Roberge-Weisse phase transition.
Using the same lattice setup as that in the present study, the quantity T RW was estimated to be approximately 1.1T c [33]. This phase transition exerts a

Re[ξ B ]
Full -1 -2 Figure 10: LYZ diagram from RHIC data at √ s NN = 19.6 GeV. "Full", "-1" and "-2" denote the same as in Fig.9 clear effect on the LYZ diagram at T ∼ 1.2T c , while no such effect appears at T T c .
The fugacity, ξ, on the unit circle stands for the pure imaginary chemical potential. The relation between the pure imaginary, zero and real chemical potential are discussed in Refs. [34,35].
In the left panel of Fig.9, we check the volume dependence by plotting 8 3 × 4 and 10 3 × 4 cases; here N max is chosen so that N max /V is approximately the same for V = 10 3 and 8 3 . In this simulation, N max /V ∼ 1.7 fm −3 , i.e., in a cube with one f m on a side, up to 3 × 1.7 quarks are included. See Ref. [36], where the authors estimate N max / √ V for obtaining reliable results as a function of the temperature based on the quark-meson model.
Using the relation, we can see the behavior of the grand partition function, Z(ξ = exp(µ/T )), for the pure imaginary chemical potential. We show this behavior in Fig.8 for several temperatures. Lee and Yang pointed out that phase transition regions correspond to zeros of the grand partition function. In Fig.8, such characteristic behaviors are observed at the Roberge-Weiss phase transition points in the pure imaginary chemical potential above T c . As a consequence, the thermo potential changes rapidly. Because our N max is finite, exact Lee-Yang zeros do not appear on the pure imaginary chemical potential. However, rapid decrease of Z(θ) is seen. Since Z nq = 0 for n q = 1, 2 (mod 3) where n q are the quark number. 2 The symmetry of leads to the relation (9). Whether zeros appear in Eq.(8) on the pure imaginary regions, i.e. on the unit circle in the complex fugacity plane, or not depends on the nature of Z n . In other words, whether the Roberge Weise phase transition occurs or not is the outcome of the dynamics. In Ref. [32], an explicit example of Z n that leads to Roberge Weiss phase transition above T c is given and its mechanism is explained.

Relativistic Heavy Ion Collisions
We next consider the LYZ diagrams obtained from RHIC data. We construct the grand partition function Eq. (1) Fig.9  and 10. To clarify the effect of a finite N max , we also calculate the LYZs while neglecting Z n for n = ±N max ("-1") and while neglecting Z n for n = ±N max and ±(N max − 1) ("-2").
Although some zeros exist on the negative real axis, they do not form a line that clearly characterizes the RW transition, which suggests that the data correspond to temperatures below T RW . Remember that the net-proton number is not an exactly conserved quantity. It is therefore very interesting to construct the Lee-Yang zeros and to study their structure for the conserved charge, such as net charge or strangeness.
In the LYZ diagrams obtained from the RHIC data, no zero appears on the positive real-ξ axis. Possible explanations for this result are that (i) there is no phase transition, but a crossover occurs at these temperatures, (ii) the systems are finite, and/or (iii) the N max values are too small. The size of the fireball produced is comparable to the QCD scale, and thus, explanation (ii) is at least partially correct. To further explore the QCD phase transition, a larger N max must be attained.
Several LYZ points appear on the unit circle for the all RHIC data. To understand the meaning of this, we calculate the LYZ diagram for the Skellam model. This is a simple probabilistic model based on the difference between two Poisson distribution variables N 1 and N 2 , and in our case Z n = I n (a). See Ref. [37]. Here I n is the modified Bessel function; the parameter a is unique, once the averages of N 1 and N 2 are given. In Fig.11, the points with a = 0.666 corresponds to the multiplicity in the case of √ s NN = 200 GeV. Two circles appear inside and outside of the unit circle, and no LYZ point on the unit circle. However, if we increase a artificially, the two circles move closer to the unit circle and several points coalesce on the unit circle. This indicates that a gross feature of the RHIC multiplicity data is configured by the probabilistic (red) a = 10. In both cases, N max is set to be 15, as in Fig. 9 origin, but the LYZ distribution includes additional information on the QCD dynamics.
Finally, we estimate where the LYZs would intercept or approach the positive real axis as the volume increases, which indicates the QCD phase transition. These zones are indicated by double -headed arrows in the inset of Fig.9. In  Fig.5, the corresponding regions are indicated by horizontal lines. If the multiplicities P n with larger n were to be measured in future experiments, the QCD phase transition could be pinpointed with more precision. Further study of the relation between the baryon and proton number distributions will improve the analysis here [38].

Concluding Remarks
A simple but important relation discussed in this paper is e.g., the fugacity expansion of the grand partition function Z with the coefficients Z n , which are the canonical partition functions. We have shown how to construct Z n from experimental net-multiplicities. In high energy heavy ion central collision, a fireball with µ and T is a good picture, at least as a global feature. Then the selection of an event with a net-multiplicity is a filter of Z n .
In lattice QCD simulation, we consider the path integral formula of the grand partition function, where det D, N f and S G are fermion determinant, the number of flavor and gluon action, respectively. We expand det D(µ) in a fugacity series. Then we get a formula (11). The canonical partition functions, Z n , depend on only the temperature T , but not on the chemical potential µ. Therefore once we extract Z n from Z(µ, T ), we can evaluate Z at different values of µ ′ . This is very important for exploring the QCD phase boundary: So far, analyses such as the moments have been done on the freeze-out points. But the freeze-out points realized in the BES experiments are in the confinement regions, and not very near to the phase transition. Indeed, Fukushima estimated the baryon density on the freeze-out line using the hadron resonance gas model, and found that the maximum of the baryon density is realized at √ s N N = 8 GeV, but its value is only slightly larger than the normal nuclear density ρ 0 [39]. Using our formula, (11), we can probe higher chemical potential regions from Z n constructed at the chemical freeze-out point.
We calculated the moments (5), and saw their behavior when increasing µ/T . We checked the effects of the finite N max .
The formula (11) makes it possible to calculate Z(µ, T ) at complex values of µ. We studied the Lee-Yang zero distribution data of both the experiments and the lattice simulations. To our knowledge, this is the first Lee-Yang zero analysis of the RHIC data. Lattice results told us that the zeros corresponding to Roberge-Weiss transitions appear above T c . We did not see such structure in the LYZ of RHIC data. This suggests that these experimental data are produced below T RW , where T RW is Roberge-Weiss transition temperature that is around 1.1T c .
There were claims that the RW transition formulated in the Euclidean space can not be interpreted as a physics object in the Minkowski space by considering its cosmological consequences [40,41]. In this paper, we discussed the relation between the RW phase transition and the high energy heavy ion experiment, which hopefully advance understanding of the QCD phase diagram 3 .
The approach investigated here is based on the simple statistical mechanics, and is easy to use for extracting information from experimental data. It works equally in analyzing experimental data and lattice QCD simulations. In the latter case, we can study the real chemical potential regions without Taylor expansions.
Several problems that should be clarified in future are • We can study finite real chemical potential regions in lattice QCD simulation and the sign problem does not appear here. A possible obstacle is the overlap problem. In our lattice Monte Carlo simulations, the gauge configurations are produced at µ = 0 as in Eq. (13). Such configurations may not overlap enough with states of large n. The canonical partition functions Z n are given in where µ I is pure imaginary chemical potential. In other words, Z(ξ) on the unit circle has whole information 4 . On this circle, det D is real, and the Monte Carlo simulation is possible. Therefore, in addition to µ = 0, more points on the unit circle may improve the overlap problem.
• Since (det D(µ)) 2 is real positive on the unit circle in ξ plane, Z(ξ) cannot be zero. Therefore Lee-Yang zeros on the unit circle are artifact of N max 5 , or the u− and d− quark contribution, det D(µ u ) det D(µ d ) gives such a effect. Therefore we need larger N max both in experiments and lattice QCD calculations.
• It is difficult to measure experimentally the net baryon multiplicity. One possible approach is to study the difference between net-proton multiplicity and net-baryon multiplicity.
• Net strangeness and net charge multiplicity are analyzed in the same way.
• The relation between the order of the RW transition and the Lee-Yang zeros should be studied more quantitatively near the end-point. The volume dependence of Lee-Yang zeros in the vicinity of a transition point depends on the order of the transition. If the transition is the second order, Lee-Yang zeros approach the RW phase points according to the volume dependence with a critical exponent.
• The Lee-Yang zeros calculated in the lattice QCD simulation at high temperature suggest the RW transition. But there might be danger to misidentify the thermal singularity with the Lee-Yang zeros. The former appears as a branch point at ξ = e −mN /T inside the unit circle and a cut to −∞ on real negative ξ axis [34,17]. In order to avoid such misidentification, it is important to check the volume dependence.
• The Lee-Yang zeros are closing to the real positive axis if the system indicates the phase transition. Therefore it is interesting the volume dependence of the Lee-Yang zeros. Such information may be obtained by varying the atomic number of beam/target.
• Recently, based on the canonical partition functions, detailed analyses of effects of N max on the Lee-Yang zeros and the RW transition were reported, using the random matrix model and the saddle point approximation [42,43]. It is valuable to repeat calculations presented in this paper with referring these analyses.
• Both experimental and lattice QCD data include errors. It is important to check whether obtained Lee-Yang zeros are stable or not against these errors, since it requires extreme care to calculate zeros of high order polynomials. In lattice QCD, we have investigated statistical errors of Lee-Yang zeros caused by statistical error of Z n by using a bootstrap analysis. We found that the Lee-Yang zeros of high temperature QCD are stable near the unit circle on the fugacity plane [32]. If errors of P n are available, we can investigate the error of Lee-Yang zeros for experimental data in a similar manner to lattice QCD. In Ref. [44], the authors addressed the numerical error of Lee-Yang zeros using a random matrix model with finite chemical potential. They calculated zeros of the grand partition function in the complex µ and observed their behavior when they change the exact polynomial coefficients, c k toc k asc k = c k (1 + R k ǫ) . Here R k are random numbers between [−1, 1]. They found that the deviation of the zeros on the real µ axis from the exact values is approximately proportional to log ǫ, and that the structure os Lee-Yang zeros is not spoiled.
A Canonical partition functions, Z n for RHIC data In this appendix, we give Z n obtained in Sec.2. Since Z −n = Z n , we show only n ≥ 0 values. The data are normalized as Z 0 = 1 at each energy. ∆Z n are estimated errors.