Connecting fluctuation measurements in heavy-ion collisions with the grand-canonical susceptibilities

We derive the relation between cumulants of a conserved charge measured in a subvolume of a thermal system and the corresponding grand-canonical susceptibilities, taking into account exact global conservation of that charge. The derivation is presented for an arbitrary equation of state, with the assumption that the subvolume is sufficiently large to be close to the thermodynamic limit. Our framework -- the subensemble acceptance procedure -- quantifies the effect of global conservation laws and is an important step toward a direct comparison between cumulants of conserved charges measured in central heavy ion collisions and theoretical calculations of grand-canonical susceptibilities, such as lattice QCD. As an example, we apply our formalism to net-baryon fluctuations at vanishing baryon chemical potentials as encountered in collisions at the LHC and RHIC.

Introduction. Studies of the QCD phase diagram are one of the focal points of current experimental heavy-ion collision programs [1]. Observables characterizing fluctuations of the QCD conserved charges -baryon number, electric charge, and strangenesshave attracted a particular attention, as these are sensitive to the finer details of the QCD equation of state and its phase structure in particular [2][3][4]. Consider, for simplicity, a case of a single conserved charge, say baryon number B, for a system in equilibrium with volume V at temperature T . The nth order scaled susceptibility χ B n is defined as a derivative of the pressure with respect to the chemical potential µ B , and it determines the cumulants, κ n [B], of the distribution of the charge B in the grand canonical ensemble (GCE). The susceptibilities, χ B n , characterize the properties of the thermal system under consideration, in particular they provide information about the possible phase changes, including remnants of the chiral criticality at vanishing chemical potential [5]. Theoretically they are calculated either using firstprinciple lattice QCD simulations [6,7], or in various effective QCD approaches [8,9]. An important question is how to relate these quantities to experimental measurements. The total net charge B does not fluctuate in the course of a heavy-ion collision, as opposed to the case of the GCE where the system can freely exchange the charge with an external heat bath. However, experimental measurements typically have limited acceptance and only cover a fraction of the total momentum space, which we subsequently assume to be characterized by a finite acceptance window in rapidity, ∆Y acc . As discussed e.g. in [4], for a sufficiently small acceptance window ∆Y acc ∆Y 4π conditions corresponding to the GCE may be imitated, i.e. effects of global charge conservation become negligible. However, in order to capture the relevant physics the acceptance window ∆Y acc must be much larger than the correlation length ∆Y cor . Consequently, ∆Y cor ∆Y acc ∆Y 4π is the minimum necessary condition for the applicability of the GCE to fluctuation measurements.
In practice the situation is more subtle. Deviations of lattice QCD calculations of the cumulants of conserved charges at T ∼ 160 MeV from the ideal hadron resonance gas (HRG) expectation do not exceed the magnitude of charge conservation effects already for an acceptance as small as ∆Y acc /∆Y 4π ∼ 0.1 [10]. Therefore, in order to capture the physics of e.g. chiral criticality the effect of charge conservation needs to be understood very well, since simply reducing the acceptance window even further risks eliminating all the non-trivial effects associated with relevant QCD dynamics [11].
In the present letter we generalize the relation (1)    A subsystem (dashed red rectangle) within a thermal system (solid black rectangle). The subsystem can exchange particles (conserved charges), shown by the circles, with the rest of the system. The filled red circles in (a) depict the particles within the subsystem, as considered in the present subensemble acceptance procedure. In contrast, the filled red circles in (b) highlight the particles from a typical configuration resulting from the binomial filter.
between the GCE susceptibilities χ B n and measured cumulants of conserved charge κ n [B] to make it valid for subsystems that are comparable in size to the total system. We will still assume that the size of the subsystem is large enough to capture the relevant physics. Further assuming strong space-momentum correlations, as is the case for LHC and top RHIC energies, the formalism presented here connects the measured cumulants with those obtained in lattice QCD over a wide range of acceptance windows.
Formalism. Consider a spatially uniform thermal system at a fixed temperature T , volume V , and total net charge, say net baryon number, B, which is described by statistical mechanics in the canonical ensemble and characterized by its canonical partition function Z(T, V, B). We pick a subsystem of a fixed volume V 1 = α V within the whole system, which can freely exchange the conserved charge B with the rest of the system (see Fig. 1). Our goal is to evaluate the cumulants κ n [B 1 ] of the distribution of charge B 1 within the subsystem [the red points in Fig. 1(a)]. Our considerations will be quite different from a binomial filter previously considered in Refs. [10,12]. The binomial filter corresponds to an independent acceptance of particles with a probability α from the entire volume V [the red points in Fig. 1(b)]. The binomial filter is appropriate for non-interacting systems, e.g. for the ideal HRG model, where there are no spatial correlations between the particles. Given a finite correlation length, however, particles will be more strongly corre-lated with their neighboring particles than with those far away. The binomial filter artificially suppresses these correlations and thus will not provide the correct results for the subvolume V 1 .
Our arguments will be based purely on statistical mechanics. Assuming the subvolume V 1 as well as the remaining volume V 2 = (1 − α)V to be large compared to correlation length ξ, V 1 ξ 3 and V 2 ξ 3 , the canonical ensemble partition function of the total system with total baryon number B is given by Here β ≡ 1 − α. The probability P (B 1 ) to find B 1 baryons in the subsystem with volume V 1 is proportional to the product of the canonical partition functions of the two subsystems: The procedure based on Eqs. (2) and (3) will be called the subensemble acceptance procedure. Note that for the ideal HRG it reduces to the binomial acceptance sampling 1 .
In the thermodynamic limit, i.e. for V → ∞, 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 : Here ρ B2 = B−B1 V −V1 is the charge density in the second subsystem andC is an irrelevant normalization 1 Strictly speaking, this is valid for the classical ideal HRG when quantum statistics effects can be neglected. This is the case especially for baryons, where, due to their large mass, corrections to baryon number cumulants arising from Fermi statistics are small at the chemical freeze-out, T 155 MeV and µ B 0.
constant. The cumulants, κ n [B 1 ], correspond to the Taylor coefficients of G B1 (t): Here we have introduced a shorthand,κ n [B 1 (t)], for the n-th derivative of generating function at arbitrary values of t, which we subsequently will refer to as tdependent cumulants. Clearly, all higher order cumulants are given as a t-derivative of the first order , which is given bỹ with the (un-normalized) t-dependent probabilitỹ In the thermodynamic limit, V → ∞,P has a sharp maximum at the mean value of B 1 , B 1 (t) [13]. The condition ∂P (B 1 ; t)/∂B 1 = 0 determines the location of this maximum resulting in an implicit relation that determines B 1 (t) : the net baryon number is uniformly distributed between the two subsystems, as it should be by construction. Therefore, The second cumulant is given by the t-derivative of which at t = 0 gives the 2nd order cumulant In order to evaluate the higher-order cumulants κ n [B 1 ] for n ≥ 3 we iteratively differentiate the t-dependent cumulantsκ n [B 1 (t)] with respect to t, starting fromκ 2 [B 1 (t)], and make use of the expression (10) for B 1 (t) . The result for the cumulants up to the 6th order is the following: In the limit α → 0 all susceptibilities, i.e. the cumulants scaled by V 1 T 3 = αV T 3 , reduce to the GCE susceptibilities, as expected, since in this limit effects of global conservation become negligible. However, as α → 0 our assumption of the subsystems being close to the thermodynamic limit breaks down, as the subsystem becomes much smaller than the correlation length, αV ξ 3 . Instead, the cumulants approach the Poisson limit [14]. In the other limit, α → 1, all cumulants of order n ≥ 2 tend to zero, reflecting the dominance of the global conservation laws and the absence of conserved charge fluctuations in the full volume.
It is instructive to consider ratios of cumulants, in which the volume V cancels. The explicit relations for the commonly used scaled variance, skewness, and kurtosis are: The modification of the scaled variance κ 2 [B 1 ]/κ 1 [B 1 ] due to global conservation laws is a multiplication of the grand canonical scaled variance by a factor (1 − α). This is similar to the binomial filter effect studied in prior works [10,12,15]. Same for the skewness κ 3 [B 1 ]/κ 2 [B 1 ], where the corresponding grand canonical ratio is multiplied by (1 − 2α). An interesting case is the kurtosis κ 4 [B 1 ]/κ 2 [B 1 ]: this ratio in the subvolume depends not only on the GCE kurtosis χ B 4 /χ B 2 but also on the GCE skewness χ B 3 /χ B 2 . If α is known, Eqs. (16)-(18) may be inverted to express the GCE cumulant ratios in terms of those of the subsystem.
Net baryon fluctuations at LHC and top RHIC energies. We apply our formalism to study the effect of baryon number conservation in view of measurements of net proton number distributions in heavy-ion collisions at the RHIC and LHC. The ALICE collaboration has published measurements of the variance of net proton distribution [16] and the analysis of higher orders up to κ 4 is in progress. In the future runs, sufficient statistics may be accumulated to extend the measurements up to the 6th order [17]. The STAR collaboration has measured the cumulants of the net proton distribution up to κ 4 [18,19], preliminary results for κ 6 are also available [20].
It should be noted that experimental measurements in heavy-ion experiments are performed in momentum space rather than in coordinate space. However, the momenta and coordinates of particles at freeze-out are correlated due to the presence of a sizable collective flow, in particular the longitudinal flow. The correlation is one-to-one in the case of a Bjorken scenario, which can be expected to be approximately realized at the highest collision energies achievable at the RHIC and LHC. In that case, the experimental momentum cuts correspond to cuts in coordinate space and our formalism is applicable. 2 In the other extreme, when no collective motion is present, cuts in momentum space do not correlate with a definite subvolume in coordinate space. In that case the binomial acceptance may be the appropriate procedure. We also note that our calculations apply to net baryon fluctuations rather than net proton ones. Experimentally, the former can be reconstructed from the latter following the binomial-like method developed in Refs. [21,22]. 2 We note that thermal smearing by ∆Y th ∼ (T /m) 1/2 somewhat dilutes the space-momentum correlation. In case of baryons, which are heavy, this effect is rather small, especially if a rapidity window of ∆Yacc 2 is considered. We expect baryon smearing to slightly shift our results for the cumulant ratios towards that obtained using the binomial filter. A detailed study of this correction is under way.
The typical chemical freeze-out temperatures, T ch ∼ 155 − 160 MeV at the LHC [23][24][25] and T ch ∼ 160 − 165 MeV at the top RHIC energies [26], are close to the pseudo-critical temperature of the QCD crossover transition determined by lattice QCD T pc 155 − 160 MeV [27,28] at µ B = 0. Also, in the vicinity of T pc lattice calculations predict a change of sign of χ B 6 , which is thought to be related to the remnants of the chiral criticality [5], although alternative explanations do also exist [29,30]. Therefore, it would be of great interest to verify the theory prediction of a negative χ B 6 experimentally. As all odd order susceptibilities vanish at µ B = 0, the relations between the higher-order cumulants  15) and (11), We then study how the cumulant ratios κ 4 /κ 2 and κ 6 /κ 2 of net baryon distribution depend on the value of the parameter α characterizing the subsystem where the fluctuations are measured. In our subensemble acceptance procedure we use as input the lattice data for χ B 4 /χ B 2 and χ B 6 /χ B 2 at T = 155 and 160 MeV from Ref. [7]. The results for the αdependence of κ 4 /κ 2 and κ 6 /κ 2 ratios calculated from Eqs. (18) and (19) are shown in Fig. 2. The solid red lines correspond to T = 160 MeV and the bands depict the error propagation of the lattice data. The dash-dotted lines correspond to T = 155 MeV. We also show the binomial acceptance results as dashed lines. These results are only valid for the classical ideal HRG which gives χ B 4 /χ B 2 = 1 and χ B 6 /χ B 2 = 1. As both the kurtosis and hyperkurtosis are symmetric with respect to α ↔ (1−α) [14] we plot our results only up to α = 0.5. In the limit α → 0 both κ 4 /κ 2 and κ 6 /κ 2 approach their GCE values. The computed values of κ 4 /κ 2 and κ 6 /κ 2 lie below the binomial acceptance baseline for all values of α, which reflects the suppression of the lattice values for χ B 4 /χ B 2 and χ B 6 /χ B 2 relative to the ideal HRG baseline. Interestingly, the difference between the ideal gas and QCD is the smallest for κ 6 /κ 2 at α = 0.5, where the effects of baryon conservation are the strongest. Actually, in the entire region 0.2 < α < 0.8 the difference is so small that it may be difficult to distinguish the true dynamics of QCD from that of an ideal HRG. Measurements in this region of α, on the other hand, may serve as a model-independent test of baryon number conservation effects. For α < 0.2, however, the measurable ratio κ 6 /κ 2 become sensitive to the equation of state, i.e. to the actual value for χ B 6 /χ B 2 . We find that a negative κ 6 /κ 2 for α 0.1 is consistent with χ B 6 /χ B 2 which is either negative or close to zero. Such a measurement would constitute a potentially unambiguous experimental signature of the QCD chiral crossover transition.
If we apply these conditions to actual experiments such as ALICE and STAR, it translates into the following: At the LHC (ALICE) with √ s NN = 5.02 TeV, the beam rapidity is y beam ln( √ s NN /m N ) 8.5 while for the top RHIC energy ( √ s NN = 200 GeV) one has y beam 5.4. Thus, α 0.1 would correspond to measurements within approximately two (one) units of rapidity for LHC (RHIC).
As discussed above, at α below a certain value, α < α lim , our formalism breaks down and the cumulants approach the Poisson limit instead. The value of α lim can be estimated. The physical volume used in lattice calculations [6,7,31,32] of χ B 2n at T 160 MeV is of order V lat = (a N σ ) 3 ∼ 10 2 fm 3 for a = 0.1 fm and N σ = 48 lattices. We can thus assume that volumes V ≥ V lat are sufficiently large to capture the relevant physics. The total volume in central collisions at the LHC can be estimated as V tot = (dV /dY ) 2 y beam , which for dV /dY ∼ 5000 fm 3 [23] and y beam ∼ 8 at the LHC yields V tot ∼ 80000 fm 3 . Therefore, α lim = V lat /V tot ∼ 10 −3 . Even an order of magnitude error in this estimate implies that α lim does not exceed 10 −2 , and thus our method is applicable for virtually the entire linear scale shown in Fig. 2. The same estimate for RHIC, where dV /dy 1500 fm 3 [26], gives α lim ∼ 7 · 10 −3 .
Summary. We presented a novel procedure to connect measurements of cumulants of conserved charge fluctuations in a finite acceptance to the grandcanonical susceptibilities, taking into account effects due to exact charge conservation. In contrast to prior works studying the ideal HRG model, our subensemble acceptance procedure works for an arbitrary equation of state, under the assumption that the acceptance is sufficiently large to reach the thermodynamic limit, and thus to capture all the relevant physics. The formalism is most suitable for central collisions of ultrarelativistic heavy-ion collisions at the highest energies where we have a strong space-momentum correlations, and it enables direct comparisons between experimental data on cumulants of conserved charges and theoretical calculations of grand-canonical susceptibilities within effective QCD theories and lattice QCD simulations. We consider our results to be particularly helpful for the ongoing experimental effort to study the QCD phase structure with fluctuation measurements. As a first application, we have studied the conditions under which a measurement of a net baryon hyperkurtosis κ 6 /κ 2 can serve as an experimental signature of the QCD chiral crossover at µ B = 0, and we found a rapidity window of ∆Y acc 2(1) at LHC (RHIC) to be the sweet spot, where the QCD dynamics is not overshadowed by baryon number conserva-tion effects.
Our framework opens a number of new avenues to explore. In the future we plan to test the limits of our approach in finite systems close to the critical point, where the correlation length becomes comparable to the system size [33]. Another extension is a simultaneous incorporation of multiple conserved charges, which will enable the analysis of the offdiagonal susceptibilities, a relevant topic in light of the corresponding measurements recently performed by the STAR collaboration at RHIC [34].