Net baryon number probability distribution near the chiral phase transition

We discuss properties of the net baryon number probability distribution near the chiral phase transition to explore the effect of critical fluctuations. Our studies are performed within Landau theory, where the coefficients of the polynomial potential are parametrized, so as to reproduce the mean-field (MF), the $Z(2)$ and $O(4)$ scaling behaviors of the cumulants of the net baryon number. We show, that in the critical region, the structure of the probability distribution changes, dependently on values of the critical exponents. In the MF approach, as well as in the $Z(2)$ universality class, the contribution of the singular part of the thermodynamic potential tends to broaden the distribution. By contrast, in the model with $O(4)$ scaling, the contribution of the singular part results in a narrower net baryon number probability distribution with a wide tail.


Introduction
Fluctuations of conserved charges reflect the critical properties of a system. In particular, in a strongly interacting medium, fluctuations of the net baryon number and of the electric charge are valuable probes of the QCD phase transition [1,2,3,4,5,6,7]. Such fluctuations may provide a signature for the conjectured chiral critical point [8,9,10], as well as for the residual criticality of the underlying O (4) transition [4,5,11,12,13,14,15], expected in QCD at small densities in the limit of massless u and d quarks [16,17].
Since fluctuations of conserved charges are experimentally accessible through measurements of their multiplicity distributions, they have been suggested as probes of the proximity of the chiral crossover to the freeze-out line in heavy-ion collisions [5,11,18]. A particular role has been a e-mail: morita@fias.uni-frankfurt.de attributed to higher order cumulants, which exhibit an enhanced sensitivity to criticality with increasing order [4,11,13,14]. In particular, the sixth-and higher-order cumulants of the net baryon number and of the electric charge vary rapidly in the transition region already at vanishing baryon chemical potential. In the chiral limit, i.e, for vanishing pion mass, O(4) scaling implies that e.g. the six order cumulant diverges to +∞ or −∞ when the critical temperature is approached from below or above, respectively. For a small non-zero pion mass, the divergence is replaced by a rapid temperature dependence, resulting in a change of the sign of the cumulants close to the chiral crossover transition [5,11,14,18].
Consequently, the data on charge fluctuations in highenergy heavy ion collisions may provide experimental evidence for the chiral transition of QCD [5,19]. First measurements of the multiplicity distribution and corresponding cumulants of the net baryon number, more precisely the net proton number, have been obtained in heavy-ion collisions at RHIC by the STAR Collaboration [20]. The data are, at least on a qualitative level, consistent with the residual criticality, owing to the underlying chiral phase transition. The final conclusion, however, requires additional studies of non-critical effects such as e.g. volume fluctuations [21], acceptance corrections [22], or an exact baryon number conservation [23].
In statistical physics, the n-th order cumulant c n (T, µ) of a conserved charge N is obtained by differentiating the thermodynamic pressure p(T, µ) with respect to the corre-sponding chemical potential µ, The pressure is related to the grand canonical partition function p = (T /V ) ln Z . The cumulants can also be expressed as polynomials in the central moments (δ N) n , where δ N = N − N . The n-th moment N n is linked to the net charge probability distribution Cumulants of the net charge (1) can exhibit singularities near the phase transition. Thus, due to criticality, the corresponding probability distribution (2) should have the characteristic shape governed by the universal properties of this transition.
The main objective of this paper is to explore the qualitative features of the probability distribution of the net baryon number near the chiral phase transition and their relation to the cumulants. Such studies are not only of theoretical, but also of phenomenological interest. In heavy ion collisions the probability distribution of conserved charges and corresponding cumulants are measured to verify the chiral phase transition or its remnant. In the absence of criticality, different moments of conserved charges, as well as particle multiplicities, are consistent with predictions of the hadron resonance gas model. The probability distribution of the net baryon number is then governed by the Skellam distribution [11]. Thus, it is desirable to verify how the critical fluctuations lead to deviations from the Skellam distribution.
In this paper we explore the influence of the chiral phase transition on the probability distribution in the framework of the Landau theory of phase transitions. We construct the thermodynamic potential in such a way, that its regular part results in the Skellam distribution of the net charge. The singular part is modeled within the Landau theory, where the coefficients of the polynomial potential are parametrized, so as to reproduce the mean-field, O(4) and Z(2) scaling behaviors of the cumulants of the net charge. We quantify the modification of the Skellam distribution in the presence of the phase transition. We show that the structure of the probability distribution changes, dependently on the values of the critical exponents.
In this article, we consider not only expected O(4) universality class, for simulating chiral QCD, but also Z(2), because this possibility has not been reliably ruled out yet [33]. Comparison with the Z(2) universality class will also help us to demonstrate significant properties of the O(4) criticality on the level of the net charge probability distribution. We consider a system in the chiral limit in which the critical behavior is most prominent, and an analytical treatment of the mean field model is feasible. The effects of an explicit breaking of the chiral symmetry on the probability distribution in the O(4) universality class were recently discussed in Ref. [34]. The paper is organized as follows: in the next Section we introduce the model for the phase transition. In Section 3, we calculate the probability distribution and corresponding cumulants near phase transition. In Section 4, we present our summary and conclusions.

Fluctuations in the Landau theory
To explore the influence of the second-order phase transition on the probability distribution of conserved charges, we consider a model based on the Landau theory of phase transitions. At finite temperature T and chemical potential µ, the effective potential density divided by T 4 ,ω(T, µ; σ ) = ω/T 4 , is introduced as a polynomial in the order parameter σ , aŝ whereω bg is a non-singular background contribution parameterized aŝ The potential (3) exhibits a second order phase transition located along the critical line T = T c (µ), determined by the condition a(T, µ) = 0. For d = π 4 /30, the characteristic energy scale of the chiral phase transition in QCD is reproduced. At µ = 0 we set the transition temperature to T c = 0.15 GeV. The precise choice of parameters does not influence our conclusions.
For a positive a(T, µ), the minimum of the Landau potential (3) is located at σ = 0, while for a(T, µ) < 0 the order parameter is non-zero, σ = ± −a(T, µ) and the symmetry is spontaneously broken. The parametrization of a(T, µ) is chosen, so that the broken symmetry phase is located below the critical line T = T c (µ). Consequently, For the coefficient a(T, µ) we employ the following parameterization For µ/T ≪ 1, Eq. (7) reduces to a(T, µ) ≃ A(T − T c ) + Bµ 2 ≡ t µ , with A > 0 and B > 0. The scaling variable, t µ , is frequently used in the literature [4,15].
The cosh(µ/T ) terms in Eqs. (4) and (7) account for the periodicity of the thermodynamic potential in an imaginary chemical potential, ω(T, iθ T ) with θ = µ I /T , which is a consequence of the U B (1) symmetry. We note, that this periodicity is identical to that of QCD thermodynamic potential [35], if µ corresponds to the baryon chemical potential. However, since we do not consider confinement, the µ is identified as the chemical potential of the elementary fermion. If one considers quarks, the thermodynamic potential has a periodicity 2π/N c in an imaginary µ. It also exhibits a more complicated temperature dependent structure than that used in Eq. (7), when coupled to gauge fields [36,37,38,39].
In the presence of critical fluctuations, the singular part of the potential scales as ω 1 ∼ |t µ | 2−α near T c . Motivated by the O(4) scaling function [14,40], we extend the mean field potential (6) by introducing the critical exponent α, The total thermodynamic potential is then given bŷ The exponent α in Eq. (8) is introduced to parametrize the scaling properties of the thermodynamic potential in the critical region. The sgn(α) is introduced to reproduce the scaling function also for α < 0. Thus, by tuning α, we can switch between O(4), Z(2) and mean-field scaling 1 . For α = 0, the singular part of ω exhibits mean-field scaling.
On the other hand, for α ≃ −0.21 (α ≃ 0.11), the Landau potential emulates the critical behavior of the O(4) (Z(2)) spin system in 3-dimensions 2 . Figure 1 shows the critical line T c (µ) of the second order transition obtained in the Landau theory, using the parametrization (7). The behavior of the first four cumulants along the line of constant temperature T /T c = 0.98 for different µ is illustrated in Fig. 2. Their critical properties can be quantified from the singular potential (8).
Close to the critical point, |T − T c |/T c ≪ 1, the singular part of the potential scales as ω 1 ∼ |t µ | 2−α . Consequently, at µ = 0, only even order cumulants of the net baryon number are finite and their singular part scales as Thus, for the O(4) exponent α ≃ −0.21, the first divergent cumulant is of the sixth order [5,14], while for the Z(2) it is of the fourth order. For µ = 0, the leading singularities of cumulants imply, that for O(4) (Z(2)) only cumulants with n ≥ 3 (n ≥ 2) diverge at the critical point.
With the mean-field critical exponent α = 0, the singular part of the cumulants remains finite near the critical point. Indeed, differentiating Eq. (6), one finds, that for T < T c , Thus, in the mean-field approach, the cumulants are finite and in general, discontinuous at the critical point. The singular behavior of the cumulants shown in Fig. 2 and their dependence on the value of the critical exponent, must be also reflected in the corresponding probability distributions P(N) of the net baryon number. In the following, we compute P(N) within the mean-field, the O(4) and the Z(2) parametrization of the thermodynamic potential, to explore the characteristic features of P(N) which are responsible for the critical behavior of the cumulants seen in Fig. 2.

The probability distribution of the net baryon number
Consider a grand canonical thermodynamic system of charged particles q and anti-particlesq at volume V , temperature T and at chemical potential µ. The latter is related to the conserved net charge N = N q − Nq.
The normalized probability distribution P(N) to find N net charges in volume V , is given in terms of the canonical Z(T,V, N) and the grand canonical Z (T,V, µ) partition function [41,42,43], where β = 1/T . The partition functions are related by the fugacity expansion, where λ = e β µ is the fugacity parameter. Thus, Z(T,V, N) is just the N-th order coefficient in the Laurent expansion of Z (T,V, µ) around λ = 0. Consequently, Z(T,V, N) can be obtained from where the integration contour C must lie in an annulus enclosing the origin in the complex λ -plane. Inside the annulus the integrand Z (T,V, µ) must be analytic (see Fig. 3).
Taking C on the unit circle, λ = e iθ , one finds the wellknown result [41,44], We note, that Eq. (16) holds only when Z (T,V, µ) is analytic on the unit circle. In general, the partition function for finite systems 3 , e.g. in lattice QCD simulations in a finite volume [46,47,48,49,50,51], is an analytic function of µ. Hence, in this case, the integral (15) is well defined and independent of the radius of the contour. However, for the partition function Z = e −βV ω , where ω =ωT 4 in Eqs. (5), (6) and (9) is approximated by the thermodynamic potential density in the thermodynamic limit, while V is kept to be finite, there are singularities in the complex µ plane, which must be properly accounted for. This approximate treatment is justified for |N| ≪ N * , where N * is a characteristic of the partition function of a finite system in the Yang-Lee theory of phase transitions 4 .
In general, to determine N * for a given system, a microscopic computation of the grand canonical partition function in a finite volume is needed. However, such a calculation cannot be carried out within the Landau model considered here. Consequently, in the following, we assume that N * is larger than the maximal N, needed in the numerical evaluation of the probability distribution P(N). We shall discuss the properties of the canonical partition function under this assumption.
In the chiral limit and for T < T c , a critical point exists on the positive real fugacity axis, λ = λ c > 1, [53]. Charge conjugation symmetry of the partition function, µ → −µ, implies, that there is a corresponding critical point at λ = 1/λ c . This symmetry is respected also by the thermodynamic potential through (4) and (7).
In the mean-field case β = 1 2 , thus the singularity is of the square root type. In addition, sinceω sing ∼ |a(T, µ)| 2 θ (−a(T, µ)) and a ∼ t µ , the singularity of the thermodynamic potential at the critical point is a discontinuity in the higher derivatives, i.e. in the cumulants c n for n ≥ 2, as seen in Fig. 2.
In the O(4) and Z(2) case, the critical points are branch points also of the thermodynamic potential,ω sing ∼ |t µ | 2−α . Consequently, the thermodynamic potential has cuts originating at the branch points. A convenient choice is to place the cuts on the real axis, between the critical points µ = µ c and µ = ∞, as well as between µ = −µ c and µ = −∞. The corresponding cuts in the fugacity are located between λ = λ c > 1 and λ = ∞ and between λ = 1/λ c and λ = 0, respectively [53]. As a result, in the thermodynamic limit, the thermodynamic potentials of different phases (in the present case, the ω 0 and ω 1 ) correspond to different Riemann surfaces connected through cuts [54].
If the grand canonical partition function Z exhibits branch singularities, then the coefficients of the Laurent expansion Z (T,V, N), depend on the integration contour in Eq. (15). In the broken-symmetry phase, the annulus is singularity free for ρ 1 = λ c > 1 and ρ 2 = 1/λ c < 1 (cf. Fig. 3). Thus, the coefficients of the Laurent expansion corresponding to the broken-symmetry phase, are obtained from (16) with Z = e −βV ω 1 . On the other hand, the partition function Z = e −βV ω 0 of the symmetric phase is singularity free outside the annulus of the broken symmetry phase. Hence, the corresponding Laurent coefficients are obtained by integrating (15) along e.g. a circular contour with radius ρ < 1/λ c or ρ > λ c .
Thus, in general, for a given N we obtain two competing partition functions. Obviously, the partition function of a thermodynamic system must be unique. The reason for this ambiguity is that we approximate ω in Z by the grand canonical thermodynamic potential density in the thermodynamic limit, while the system is kept in a finite volume V . If V in the thermodynamic potential density is kept finite when the canonical partition function is calculated, the uniqueness of the canonical partition function is restored, even after taking the thermodynamic limit. This is because, the solution corresponding to a larger value of ω is suppressed.
In the following, we consider only the canonical partition function obtained from the grand canonical one, using Eq. (15) for µ < µ c , i.e. for the broken symmetry phase. The relation between the probability distribution P(N) and the canonical partition function (13) implies, that the structure of P(N) is entirely governed by the properties of Z(T,V, N).

Probability distribution of the nonsingular potential
The probability distribution corresponding to the nonsingular part of the thermodynamic potential density ω 0 (4) can be computed analytically. Indeed, using the generating function of the modified Bessel function I n (x), one can directly expand the grand canonical partition function Z = e −βV ω 0 , without passing through the integral representation (15), and obtain [41,55] The probability distribution from the nonsingular part of the Landau potential (5) is then obtained from Eq. (13) as the Skellam function, where Ω 0 = V ω 0 = −2dV T 4 cosh(µ/T ) is the corresponding thermodynamic potential.

Probability distribution of the mean-field potential
For the Landau potential (6) with the mean-field exponent α = 0, the canonical partition function Z(T,V, N) can also be obtained analytically, using the same procedure as was discussed above. The partition function for the singular case reads By expanding the argument of the exponential with t = 1 − T /T c and using Eqs. (17) and (16) one finds The probability distribution P MF (N) can be then obtained from Eqs. (20), (22) and (13), as Figure 4 shows the resulting probability distributions calculated at µ = 0 for V = 30 fm 3 . The mean-field probability distribution P MF (N) is broader than P NS (N) for |N| > 30. This feature is particularly evident in the ratio P MF (N)/P NS (N), shown in Fig. 5. Thus, for the Landau potential with α = 0, the criticality is expressed in the probability distribution as broadening of the regular Skellam function. 5

Probability distribution of the O(4) and Z(2) potentials
With the mean-field exponent α = 0, the net charge probability distributions were calculated exactly. However, such an exact derivation cannot be applied for the singular part with α = 0. In this case, one can only use the numerical or some approximate analytical methods to obtain the canonical partition function and the corresponding probability distribution P α (N). The most direct approach to get Z(T,V, N), is the numerical integration of Eq. (16). There are, however, limitations of this method, in particular at large N, due to the oscillatory structure of the integrand. We adopt the qawo subroutine in QUADPACK package library for computing oscillatory integrals and check convergence within 1% accuracy.
To test the numerical method, we first compare the probability distribution P MF (N) obtained analytically in Eq. (23) with the numerical integration of Eq. (16). The numerical integration reproduces the analytical result up to |N| ≃ 60, at which P(N) ∼ 10 −12 , but fails for larger N due to the round-off error. Thus, to calculate P α (N) for the model with nontrivial α, one can trust the numerical methods only up to some finite value of |N|, where P(N) ∼ 10 −12 .
The numerical results for P α (N) are shown in Fig. 4. In Fig. 5, the same probability distributions are normalized to the non-singular distribution P NS (N), given in Eq. (19). For α = −0.21 the ratio P α (N)/P NS (N) differs considerably from that of the mean-field model. The ratio P MF (N)/P NS (N) increases with |N| for |N| > 30 and exhibits oscillations about unity for |N|. While the ratio P O (4) (N)/P NS (N) also oscillates for small |N|, albeit with the opposite phase, it then decreases up to an intermediate value of |N| ≃ 40, and then increases sharply. This indicates that for the critical exponent α ≃ −0.21, the probability distribution, relative to the Skellam distribution, is narrower for the intermediate, and wider for larger values of |N|. With the Z(2) exponent α ≃ 0.11, the ratio is very similar to the MF one for |N| < 40, while at larger |N| it shows a much sharper increase than either the MF and O(4) distribution. These properties are reflected also in the critical behavior of the cumulants.
The canonical partition function Z(T,V, N) can also be computed by applying the method of the steepest descent to the integral in Eq. (15). The region of applicability of this method is found, however, to be limited to quite small N ∼ 10, for the same set of parameters as in Fig. 4, because of the analytic structure of the integrand.

Reconstructing cumulants from the probability distributions
We have argued, that the structure of the net charge probability distribution is sensitive to the phase transition and its properties. To verify if the observed characteristic structures of P(N) are entirely due to criticality, one would need to check, whether the cumulants reconstructed from these distributions have properties expected near the critical point. We calculate different cumulants with exact results, except for some deviations in the immediate neighborhood of the phase transition, at a non-vanishing chemical potential. In the model with the O(4) exponent, the P α (N) calculated numerically up to |N| ≤ 67, at µ = 0 and T /T c = 0.99, reproduces the exact results of cumulants in the vicinity of the critical point, as seen in Fig. 6. This is not only the case for the non-critical second-and fourth-order cumulants, but also for the critical sixth-order cumulant, which diverges at the critical point. Similarly, one observes that the divergent behavior of c 4 and c 6 for the Z(2) exponent is also well reproduced by the corresponding probability distribution. 6 6 We note, that at finite chemical potential a multiplicative factor e β µN in Eq. (13) enhances the significance of large N-contributions to P(N). This causes difficulties to reproduce the critical structure of a higher order cumulants at finite µ directly from the numerically calculated probability distribution.
Since both the O(4) and Z(2) models exhibit the positively divergent c 6 at µ = 0, the long tail in the probability distribution can be attributed to this divergence. On the other hand, the difference of the O(4) and Z(2) probability distributions up to |N| ∼ 40 shown in Fig. 5 can be understood as a consequence of the different sign of the critical exponent. In the chiral limit, the contribution of the singular part of the thermodynamic potential to the fourth order cumulant is negative and vanishes at the critical temperature [14]. This means that, for α < 0, c 4 is smaller than the Skellam counterpart for T = T c . This accounts for the narrowing around |N| ∼ 40 shown in Fig. 5.
These results confirm our main conclusion, that the characteristic behavior of the net charge probability distributions in the O(4) and Z(2) symmetric models near the critical temperature, is a reflection of the corresponding second order phase transition.

Conclusions
We have explored the influence of the second order chiral phase transition on the properties of the probability distribution P(N) of the net baryon number.
We have modeled the critical behavior within the Landau theory of phase transitions, where the singular part of an effective thermodynamic potential is expressed as a polynomial in the order parameter. The coefficients of the potential were parameterized, so as to reproduce the mean-field, the O(4) and the Z(2) scaling behavior of cumulants of the net charge. The structure of the regular part of the potential was adopted from the hadron resonance gas model to reproduce the Skellam function for the net-baryon number probability distribution.
In the mean-field approach, the probability distribution was computed analytically. In the model which reproduces the O(4) and Z(2) scaling of the net charge fluctuations, the probability distribution was obtained numerically.
We have shown, that the shape of the distribution depends on the universality class of the second-order phase transition. In particular, with the mean field and Z(2) critical exponents, the probability distribution P(N) is characterized by a somewhat sharper peak and broader tail, compared to the non-singular Skellam function. On the other hand, in the O(4) model, the distribution is less peaked but narrower for intermediate values of |N|, followed by a wider tail for large |N|.
We have explicitly demonstrated, at a vanishing chemical potential, that the above behaviors of the net baryon number probability distributions indeed capture the critical properties of a system with the O(4) and Z (2) scaling. This has been done by reconstructing the critical structure of the higher order cumulants directly from these distributions.
Although our model approach does not provide a quantitative description of the probability distribution relevant for QCD and heavy ion phenomenology, the qualitative aspects of our results on the influence of criticality on the probability distribution and its properties are nevertheless of general validity in the chiral limit. Effects of an explicit breaking of the chiral symmetry, together with a more consistent treatment of the critical fluctuations, have been recently examined in [34]. These results indicate, that a direct comparison of measured probability distribution of conserved charges obtained in heavy-ion collisions with a non-critical distribution of an ideal gas, e.g. the Skellam function, can provide pertinent information on the underlying chiral phase transition of QCD and its origin.