Stable Yang-Lee zeros in truncated fugacity series from net-baryon number multiplicity distribution
aa r X i v : . [ h e p - ph ] N ov YITP-15-44, RIKEN-QHP-190
Stable Yang-Lee zeros in truncated fugacity series from net-baryon numbermultiplicity distribution
Kenji Morita
1, 2, 3, ∗ and Atsushi Nakamura
4, 5 Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Frankfurt Institute for Advanced Studies,Ruth-Moufang-Strasse 1, D-60438, Frankfurt am Main, Germany Institute of Theoretical Physics, University of Wroclaw, PL-50204 Wroc law, Poland Research Center for Nuclear Physics (RCNP), Osaka University, Ibaraki, Osaka, 567-0047, Japan Theoretical Research Division, Nishina Center, RIKEN, Wako, Saitama 351-0198, Japan (Dated: September 1, 2018)We investigate Yang-Lee zeros of grand partition functions as truncated fugacity polynomials ofwhich coefficients are given by the canonical partition functions Z ( T, V, N ) up to N ≤ N max . Such apartition function can be inevitably obtained from the net-baryon number multiplicity distributionin relativistic heavy ion collisions, where the number of the event beyond N max has insufficientstatistics, as well as canonical approaches in lattice QCD. We use a chiral random matrix model asa solvable model for chiral phase transition in QCD and show that the closest edge of the distributionto real chemical potential axis is stable against cutting the tail of the multiplicity distribution. Thesimilar behavior is also found in lattice QCD at finite temperature for Roberge-Weiss transition.In contrast, such a stability is found to be absent in the Skellam distribution which does not havephase transition. We compare the number of N max to obtain the stable Yang-Lee zeros with thoseof critical higher order cumulants. PACS numbers: 12.38.Gc, 12.38.Mh, 25.75.Nq, 25.75.Gz
I. INTRODUCTION
Phase transition in quantum chromodynamics (QCD)is one of the central subjects in high energy nuclearphysics both theoretically and experimentally. Firstprinciple lattice QCD (LQCD) calculations have shownthat the transition from quark-gluon plasma (QGP) tohadronic matter is of crossover type at physical quarkmasses [1], in which order parameters and thermody-namic quantities change smoothly as functions of tem-perature. At finite baryon density, one expects that thenature of the transition can change. Unfortunately, lit-tle is known about the state of matter at high baryondensity from LQCD calculations because of the difficultyin numerical simulation at finite baryon chemical poten-tial µ [2, 3]. Various approximation methods appliedso far seem to work only µ < T or a small volume orheavy quark mass region. Nevertheless, effective modelswhich implement relevant symmetries in QCD and large N c studies have shown that rich phase structure exists inhigh density [4, 5]. In particular, if there is a first orderphase transition at T = 0 and large µ , there must be acritical point (CP) at which the first order phase tran-sition line terminates and the transition becomes secondorder. Existence of CP is supported by many chiral ef-fective models [6], but its location depends on the detailof the models [7].Stimulated by these theoretical results, the first beamenergy scan program at Relativistic Heavy Ion Collider(RHIC) has been carried out in search for the CP [8, 9]. ∗ [email protected] Since lower colliding energies leaves the incident nucle-ons in the central region, one expects to explore higherbaryon density region at lower energies. There are anumber of observables which might have potential toindicate the transition from QGP to hadronic matter.Among them, event-by-event fluctuations of conservedcharges are intimately connected to critical behavior as-sociated with the phase transition [10–13]. Measure-ments of the net-proton number fluctuations as a proxyof the net-baryon one [14] and net electric charges havebeen presented for Au+Au collisions at various energiesfrom √ s NN = 7 . T and baryonic chemical potential µ .Through systematic analyses of the location of ( T, µ ) cor-responding to each colliding systems [20], one can mapexperimental measurements for property of the matteron T − µ plane. Furthermore, recent LQCD results atphysical quark masses indicate that the crossover regioncoincides with the chemical freeze-out at least for µ < T [21, 22]. One may look for remnant of the chiral critical-ity in the crossover region, originating from second orderphase transition in the vanishing quark mass [23–25].Property of the transition can be characterized by be-havior of fluctuations of conserved charges as well as anorder parameter and its fluctuations [26, 27]. In the caseof the chiral phase transition in QCD, the chiral orderparameter or quark-antiquark condensate h ¯ qq i couplesto quarks carrying the baryon number and the electriccharge. Thus, the second order chiral phase transitionin the chiral limit at finite temperature is characterizedby not only divergent fluctuation of the order parame-ter but also higher order cumulants of the net baryonnumber and the net electric charge [28]. The divergenceof the conserved charge fluctuations is governed by thecritical exponent of the specific heat α which dependson the universality class QCD belongs to. Although itis not completely determined yet [29, 30], recent simula-tions [31] indicate O (4) in the three dimensions, as con-jectured by Pisarski and Wilczek [32]. In this case, thefirst divergent cumulant appears at the sixth order. Atfinite but small quark masses, the divergence is replacedby sign change, owing to the property of the universal O (4) scaling function [24, 33]. At nonzero net baryonnumber density, the divergence at 2 n -th order cumulantsin the chiral limit appears at n -th order one. The tricrit-ical point in the chiral limit becomes the CP, where thesecond order cumulants diverge [26, 34, 35]. The chemicalfreeze-out line may locate at lower temperature than thechiral phase boundary [36] such that measured fluctua-tions might not reflect those at the phase transition [37].Nevertheless, the existence of the CP is accompanied byanomalous behavior of the cumulants such as negativethird and fourth order cumulants around the CP [38–40]and may lead to non-monotonic behavior of the higherorder cumulants as functions of √ s NN . Indeed, the mea-sured net-proton number cumulants in [16] seems to fol-low this expectation, although still inconclusive due touncertainty.The measurement of the cumulants is based on event-by-event multiplicity distribution. Once the fluctuationsare regarded as those of the grand canonical ensemble,the multiplicity distribution can be identified with un-normalized probability distribution.While the cumulants are expressed by central momentsof the probability distribution, it is convenient for theo-retical studies to compute them by differentiating thethermodynamic pressure with respect to chemical poten-tials. Recently, one of the authors (K.M.) investigatedthe probability distribution of the net baryon number inmodels with phase transitions [25, 41, 42]. It turns outthat sufficient information on the tail in the probabilitydistribution is responsible for the critical behavior of thehigher order cumulants and that the remnant of the O (4)criticality can be characterized by narrower tail than thecorresponding reference distribution. In the probabilitydistribution, such information on the phase transition isencoded in the N dependence of the canonical partitionfunction Z ( T, V, N ).Since the grand partition function is more straight-forward in relativistic quantum field theories where thenumber of particles are not definite, computations ofthe canonical partition function are not generally easy.In Ref. [43], Hasenfratz and Toussaint proposed thatthe canonical partition function, Z ( T, V, N ), is calcu-lated through the Fourier transformation of the grandcanonical partition function, Z ( T, V, µ ), evaluated at pureimaginary µ . The difficulty associated with the complexfermion determinant is replaced by the highly oscillating integral which requires extraordinary numerical precision[41, 44, 45].Nevertheless, the probability distribution gives furtherinsights into property of the system including phase tran-sitions.In Ref. [46], one of the authors (A.N.) pointed outthat one can extract the fugacity parameter λ = e µ/T atthe chemical freeze-out and construct Z ( T, V, N ) for thenet baryon number without any assumption on the prop-erty of equilibrium P ( N ). Furthermore, once Z ( T, V, N )is known, one can obtain the grand partition function Z ( T, V, µ ) as a series of fugacity. This enables us to applyYang-Lee theory for the phase transition [47, 48](For re-cent reviews, see, e.g., [49, 50]), in which zeros of the par-tition function give information on the thermodynamicproperty of the system. The zeros of the partition func-tion are distributed on a line in the complex plane of anexternal parameter and its density grows up with the sys-tem volume, then finally coalesce into the line in the ther-modynamic limit. This property leads us, in principle,to obtain the location and order of the phase transitionfrom the distribution of the zeros. Even in the absenceof the phase transition, the zeros accumulated on theedge of the distribution exhibit singular behavior. Thissingularity, known as Yang-Lee edge singularity [51, 52],can be regarded as a CP in the complex plane and givesinfluence on the thermodynamics on the real axis [53].In both experiments and the canonical approach inLQCD, Z ( T, V, N ) at large N requires such high statis-tics that obtained information is limited to some finite N ,thus one has to truncate the fugacity polynomial there inreconstructing the grand partition function (See [54, 55]for recent LQCD calculations). It is not a priori clearwhether one can obtain the correct information on thephase transition from such a truncated partition func-tion. The purpose of this paper is to clarify this point.We employ a solvable model for the chiral phase transi-tion in QCD. We present the Yang-Lee zeros in a chiralrandom matrix model, both for the exact grand partitionfunction and for the reconstructed one as a truncated fu-gacity series with the canonical partition function beingthe coefficients. We discuss effects of the truncation onthe distribution of the Yang-Lee zeros and compare itwith the spurious zeros of the Skellam partition function,originated from the truncation.In the next section, we briefly summarize the gen-eral relation among the probability distribution, partitionfunctions, and Yang-Lee zeros. A chiral random matrixmodel and its Yang-Lee zeros are presented in Sec. III.We demonstrate differences of the truncation effects onthe Yang-Lee zeros between the models with and withoutphase transition in Sec. IV. Implications for heavy ion ex-periments are discussed in Sec. V. Section VI is devotedto concluding remarks. Detailed expressions for partitionfunctions in the chiral random matrix model are given inthe Appendix. II. GENERAL FRAMEWORK
We start from experimentally measured data of net-baryon number multiplicity distribution P ( N ), where N is the net-baryon number. In real experiments one mea-sures the net-proton number ∆ N p = N p − N ¯ p . In princi-ple one can reconstruct P ( N ) from P (∆ N p ), P ( N p ) and P ( N ¯ p ) [56]. In this study we entirely assume the isospininvariance and regard P (∆ N p ) as a proxy of P ( N ). Theshape of the distribution depends on the colliding en-ergies, centrality etc. The net-baryon number can takeany value as long as it can be packed within the systemvolume. Owing to limited statistics, however, we do notobserve such states that have too large N far from itsmean value M . Thus, we define the possible minimumand maximum of N min and N max as P ( N < N min ) = 0 P ( N > N max ) = 0 . (1)In thermal equilibrium, probability distribution of thenet-baryon number in the grand canonical ensemblereads, for the fugacity factor λ = e µ/T , P ( T, V, N, µ ) = Z ( T, V, N ) λ N Z ( T, V, µ ) , (2)where Z ( T, V, µ ) is the grand partition function Z ( T, V, µ ) = Tr[ e − ( ˆ H − µ ˆ N ) /T ] (3)and Z ( T, V, N ) is the canonical partition function Z ( T, V, N ) = Tr[ e − ˆ H/T δ ˆ N,N ] . (4)Assuming the measured multiplicity distribution is theequilibrium one, one finds N dependence of P ( N ) comesfrom Z ( T, V, N ) λ N . Using the charge-conjugate symme-try Z ( T, V, − N ) = Z ( T, V, N ), one can determine λ from P ( N ) and obtain the canonical partition function [46] Z ( T, V, N ) = P ( N ) λ − N . (5)Because of limited range of N (1), the canonicalpartition function (4) can be obtained for N ∈ [max(0 , N min ) , N max ]. For the energy scan range in RHICexperiments, N min < i.e. , there are a few events inwhich more anti-protons are observed than protons, ex-cept for √ s NN = 7 . N min = 1 [16]. Notethat we need only N ≥ − N max ≤ N ≤ N max .This limitation in N also applies to theoretical ap-proaches such as model studies [25, 41, 42] and latticesimulations [54, 55]. In the former, the canonical par-tition functions have been calculated using a projectionformula Z ( T, V, N ) = 12 πi I dλ Z ( T, V, λ ) λ N +1 . (6) where integration contour C in complex λ plane can bearbitrary, but it is convenient to take the unit circle λ = e iθ . Then the formula becomes Z ( T, V, N ) = 12 π Z π − π dθ cos( N θ ) Z ( T, V, θ ) (7)where θ is related to an imaginary chemical potential − iµ/T . In Refs. [25, 41, 42], thermodynamic potentialΩ = − pV in Landau theory [41] and in chiral quark-meson model [25, 42] was used through Z ( T, V, θ ) = e − Ω( T,θ ) /T . Owing to the rapid oscillation in large N ,it turned out that the numerical integration in doubleprecision works up to P ( N ) ≃ − .In lattice QCD simulations, two approaches can pro-vide the canonical partition functions, i.e., (i) the fugac-ity expansion of the fermion determinant [57] and (ii)Hasenfratz and Toussaint method [43]. In the fugacityexpansion, we must diagonalize a matrix whose rank isproportional to the lattice spacial volume. This requireslarge computational resource and currently one cannotgo to simulations on large lattices. In the method (ii), as N increases, more accuracy is needed, and consequently N cannot go to very large.Once the canonical partition function Z ( T, V, N ) is ob-tained, one can also construct a truncated grand canoni-cal partition function as a series in λ Z tr ( T, V, λ ; N max ) = N max X N = − N max Z ( T, V, N ) λ N . (8)Owing to the truncation of the series at N = − N max and N max , this partition function is an approximationof the exact partition function which could be obtainedif one can take N max = N ∗ with N ∗ being the num-ber of net-baryons fulfilling the system volume [47]. Forlattice QCD at finite temperature with N s × N t lattice, N ∗ = 2 N s [57] . Thus, one needs to establish relations ofphysical quantities obtained from the truncated partitionfunction (8) with those from the exact partition function.As seen in the summation running from − N max , relativis-tic partition functions contain negative powers of λ . Thesuppression of high N contribution to Z cannot be real-ized by small λ . One needs to know large N behavior in Z ( T, V, N ).Similar studies on higher order cumulants of the netbaryon number have been carried out in Ref. [42], inwhich sufficiently large N max depending on the order ofthe cumulants is shown to be necessary to obtain a cor-rect value of the cumulants. In this paper, we focus onYang-Lee zeros for the baryon chemical potential.The zeros of the partition function in complex chemicalpotential plane can be obtained by solving an equation Z ( T, V, µ ) = 0 , (9) In Ref. [57], N ∗ was derived for the quark fugacity series as N ∗ q = 2 N c N s for complex µ . Owing to the negative powers of the fu-gacity, the equation is a polynomial one in λ with order2 N ∗ . For the truncated partition function Z tr , one needsto solve Z tr ( T, V, µ ; N max ) = 0 (10)of which the order of the polynomial is 2 N max .In the exact case, the roots have both real and imag-inary part and its distribution in the complex chemicalpotential plane is expected to form a line, which crossesthe real axis at the transition point in the thermodynamiclimit. The behavior of the distribution depends on thenature of the phase transition. In Ref. [58], the behaviorof the Yang-Lee zeros around the CP was studied by us-ing a chiral random matrix model. The singularity asso-ciated with the CP appears as a branch point in complex µ plane and its property is shown to be connected withthe universality. In lattice QCD, the phase transitionbetween different Z ( N c ) sector in the deconfined phase,Roberge-Weiss phase transition, has been recently ana-lyzed from a view point of Yang-Lee zeros [55]. In thiswork, we use a chiral random matrix model similar toused in Ref. [58] but with an extension to periodic prop-erty in imaginary chemical potential as it is necessary tohave integer net baryon number.The partition function is written as a polynomial in λ = e µ/T . Since a complex root λ is accompanied withits conjugate ¯ λ and the charge conjugate symmetry im-plies 1 /λ and 1 / ¯ λ are also roots, only the roots locatedin the first quadrant of the complex µ plane are inde-pendent. In practice, it is convenient to use Joukowskitransformation ω = λ + 1 /λ (= 2 cosh βµ ) and reorganizethe series in terms of ω to reduce the number of roots tosearch for. Using a property of Chebychev polynomial T n (cosh x ) = cosh( nx ), one finds λ N + 1 λ N = 2 cosh( βµN ) = 2 T N (cosh βµ )= 2 T N ( ω/ . (11)Then Eq. (8) reduces to a series expression containingonly positive powers. After expanding the Chebychevpolynomial by Eq. (A7), one finds, Z tr ( T, V, ω ; N max )= Z ( T, V, N = 0) + N max X n =1 nZ ( T, V, n ) × [ n/ X k =0 ( − k ( n − k − k !( n − k )! ω n − k . (12)This formula could be also useful to compare a relativisticsystem with nonrelativistic ones. The roots of ω space iseasily converted into those in λ and µ plane as µT = ± cosh − ω µ and λ plane. III. CHIRAL RANDOM MATRIX MODEL
In this section, we introduce a chiral random matrixmodel which is an effective model for the spontaneouschiral symmetry breaking in QCD. Since this model is an-alytically solvable in the chiral and thermodynamic limit[59] and analytic expression for the partition function infinite volume is known [58], we find this model as themost suitable one for the present purpose. An apparentshortcoming of the model for applying to the net baryonnumber probability distribution is lack of periodicity inimaginary chemical potential, which is a consequence of U (1) B symmetry. Thus, we first extend the model to ex-hibit the appropriate periodicity and the phase structurein the imaginary baryon chemical potential.In QCD, the partition function has a periodicity 2 π/N c in the imaginary quark chemical potential θ q = θ/
3, thus2 π in the baryon number. LQCD simulations have shownthat there is no phase transition in imaginary baryonchemical potential at temperatures below chiral crossovertemperature and thermodynamic quantities smoothly be-have as ∼ cos θ [60–62]. This fact combined with Eq. (7)implies that the phase transition at large baryon numberdensity is encoded in higher Fourier coefficients of thesmoothly oscillating function. A. Partition function and thermodynamics
We start with a partition function of the chiral randommatrix model with N s sites given in [59] Z RM = Z D X exp (cid:18) − N s σ Tr XX † (cid:19) det N f ( D + m ) (14)where σ denotes the variance of the random matrix X which has N s × N s dimension and D is the 2 N s × N s matrix approximating the Dirac operator. At T = µ = m = 0, σ is the only dimensionful parameter. We use it asa unit of mass in the model and put σ = 1 in expressionsbelow.The Dirac operator takes the form D = (cid:18) iX + iCiX † + iC (cid:19) . (15)The matrix C describes the effect of temperature andchemical potential. In Ref. [59], it was chosen as C k = aπT + bµiN c (16)for one half of eigenvalues and C k = − aπT + bµiN c (17)for the other half , with a and b being the dimensionless Note that µ is the chemical potential of the baryon number, thus µ/N c stands for that of the quark number. parameters.The linear ansatz for the matrix C (16)-(17) accountsfor the fact that there are the two smallest Matsubarafrequencies ± πT . This model does not have any thermaldistribution which gives the fugacity factor e µ/T thus northe periodicity in imaginary chemical potential, since itappears as a result of summation over the Matsubarafrequencies. In order to make the partition function pe-riodic, we perform a following replacement bN c µ + iπaT = πaT (cid:18) i + baπN c µT (cid:19) (18) → πaT (cid:18) i + baπN c µ T (cid:19) (19)which gives a periodicity 2 πT in µ I to the partition func-tion. Compared to the original linear ansatz, this replacedoes not change anything at µ = 0 but alters the phasestructure at µ > T .The phase structure of the model is easily evaluatedby taking N s → ∞ limit. Introducing an auxiliary N f × N f complex matrix field φ and performing the Gaussianintegration with respect to X , one obtains the partitionfunction [63] Z RM = Z D φ exp[ − N s Ω( φ )] (20)where Ω( φ ) stands for the effective potential. Then thepartition function can be determined by the minimum ofthe potential, which is evaluated at the saddle point φ of the integrand: ∂ Ω( φ ) ∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ = φ = 0 , (21)and lim N s →∞ N s ln Z RM = − min φ Ω( φ ) . (22)The saddle point φ is related to the chiral condensatethrough h ¯ ψψ i = 1 N f V ∂ ln Z RM ∂m , (23)= 1 N f V N s σ φ (24)where the four dimensional volume V corresponds to N s such that N s represents the typical number of the instan-ton (or anti-instanton) in V . For real m , one expects φ is a real matrix proportional to the unit matrix. There-fore, the saddle point can be obtained by solving (21) forthe potentialΩ( φ ) /N f = φ −
12 ln (cid:26)(cid:20) ( φ + m ) − ˜ T (cid:16) A sinh µ T + i (cid:17) (cid:21)(cid:20) ( φ + m ) − ˜ T (cid:16) A sinh µ T − i (cid:17) (cid:21)(cid:27) (25) T / T c µ /T c m=0, 2nd orderm=0, 1st orderm=0.05, 1st orderTCPCP FIG. 1. Phase diagram of the chiral random matrix modelwith periodicity in the imaginary baryonic chemical potential. where ˜ T ≡ πaT (26)and A ≡ baπN c . (27)In the chiral limit m = 0, one finds that φ = 1 at T = µ = 0 and a second order phase transition occurs at T = 1 / ( πa ) and µ = 0, where φ continuously approachesto zero. Thus, φ can be regarded as an order parameterof the chiral phase transition.The parameters in the model, σ , a and b , are deter-mined as follows. The only dimensionful parameter σ is estimated to be σ ∼ h ¯ ψψ i ∼ − at T = µ = 0. Since T c = 1 / ( πa )at µ = 0, putting T c = 160 MeV yields a = 0 .
2. Theremaining parameter b connects the model to the den-sity scale. With the linear ansatz for C (16)-(17), onefinds the first order phase transition at T = 0 and bµ/N c = 0 . b = 0 .
13, corresponding to the first order transitionpoint at µ c ≃ m = 0 . T > T and µ < µ , where T ≃ . T c and µ ≃ .
504 is the loca-tion of the tricritical point (TCP). Below T , there is thefirst order phase transition line. At finite quark mass,the second order line is replaced by smooth crossoverand TCP becomes CP with slightly decreased temper-ature and increased chemical potential, T CP = 0 . T c π - π π π φ θ m=0.05, T/T c =0.675m=0, T/T c =0.7m=0.05, T/T c =0.98m=0, T/T c =0.98 FIG. 2. Behavior of order parameter φ of a periodic chiralrandom matrix model in imaginary chemical potential θ = µ I /T at µ R = 0. and µ CP = 4 .
72, respectively. While these structures arethe same as those in Refs. [58, 59], the apparent singular-ity at T = 0 in the periodic parametrization significantlymodifies the phase boundary at low temperature. Westress that our purpose in this paper is to explore theproperty of partition function zeros rather than determi-nation of the phase structure.With the parameter set for a, b and σ , we find that this form also gives a reasonable thermodynamic quanti-ties at imaginary chemical potential. Figure 2 displaysthe behavior of the order parameter φ in the imaginarybaryonic chemical potential θ = µ I /T . One sees thatour parameterization (19) gives the correct periodicity2 π and expected temperature dependence such as largeramplitude at higher temperature below T c [64, 65]. Ow-ing to lack of a Z (3) sector such as the Polyakov loopbackground, this model does not exhibit the Roberge-Weiss phase transition [66] at high T .At finite N s , the partition function can be expressedas [58] Z RM = N s / X k ,k =0 (cid:18) N s / k (cid:19)(cid:18) N s / k (cid:19) ( N − k − k )! × F ( k + k − N s ; 1; − m N s )( − N s ˜ T ) k + k × (cid:16) i + A sinh µ T (cid:17) k (cid:16) i − A sinh µ T (cid:17) k . (28)where an irrelevant constant factor is ignored and F ( a, b ; x ) denotes the confluent hypergeometric func-tion. One may directly obtain zeros of this partitionfunction, but one needs to expand Z in a series of thefugacity λ to examine effects of tails in the probabilitydistribution function. We put the details in the AppendixA and write down only the result for the canonical par-tition function, for δ ≡ | k − k | , Z ( T, N s , N ) = N s / X k ,k =0 (cid:18) N s / k (cid:19)(cid:18) N s / k (cid:19) ( N s − k − k )! F ( k + k − N s ; 1; − m N s )( − N s ˜ T A / k + k × k + k X k =0 (cid:18) k + k k (cid:19) (cid:20) − A ) A (cid:21) k + k − k (cid:18) k k − N (cid:19) k = k δ δ X k =0 (cid:18) − A (cid:19) k ( δ + k − δ − k )!(2 k )! k + k − k X k =0 (cid:18) k + k − k k (cid:19) (cid:20) − A ) A (cid:21) k + k − k − k (cid:18) k k − N (cid:19) k = k (29) B. Phase boundary and Yang-Lee zeros
We compute the Yang-Lee zeros for the truncated par-tition function (12) with the canonical partition func-tion of the chiral random matrix model (29). Taking Note that in Refs. [58, 59] the coefficients in the temperatureand chemical potential are absorbed into T and µ . While thequalitative phase structure does not depend on the parametersin the linear ansatz, it does so when one employs the periodicparametrization (19). N max = N s in Eq. (12), one recovers the exact grandpartition function (28). The computation of the zerosrequires a special care in numerical digits as cautioned inliterature [46, 67]. We perform the calculations in 50-300digits utilizing FMLIB package [68] in fortran 90.Figure 3 shows the distribution of the Yang-Lee zero ofthe periodic chiral random matrix model in the complex ω plane, for m = 0 and at T /T c = 0 .
99. The distributionof the zeros is symmetric with respect to the horizontalaxis because the partition function is an even order poly-nomial of ω and a root has its complex conjugate. Thesolid line in Fig. 3 stands for Stokes boundary, which can -300-200-100 0 100 200 300-700 -600 -500 -400 -300 -200 -100 0 100 I m ω Re ω m=0, T/T c =0.99 ω plane N s =2060100Stokes BoundaryBranch Point FIG. 3. Yang-Lee zeros of the periodic random matrix modelin complex ω plane for m = 0 at T /T c = 0 .
99. Open symbolsstand for the zeros in different N s and solid line indicates theStokes boundary. Branch point is denoted by closed circles. be regarded as an extension of the phase boundary to acomplex chemical potential plane. In the thermodynamiclimit N s → ∞ , it satisfiesRe ∂ Ω( φ ) ∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ = φ ! > , (30)ReΩ( φ = φ , ) = ReΩ( φ = φ , ) , (31)where the first condition ensures the well-defined parti-tion function at the saddle point of integrand in Eq. (20)and the second condition denotes the continuity of thereal part of pressure at the boundary [58]. φ , and φ , stand for the two out of five solutions of the gap equation ∂ Ω /∂φ = 0 and give the minimum of ReΩ in both sidesof the boundary, respectively. The density of the zerosincreases with N s and turns into the cut which consti-tutes the Stokes boundary in the thermodynamic limit.This is clearly seen in Fig. 3. There are two branch pointslocated on the real axis. Since ω = λ +1 /λ > µ ,the one at Re ω = ω c = 3 . > µ , while the otherone, Re ω = − . µ/T = π .The Stokes boundary exhibits a closed curve, reflectingthe periodicity in imaginary µ and existence of the phaseboundary at real µ axis and Im µ/T = π .The phase structure can be more intuitively under-stood by going to complex µ plane. Figure 4-left dis-plays the distribution of the same Yang-Lee zeros as inFig. 3, but the zeros in Re µ < µ → − µ . The branching point on the hor-izontal axis indicates the second order phase transition point. The Stokes boundary extends to both direction inimaginary µ and ends up at the other branch point. Notethat the branch points at Im µ/T = − π and π are essen-tially the same because of the periodicity. We refer to[69, 70] for behavior of the order parameter in complex µ plane and related topics. The zeros distribute alongthe boundary and becomes more dense for large N s , butdistance to the real axis is not so close for these valuesof N s . The behavior of the density of the zeros is relatedto a property of the thermodynamic potential which canbe described by an analogy to electrostatics [58]. In thiscase, ReΩ can be regarded as the electrostatic potentialon the (Re µ/T, Im µ/T ) plane and the normal componentof the electric field E = −∇ (ReΩ) to the Stokes bound-ary has a discontinuity of which amount is proportionalto the density of the zeros. We confirmed that in thismodel these discontinuities at large Re µ , where the zerosare dense, are much larger than those at small µ , follow-ing the expectation. Although the density of the zerosfar from the branching point is a model-dependent fea-ture dependent on the shape of the Stokes boundary, itis governed by the universality near the branch point onthe real axis as pointed out in Ref. [58].Effects of the finite but small quark mass can be seenin the right panel of Fig. 4 where the distribution ofYang-Lee zeros for m = 0 .
05 at the same temperatureis displayed. Owing the explicit chiral symmetry break-ing, the phase transition becomes a crossover such thatthe branch point on the real axis moves to above. Asa result, there are two branch points of which are com-plex conjugate each other. The same thing occurs alsoto the branch point at Im µ/T = ± π . Here we emphasizethat these complex singularities are, albeit unphysical,indicating existence of a chiral phase transition in thechiral limit. These are also known as Yang-Lee edge sin-gularities [52]. The critical point at finite density (SeeFig. 1) is realized by coalescence of the branch pointsclose to real µ axis when temperature is decreased [53].As seen in Fig. 4-right, the Yang-Lee zeros are fairly onthe boundary line and exhibit expected behaviors. IV. YANG-LEE ZEROS FROM TRUNCATEDPARTITION FUNCTIONS
As described in Sec. II, the connection of net baryonnumber multiplicity distribution (2) with the recon-structed grand partition function (8) could potentiallyenables us to extract the Yang-Lee zeros from experi-mental data. Since the results presented in the previoussection correspond to N max = N ∗ , i.e., no informationon the exact partition function is lost, we need to evalu-ate whether one can obtain the correct distribution of theYang-Lee zeros when the fugacity expansion is truncated.Furthermore, even if one starts from a partition functionwhich does not exhibit any phase transition such as anideal Boltzmann gas, the truncation produces the zerosof partition function because it is a polynomial of order - π - π /20 π /2 π -2 0 2 4 6 8 I m µ / T Re µ /T m=0, T=0.99T c µ plane N s =2060100Stokes BoundaryBranch Point - π - π /20 π /2 π -2 0 2 4 6 8 I m µ / T Re µ /Tm=0.05, T=0.99T c µ plane N s =2060100Stokes BoundaryBranch Point FIG. 4. Yang-Lee zeros of the periodic random matrix model in complex µ plane. The left panel corresponds to the case ofFig. 3. Right panel shows the case with a finite but small quark mass, m = 0 .
05 at the same temperature. N max . In this section we investigate in detail the effects ofthe truncation on the distribution of the Yang-Lee zeros. A. Random matrix model
Figures 5 and 6 display the distribution of the Yang-Lee zeros from the truncated partition function of the pe-riodic chiral random matrix model for various N max and m = 0 . Hereafter we set N s = 60. We confirmed thefollowing results does not depend on the specific choice of N s . We plot only the first quadrant in complex µ planeaccording to the symmetry structure of the distribution.The left panel in Fig. 5 shows the case of T = T c , atwhich transition is of crossover type as seen in the branchpoint at (Re µ/T, Im µ/T ) ≃ (3 , π/ N max = 60 = N s , the zeros are located on the Stokes boundary (dashedline). Reducing N max by one, i.e., removing Z ( N = 60)from the series, one sees a drastic change in the distri-bution. The distribution of the zeros at large Re µ andIm µ splits into the two lines, but the rest of the zerosremains unchanged. Further reduction of N max substan-tially modifies the distribution such that the splittingoccurs closer to the edge closer to the real µ axis. Nev-ertheless, up to N max = 21, the edge of the distributionwhich is the closest Yang-Lee zero to the real µ axis re-mains the same. Beyond N max = 20, the distribution nolonger holds the information on the exact Yang-Lee zerosthus the apparent relation to the phase boundary is lost. Note that T c is defined for m = 0. Thus it is slightly lower thanthe chiral crossover temperature for m = 0 . The behavior with respect to changing N max does notdepend on temperature or corresponding phase transi-tion. In the right panel of Fig. 5, we plot the resultof the same analysis for T = T CP = 0 . T c where thebranch point appears on the real axis, indicating the crit-ical point. Reflecting the location of the branch point,the edge of the distribution also become closer to thereal axis compared to the crossover case. The edge isstable against decreasing N max down to N max = 19,then it starts to deviate slowly when decreased further.This is so also in the case of first order phase transition( T = 0 . T c ) depicted in Fig. 6. The branch point is hid-den in unphysical Riemann sheets [70] and the edge isvery close to the real axis.We also note that there is always a zero at Im µ/T = π when N max is odd. These zeros look special since itcorresponds to negative real axis in both complex λ and ω plane. However, this is a mathematical consequencebecause in this case the truncated partition function (12)is an odd order polynomial, thus it has at least one realroot. As seen in Figs. 5 and 6, it becomes the edge of oneof the lines bifurcating from the exact Yang-Lee zeros.These results indicate the stability of the edge doesnot depend on the detail of the phase structure, althoughthe location of the edge seems to be connected with theshape of the Stokes boundary which is model dependentthrough the µ dependence of the partition function. Inparticular, the present results are obtained by employingthe periodicity (19) in the random matrix model whichdoes not correctly take into account degrees of freedomwith baryon charges [59]. We note that this modificationcauses unphysical behavior in thermodynamics, such asnegative Z ( T, V, N ) at some small N at low T , which π /4 π /23 π /4 π
0 1 2 3 4 5 6 7 8 I m µ / T Re µ /T m=0.05, T=T c Closest Zero N max =60594836212016 Branch Point π /4 π /23 π /4 π
0 1 2 3 4 5 6 7 8 I m µ / T Re µ /T m=0.05, T=T CP Closest ZeroBranch PointN max =60594836191612
FIG. 5. Distribution of the Yang-Lee zeros from the truncated partition function of the periodic chiral random matrix modelfor m = 0 .
05. Left and right panels stand for T = T c and T = T CP , respectively. π /4 π /23 π /4 π
0 1 2 3 4 5 6 7 8 I m µ / T Re µ /Tm=0.05, T=0.6T c Closest ZeroBranch PointN max =605948362316
FIG. 6. Same as Fig. 5, but for T = 0 . T c presumably reflect the unusual curvature of the phaseboundary in Fig. 1. Therefore, we note that the shapeof the distribution itself might not be relevant for realis-tic situations. Nevertheless, below we shall see that thestability of the edge is specific to the case with a phasetransition. π /2 π π /22 π
0 0.5 1 1.5 2 2.5 3 I m µ / T Re µ /T β =1.89 (T/T c ∼ RW Transition
Lattice QCDN max =17151311
FIG. 7. Distribution of Yang-Lee zeros for a lattice QCD data[54].
B. Lattice QCD
Figure 7 displays distribution of the Yang-Lee zerosabove T c calculated in lattice QCD simulation via thecanonical method [54]. While calculations in the confine-ment phase is still numerically difficult thus we do notsee clear indications of a phase transition at low T , theRoberge-Weiss (RW) transition [66] provides us a well-0defined phase transition in high temperature quark-gluonplasma phase, though at imaginary chemical potential.In this figure, the data are calculated on 8 × β = 1 .
89 which corresponds to
T /T c ≃ .
94. A moredetailed analysis in lattice QCD with different lattice se-tups can be found in Ref. [55]. Since quark mass is heavy,the calculation is not relevant for chiral phase transition.The RW transition is regarded as a transition from one Z (3) sector to another one when single quarks can beexcited owing to deconfinement and is known to exhibita first order phase transition at Im µ q /T = ± π/
3. Interms of baryon chemical potential, the transition linesreduce to Im µ/T = ± π which is shown as a dotted line inFig. 7. A brief explanation of the Roberge-Weiss phaseboundary can be found in Appendix B. Since it is hardto compute Z ( N ) near N = N ∗ , the canonical approachin lattice QCD lacks large N contribution when one con-structs the truncated partition function (12). One seesthat in Fig. 7 the behavior of the distribution of the Yang-Lee zeros against changing N max is similar to that of therandom matrix model, despite the completely differentorigin of the phase transition. Therefore, we expect thesimilar splitting behavior of the distribution also appear-ing in Refs. [46, 55] is also due to the truncation effect.Indeed, N s and N max dependence of the Yang-Lee zeroshown in Ref. [55] agrees with the truncation effects dis-cussed here. We expect that the bifurcation of the zerostarts at large Re µ by improving the fugacity expansion,but one needs to take N max = N ∗ to completely producethe Yang-Lee zero along the transition line. In the RWtransition where the transition point at Im µ/T = ± π , theedge of the distribution is the closest zero to the imagi-nary axis. One sees that this point is also stable againstchanging N max . This fact suggests that the stability ofthe edge is not specific to the random matrix model orchiral phase transition but might be a general propertyof the distribution when partition function is truncated. C. Skellam distribution
Finally we examine a model without phase transitionin order to check whether the stability of the edge is spe-cific to phase transition or not. We employ the Skellamdistribution [71] of which probability distribution of thenet baryon number N is given by P S ( N ) = (cid:18) N B N ¯ B (cid:19) N/ I N (2 p N B N ¯ B ) e − ( N B + N ¯ B ) (32)where N B and N ¯ B denote the thermal averages of thenumbers of baryons and anti-baryons, respectively. Themean M and variance σ of the distribution are givenby M = N B − N ¯ B and σ = N B + N ¯ B , respectively.For N B = N ¯ B , The distribution becomes symmetric andthe argument of the modified Bessel function I N ( x ) isreduced to 2 N B = σ . This distribution can be derivedfrom non-interacting Boltzmann gas [18], thus the canon- ical and grand canonical partition functions read Z ( N ) = I N ( σ ) (33) Z ( λ ) = exp (cid:20) σ (cid:18) λ + 1 λ (cid:19)(cid:21) (34)where the temperature and volume dependence is en-coded in σ . Obviously the grand partition function (34)does not have any roots thus no phase transition exists.When one constructs the truncated grand partition func-tion (8) from the canonical partition function (33), how-ever, there exist complex roots. Consequently, one mightsee these spurious zeros even if the system does not haveany phase transition, when one constructs the partitionfunction through the fugacity expansion.Here we investigate such spurious zeros from the Skel-lam partition function (33) such that it has the samevariance with the random matrix model at N s = 60, T = T c and m = 0 .
05 of which distribution of the Yang-Lee zeros is displayed in Fig. 5. Since the information onthe phase transition is encoded in the tail of the prob-ability distribution P ( N ), the Skellam distribution withthe same variance serves a useful reference distribution[25]. The probability distribution of the random matrixmodel and corresponding Skellam distribution are shownin Fig. 8-left. Both distributions almost agree for small N , according to the same σ , but the deviation appearsin the tail of the distribution with tiny probability.Figure 8-right shows the distribution of zeros of thetruncated partition function for the Skellam distributionwith σ = 0 . N max ≥
70 which is similar to those in the randommatrix model, the distributions consist almost parallellines moving to large real µ direction as N max increases.This behavior reflects the fact that all the zeros go awayto infinity as N max → ∞ since the exact grand partitionfunction does not have roots. Remarkably, the edges ofthe distributions also move together with the rest of ze-ros, in contrast to the random matrix model and latticeQCD. Furthermore, one notes that the distributions for N max ≤
20 in the random matrix model, shown in Fig. 5,resemble those from the Skellam distribution. This ob-servation indicates that the stability of the edge against N max is a consequence of the existence of phase transi-tion and information on the phase transition is lost fortoo small N max . V. DISCUSSIONA. Comparison with N max for cumulants In the previous section, we have shown that the edgeof the distribution of the Yang-Lee zeros remains un-changed when the tail part of the canonical partitionfunction is missing. In practice, this property gives im-plications for necessary statistics in heavy ion studies ofthe net-baryon number fluctuations and in lattice QCD1 -180 -160 -140 -120 -100 -80 -60 -40 -20 -80 -60 -40 -20 0 20 40 60 80 P ( N ) N σ =0.365Random matrix: N s =60Skellam (Same σ ) 0 π /4 π /23 π /4 π
0 1 2 3 4 5 6 7 8 I m µ / T Re µ /TSkellam, σ =0.365N max =8070605040302010 FIG. 8. Left: P ( N ) for the random matrix model at T = T c and m = 0 .
05, and the corresponding Skellam distribution. Right:Distribution of Lee-Yang zeros for a truncated Skellam partition function.TABLE I. N max necessary for reconstructing i -th order cu-mulants and edge of the Yang-Lee zeros from Z ( N ) in therandom matrix model at T = T c and m = 0 . N s N (2)max N (4)max N (6)max N max
60 3 4 6 2180 3 5 6 26100 4 5 7 30 calculations. Since the sufficient N max to see the stableedge depends on the system volume, here we compareit with corresponding N ( i )max for i -th order cumulants.Here we consider only even order ones for net-baryonnumber at µ = 0, since we are looking at Z ( T, V, N )rather than P ( N ) which becomes asymmetric with re-spect to N at µ >
0. Thus the first central moment δN = N − h N i = N . The second, fourth and sixth ordercumulants c n ( n = 2 , ,
6) read c = h ( δN ) i (35) c = h ( δN ) i − h ( δN ) i (36) c = h ( δN ) i − h ( δN ) ih ( δN ) i + 30 h ( δN ) i (37)The property of the higher order cumulants of net-baryon number probability distribution for changing N max was studied in Ref. [42] by using a chiral quark-meson model. For sufficiently large volume, it was shownthat N ( i )max for the cumulants approximately scale with √ V .We summarize the values of each N max in Table I.The calculations are done for T = T c and m = 0 .
05 inthe random matrix model. Owing to the narrow Z ( N ), even the sixth order cumulant for N s = 100 only re-quires N (6)max = 7, i.e, − ≤ N ≤ Z ( T, V, N ), while the edge of Yang-Lee zeros demands N max = 21. The small N (6)max implies the system vol-ume is not large enough to exhibit √ V scaling regime ofthe cumulants. This can be understood from the smallvalue of σ in the random matrix model calculations. Therapid decay of P ( N ) give a rather weak dependence of N max for the higher order cumulants. For a sufficientlylarge volume, one expects that P ( N ) resembles Gaussiannear the peak, while the probability distribution (Fig. 8)has a sharp peak. Thus we cannot assess the value of N max needed in a realistic situation relevant for heavyion collisions. Moreover, the baryon number carried inthis model is not a physical one, as mentioned above. Allwe can say is that one may need much more statistic thanhigher order cumulants. B. Skellam distribution for large volume
In the fluctuation measurements at RHIC, observed P ( N ) can be well described by the Skellam distributionand deviation from the Skellam distribution exists in thetail, resulting in higher order cumulants different fromthe Skellam case. The obtained variance reaches σ ∼ √ s NN = 7 . σ = 14 .
89 and M = 14 .
41 with available bin2 π /4 π /23 π /4 π
0 0.5 1 1.5 2 2.5 3 I m µ / T Re µ /TSkellam, σ =14.89N max =8070605040302010 FIG. 9. Distribution of the zeros for the Skellam partitionfunction for σ = 14 . from N = 0 to N = 34. As mentioned in Sec. II, one canconstruct Z ( N ) from − ≤ N ≤
34 according to chargeconjugation symmetry [46]. In the Skellam distributionfor σ = 14 .
89, we find that N ( i )max = 13 ,
20 and 26 forsecond, fourth and sixth order cumulant, respectively.Note that these N ( i )max apply to the cumulants at µ = 0.The data does not have enough statistic for the sixthorder cumulant at freeze-out µ .We plot the distribution of the Yang-Lee zeros for theconstructed Skellam distribution with σ = 14 .
89. Thebasic feature is the same as the small σ case (Fig. 8).The line of zeros moves toward infinity as N max increases.One sees that some zeros below N = 30 appear on theimaginary µ axis, which corresponding to the unit circlein complex fugacity plane. This means that, for a largevolume case, the zeros can appear on the imaginary axiswhen the tail of the Z ( N ) is not provided. In the Skellamdistribution, these zeros can be obtained directly by look-ing at the truncated partition function on the imaginary µ axis, for λ = e iθ , Z trSkellam ( θ ) = 2 N max X N = − N max I N ( σ ) cos N θ, (38)which converges into Eq. (34) with oscillations giving ze-ros on imaginary µ . The data for N = 0 and N = 34 have only 1 event. VI. CONCLUDING REMARKS
In this paper, we present analyses on partition func-tion zeros which can be obtained from a truncated seriesof the fugacity expansion. By solving an extended chi-ral random matrix model which has a periodicity in theimaginary chemical potential, we compare the exact lo-cation of the Yang-Lee zeros and those obtained fromthe truncated series. We found that the edge of the dis-tribution of the zeros is insensitive to the truncation ofhigher order terms in the fugacity expansion to some de-gree. We found the similar behavior in lattice QCD athigh temperature in the context of the Roberge-Weissphase transition. This observation indicates that thosehigher order terms may have limited influences in searchfor the location of the phase boundary in lattice QCDcalculations and heavy ion experiments. Although thedistribution of zeros exist in systems without phase tran-sition, due to the truncation, the zeros closest to the real µ axis are stable against truncation if the system has aphase transition or crossover. The spurious zeros in theSkellam distribution moves toward infinity against thetruncation. Therefore, one can distinguish whether thedistribution is related to the phase transition or not bylooking at the stability of the edge of the distributionagainst the truncation.Although the information on the Stokes boundary islost in the case of too small N max , we expect that itdoes not mean that all the relevant information on thephase transition gets lost in the truncated partition func-tion. This expectation follows from the fact the sixth andhigher order cumulants at µ = 0 should be influenced bythe phase transition and the truncated series is still ableto reproduce them. It would be interesting to see howthe distributions of the zeros in the small N max cases inFigs. 5-6 are related to the remnant of the phase transi-tion.The order of the truncation in the fugacity series to ob-tain the stable edge of the Yang-Lee zeros, N max , is nev-ertheless found to be much larger than those for higherorder cumulants. We cannot make a quantitative assess-ment on realistic values for heavy ion experiments due tothe lack of connection in the model to the real world. Wehope that such an estimate becomes feasible in the nearfuture. ACKNOWLEDGMENTS
The authors would like to thank X. Luo and N. Xufor providing numerical data of the STAR collabora-tion. They would like to gratefully thank B. Frimanand K. Redlich for fruitful discussion and continuous en-couragement. They acknowledge stimulating discussionswith Ph. de Forcrand, F. Karsch, J. Knoll, V. Koch,K. Nagata and J. Wambach. This work was sup-ported by the Grants-in-Aid for Scientific Research onInnovative Areas from MEXT (No. 24105008), by the3Grants-in-Aid for Scientific Research from JSPS (No.15H03663, No.26610072), HIC for FAIR, and the Pol-ish Science Foundation (NCN), under Maestro grant2013/10/A/ST2/00106.
Appendix A: Derivation of the canonical partitionfunction in a chiral random matrix model
In the following, we derive an analytic expression for Z ( T, N s , N ) from Eq. (28). First we rewrite the µ depen-dent part in terms of the fugacity λ = e µ/T .SinceRe (cid:20)(cid:16) A sinh µ T + i (cid:17) k (cid:16) A sinh µ T − i (cid:17) k (cid:21) (A1)= (cid:16) A sinh µ T + 1 (cid:17) k + k cos[2( k − k ) φ ] (A2)where tan φ = (cid:16) A sinh µ T (cid:17) − , (A3)and imaginary part vanishes after summation over k and k , using the Chebychev polynomial T k − k (cos 2 φ ) = cos[2( k − k ) φ ] (A4)andcos 2 φ = (cid:16) A sinh µ T − (cid:17) / (cid:16) A sinh µ T + 1 (cid:17) , (A5)we have the partition function Z as Z = N s / X k ,k =0 (cid:18) N s / k (cid:19)(cid:18) N s / k (cid:19) ( N s − k − k )! × F ( k + k − N s ; 1; − m N s )( − N s ˜ T A / k + k × (cid:18) λ + 1 λ + 2(2 − A ) A (cid:19) k + k × T k − k (cid:18) λ + λ − − A + 2) /A λ + λ − − A − /A (cid:19) . (A6)Expanding the Chebychev polynomial by the followingexpression T n ( x ) = n = 0 n n X k =0 ( − k ( n + k − n − k )!(2 k )! (1 − x ) k n ≥ , (A7)and using binomial expansion in the third line of (A6), wecan express Z in terms of λ + λ − . We obtain Eq. (29) byapplying the projection (6). Note that maximum powerof λ is given by N s . Appendix B: Roberge-Weiss transition as a thermalcut
In this appendix, we give a brief explanation of the cutarising from the Fermi distribution function and apply itto the Roberge-Weiss transition in QCD.
1. Thermal cut in free Fermi gas
The thermodynamic potential of the free Fermi gas isgiven byΩ f ∼ − Z d p (2 π ) ln[1 + e − β ( E p − µ ) ] + ( µ → − µ ) (B1)where E p = p p + m . When the chemical potential µ has an imaginary part, µ I = θT , the imagary part givesa phase in front of the Boltzmann factor:1 + e − β ( E p − µ ) = 1 + e iθ e − β ( E p − µ R ) (B2)where µ = µ R + iµ I . Therefore, for θ = ± π , the phasegives − θ = ± π and m ≤ µ R < ∞ . The anti-particle term also gives the cut symmetric with respectto the imaginary axis. In Ref. [73], it is pointed out thatthe branch point singularity limits the convergence ra-dius when one tries to analytically continue the resultsin the imaginary chemical potential to the real one. Sincethis cut originates from the Fermi distribution, the sameanalytic structure appears in chiral models with fermions[69].
2. Roberge-Weiss transition
In QCD at high temperature, quarks are deconfinedand have a light mass owing to chiral restoration. Sincethe deconfinement can be expressed as a breaking of Z ( N c ) symmetry, it is useful to resort to chiral effec-tive models with the Polyakov loop background [74–76]which successfully describe the Roberge-Weiss transition[64, 65, 77]. Then, the relevant leading single quark con-tribution to the thermodynamic potential readsΩ q ¯ q ∼ − Z d p (2 π ) ln[1 + 3Φ e − β ( E p − µ q ) ]+( µ q → − µ q , Φ → ¯Φ) (B3)where µ q = µ/ ϕ . Onemay express Φ = | Φ | e iϕ . Then the thermodynamic con-tribution becomes1 + 3Φ e − β ( E p − µ q ) = 1 + 3 | Φ | e i ( ϕ − θ q ) e − β ( E q − µ q,R ) . (B4)4The phase of the Polyakov loop ϕ varies as a func-tion of the imaginary quark chemical potential θ q . TheRoberge-Weiss transition at θ q = π/ Z (3) sector with ϕ = 0 to an-other one ϕ = − π/ ϕ and θ q gives the prefactor − | Φ | ∼ µ q,R = 0. This feature gives the cut drawn asRW transition line in Fig. 7. A derivation based on theGaussian P ( N ) can be found in Ref. [55]. In the confine-ment phase where | Φ | ∼
0, this term is suppressed andthe thermal cut from the quark does not appear. [1] Y. Aoki, G. Endr¨odi, Z. Fodor, S. D. Katz, and K. K.Szab´o, “The order of the quantum chromodynamicstransition predicted by the standard model of particlephysics,” Nature (London) , 675 (2006).[2] S. Muroya, A. Nakamura, C. Nonaka, and T. Takaishi,“Lattice QCD at finite density -An introductory review-,” Prog. Theor. Phys. , 615 (2003).[3] P. de Forcrand, “Simulating QCD at finite density,” Proc.Sci.
LAT2009 , 010 (2009).[4] K. Fukushima and T. Hatsuda, “The phase diagram ofdense QCD,” Rep. Prog. Phys. , 014001 (2011).[5] K. Fukushima and C. Sasaki, “The phase diagram of nu-clear and quark matter at high baryon density,” Prog.Part. Nucl. Phys. , 99 (2013).[6] M. Asakawa and K. Yazaki, “Chiral restoration at finitedensity and temperature,” Nucl. Phys. A504 , 668 (1989).[7] M. Stephanov, “QCD phase diagram and the criticalpoint,” Prog. Theor. Phys. Suppl. , 139 (2004).[8] N. Xu (STAR Collaboration), “An overview of STARexperimental results,” Nucl. Phys.
A931 , 1 (2014).[9] R. A. Soltz (PHENIX Collaboration), “PHENIX beanenergy scan results,” Nucl. Phys.
A931 , 780 (2014).[10] M. Stephanov, K. Rajagopal, and E. Shuryak, “Signa-tures of the tricritical point in QCD,” Phys. Rev. Lett. , 4816 (1998).[11] M. Stephanov, K. Rajagopal, and E. Shuryak, “Event-by-event fluctuations in heavy ion collisions and the QCDcritical point,” Phys. Rev. D , 114028 (1999).[12] M. Asakawa, U. W. Heinz, and B. M¨uller, “Fluctuationprobes of quark deconfinement,” Phys. Rev. Lett. ,2072 (2000).[13] S. Jeon and V. Koch, “Charged particle ratio fluctuationas a signal for QGP,” Phys. Rev. Lett. , 2076 (2000).[14] Y. Hatta and M. A. Stephanov, “Proton-number fluctu-ation as a signal of the QCD critical end point,” Phys.Rev. Lett. , 102003 (2003).[15] M. M. Aggarwal et al. (STAR Collaboration), “Highermoments of net proton multiplicity distributions atRHIC,” Phys. Rev. Lett. , 022302 (2010).[16] L. Adamczyk et al. (STAR Collaboration), “Energy de-pendence of moments of net-proton multiplicity distribu-tions at RHIC,” Phys. Rev. Lett. , 032302 (2014).[17] L. Adamczyk et al. (STAR Collaboration), “Beam en-ergy dependence of moments of the net-charge multiplic-ity distributions in Au+Au collisions at RHIC,” Phys.Rev. Lett. , 092301 (2014).[18] P. Braun-Munzinger, K. Redlich, and J. Stachel, “Par-ticle production in heavy ion collisions,” in Quark-GluonPlasma 3 , edited by R. C. Hwa and X. N. Wang (WorldScientific, 2004) p. 491.[19] A. Andronic, P. Braun-Munzinger, and J. Stachel,“Hadron production in central nucleus-nucleus collisions at chemical freeze-out,” Nucl. Phys.
A772 , 167 (2006).[20] J. Cleymans, H. Oeschler, K. Redlich, and S. Wheaton,“Comparison of chemical freeze-out criteria in heavy-ioncollisions,” Phys. Rev. C , 034905 (2006).[21] C. R. Allton, M. D¨oring, S. Ejiri, S. J. Hands, O. Kacz-marek, F. Karsch, E. Laermann, and K. Redlich, “Ther-modynamics of two flavor QCD to sixth order in quarkchemical potential,” Phys. Rev. D , 054508 (2005).[22] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao,S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner,and W. Unger, “Phase boundary for the chiral transi-tion in (2+1)-flavor QCD at small values of the chemicalpotential,” Phys. Rev. D , 014504 (2011).[23] F. Karsch and K. Redlich, “Probing freeze-out conditionin heavy ion collisions with moments of charge fluctua-tions,” Phys. Lett. B , 136 (2011).[24] B. Friman, F. Karsch, K. Redlich, and V. Skokov, “Fluc-tuations as probe of the QCD phase transition and freeze-out in heavy ion collisions at LHC and RHIC,” Eur. Phys.J. C , 1694 (2011).[25] K. Morita, B. Friman, and K. Redlich, “Criticality thenet-baryon number probability distribution at finite den-sity,” Phys. Lett. B , 178 (2015), arXiv:1402.5982v1.[26] Y. Hatta and Y. Ikeda, “Universality, the QCD criticaland tricritical point, and the quark number susceptibil-ity,” Phys. Rev. D , 014028 (2003).[27] V. Koch, “Hadronic fluctuations and correlations,”arXiv:0810.2520.[28] M. A. Stephanov, “Non-gaussian fluctuations near theQCD critical point,” Phys. Rev. Lett. , 032301(2009).[29] S. Aoki, H. Fukaya, and Y. Taniguchi, “Chiral symmetryrestoration, the eigenvalue density of the Dirac operator,and the axial U (1) anomaly at finite temperature,” Phys.Rev. D , 114512 (2012).[30] T. Bhattacharya, M. I. Buchoff, N. H. Christ, H. T. Ding,R. Gupta, C. Jung, F. Karsch, Z. Lin, R. D. Mawhinney,G. McGlynn, S. Mukherjee, D. Murphy, P. Petreczky,D. Renfrew, C. Schroeder, R. A. Soltz, P. M. Vranas, andH. Yin, “QCD phase transition with chiral quarks andphysical quark masses,” Phys. Rev. Lett. , 082001(2014).[31] S. Ejiri, F. Karsch, E. Laermann, C. Miao, S. Mukherjee,P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger,“Magnetic equation of state in (2+1)-flavor QCD,” Phys.Rev. D , 094505 (2009).[32] R. D. Pisarski and F. Wilczek, “Remarks on the chiraltransition in chromodynamics,” Phys. Rev. D , 338(1984).[33] J. Engels and F. Karsch, “Scaling functions of the freeenergy density and its derivatives for the 3 d O (4) model,”Phys. Rev. D , 094506 (2012). [34] C. Sasaki, B. Friman, and K. Redlich, “Quark numberfluctuations in a chiral model at finite baryon chemicalpotential,” Phys. Rev. D , 054026 (2007).[35] C. Sasaki, B. Friman, and K. Redlich, “Susceptibilitiesand the phase structure of a chiral model with Polyakovloops,” Phys. Rev. D , 074013 (2007).[36] S. Floerchinger and C. Wetterich, “Chemical freeze-outin heavy ion collisions at large baryon densities,” Nucl.Phys. A890-891 , 11 (2012).[37] K. Fukushima, “Hadron resonance gas and mean-field nu-clear matter for baryon number fluctuations,” Phys. Rev.C , 044910 (2015).[38] M. Asakawa, S. Ejiri, and M. Kitazawa, “Third momentsof conserved charges as probes of QCD phase structure,”Phys. Rev. Lett. , 262301 (2009).[39] M. A. Stephanov, “Sign of kurtosis near the QCD criticalpoint,” Phys. Rev. Lett. , 052301 (2011).[40] V. Skokov, B. Friman, and K. Redlich, “Quark numberfluctuations in the Polyakov loop-extended quark-mesonmodel at finite baryon density,” Phys. Rev. C , 054904(2011).[41] K. Morita, V. Skokov, B. Friman, and K. Redlich, “Netbaryon number probability distribution near chiral phasetransition,” Eur. Phys. J. C , 2706 (2014).[42] K. Morita, B. Friman, K. Redlich, and V. Skokov, “Netquark number probability distribution near the chiralcrossover transition,” Phys. Rev. C , 034903 (2013).[43] A. Hasenfratz and D. Toussaint, “Canonical ensemblesand nonzero density quantum chromodynamics,” Nucl.Phys. B371 , 539 (1992).[44] R. Fukuda, A. Nakamura, and S. Oka, “Canonical ap-proach to finite density QCD with multiple precisioncomputation,” (2015), arXiv:1504.06351 [hep-lat].[45] A. Nakamura, S. Oka, and Y. Taniguchi, “Canonicalapproach to finite density QCD with winding numberexpansion,” (2015), arXiv:1504.04096 [hep-lat].[46] A. Nakamura and K. Nagata, “Probing QCD phasestructure by baryon multiplicity distribution,” (2013),arXiv:1305.0760.[47] C. N. Yang and T. D. Lee, “Statistical theory of equationsof state and phase transitions. I. theory of condensation,”Phys. Rev. , 404 (1952).[48] T. D. Lee and C. N. Yang, “Statistical theory of equationsof state and phase transitions. II. lattice gas and isingmodel,” Phys. Rev. , 410 (1952).[49] R. A. Blythe and M. R. Evans, “The Lee-Yang theoryof equilibrium and nonequilibrium phase transitions,”Brazilian Journal of Physics , 464 (2003).[50] I. Bena, M. Droz, and A. Lipowski, “Statistical mechan-ics of equilibrium and nonequilibrium phase transitions:the Yang-Lee formalism,” Int. J. Mod. Phys. B , 4269(2005).[51] P. J. Kortmann and R. B. Griffiths, “Density of zeros onthe Lee-Yang circle for two ising ferromagnets,” Phys.Rev. Lett. , 1439 (1971).[52] M. E. Fisher, “Yang-Lee edge singularity and φ fieldtheory,” Phys. Rev. Lett. , 1610 (1978).[53] S. Ejiri, Shinno Y, and H. Yoneyama, “Complex sin-gularities around QCD critical point at finite densities,”Prog. Theor. Exp. Phys. , 083B02 (2014).[54] K. Nagata, S. Motoki, Y. Nakagawa, A. Nakamura, andT. Saito (XQCD-J Collaboration), “Towards extremelydense matter on the lattice,” Prog. Theor. Exp. Phys. , 01A103 (2012). [55] K. Nagata, K. Kashiwa, A. Nakamura, and S. M.Nishigaki, “Lee-Yang zero distribution of high tempera-ture QCD and Roberge-Weiss phase transition,” (2014),arXiv:1410.0783 [hep-lat].[56] M. Kitazawa and M. Asakawa, “Revealing baryon num-ber fluctuations from proton number fluctuations in rel-ativistic heavy ion collisions,” Phys. Rev. C , 021901(2012).[57] K. Nagata and A. Nakamura, “EoS of finite density QCDwith Wilson fermisions by multi-parameter reweightingand Taylor expansion,” JHEP , 092 (2012).[58] M. Stephanov, “QCD critical point and complex chemicalpotential singularities,” Phys. Rev. D , 094508 (2006).[59] M. A. Halasz, A. D. Jackson, R. E. Shrock, M. A.Stephanov, and J. J. M. Verbaarschot, “Phase diagramof QCD,” Phys. Rev. D , 096007 (1998).[60] P. de Forcrand and O. Philipsen, “The QCD phase tran-sition for small densities from imaginary chemical poten-tial,” Nucl. Phys. B642 , 290 (2002).[61] P. de Forcrand and O. Philipsen, “The chiral critical lineof n f = 2 + 1 QCD at zero and non-zero baryon density,”JHEP , 077 (2007).[62] M. D’Elia and M. P. Lombardo, “Finite density QCDvia an imaginary chemical potential,” Phys. Rev. D ,014505 (2003).[63] J. B. Kogut and M. A. Stephanov, The Phases of Quan-tum Chromodynamics: From Confinenment to ExtremeEnvironments (Cambridge University Press, 2004).[64] Y. Sakai, K. Kashiwa, H. Kouno, and M. Yahiro,“Polyakov loop extended Nambu-Jona-Lasinio modelwith imaginary chemical potential,” Phys. Rev. D ,051901(R) (2008).[65] K. Morita, V. Skokov, B. Friman, and K. Redlich,“Probing deconfinement in a chiral effective model withPolyakov loop at imaginary chemical potential,” Phys.Rev. D , 076009 (2011), 1107.2273.[66] A. Roberge and N. Weiss, “Gauge theories with imagi-nary chemical potential and the phases of QCD,” Nucl.Phys. B275 , 734 (1986).[67] M. A. Halasz, A. D. Jackson, and J. J. M. Verbaarschot,“Yang-Lee zeros of a random matrix model for QCD atfinite density,” Phys. Lett. B , 293 (1997).[68] D. M. Smith, http://myweb.lmu.edu/dmsmith/FMLIB.html .[69] V. Skokov, K. Morita, and B. Friman, “Mapping outthe phase diagram of strongly interacting matter,” Phys.Rev. D , 071502(R) (2011).[70] B. Friman, “Phase transitions at finite density,” Acta.Phys. Pol. B (Proc. Suppl.) , 707 (2012).[71] J. G. Skellam, “The frequency distribution of the differ-ence between two poisson variates belonging to differentpopulations,” Journal of the Royal Statistical Society Se-ries A , 3 (1946).[72] P. Braun-Munzinger, B. Friman, F. Karsch, K. Redlich,and V. Skokov, “Net-proton probability distribution inheavy ion collisions,” Phys. Rev. C , 064911 (2011).[73] F. Karbstein and M. Thies, “How to get from imagi-nary to real chemical potential,” Phys. Rev. D , 025003(2007).[74] K. Fukushima, “Chiral effective model with the Polyakovloop,” Phys. Lett. B , 277 (2004).[75] C. Ratti, M. A. Thaler, and W. Weise, “Phases of QCD:Lattice thermodynamics and a field theoretical model,”Phys. Rev. D , 014019 (2006).[76] B. J. Schaefer, M. Wagner, and J. Wambach, “Thermo- dynamics of (2+1)-flavor QCD: confronting models withlattice studies,” Phys. Rev. D , 074013 (2010).[77] K. Morita, V. Skokov, B. Friman, and K. Redlich, “Roleof mesonic fluctuations in the Polyakov loop extended quark-meson model at imaginary chemical potential,”Phys. Rev. D84