Critical point fluctuations: Finite size and global charge conservation effects

We investigate simultaneous effects of finite system size and global charge conservation on thermal fluctuations in the vicinity of a critical point. For that we consider a finite interacting system which exchanges particles with a finite reservoir (thermostat), comprising a statistical ensemble that is distinct from the common canonical and grand canonical ensembles. As a particular example the van der Waals model is used. The global charge conservation effects strongly influence the cumulants of particle number distribution when the system size is comparable to that of the reservoir. If the system size is large enough to capture all the physics associated with the interactions, the global charge conservation effects can be accurately described and corrected for analytically, within a recently developed subensemble acceptance method. The finite size effects start to play a significant role when the correlation length grows large due to proximity of the critical point or when the system is small enough to be comparable to an eigenvolume of an individual particle. We discuss our results in the context of fluctuation measurements in heavy-ion collisions.


I. INTRODUCTION
The structure of the phase diagram of QCD matter is one of the most interesting unsolved problems in physics. Within phenomenological statistical models as well as in lattice QCD simulations mainly the grand canonical ensemble (GCE) is used. The critical behaviour is then probed by the statistical fluctuations of conserved charges [1][2][3][4][5][6]. Useful measures of these fluctuations are the scaled variance ω, as well as the (normalized) skewness Sσ and kurtosis κσ 2 . For example, for the net baryon number B they are defined as the following, where ... denotes the GCE averaging and ∆B ≡ B − B . These quantities can also be expressed through baryon number cumulants κ n : The GCE cumulants are calculated as the partial derivatives of the system pressure p with respect to a corresponding chemical potential µ: Here V and T are the system volume and temperature, respectively. The ratios of cumulants in Eq. (3) are intensive (size-independent) measures in the GCE.
The GCE cumulants evaluated in effective QCD models can be directly compared with lattice QCD predictions, a procedure often used for testing and constraining various models and approaches [7][8][9][10][11][12][13][14][15]. On the other hand, a comparison of theoretical predictions with the event-by-event fluctuation measurements in relativistic heavy ion collisions looks rather challenging. In the GCE the system of volume V may exchange particles (and conserved charges) with a reservoir (thermostat) of volume V 0 − V . In the total volume V 0 the conserved charge is strictly fixed. Thus, volume V 0 corresponds to a canonical ensemble (CE). To reach the GCE conditions inside the volume V one has to require V /V 0 1. And while direct comparisons of the GCE cumulants with experimental data are commonplace in the literature [11,[16][17][18], it is clear that the global charge conservation will arXiv:2004.14358v1 [hep-ph] 29 Apr 2020 influence to some extent the conserved charge distribution measured in experiment, making it different from the GCE baseline. Studies based on the ideal hadron gas model indeed show that higher-order cumulants of baryon number are strongly affected by the global conservation [19]. In addition, the volume V should also be large enough to take into account all relevant physical effects due to particle interactions. If both V 0 and V are large enough, one can derive analytically modifications which come from global conservation laws. This has been shown in a recent paper [20] and will be discussed later in the present study.
In high energy nucleus-nucleus collision experiments not all final particles are measured on an event-by-event basis. Within a statistical approach the subset of measured particles can be treated as a subsystem with finite volume V , whereas non-detected particles play the role of the finite reservoir (thermostat). In this situation, the effects of exact charge conservation on fluctuations are usually modeled by a binomial acceptance correction procedure [19,[21][22][23][24]. This procedure assumes that the probability to be measured is the same for each particle of a given type and it is independent of any inter-particle correlations. As will be seen below, the binomial acceptance procedure can be justified only for a statistical system of classical non-interacting particles.
Previously, the finite size effects (without conservation law effects) for the first order liquid-gas phase transition were discussed in Refs. [25][26][27][28]. The effects of finite particle number sampling on baryon number fluctuations were studied in Ref. [29] within fluid dynamical simulations. A Monte Carlo procedure allowing to sample particle multiplicities in the presence of excluded volume effects was developed in Ref. [30].
The size of the considered system becomes especially important in the vicinity of the critical point (CP). The CP as the end point of the first order phase transition exists as a universal feature of all molecular systems. At the CP the intensive fluctuation measures become singular in the thermodynamic limit V → ∞. These infinite values evidently cannot appear in a finite V . Therefore, both the charge conservation and finite size effects can be equally important in the vicinity of the CP.
It was demonstrated [12] that the nuclear CP, i.e., the end point of the liquid-gas transition in the system of interacting nucleons at small T and large µ B , affects the susceptibilities of conserved charges even at µ B = 0 and large T , and limits the radius of convergence of Taylor expansion in µ B /T at µ B = 0 [31]. The sought-after hypothetical chiral QCD CP is expected to produce strong signals in high-order fluctuation measures [2,4,32,33]. It is possible that conserved charge susceptibilites are determined by a complex interplay of the chiral and liquidgas phase transitions in certain regions of the phase dia-gram [34,35].
Our paper presents a first step to study both the finite size and global charge conservation effects in the vicinity of a CP. To give a specific example we consider a classical statistical system of interacting nucleons described by the van der Waals (vdW) model. This model was previously applied to nuclear matter considered as a system of interacting nucleons in Ref. [36]. The production of antibaryons will be neglected. In this situation the number of nucleons becomes a conserved charge. We will not consider the mixed phase region at T < T c in the present study, and will focus our studies at (super)critical temperatures, T ≥ T c .
Our considerations will be based purely on equilibrium statistical mechanics. The non-equilibrium effects in heavy-ion collisions are certainly important, especially in the vicinity of the CP, and the dynamical theory of critical fluctuations is under development [37][38][39][40][41]. We plan to incorporate the non-equilibrium effects in future works.
The paper is organized as follows. Section II presents the main properties of the vdW model in the thermodynamic limit. Section III describes the model results for the particle number fluctuations in the finite systems. The summary in Sec. IV closes the paper.

II. VAN DER WAALS MODEL
The CE partition function, Z ce , for the vdW system of classical particles can be written as [42] Z ce (N, where N , V , and T are respectively the number of particles, volume, and the temperature of the system, while a > 0 and b > 0 are the vdW interaction parameters. The a parameter regulates the attraction, while b corresponds to a repulsion between particles via the excluded volume effects. The function ϕ is given as where g and m are respectively the degeneracy factor and the mass of the particles, and K 2 is the modified Bessel function of the second kind.
The system pressure in the CE is calculated as where n ≡ N/V . The CP is defined by the conditions [42,43] ∂p ∂n T = 0 , ∂ 2 p ∂n 2 which gives Introducing the reduced variables T = T /T c , n = n/n c , and p = p/p c one can rewrite the vdW equation (7) 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. [42]).
To calculate the particle number fluctuation measures one usually transforms the CE description into the GCE one. This requires to introduce a reservoir and to take the thermodynamic limit with V → ∞. For the vdW equation of state these steps were done for the first time in Ref. [44]. In the vdW model the particle number fluctuation measures with ∆N ≡ N − N , were calculated analytically in the GCE in the thermodynamic limit V → ∞ [6,44]: The GCE fluctuation measures (13)-(15) are presented in Fig. 1. All three of them exhibit singular behavior at the CP. While the ω gce (13) tends to +∞ at the CP, the Sσ gce (14) and κσ 2 gce (15) have a richer structure in a vicinity of the CP. They can tend to +∞, −∞, or 0 depending on the path of the approach to the CP. Introducing quantities ρ = n − 1 and τ = T − 1 one finds  13) and (15) in the ( n, T ) plane. The Poisson limit with ωgce = Sσgce = κσ 2 gce = 1 corresponds to those regions of the thermodynamic plane where the vdW interactions become negligible. The two points on the phase diagram marked by crosses are analyzed in detail in the present paper.
at τ 1 and ρ 1: While the CP signals of ω fade out as one moves away from the CP in the phase diagram, they remain stronger in the higher-order fluctuation measures, Sσ and κσ 2 , even far away from the CP [45]. Note that in the classical ideal gas case, a = b = 0, all fluctuation measures in Eqs. (13)-(15) are reduced to ω gce = Sσ gce = κσ 2 gce = 1 which corresponds to the Poisson N -distribution. The general features of the GCE fluctuations presented in this section, especially those connected with the CP remain the same for all models from the mean-field universality class, to which the vdW model belongs to (see e.g. Ref. [46]).

III. FLUCTUATIONS IN A SUBENSEMBLE
Let us partition a finite system of volume V 0 into a subsystem of volume V < V 0 and another subsystema reservoir -of volume V 0 − V . We assume that both subsystems can exchange particles, but the total number of particles N 0 in the whole system is fixed. The corresponding statistical ensemble will be referred as a subensemble, distinguishing it from both the CE and GCE. We neglect all interactions at the interface, i.e. between all particles from different subsystems. The partition function of the system in volume V can then be written as [20,42]: The probability to find N particles in the volume V takes the form The mean value N and the central moments (∆N ) k with k = 2, 3, . . . in the subensemble are calculated with the probability distribution (20): In our example of the vdW model, the CE partition functions Z ce in Eq. (19) are given by Eq. (5). Introducing variables the partition function (19) is written as Here and it cancels out in the probability distribution (20).
The minimal and maximal numbers of particles in the subensemble, (5), which is due to the excluded volumes. Here ... is a floor function. N min and N max can also be rewritten as: The moments (21) are independent of particle's degeneracy g and mass m, since they only enter Eq. (23) through the common factor η 1 . In the following we explore the behaviour of fluctuations in the subensemble for different values of x and N 0 .

A. Charge conservation effects
As a first specific case, we consider the thermodynamic limit, V 0 → ∞, at 0 < x < 1. Thus, both V 0 → ∞ and V → ∞, but the values of x = V /V 0 remain finite.
In this case we follow a recently developed subensemble acceptance procedure [20]. It allows to obtain the cumulants of particle number distribution in the subensemble in terms of the corresponding GCE cumulants and the volume fraction x, quantifying the corrections to the GCE cumulants because of the global conservation of particle number. One obtains (see Ref. [20] for the derivation details): where κ gce n ∝ V 0 is the n-th cumulant in the GCE, and Note that ξ n (x) correspond to the n-th cumulant of the Bernoulli distribution, p x (l) = x l (1 − x) 1−l for l = 0, 1. Similarly, higher order κ n and ξ n cumulants can be obtained. Using Eqs. (29) and (32) one finds the scaled variance, skewness, and kurtosis: Equations (33) It is instructive to consider the limit of an ideal classical gas. This limit is recovered for a = 0 and b = 0. In this case, κ gce n = N 0 for all n = 1, 2, . . ., and Eqs. (29) reduce to κ id n = ξ n N 0 . Therefore, one obtains the following for the classical ideal gas: The fluctuation measures (39)-(41) are presented in Fig. 2. Equations (39)-(41) coincide with those obtained after the binomial acceptance correction procedure (see Ref. [24] for details). The binomial acceptance procedure is suitable for describing the global charge conservation 2 Note that we still assume here that the volume V is large enough to neglect the finite-size effects, no matter how small x is. effects in non-interacting systems. The applicability of the binomial acceptance, however, does not extend to interacting systems, the vdW model in particular.

B. General case of a finite reservoir
In this subsection we consider the general case when both the system and reservoir are finite. Thus, both the finite size and conservation law effects are present. Figure 3 shows examples of the subensemble particle number fluctuations calculated accordingly to the general Eqs.  Fig. 3) are considered. These two points are marked by crosses in Fig. 1. The choice of these specific points for the illustration is due to the following reasons. First, these two locations correspond to rather different GCE values for the fluctuation measures which are shown by full red circles in Fig. 3. The deviations from the ideal gas limit in both cases are large. Second, these two points in the n-T plane are in different proximities to the CP at ( n = 1, T = 1).
The thermodynamic limit results given by Eqs. (33) and (35) are represented by red lines in Fig. 3. The points x → 0 on these red lines correspond to the GCE values (13)- (15). They are shown in Fig. 3 by full red circles. The x-dependence according to Eqs. (33) and (35), shown by the red lines in Fig. 3, reflects the global N 0 conservation. The comparison of these lines with those in Fig. 2 shows a strong sensitivity of the skewness and kurtosis to the presence of interactions between particles. At finite x values the effects of the N 0 -conservation keep being significant even in the thermodynamic limit N 0 → ∞. At finite N 0 there are additional finite size effects. How large N 0 should be to approach the thermodynamic limit shown by red lines in Fig. 3 with a certain accuracy? This depends on both the proximity of the point ( n, T ) to the CP on the phase diagram and on the numerical value of x. The closer the system is to the CP, the larger are the finite-size effects. This evidently reflects the growth of the correlation length as one approaches the CP, which is known to become of a macroscopic magnitude at the CP.
The magnitude of the finite size effects at a fixed N 0 is minimal at x = 1/2, as seen from Fig. 3. This is because both volumes V = V 0 /2 and V 0 − V = V 0 /2 are relatively large in this case. Thus, to minimize the finite size effects in the event-by-event fluctuation data it may be worthwhile to aim for an acceptance which encompasses close to 50% of all final particles on average. The effects from the global charge conservation are not small in this case.
However, they can be estimated (and then corrected for) using the formulas of the subensemble acceptance procedure, Eqs. (33)- (35). It should be noted, however, that the skewness goes to zero at x = 1/2, as this quantity is an asymmetric function of x in the interval [0, 1]. It would therefore be necessary to consider acceptance away from x = 1/2 for this quantity.
The thermodynamic limit can be reached also at smaller x-values. The smaller x-values would however require the larger N 0 to reach the same level accuracy with respect to the finite-size effects. Let us consider for example the lines with N 0 = 400, shown in the right panel in Fig. 3 for the phase diagram point n = 1.5, T = 1.5. 3 To have ω, Sσ, and κσ 2 deviate from their thermodynamic limits by no more than 10% one has to take, respectively, x 0.05, x 0.10 and x 0.15. This numerical example as well as the general trend of the data presented in Fig. 3 demonstrate an important conclusion: For the same ( n, T ) point, the proximity to the thermodynamic limit behavior is different for different fluctuation measures. To reach the same proximity for the higher moments of particle number distribution one needs a larger system (larger N 0 ) and/or larger experimental acceptance (larger x).
The finite size effects become stronger in the vicinity of the CP. For example, the results at n = 1, T = 1.2 are shown in the left panel of Fig. 3. This n, T point is closer to the CP. To reach the proximity to the thermodynamic limits shown by red lines the values of N 0 and/or x must be higher than those in the right panel.
When x is close to 0 or 1 some "oscillations" in the xdependence are visible for moderate values of N 0 . This is connected to the excluded volume restrictions when only few finite-sized particles can fit in the volume V . We explore the finite size effects specifically in the next subsection.

C. Finite size effects
In this subsection we discuss the thermodynamic limit, V 0 → ∞, for finite values of volume V . This corresponds to x = V /V 0 → 0 as V 0 → ∞ simultaneously. In this case, the free energy, F = − T ln Z ce , of the reservoir with N 0 −N particles in the volume V 0 −V can be written as The partition function (19) can then be expressed as where (∂F/∂N 0 ) V0,T = µ 0 , is the chemical potential of the reservoir and N max = V /b . Equation (43) includes the finite size effects because of the finite value of N max . This finite size restriction is not very important at the regions of the phase diagram located far away from the CP. It is however crucial at the CP when the intensive fluctuation measures become divergent. In the thermodynamic limit V → ∞, Eq. (43) leads to N max → ∞ and the N -fluctuation measures in the subensemble approach their GCE values. Their behavior was discussed in Sec. II. In the vdW model, µ 0 is calculated as [44]: The probability distribution (20) is then calculated as where This agrees with the result of Ref. [27] where particle number distributions for finite vdW system were calculated in the (T ,µ) plane. The results for the scaled variance ω in the subensemble with the probability distribution (46) are presented in Fig. 4. The finite size effects disappear with increasing N max = V /b. To approach the GCE limit with the same accuracy larger V -values are needed the closer n and T are to their critical values. At the CP point n = 1 and T = 1 the scaled variance in the subensemble behaves as ω ∼ √ V and is divergent at V → ∞.  The "oscillatory" behavior of the fluctuations in the subensemble at small volumes V is observed in Fig. 3 at x ∼ 0 and x ∼ 1. This will be illustrated now on example of the scaled variance ω. When the volume V is so small that only one particle can fit in, b < V < 2b, the partition function (19) of the subensemble is a sum of only two terms with N = 0 and N = 1. In this case, one obtains N k = N for k = 1, 2 and, thus, At V /b → 1 + 0 one finds N → 0. Thus, ω[N ] → 1 which is in agreement with the Poisson distribution. We demonstrate the excluded volume threshold effects by depicting in Fig. 5 the scaled variance ω as a function of V /b in the region of small V /b. The vertical dashed lines show the thresholds of the system volume V /b at 1, 2, 3 and 4 particle level. One sees that the excluded volume threshold effects for ω[N ] are substantial at V /b 5 which corresponds to x 5(bn)/N 0 .

IV. SUMMARY
We investigated particle number fluctuations in an interacting thermal subsystem, taking into account effects associated with the global conservation of particle number (conserved charge) and finite system size. The total number of particles N 0 (total conserved charge) in the whole volume is fixed, in analogy to the final state (net) baryons in heavy-ion collisions, and treated in the canonical ensemble. The fluctuations of particle number N in a subvolume (acceptance) V < V 0 are described by a statistical ensemble which is distinct from both the canonical and grand-canonical ensembles.
The specific calculations have been performed for the van der Waals (vdW) equation of state, which contains a first-order phase transition and a critical point. The supercritical temperatures have been considered. Due to the universality of the critical behavior, we expect our results to reflect generic features of fluctuations near a critical point of a first-order phase transition in the presence of global charge conservation effects.
The global charge conservation influences the fluctuations at any finite value of the subvolume fraction x ≡ V /V 0 . In the thermodynamic limit, N 0 → ∞, these effects are in agreement with the recently developed subensemble acceptance procedure [20] and thus can be corrected for analytically.
In a more general case of a finite N 0 and finite x, both the finite size and global charge conservation effects simultaneously influence the fluctuation measures. The finite size effects at a fixed value of N 0 are the smallest at x = 1/2, where the two subsystems are both large. The magnitude of the finite size effects depends on the proximity of the critical point: the closer the system is to the critical point, the larger are the finite size effects. This can be understood due to the growth of the correlation length and, correspondingly, fluctuations in the vicinity of the critical point, which become comparable to the total system size.
Threshold effects are observed for very small volumes, V b, when only few finite-sized particles fit into the volume. An oscillatory behavior is observed, associated with the opening of new channels at the thresholds.
The following strategy may be adopted for extracting the GCE values of the cumulant ratios in relativistic heavy-ion collisions: 1) The behavior of ω and Sσ of the fluctuations of a conserved charge should be analyzed within several different acceptances (which corresponds to different x values). If the linear x-dependence of ω and Sσ is established, it can be considered as a signal of approaching the thermodynamic limit (see Fig. 3). Linear fits can then be performed to extract ω gce and Sσ gce .
2) The finite-size effects have a stronger influence on the kurtosis κσ 2 compared to ω and Sσ. As the finitesize effects are the smallest at x = 1/2, it is advisable to measure κσ 2 in an acceptance as close to x = 1/2 as possible. One can then extract κσ 2 gce from experimentally measured κσ 2 and the previously reconstructed Sσ gce using Eq. (35).
It should be noted that our analysis is based on an idealized picture of a homogeneous system in statistical equilibrium. The analysis does not incorporate the various dynamical effects present in relativistic heavyion collision experiments, as well as detector limitations. Moreover, measurements in heavy-ion experiments are performed in the momentum space rather than in the coordinate space. The degree of correlation between momenta and coordinates of particles at freeze-out depends on the collective flow, for example, the longitudinal flow. In future works we plan to include the influence of dynamical effects as well as to address the connection between the system's separation in coordinate space with the corresponding separation in the momentum space. We also plan to extend our approach to fully relativistic systems with multiple conserved charges, as is appropriate for relativistic heavy-ion collisions. der contract number DE-AC02-05CH11231. J.S. thanks the Samson AG and the BMBF through the ErUM-Data project for funding. This work was supported by the DAAD through a PPP exchange grant. Computational resources were provided by the Frankfurt Center for Sci-entific Computing (Goethe-HLR). H.St. acknowledges the support through the Judah M. Eisenberg Laureatus Chair by Goethe University and the Walter Greiner Gesellschaft, Frankfurt.