Higher order conserved charge fluctuations inside the mixed phase

General formulas are presented for higher order cumulants of the conserved charge statistical fluctuations inside the mixed phase. As a particular example the van der Waals model in the grand canonical ensemble is used. The higher order measures of the conserved charge fluctuations up to the hyperkurtosis are calculated in a vicinity of the critical point (CP). The analysis includes both the mixed phase region and the pure phases on the phase diagram. It is shown that even-order fluctuation measures, e.g. scaled variance, kurtosis, and hyperkurtosis, have only positive values in the mixed phase, and go to infinity at the CP. For odd-order measures, such as skewness and hyperskewness, the regions of positive and negative values are found near the left and right binodals, respectively. The obtained results are discussed in a context of the event-by-event fluctuation measurements in heavy-ion collisions.


I. INTRODUCTION
The structure of the QCD phase diagram is one of most interesting unsolved problems in physics. Statistical fluctuations of conserved charges are regarded to be sensitive probes of the critical behavior in strongly interacting matter [1][2][3][4][5][6]. The fluctuations can be quantified in terms of cumulants (susceptibilities) of the conserved charge distribution. Without the loss of generality, we will specifically refer to net baryon number B throughout this work. Cumulant of order j can be written as follows: , j = 1, 2, . . . , where ... denotes the ensemble average. Cumulants can be expressed in terms of the moments B r explicitly [7] κ j = j k=1 (−1) k−1 (k − 1)!B j,k ( B , . . . , B j−k+1 ). (2) Here B j,k are partial exponential Bell polynomials. It useful to consider ratios of cumulants, as such quantities are intensive, i.e. they are size-indepedent in the thermodynamic limit. The most familiar such measures are scaled variance ω, skewness Sσ, and kurtosis κσ 2 (see e.g. Ref. [8]): Higher order measures such as hyperskewness κ 5 /κ 2 and hyperkurtosis κ 6 /κ 2 are also used.
The statistical fluctuations are sensitive to presence of a first order phase transition (FOPT). The endpoint of the FOPT is the critical point (CP), where the cumulant ratios exhibit singular behavior. At the end point of the FOPT, denoted as in the thermodynamic limit V → ∞. The larger the order of a cumulant is, the stronger is its sensitivity to critical phenomena.
At least two FOPTs are relevant for the QCD phase diagram: (i) the nuclear liquid-gas transition at small temperature T and large baryon chemical potential µ well established both theoretically [9][10][11][12][13][14][15][16] and experimentally [17][18][19] and (ii) the hypothetical first-order chiral phase transition at finite baryon densities [2,4,20,21]. Both transitions are expected to influence the baryon number fluctuations. Model calculations suggest that the behavior of baryon number cumulants in certain regions of the phase diagram is determined by a complex interplay of the chiral and liquid-gas phase transitions [22,23].
Large attention has been given to the structure of higher order measures of fluctuations of conserved charges at supercritical temperatures and in pure phases (see e.g. Refs. [2,4,[20][21][22][23][24][25]). On the other hand, less attention has been paid to the mixed phase. Nevertheless, it is feasible that system created in relativistic nucleus-nucleus collisions can enter the mixed phase of a FOPT under certain conditions. This is especially relevant in view of the plans of HADES collaboration at Helmholtzzentrum für Schwerionenforschung (GSI) to arXiv:2011.06420v1 [hep-ph] 12 Nov 2020 measure the higher order net-proton and net-charge fluctuations in central Au+Au reactions at collision energies E lab = 0.2 − 1.0 AGeV to probe the liquid-gas FOPT region [26]. The freeze-out of the expanding system created in collisions at these energies may well take place in the mixed phase of the nuclear liquid-gas FOPT. From a theoretical point of view, it is convenient to study the statistical fluctuations using the grand canonical ensemble (GCE). In the GCE, the cumulants are determined by partial derivatives of the pressure p with respect to a corresponding chemical potential µ: Here V and T are the system volume and temperature, respectively. In (µ, T )-variables, the FOPT is a line at subcritical temperatures (T < T c ) that ends at the CP. Each point on this line corresponds to a coexistence of two phases: a diluted 'gas' phase with density n 1 and a dense 'liquid' phase with density n 2 . The total baryon density n is a superposition of the gaseous and liquid phase densities, and lies anywhere in the range n ∈ [n 1 , n 2 ]. For this reason it is more appropriate to study the mixed phase phenomena using density-temperature variables (n, T ) instead. Figure 1 depicts a typical phase diagram for a system with the liquid-gas FOPT calculated within the van der Waals (vdW) model (see Sec. III). A large fraction of the (n, T ) plane at T < T c corresponds to the mixed phase. At each point of the mixed phase the pressures of the first and second phases are equal, p 1 (T, µ) = p 2 (T, µ). This is a manifestation of the so-called Gibbs equilibrium condition for the FOPT. However, the T -and µ-derivatives of the functions p 1 and p 2 are different. Therefore, the statistical fluctuations of conserved charges given by Eq. (4) differ between the first and second phases.
In this paper, we present a general formalism to calculate the GCE conserved charge cumulants in the mixed phase. The formalism, presented in Sec. II, takes into account the statistical fluctuations in each of the two phases that comprise the mixed phases, as well as fluctuations in the volume fractions occupied by each of the phases. The formalism is then applied to describe the behavior of cumulants up to sixth order in the mixed phase of a vdW fluid (Sec. III). Summary in Sec. IV closes the paper.

II. GRAND-CANONICAL FLUCTUATIONS IN THE MIXED PHASE
The total system volume V is partitioned in the mixed phase into volumes V 1 = xV and V = yV occpupied by the first and second phases, respectively. Here y ≡ 1 − x.
The r-th moment of conserved charge distribution is the following: Here ρ 1 ≡ B 1 /V 1 and ρ 2 = B 2 /V 2 are the baryon densities in the first and second phase, respectively, and . . . corresponds to the GCE averaging. The fluctuating quantities are the densities ρ 1 , ρ 2 and the volume fraction x, whereas the total volume V is fixed. Following Refs. [27,28] we assume that the fluctuations of all these quantities are independent in the thermodynamic limit, i.e. ρ l 1 ρ m 2 x n = ρ l 1 ρ m 2 x n for any non-negative integers l, m, and n.
It is instructive to start with the first moment, r = 1. Equation (5) in this case reduces to where x 0 = x is the mean volume fraction occupied by the first phase, y 0 ≡ 1 − x 0 , and n 1 = ρ 1 , n 2 = ρ 2 are the mean densities in the first and second phases, respectively. B defines the mean baryon density n in the system: n ≡ B /V . Equation (6) defines x 0 in terms of the mean densities: To obtain all other cumulants one substitutes Eq. (5) into (2). The first three cumulants read Here κ j,1 (κ j,2 ) is a jth order cumulant of the B 1 (B 2 ) fluctuations in the first (second) phase. These cumulants describe fluctuations in a pure phase, thus, they should be calculated according to Eq. (4). κ j,x is the jth order cumulant of the volume fraction x distribution. It is expressed in terms of cumulants of a subvolume In the thermodynamic limit, V → ∞, all cumulants of extensive quantities are proportional to the system volume V : i.e. κ j,1(2) ∼ V and Leaving in Eqs. (8)-(10) only the terms that are linear in V one obtains the following expressions for the cumulants in the thermodynamic limit: Cumulants κ j,x of the volume fraction parameter x distribution can be expressed in the thermodynamic limit in terms of the GCE cumulants κ j,1 and κ j,2 of the two phases. Details of this calculation are given in the Appendix A. Note that the total cumulants reduce to the sum of cumulants of fluctuations in the two phases if the xfluctuations are neglected, i.e. for x ≡ x 0 one has κ j κ j,1 + κ j,2 , j = 1, 2, . . . .
Note that Eq. (13) is also valid away from the thermodynamic limit, i.e. at finite values of the system volume V . It is also instructive to rewrite Eqs. (11) and (12) in terms of susceptibilities, χ j ≡ κ j /(V T 3 ), which are the intensive measures of particle number fluctuations. One obtains The susceptibility χ j of particle number fluctuations in the mixed phase corresponds to a linear combination of the pure phase susceptibilities at the left and right binodals plus the contribution from the x-fluctuations.
Equations (11) and (12) [as well as (14), (15)] are model-independent expressions describing the GCE fluctuations of a conserved charge in the mixed phase of a FOPT in the thermodynamic limit. Model dependence will enter only through explicit form of the cumulants κ j,1 and κ j,2 .

III. VAN DER WAALS FLUID
In this section we illustrate the general formalism introduced in the previous section using the vdW equation of state describing Maxwell-Boltzmann interacting particles. Here we neglect the antiparticles, therefore, the number of particles plays the role of the conserved charge.
The system pressure of a vdW fluid in a pure phase reads where a > 0 and b > 0 are the model parameters describing the attractive and repulsive interactions, respectively. The CP is defined by conditions [29,30] ∂p which give Introducing reduced variables T = T /T c , n = n/n c , and p = p/p c one can rewrite the vdW equation (16) in a universal form which is independent of the specific numerical values of the interaction parameters a and b. This is a particular case of the principle of the corresponding states (see, e.g. Ref. [29]). The phase diagram of the vdW fluid in the (n, T ) plane is presented in Fig. 1.
In the GCE, the vdW model particle number density can be written as follows [31]: where n id (T, µ) is the ideal gas density in the GCE. The cumulants of vdW model particle number distribution in the GCE can be calculated up to a desired order by iteratively differentiating Eq. (20) with respect to the chemical potential µ (see Ref. [32] for the technical details). The resulting expressions for the scaled variance, skewness, and kurtosis for the case of pure phases read [6,31]: Following the same procedure we also calculate the GCE hyperskewness κ 5 /κ 2 and hyperkurtosis κ 6 /κ 2 . As the resulting expressions are very lengthy, we do not list them here. Note that Eqs. (21)-(23) describe cumulants of the GCE particle number distribution in pure phases, i.e. at all densities at T > T c and outside the mixed phase region at T < T c , as well as in metastable phases. Calculation of fluctuations in the mixed phase, however, requires the use of Eqs. (11), (12). The boundaries of the mixed phase -the left and right binodals n 1 ( T ) and n 2 ( T ) -are defined by the Gibbs equilibrium conditions: µ( T , n 1 ) = µ( T , n 2 ) and p( T , n 1 ) = p( T , n 2 ). In the case of vdW model, these conditions lead to 9 4 The binodals are depicted in Figs The calculations show that the x-fluctuations do not affect ω, Sσ, and κσ 2 significantly. The only exception is a close proximity to the right binodal, where x 1 (see Appendix A). The differences between Eqs. (13) and (12) would be barely visible in Fig. 2. Thus, ω, Sσ, and κσ 2 calculated using Eq. (12) are not presented in Fig. 2.
As seen from Fig. 2 (a), the scaled variance ω → +∞ at the CP, both inside and outside the mixed phase. This is in agreement with Ref. [27]. A behavior of Sσ shown in Fig. 2 (b) is also similar inside and outside the mixed phase. However, the behavior of κσ 2 inside and outside the mixed phase differs drastically. This is shown in Fig. 2 (c). Outside the mixed phase, κσ 2 attains large positive or negative values in vicinity of the CP, it approaches +∞, −∞, or 0 at T → T c , n → n c , i.e., its value depends on the path of approach toward the CP in the (n, T ) plane. Inside the mixed phase, κσ 2 is only positive, with large values attained in vicinity of the CP. κσ 2 tends to +∞ when approaching the CP from inside the mixed phase. Negative values of κσ 2 are observed at T > T c only. At T < T c the values of κσ 2 are only positive, this is a reflection of the fact that κσ 2 in pure phases outside the mixed phase is positive at subcritical temperatures [27]. It thus follows from Eq. (13) that κσ 2 in the mixed phase is positive as well. Large negative values of κσ 2 near the CP are only possible in a small region outside the mixed phase at supercritical temperatures T T c The same qualitative structure of kurtosis in pure phases around the CP is also present in the Ising model [21]. Our conclusions regarding the behavior of κσ 2 in and around the mixed phase region near the CP thus applies to the Ising model as well.
Similar arguments are applicable for the hyperkurtosis shown in Figs. 3 (c) and (d) and higher fluctuation measures of even order. For instance, the structure of hyperkurtosis is more involved in comparison with κσ 2 . The band of large positive values at n ≈ n c is surrounded by two bands of large negative values. However, close to binodals κ 6 /κ 2 becomes positive again (in this aspect κ 6 /κ 2 is similar to κσ 2 ). As a consequence, the hyperkurotisis is positive in the whole mixed phase region in accordance with Eq. (13). Thus, we conclude that hyperkurtosis is mostly positive in a vicinity of the CP if the mixed phase region is taken into account. Negative values of κ 6 /κ 2 appear within two relatively narrow bands at supercritical temperatures, T > T c , on the (n, T ) phase diagram.
In general, the structure of the fluctuation measures outside the mixed phase becomes increasingly complex with an increase of their order. This is not the case inside the mixed phase region. Only positive values of the even-order measures are found in a vicinity of the CP inside the mixed phase. For the odd-order fluctuation measures, e.g. hyperskewness, the mixed phase is split into two regions: (i) the first region has positive values and generally corresponds to lower densities close to the left (gaseous) binodal and (ii) the second region with negative values of odd order cumulants at higher densities close to the right (liquid) binodal. This is shown for the case of hyperskewness in Figs

IV. SUMMARY
In the present work we determined the thermal (grandcanonical) fluctuations of a conserved charge inside the mixed phase of a first-order phase transition. As opposed to fluctuations in pure phases, where they are determined solely by the equilibrium properties of that single phase, in the mixed phase the cumulants receive contributions from fluctuations in both the gaseous and the liquid phases, as well from the fluctuations in the volume fractions occupied by the two phases. Our main result here is given by Eqs. (11) and (12), which express the grand-canonical conserved charge cumulants for a point inside the phase coexistence region in the thermodynamic limit for any equation of state with a first-order phase transition.
A mixed phase cumulant κ j of order j reduces to a sum of the corresponding pure phase cumulants κ j,1 and κ j,2 from each of the two phases plus a contribution from the fluctuations of the relative volume fraction x occupied by the gaseous phase, the cumulants of the x-distribution are denoted as κ j,x . The cumulants κ j,1 (2) should be calculated in the standard way -as derivatives of the grand potential with respect to the chemical potential on the left (right) side of the phase coexistence line in the (µ, T ) plane. The x-fluctuation cumulants κ j,x can be expressed in terms of the cumulants κ j,1 (2) , although their calculation can be quite involved (see Appendix A for explicit results up to fourth order). We do observe, however, that the effects of the x-fluctuations are found to be negligible for cumulants up to fourth order. For the fifth and the sixth orders notable contributions of the xfluctuations appear in a vicinity of the right binodal 1 . Therefore, the simple approximate relation (13) can be used in most practical applications. This also implies that fluctuations in the mixed phase are mainly determined by the intrinsic properties of the two coexisting phases.
To illustrate our results more explicitly, we used the equation of state of a van der Waals fluid. The cumulant ratios of conserved charge fluctuations such as scaled variance κ 2 /κ 1 , skewness κ 3 /κ 2 , kurtosis κ 4 /κ 2 , hyperskewness κ 5 /κ 2 , and hyperkurtosis κ 6 /κ 2 were calculated both outside and inside the mixed phase region. Outside the mixed phase we reproduce the earlier results from the literature [21,27,33,34], where cumulant ratios exhibit increasingly involved structures as the order is increased. Inside the mixed phase, on the other hand, the structure of cumulants is simpler. The even-order cumulants are predominantly positive while odd-order cumulants can have either sign, generally attaining positive (negative) values close to the left (right) binodal. In particular, we conclude that negative values of the kurtosis are consistent with a crossover region just above the critical point, as discussed in Ref. [21], but not with any region inside the mixed phase of a FOPT.
The obtained results are relevant in the context of heavy-ion collisions, which create a strongly interacting fluid that may pass through a mixed phase of a FOPT at finite baryon density. Such a scenario can be probed by event-by-event fluctuation measurements. In particular, one can consider the well-established nuclear liquid-gas phase transition. A beam energy scan of different colliding ions in a sub-GeV collision energy regime will probe the phase diagram in the vicinity of the nuclear LGPT.
Such a program can be performed by the HADES experiment at GSI. If the specific non-monotonic behaviour of high-order cumulants, as calculated here in the grandcanonical limit, is observed in measurements of higherorder fluctuations, this can be interpreted as a signal of the proximity to the CP.
We note, however, that the results are not directly suitable for quantitative conclusions. The measurements in heavy-ion collisions are performed in coordinate rather than momentum space. Also, even in the scenario where thermalization of fluctuations is achieved in heavy-ion collisions, the measurements will still be affected by global charge conservation, finite system-size, and volume fluctuation effects, as has been discussed in the literature [35][36][37][38]. To address the effects of global conservation in pure phases a subensemble acceptance method has recently been introduced by us in Refs. [39][40][41]. In the future we plan extend this method to address global charge conservation influence on conserved charge fluctuations in the mixed phase. Another interesting avenue is the so-called strongly intensive fluctuation measures [36,42]. These quantities are designed to cancel out the geometric effects of the total volume fluctuations, but they are ex-pected to be sensitive to the critical point of a FOPT [43], thus it is of interest to elucidate their behavior in the mixed phase. As stated in Sec. II, we assume that in thermodynamic limit, i.e. for V → ∞, correlations between baryon numbers B 1 ,B 2 and x can be neglected. Thus, in the following calculations we will use the canonical ensemble with B 1 , B 2 fixed to their mean values B 1 = B 1 , B 2 = B 2 . Furthermore, in thermodynamic limit the surface terms between coexistent phases can be neglected and, in similarity with Refs. [39,41], the partition function can be presented as a product of partition functions of the two subsystems. Then the probability P (x) is proportional to the product of the canonical partition functions of the first and second phase: Here y ≡ 1 − x, B 1 = V x 0 n 1 , B 2 = V y 0 n 2 , and x 0 is given by Eq. (7). In the thermodynamic limit, the above results can be generalized, since in this case the canonical partition function can be expressed through the volume-independent free energy density f : Z(T, V, B) = exp − V T f (T, n) with n ≡ B/V being the conserved baryon density. Thus, To evaluate κ j,x we introduce the cumulant generating function ψ x (t): HereC is an irrelevant normalization constant. The cumulants, κ j,x , correspond to the Taylor coefficients of ψ x (t): Here we have introduced a shorthand,κ j,x (t), for the n-th derivative of generating function at arbitrary values of t, which we subsequently will refer to as t-dependent cumulants. Clearly, all higher order cumulants are given as a t-derivative of the first order t-dependent cumulant,κ 1,x (t), which is given bỹ with the (un-normalized) t-dependent probabilitỹ One can check that in the thermodynamic limit, V → ∞,P has a sharp maximum at the mean value of x, x(t) . The condition ∂P (x; t)/∂x = 0 determines the location of this maximum resulting in an implicit relation that determines x(t) : Here we also used a thermodynamic relation The solution to Eq. (A7) at t = 0 is x(t = 0) = x 0 , as should be by construction.