Particlization of an interacting hadron resonance gas with global conservation laws for event-by-event fluctuations in heavy-ion collisions

We revisit the problem of particlization of a QCD fluid into hadrons and resonances at the end of the fluid dynamical stage in relativistic heavy-ion collisions in a context of fluctuation measurements. The existing methods sample an ideal hadron resonance gas, therefore, they do not capture the non-Poissonian nature of the grand-canonical fluctuations, expected due to QCD dynamics such as the chiral transition or QCD critical point. We address the issue by partitioning the particlization hypersurface into locally grand-canonical fireballs populating the space-time rapidity axis that are constrained by global conservation laws. The procedure allows to quantify the effect of global conservation laws, volume fluctuations, thermal smearing and resonance decays on fluctuation measurements in various rapidity acceptances, and can be used in fluid dynamical simulations of heavy-ion collisions. As a first application, we study event-by-event fluctuations in heavy-ion collisions at the LHC using an excluded volume hadron resonance gas model matched to lattice QCD susceptibilities, with a focus on (pseudo)rapidity acceptance dependence of net baryon, net proton, and net charge cumulants. We point out large differences between net proton and net baryon cumulant ratios that make direct comparisons between the two unjustified. We observe that the existing experimental data on net-charge fluctuations at the LHC shows a strong suppression relative to a hadronic description.


I. INTRODUCTION
Event-by-event fluctuations in relativistic heavy-ion collisions have long been considered sensitive experimental probes of the QCD phase structure [1][2][3][4]. At the highest collision energies achievable at the LHC and RHIC they can be used to analyze the QCD chiral crossover transition at small baryon densities [5]. The equilibrium fluctuations of the QCD conserved charges in the grandcanonical ensemble have been computed at µ B = 0 from first principles, via lattice gauge theory simulations [6,7]. An appropriately performed comparison between experimental measurements and lattice QCD predictions can, in principle, establish whether a locally equilibrated QCD matter is indeed created in experiment. At lower collision energies, the fluctuations are used in the experimental search for the hypothetical QCD critical point and the first-order phase transition at finite baryon density. This is motivated by the fact that fluctuations, in particular the net proton cumulants of higher order, are increasingly sensitive to the proximity of the critical point [8,9]. The corresponding measurements are in the focus of several experimental programs, including beam energy scans performed at RHIC [10,11] and CERN-SPS [12]. The experimental data in the literature includes second order cumulants, both diagonal [13][14][15][16] and off-diagonal [17][18][19], as well as higher-order fluctuation measures [11,[20][21][22].
A proper theoretical modeling is crucial for interpreting the experimental data. It is not uncommon in the literature to directly compare the theoretical fluctuations evaluated in the grand-canonical ensemble with experimental measurements [23][24][25][26][27][28][29][30][31][32]. Such comparisons, however, have several important drawbacks. For one thing, the experimental measurements are performed in momentum space whereas the theoretical approaches operate in configuration space. Cuts in the momentum space may be identified with the coordinate space if strong space-momentum correlations are present, for instance due to Bjorken flow, but even in this case a degree of smearing will be present because of the thermal motion [33,34]. Event-by-event fluctuations, especially the high-order cumulants, are strongly affected by global conservation laws [35][36][37], requiring large corrections to the grand-canonical distributions. Other mechanisms include volume fluctuations [38][39][40], finite system size [41], as well as non-equilibrium dynamics such as memory effects [42] or hadronic phase evolution [43]. Proper modeling of these effects is thus required for analyzing the experimental data quantitatively.
The standard approach to describe the evolution of strongly interacting QCD matter created in heavy-ion collisions is relativistic fluid dynamics [44,45]. The hydrodynamic description terminates at a so-called particlization stage [46], where the QCD fluid is transformed into an expanding gas of hadrons and resonances. This picture forms the basis of the hybrid models of heavy-ion collisions [47,48] and it works quite well in describing the spectra and flow of bulk hadrons measured in a broad range of collision energies [49][50][51][52].
Event-by-event fluctuations of hadron yields, on the other hand, are seldom analyzed in the hydro picture. The yields of hadrons and resonances are usually sampled in each fluid element from a Poisson distribution. Because the Poisson distribution is additive, this means that the yields of all hadron species in the full space follow the Poisson distribution as well. This picture corresponds to the multiplicity distribution of an ideal arXiv:2012.09954v1 [hep-ph] 17 Dec 2020 Maxwell-Boltzmann hadron resonance gas (HRG) in the grand-canonical ensemble. Most hydro simulations use this type of sampling [53][54][55][56]. More advanced procedures incorporate exact conservation of the QCD conserved charges and/or energy-momentum [57][58][59][60][61], however, these procedures are still restricted to the equation of state of an ideal HRG. The existing methods, therefore, are not suitable to analyze the fluctuation signals of any effect that goes beyond the physics of an ideal hadron gas.
Interacting HRG models, on the other hand, offer more flexibility. For instance, an HRG model with excluded volume corrections can describe the lattice QCD cumulants of net baryon distribution in vicinity of the chemical freeze-out at µ B = 0 [62,63], which the ideal HRG model cannot. Another example is HRG model with van der Waals interactions, which captures the physics of nuclear liquid-gas transition at large µ B [31,64]. It is the purpose of this work to formulate a particlization routine appropriate to describe event-by-event fluctuations encoded in the equation of state of such interacting HRG models.
The paper is organized as follows. In Sec. II we introduce a method for sampling an interacting HRG at particlization stage of heavy-ion collisions that we call subensemble sampler. Sec. III describes the technical details of sampling an excluded volume HRG model that we study this work as an example. In Sec. IV the subensemble sampler is applied for the description of net baryon and net proton fluctuations in heavy-ion collisions at LHC energied. Discussion and summary in V close the article.

II. SUBENSEMBLE SAMPLER
Consider the particlization stage of heavy-ion collisions at the end of the ideal hydrodynamic evolution. This stage is characterized by a hypersurface σ(x), where the space-time coordinate x is taken in the Milne basis, x = (τ, r x , r y , η s ). Here τ = t 2 − r 2 z and η s = 1 2 ln t+rz t−rz are the longitudinal proper time and space-time rapidity, respectively, r x , r y , and r z are the Cartesian coordinates. The QCD matter is assumed in local thermodynamic equilibrium at each point x on this hypersurface. 1 As the fluid is converted into hadrons at this stage, the equation of state is described by hadron and resonance degrees of freedom, i.e. this has to be a variant of the hadron resonance gas model matched to the actual QCD equation of state at each point on the hypersurface. Let us denote Z HRG (T, V, µ) as the grand partition function of a hadron resonance gas at temperature T , volume V , and chemical potentials µ = (µ B , µ Q , µ S ), 1 In a more general case the deviations from local equilibrium are described using viscous corrections. and P hrg ({N i } f i=1 ; T, V, µ) as the corresponding multiplicity distribution for all hadron species. Here f is the number of different hadron species. In case of the commonly used ideal HRG model P hrg has a form of a multi-Poisson distribution where the Poisson means correspond to the mean multiplicities of primordial hadrons and resonances. Most particlization routines work with the multi-Poisson distribution of the ideal HRG model. However, P hrg will differ from the multi-Poisson distribution in a more general case of a non-ideal HRG. Thus, in the present work we generalize the particlization routine for arbitrary hadron multiplicity distributions.

A. Uniform fireball
Let us first consider a case of the grand-canonical ensemble, where the global conservations laws are enforced on average. Later on we will relax this assumption to incorporate exact global conservation.
If we further assume for the time being that the intensive thermal parameters T , µ B , µ Q , and µ S are the same across the entire fireball, and the partition function of the entire system coincides with the grand partition function Z HRG of a uniform HRG: Here with Z HRG (T, V, Q) being the canonical partition function of the HRG model with Q = (B, Q, S), and is the effective system volume at particlization. The single-particle momentum distribution function is given by the Cooper-Frye formula [65]: Here f i (x, p) is the single-particle distribution function. In the following we neglect quantum statistics and viscous corrections but take into account the possibility of interactions between hadrons. We assume that the distribution function takes the following general form 2 . . . 1 , 1 , 1 2 , 2 , 2 , , Figure 1. A schematic view of the partition of the space-time rapidity axis at particlization into N locally grand-canonical subvolumes, each characterized by values of the local temperature Tj, the chemical potential µj, and the volume Vj. Here is the flow velocity profile, and λ int i (T, µ) is a correction factor which describes deviations from the ideal gas distribution function induced by interactions. The explicit form of this factor depends on the interacting HRG model under consideration. The mean particle number N i is obtained by integrating Eq. (4) over the momenta: The full space hadron multiplicity distribution is given by the multiplicity distribution of the grand-canonical HRG:

B. Partition in rapidities
Let us now split the hypersurface into s slices along the space-time rapidity axis (see Fig. 1). The boundaries of each slices are η min j and η max j > η min j . Furthermore, one has η min j = η max j−1 for j > 1, and η min 1 = −η max and η max s = η max , where η max is the global maximum value of the space-time rapidity. One could, for instance, identify η max with the beam rapidity.
The subvolume characterizing the physical size of slice j is The key assumption in the following is that each subvolume V j is sufficiently large for it to be in the thermodynamic limit. Or in other words, V j ξ 3 for each i where ξ is any relevant correlation length. If that is the case, one can neglect the surface effects, namely the interactions between particles from different subvolumes. Mathematically speaking, this implies a scaling Z HRG (T, V j , µ) ∼ e Vj [or, equivalently, ln Z HRG (T, V j , µ) ∼ V j ] for V j ξ 3 . Also, the total partition function factorizes into a product of partition functions for each of the subvolumes: The form of Eq. (11) allows us to relax the assumption of the constancy of thermal parameters. Let us now assume that the intensive thermal parameters depend on the space-time rapidity η s . This implies that each of the rapidity slices is characterized by its own set of values of the thermal parameters, i.e. in Eqs. (10), (11) one has T → T i and µ → µ i : Let us denote the hadron multiplicities in a subvol- The multiplicity distribution N j is given by the corresponding multiplicity distribution of the HRG model with thermal parameters of the given subvolume, i.e. P gce N j = P hrg (N j ; T j , V j , µ j ). Due to the fact that we neglected all correlations between particles from the different subvolumes, the multiplicity distribution ofN j is independent of the multiplicity distributions in all other subvolumes. The probability distribution for multiplicities N j s j=1 across all subvolumes thus factorizes as follows: The factorization in Eq. (14) will no longer hold once we introduce exact global conservation of conserved charges. The momentum distribution of hadron species i emitted from a rapidity slice j reads with Here µ j,i = b i µ B,j + q i µ Q,j + s i µ S,j .

C. Exact global conservation laws
Let us now incorporate the effect of exact global conservation of conserved charges. As we work in the thermodynamic limit, V j ξ 3 , the exact conservation will not affect the mean multiplicities due to the thermodynamic equivalence of statistical ensembles. However, as the thermodynamic equivalence does not extend to fluctuations, the fluctuation observables will be affected by the exact conservation, no matter how large the system is.
The total values of the globally conserved baryon number, electric charge, and strangeness coincide with the GCE mean values due to the thermodynamic equivalence of ensembles: To enforce the global conservation laws on the level of multiplicity distributions one has to project out all microstates that violate the global conservations laws from the grand-canonical partition function. This is achieved by introducing a Kronecker delta into the grand partition function (12) of the entire system: The presence of the delta function in Eq. (18) breaks the factorization of multiplicity distributions in different rapidity slices. The joint multiplicity distribution reads Here q i = (b i , q i , s i ) is a vector of conserved charge values carried by hadron species j.

D. Sampling the multiplicity distribution
Here we present a general method for sampling the joint multiplicity distribution Eq. (19) of hadron numbers in all the subsystems. The method is based on rejection sampling and it assumes that it is known how to sample the multiplicity distribution of the grand-canonical variant of the HRG model used. To generate a configuration from the distribution (19) 1. SampleN j for j = 1 . . . s independently for each subsystem from the grand-canonical variant of an interacting HRG model characterizing each subsystem.

Compute
s k=1 Q k via Eq. (20). Accept the configuration if Q tot = s k=1 Q k , or go back to step 1 otherwise.
The method is general in the sense that it does not assume anything about the specific HRG model used. It will work both for an ideal and interacting HRG. It should be noted, however, that the algorithm may become inefficient if the acceptance rate in step 2 becomes low. This can happen for large systems and multiple conserved charges. More efficient algorithms can be devised for specific versions of the HRG model, see e.g. a multistep method of Becattini and Ferroni in Ref. [58]. We do employ this method in our Monte Carlo simulations in Sec. IV.

E. Thermal smearing
The algorithm in the previous section allows to sample hadron multiplicity distributions differentially in spacetime rapidity. The experiments, however, perform measurements in momentum rather than coordinate space, therefore, a transition to momentum space is necessary. In some cases, such as the Bjorken flow scenario at the highest collision energies, it is possible to identify the space-time rapidity η s with the momentum rapidity Y , allowing to study rapidity-dependent hadron distributions without the transition to the momentum space. Even in this case, however, a degree of smearing between η s and Y is present due to thermal motion. The boost invariance breaks down at lower collision energies and the problem of space-momentum correlations becomes even more severe. For these reasons it is necessary to assign each of the hadrons a 3-momentum. Furthermore, if a subsequent afterburner stage is to be included into the modeling, one has to generate both the spatial and momentum coordinates for each hadron.
The procedure to generate the momenta of all the hadrons is fairly straightforward. Once the multiplicity distributions N j s j=1 for all the rapidity slices have been sampled, the coordinates and momenta of all the hadrons can be generated through the standard Cooper-Frye momentum sampling, applied independently to each hadron in each of the rapidity slices. Several implementations for this task are available, see e.g. [53,54,66]. The sampled hadrons should then be provided as input into a hadronic afterburner like UrQMD [67,68] or SMASH [69], if one is used, or a cascade of resonance decays performed to obtain the final state particles that are measured experimentally. The comparison with data can then be done in the standard way, by computing the observables in a given acceptance as statistical averages.

III. EXCLUDED VOLUME MODEL FOR NET BARYON FLUCTUATIONS
To illustrate the developed formalism we shall apply it to net proton and net baryon fluctuations in heavyion collisions at energies reachable at LHC and RHIC. In this section we describe the motivation and the technical details behind an excluded volume HRG model that we use for the analysis. A reader interested only in the final heavy-ion results may skip to Sec. IV where these are presented and discussed.
The typical chemical freeze-out temperatures, T ch ∼ 155 − 160 MeV at the LHC [70][71][72] and T ch ∼ 160 − 165 MeV at the top RHIC energies [73], are close to the pseudo-critical temperature of the QCD crossover transition determined by lattice QCD T pc 155 − 160 MeV [74,75] at µ B = 0. Lattice QCD predicts that the high-order net baryon cumulants, namely the kurtosis χ B 4 /χ B 2 and the hyperkurtosis χ B 6 /χ B 2 ratios deviate significantly from the Skellam distribution baseline of the ideal HRG model, where these ratios are equal to unity. The hyperkurtosis in particular turns negative around T pc which is thought to be related to the remnants of the chiral criticality [5] at vanishing light quark masses. It would certainly be of great interest to verify this theory prediction of a negative χ B 6 experimentally, which may serve as an experimental evidence for the chiral crossover transition. The measurement of higher-order net proton fluctuations is planned in future runs at the LHC [76].
In our previous work [77] we studied this question analytically, in the framework of the subensemble acceptance method (SAM). There, the sensitivity of measurements to the equation of state was predicted to be not overshadowed if the measurements are performed in acceptance spanning 1-2 units of rapidity. However, the entire argument in [77] has been done in the configuration space, relying on perfect momentum-space correlations due to Bjorken flow. Here we would like to determine how the results will be distorted by the thermal smearing and resonance decays.
To apply the formalism of Sec. II we need to employ an interacting HRG model that matches the lattice QCD equation of state and be able to sample the grandcanonical multiplicity distribution of such a model. Here we take an HRG model with excluded volume interactions in the baryonic sector -the EV-HRG model -which was formulated in Refs. [62,64] and shown to describe well the lattice data on the diagonal net-baryon susceptibilities at µ B = 0 at temperatures up to and even slightly above T pc .

A. Single-component EV model
Before discussing the full model let us first consider a single-component excluded volume model in order to introduce the multiplicity sampling procedure. The grand partition function at fixed temperature T , volume V , and chemical potential µ reads Here is an ideal gas density of particle species with degeneracy d and mass m at vanishing chemical potential. K 2 is the modified Bessel function of the second kind. Equation (21) defines the multiplicity distribution of the EV model, giving the following (unnormalized) probability functioñ In the thermodynamic limit, N → ∞, the particle density n ev (T, µ) = N ev /V is determined by the maximum term in Eq. (21). MaximizingP ev with respect to N gives a transcendental equation defining n ev (T, µ): The solution to Eq. (24) is given in terms of the Lambert W function (see Ref. [78] for details): or The pressure reads

Dimensionless form
In the EV model it is possible to replace the three thermal parameters (T, V, µ) and the excluded volume parameter b by two dimensionless quantities, namely a reduced volumeṼ ≡ V /b and a parameter κ ≡ b φ(T ) e µ/T that characterizes the strength of repulsive interactions. The probability distribution (23) then takes the form The mean particle number reads This reduced form implies that the multiplicity distribution is fully specified if the values of parametersṼ and κ are known.

Cumulants of particle number distribution
Cumulants of the particle number distribution in the EV model can be evaluated from the probability distribution function (28). The nth moment reads The sums over N are finite due to the presence of the θ function in Eq. (28). Thus, for finiteṼ , they can be carried out explicitly. The cumulants can be expressed in terms of the moments as Here B n,k are the partial Bell polynomials. Explicit expressions for κ n [N ] can be obtained in the thermodynamic limit,Ṽ → ∞. This is achieved through the cumulant generating function The t-dependent mean value N ev (t) is obtained from Eq. (29) by a substitution κ → κ e t : Equation (33) corresponds to the first cumulant. The higher-order cumulants are obtained by differentiating N (t) ev with respect to t. The results up to fourth order read It follows that all cumulant ratios in the EV model depend exclusively on the value of a single parameter κ in the thermodynamic limit.
B. Sampling the excluded volume model To sample particle numbers from the probability distribution (23) of the EV model we will use a rejection sampling technique. First we sample N from an auxiliary envelope distributionP aux (N ; T, V, µ), which we take to be a Poisson distribution centered around N ev : Here N ev ≡ n ev (T, µ) V with n ev defined by Eq. (26). The theta function ensures that the packing limit is not violated, i.e. if for a value N sampled from the Poisson distribution one has V − bN < 0, this value is rejected.
To correct for the difference betweenP aux (N ) and P ev (N ) we apply rejection sampling for each value of N sampled fromP aux . First, we define a weight factor w(N ) as the ratio between the true and auxiliary multiplicity distributions: Here n ≡ N/V . The number N sampled fromP aux (N ) shall be accepted if η < w(N )/w max where w max is the maximum possible value of w(N ) and η is a random number uniformly distributed in an interval [0, 1].
To determine w max ≡ w(N max ) let us rewrite Eq. (39) as where we used Eq. (24). N max is determined from ∂ w(N )/∂N = 0. One obtains an equation The solution to the above equation is n max = n ev , i.e. the weight is maximized at the mean value of N in the thermodynamic limit 3 . w max reads In numerical calculations it is more convenient to work directly with normalized weights: The reduced weight in terms dimensionless variablesṼ and κ reads where N ev is given by Eq. (29). The sampling procedure described here is similar to the Monte Carlo EV model analysis performed in Ref. [79], with one distinction. In Ref. [79] an importance sampling technique was employed, where each generated event was accepted with a weightw(N ). Here, instead, all accepted events have the same weight, but their sampling involves an additional rejection step with respect to the weights w(N ).

Testing the sampling procedure
To test the sampling procedure described above we take κ = 0.03 and perform Monte Carlo sampling for different values ofṼ . The choice κ = 0.03 is motivated by the fact that this value is obtained in the EV-HRG model with baryonic excluded volume b = 1 fm 3 at T = 160 MeV and µ B = 0 [62], therefore, the exercise approximately corresponds to sampling the baryon multiplicity distribution in the vicinity of the QCD chiral crossover transition where the EV-HRG model approximates well the QCD cumulants of the net baryon distribution.
We sample 10 8 numbers at each value ofṼ and calculate cumulants of the resulting particle number distribution up to κ 6 . Figure 2 depicts the resultingṼdependence of the scaled variance κ 2 /κ 1 , kurtosis κ 4 /κ 2 , and hyperkurtosis κ 6 /κ 2 (symbols). The solid lines correspond to an analytic calculation of these ratios via a direct summation over all probabilities [Eqs. (30), (31)]. The Monte Carlo calculations agree with the analytic expectations at all studied values ofṼ , validating the sampling method. Figure 2 allows also to establish when the condition V ξ 3 is reached. This is signalled by the approach of the cumulant ratios to their expected values in the thermodynamic limit [Eqs. (34)- (37)], shown in Fig. 2 by the horizontal dashed lines. Cumulant of a higher order generally requires larger values ofṼ to reach the thermodynamic limit, reflecting the fact that higher cumulants are more sensitive to the correlation length ξ. We observe that cumulant ratios up sixth order are within few percent or less of the thermodynamic limit forṼ 20. The cumulants then scale linearly with the volume for larger values ofṼ . The valueṼ 20 thus establishes a lower bound on the physical volume of a single rapidity slice for the subensemble sampler in Sec. II to be applicable.

C. EV-HRG model
Having established the baryon multiplicity sampling procedure in a single-component case, we now turn to the full model. Quantitative applications to heavy-ion fluctuation observables require an equation of state with hadron and resonance degrees of freedom matched to first-principle lattice QCD equation of state. For the purposes of net baryon and net proton fluctuations studied here we employ a variant of an excluded volume hadron resonance gas (EV-HRG) model introduced in Refs. [64,80]. The repulsive EV interactions are introduced for all baryon-baryon and antibaryon-antibaryon pairs in the EV-HRG model, with a common value b of the EV parameter for all these pairs.
The pressure in the EV-HRG model is partitioned into a sum of meson, baryon and antibaryon contributions Here n id M and n id B(B) correspond to cumulative number densities of mesons and (anti)baryons in the ideal HRG limit (b → 0): Here µ i = q i · µ is the chemical potential of particle species i. The expression (47) can be rewritten in terms of the Lambert W function in close to analogy to Eq. (23) of the single-component EV model: The particle number densities of individual hadrons species are calculated as derivatives of the pressure with respect to the corresponding chemical potential n ev i = ∂p/∂µ i . The mean multiplicities N i ev ≡ V n ev i in the grand-canonical EV-HRG model read Here The mean multiplicities of mesons coincide with the ideal HRG model baseline. The multiplicities of (anti)baryons, on the other hand, are suppressed relative to ideal HRG due to EV interactions. This is quantified by a factor in the r.h.s of Eq. (52). For κ B 0.03, a value corresponding to T = 160 MeV and µ = 0 (see below), the yields of baryons are suppressed by about 5%.
Equations (51) and (52) define the factor λ int i (T, µ) entering the single-particle distribution functions f i (x, p) for particle species i in the Cooper-Frye formula, Eqs. (4) and (15). For mesons, i ∈ M , one has λ i (T, µ) = 1. For (anti)baryons The EV-HRG model has been studied in Refs. [62,63] in the context of lattice QCD results on diagonal net baryon susceptibilities and Fourier coefficients of net baryon density at imaginary chemical potentials. Reasonable description of these observables at temperatures close to T pc has been obtained for b = 1 fm 3 , corresponding to κ B 0.03. We employ this value of b in the present analysis. Figure 3 depicts the temperature dependence of kurtosis κ 4 /κ 2 and hyperkurtosis κ 6 /κ 2 of net baryon fluctuations at vanishing temperatures. The calculations are compared with the lattice data of Wuppertal-Budapest (blue bands) [81] and HotQCD (green bands and symbols) [82] collaborations. The model is in quantitative agreement with the lattice data for these two quantities up to T 180 MeV. This implies that netbaryon distribution of the EV-HRG model in this temperature range closely resembles that of QCD, at least on the level of sixth leading cumulants. And while this does not necessarily imply that EV interactions is the correct physical mechanism behind the behavior of net baryon susceptibilities, we view the EV-HRG model to be an appropriate tool for the purpose of analysis net baryon and net proton cumulants in heavy-ion collisions.
The sampling procedure in Sec. III B can be generalized for the EV-HRG model that has multiple hadron components. We note that the system in the EV-HRG model is partitioned into three independent subsystems, mesons, baryons, and antibaryons, see Eq. (45). Therefore, the sampling of the grand-canonical multiplicities proceeds independently for each of the three subsystems. The multiplicities of the non-interacting mesons are sampled from the Poisson distribution, in the same manner as in the ideal HRG. The joint probability distribution of numbers of all the baryon species, on the other hand, reads where, as before,Ṽ = V /b and The auxiliary envelope distribution for the sampling is a cut multi-Poisson distribution: The theta function is introduced to avoid exceeding the packing limit.
Finally, the normalized weight for the rejection sam- µ B = 0 Figure 3. Temperature dependence of net baryon susceptibility ratios χ B 4 /χ B 2 (left) and χ B 6 /χ B 2 (right) evaluated at µB = 0 in the EV-HRG model. The blue and green bands/symbols depicts lattice QCD results of Wuppertal-Budapest [81] and HotQCD [82] collaborations, respectively. pling step reads Here n i ≡ N i /V and n ev i ≡ N i ev /V . The algorithm for sampling the multiplicity distribution of the EV-HRG model is the following: 1. Sample the multiplicities {N i } of all baryons from the cut multi-Poisson distribution (57).
2. Generate a random number η from the uniform distribution on the unit interval (0, 1). If η < w B ({N i }), go to the next step. Otherwise, go back to step 1.
3. Repeat steps 1-2 in the same fashion to sample the multiplicities of antibaryons.
4. Sample multiplicities of mesons from the multi-Poisson distribution of the ideal HRG model.
The procedure for generating the multiplicity distribution in the EV-HRG model in various rapidity slices that are constrained by global conservation of conserved charges, as described in Secs. II and III, is implemented in an extended version of the Thermal-FIST package [83]. We use this package in all our calculations.

IV. FLUCTUATIONS IN HEAVY-ION COLLISIONS AT LHC ENERGIES
A. The setup We apply our formalism to study the rapidity acceptance dependence of fluctuation observables in heavy-ion collisions. To proceed we need to specify the partition of the space-time rapidity η s axis into fireballs as well as the η s dependence of thermal parameters and volume.
Let us consider Pb-Pb collisions at √ s NN = 2.76 TeV.
At midrapidity the chemical freeze-out is characterized by vanishing chemical potentials, temperature values T 155 − 160 MeV and freeze-out volume per rapidity unit dV /dy ∼ 4000 − 5000 fm 3 [70,84]. The simplest possibility then is to assume boost invariance across the entire space-time rapidity range. In this scenario, the mean total number of particles of given kind in full space, say charged multiplicity N 4π ch or number of (anti)baryons N 4π B(B) , is then simply given by multiplying the rapidity density at y = 0 by the total (space-time) rapidity coverage −η max for the charged multiplicity. The question is how to determine η max s . One possibility is to equate this quantity to the beam rapid- However, such an estimate is too crude and will overestimate the actual N 4π ch . The rapidity density of charged multiplicity measured at √ s NN = 2.76 TeV by the ALICE collaboration [85] is consistent with a Bjorken plateau only in a rapidity range |y| 2, whereas at higher rapidities dN ch /dy drops. The entire measured rapidity dependence of dN ch /dy is described well by Gaussian with a width σ = 3.86 ± 0.05 [85]. We can use this fact to relate N 4π ch and dN ch /dy| y=0 in an empirical way: Here the error is due to the uncertainty in the value of σ. Comparing Eq. (60) with (59) one obtains We shall use a value η max s = 4.8 for Pb-Pb collisions at √ s NN = 2.76 TeV in the following. We take T = 160 MeV and µ = 0 for all rapidities. With this choice the model accurately reproduces the rapidity densities of various hadron species at |y| 2, where the Bjorken plateau is observed in the data [85], and provides an accurate estimate of the total hadron multiplicities in full phase space. As we have assumed boost invariance across the entire space-time rapidity range, the model does not describe rapidity distributions at |y| 2 and thus should not be applied to calculate observables at large rapidities. However, given the fact that the model does reproduce the 4π charged multiplicity, it is suitable to describe the influence of global conservation laws on observables computed around midrapidity, |y| 2. In the following we focus on these regions around midrapidity. In a more general study the assumption of boost invariance can be relaxed to incorporate a more accurate description of the forward-backward rapidity regions.
Our model yields a vanishing total net baryon number in the full space. Essentially, this means that we neglect baryons from the fragmentation regions. This is similar to a recent study [86] performed in the framework of the ideal HRG model. There it was estimated that the effect of fragmentation baryons at the LHC does not exceed 6% for the sixth order net proton cumulant. We therefore expect the possible influence of the fragmentation region baryons on our results to be small.
We partition the space-time rapidity axis uniformly into slices of width ∆η s = 0.1. With η max s = 4.8 this implies a total of 96 slices. The volume of a single slice in 5% most central collisions is V i = dV /dy ∆η s 400 fm 3 . This value is sufficiently large to ensure that the thermodynamic limit is reached in each of the subvolumes and thus the requirements for the validity of the sampling procedure described in Sec. II satisfied. This also implies that all intensive quantities, such as cumulant ratios, are independent of the value of V i in this regime, i.e. V i can be scaled up and down as long as V i ξ 3 . This feature is very useful for the Monte Carlo sampling procedure. Indeed, as the statistical error in higher-order cumulants increases with the volume, this error can be minimized by choosing the volume as small as possible. According to Fig. 2, a value V i = 20 fm 3 is sufficiently large to capture all the relevant physics for cumulants up to sixth order. For this reason we take V i = 20 fm 3 in our Monte Carlo simulations and then linearly scale up the resulting cumulants to match the volume V i = 400 fm 3 in 0-5% Pb-Pb collisions.
We take the EV-HRG model with b = 1 fm 3 . As discussed in Sec. III B, this model provides a reasonable description of high-order net baryon susceptibilities from lattice QCD. The grand-canonical distribution of hadron multiplicities can be efficiently sampled following the rejection sampling based algorithm described in Sec. III B. We take T i = 160 MeV and vanishing chemical potentials, µ i = 0, uniformly for all subvolumes along the rapidity axis.
For the net baryon cumulants we shall take into account only the exact conservation of baryon number, which is exactly vanishing, B = 0, in all events. In principle, one should also take into account the exact conservation of electric charge and strangeness. However, as discussed in Refs. [77,87], the influence of these conserved charges on net baryon cumulants is negligible at LHC energies. Their influence on net proton cumulants is more sizable [87] but still expected to be subleading compared to baryon number conservation. Neglecting the exact conservation of electric charge and strangeness allows to significantly speed up the Monte Carlo event generation, as this strongly reduces the rejection rate associated with exact conservation of multiple conserved charges and allows to gather enough statistics within a reasonable time period to accurately evaluate cumulants up to sixth order. We do analyze the influence of electric charge and strangeness conservations on 2nd order cumulants of various net-particle distributions in Sec. IV F Once the joint hadron multiplicity distribution from all the subvolumes has been sampled, we generate the hadron momenta, independently for each hadron. To that end we employ the blast-wave model [88], which provides a reasonable description of bulk particle's p T spectra at LHC [89]. The model corresponds to a particlization of a cylindrically shaped fireball ( r 2 x + r 2 y ≤ r max ), at a constant value of the longitudinal proper time τ = τ 0 . The longitudinal collective motion obeys the Bjorken scaling while the radial velocity scales with the transverse radius, β r ∝ r n ⊥ . This corresponds to a flow profile u µ (x) = (cosh ρ cosh η s , sinh ρ e ⊥ , cosh ρ sinh η s ), where ρ = tanh −1 β r , and β r = β s ζ n is the transverse flow velocity profile. Here ζ ≡ r ⊥ /r max is a normalized transverse radius. The momentum distribution of hadron species with mass m emerging from a jth space-time rapidity subvolume is given by Here m T = p 2 T + m 2 is the transverse mass, y = 1 2 log ω p − p z ω p + p z is the longitudinal rapidity, and I 0 is a modified Bessel function.
The sampling of momenta from the distribution (62) is readily implemented in the Thermal-FIST package that we employ. We are only left with specifying the values of the blast-wave model parameters β s and n. For this purpose we make use of the result of a recent study [90], where the blast-wave model was fitted to experimental data of the ALICE collaboration with account for modification of p T spectra due to resonance decays. For 5% most central Pb-Pb collisions at √ s NN = 2.76 TeV one has β s = 0.77 and n = 0.36, which gives a reasonable description of bulk hadron p T spectra 4 . One should note that Ref. [90] has extracted a temperature value of T = 149 MeV from the p T spectra fits rather rather than the T = 160 MeV value that we use here for fluctuations. However, the T = 160 MeV value shows a similarly good agreement of the blastwave model proton p T spectrum with the data, as the one shown in [90] for T = 149 MeV. Figure 4 compares the shape of the p T spectrum of protons as observed in the data (red symbols) [89] and predicted by the blastwave model [Eq. (62)] with T = 160 MeV, β s = 0.77, and n = 0.36. The dashed line in Fig. 4 corresponds to blastwave model spectrum which includes the modification of the proton p T spectrum due to resonance decays. This effect, computed here via Monte Carlo simulations of decays, only slightly modifies the momentum distribution.
In the final step of the Monte Carlo event generation procedure we perform all strong and electromagnetic decays until only stable hadrons are left. We generate 10 10 events in total 5 and study the rapidity dependence of various fluctuation observables. As our analysis only concerns the baryons, to speed-up the Monte Carlo procedure we omit all the primordial mesonic species (step 4 in the algorithm of Sec. III B), as these do not affect the behavior of (anti)baryons in any way within the EV-HRG model that we use.

B. Rapidity acceptance dependence of net baryon cumulants
We start with the rapidity acceptance dependence of net baryon number cumulants. First, we look at the second cumulant of net baryon fluctuations normalized by the Skellam distribution baseline, κ 2 [B −B]/ B +B . This type of ratio has been extensively studied at LHC energies by the ALICE collaboration [15] for net protons. This ratio equals unity for the case of a grand-canonical ideal HRG model at any temperature and chemical potentials. The ratio, however, does exhibit small deviations from unity in the EV-HRG model that we use. For instance, at T = 160 MeV and µ = 0 the grand-canonical value reads We note that it is currently challenging to directly compute κ 2 [B −B]/ B +B in lattice QCD, as the denominator B +B is not a conserved quantity. Given the good agreement of the EV-HRG model with lattice QCD for the higher-order cumulants, however, we expect QCD to have a similar value to the one given by Eq. (63). An interesting question now is to determine if and how the grand-canonical value in Eq. (63) is reflected in heavy-ion data.
The top panel of Fig. 5 depicts the rapidity acceptance ∆Y acc dependence of κ 2 [B−B]/ B+B that results from the Monte Carlo sampling within the subensemble sampler. Here the acceptance is centered at midrapidity, i.e. particles with rapidity |y| < ∆Y acc /2 are accepted. The red symbols depict the full result which includes the distortion of hadron momenta due to thermal smearing at particlization and subsequent resonance decays. The black symbols, on the other hand, correspond to the case when these effects are neglected, i.e. the final kinematical rapidity is taken to be equal to the space-time rapidity coordinate at particlization, y ≡ η s . Comparing the two allows to establish the effect of thermal smearing and resonance decays. We observe that the Monte Carlo results (top), κ4/κ2 (middle), and κ6/κ2 (bottom) of net baryon distribution in 0-5% central Pb-Pb collisions at the LHC in an excluded volume HRG model matched to lattice QCD. The symbols depict the results of the Monte Carlo event generator, the full black squares correspond to neglecting the momentum smearing, the open red triangles include the thermal smearing at particlization, and the full red circles incorporate the smearing due to both the thermal motion and resonance decays. The dashed black lines correspond to the predictions of the SAM framework [77]. The solid red lines correspond to adding a Gaussian rapidity smearing on top of the SAM. The dashed blue lines correspond to the binomial acceptance, which describes the effects of baryon number conservation in the ideal HRG model limit.
in the no-smearing case agree with the analytic expectations of the SAM (black lines). The SAM baseline for κ 2 [B −B]/ B +B is given by [77,87] Here α is a fraction of the total volume which corresponds to the acceptance |η S | < ∆Y acc /2: The agreement of the Monte Carlo points with the SAM is the expected result and serves as a validation of the sampling procedure. Notable differences between the red (momentum rapidity) and black (space-time rapidity) points in Fig. 5 appear when the acceptance is sufficiently small, ∆Y acc 1. This is a consequence of the dilution of momentumspace correlations due to thermal motion.
For a very small acceptance, ∆Y acc 1, the results converge to the baseline given by the binomial distribution, The additional momentum smearing due to decays of baryonic resonances is virtually negligible, being completely overshadowed by the thermal smearing. This is true not only for the variance, but also for the kurtosis and hyperkurtosis, as seen by comparing the red points (thermal smearing + resonance decays) with the open red triangles (thermal smearing only) in all three panels of Fig. 5. To understand this behavior one can consider e.g. decays ∆ → N π. In such a decay the released momentum is split evenly between the two decay products in the resonance center-of-mass frame. This leads to a larger velocity (rapidity) smearing of the lighter decay product -the pion -whereas the velocity (rapidity) of nucleon is less affected. We conclude that the smearing of baryon fluctuations due to resonance decays can be safely neglected. Note that this statement does not extend to (net-)particle fluctuations involving lighter hadrons such as pions or kaons. There the effect of resonance decays should be carefully taken into account.
In the Appendix we develop a simplified analytic model to take into account the momentum smearing in net baryon cumulants. There we assume that the shift in kinematical rapidity relative to the space-time rapidity is described for all baryon species by a Gaussian smearing. The red lines in Fig. 5 exhibit the results of such a simplified calculation. For a Gaussian width of σ y = 0.3 the simplified calculations agree very well with the full Monte Carlo results. Therefore, this model can be used to predict the rapidity dependence of p T -integrated net baryon cumulants without invoking the time-consuming Monte Carlo event generator.
Net-baryon fluctuations in a sufficiently large rapidity acceptance ∆Y acc 1 are accurately described by the analytical SAM baseline (64). This conclusion is important, because the SAM makes the connection between the grand-canonical susceptibilities and cumulants constrained by global conservation laws without any additional assumptions regarding the equation of state. In our previous work [77] where the SAM is introduced, we argued that the SAM is reliable for rapidity acceptances ∆Y acc 1, where the distortion due to thermal smearing is expected to be subleading. The results obtained in the present work using the EV-HRG model explicitly confirm this. The grand-canonical κ 2 [B −B]/ B +B ev,gce can therefore be extracted from data by fitting the αdependence of net-baryon fluctuations measured in sufficiently large rapidity acceptance via Eq. (64).
We turn now to the kurtosis of net baryon fluctuations, Here β ≡ 1 − α. At LHC energies one has χ B 3 /χ B 2 = 0, thus, the second term in Eq. (67) does not contribute.
With thermal smearing and resonance decays included, the kurtosis deviates from the SAM baseline for ∆Y acc 1 and for ∆Y acc 1 tends to the binomial distribution baseline, which at the LHC energies reads . Lattice QCD predicts a sign change of the grand-canonical hyperkurtosis at µ = 0 in the vicinity of the pseudocritical temperature (Fig. 3). This qualitative feature is thought to be a signature of the QCD chiral crossover transition [39]. Therefore, a corresponding measurement of κ 6 [B −B]/κ 2 [B −B] in heavyion collisions at the LHC can potentially serve as the first experimental signature of that transition. The EV-HRG model reproduces the available lattice QCD data for χ B 6 /χ B 2 and gives the following value at T = 160 MeV: This agrees within errors with the continuum estimate of the Wuppertal-Budapest collaboration, χ B 6 /χ B 2 = −0.26 ± 0.17 [81] as well as with N τ = 8 results of the HotQCD collaboration [82] shown in Fig. 3.
The lower panel of Fig. 5 shows the rapidity acceptance dependence of the hyperkurtosis. In the absence of momentum smearing, the Monte Carlo results are described by the analytical SAM baseline, which for LHC energies, i.e. for µ = 0, reads [77] The hyperkurtosis, in the absence of momentum smearing, is sensitive to the grand-canonical value (68) in acceptances up to ∆Y acc 1.5. For larger acceptances baryon conservation dominates, making it difficult to disentangle between the EV-HRG model and the binomial baseline, given by κ 6 3αβ). This was already pointed out in our previous study [77]. The thermal smearing distorts the signal at small acceptances, ∆Y acc 0.5, where the hyperkurtosis is closer to the binomial distribution baseline than it is to the SAM. At 0.5 ∆Y acc 1.5, on the other hand, κ 6 [B −B]/κ 2 [B −B] is overshadowed neither by the thermal smearing nor by the baryon number conservation. We, therefore, argue that a measurement of a hyperkurtosis, which is negative over this entire range may be interpreted as a signal of the chiral crossover 6 .

C. Net baryon vs net proton fluctuations
Our discussion has so far been restricted to cumulants of net baryon distribution. However, experiments typically cannot measure all baryons, in particular the measurement of neutrons is extremely challenging. For this reason one usually uses net protons as a proxy for net baryons. It is natural to expect net protons to carry at least some information about net baryon fluctuations. In fact, as shown by Kitazawa and Asakawa [92,93], under  Rapidity acceptance dependence of net baryon (black squares) and net proton (blue circles) cumulant ratios κ2/κ Skellam 2 (top), κ4/κ2 (middle), and κ6/κ2 (bottom) in 0-5% central Pb-Pb collisions at the LHC in an excluded volume HRG model matched to lattice QCD. The open blue diamonds correspond to net proton cumulants evaluated from net baryon cumulants using a binomial-like method of Kitazawa and Asakawa [92,93]. The black lines correspond to the analytical predictions of the SAM framework with (solid) and without (dashed) Gaussian rapidity smearing.
the assumption of isospin randomization at late stages of heavy-ion collisions, one can reconstruct the cumulants net baryon distribution from the measured factorial moments of proton and antiproton distributions.
However, these considerations do not imply that ratios of proton cumulants can be used directly in place of the corresponding ratios of baryon cumulants, something which has nevertheless been employed in a number of works in the literature [28,29,94]. The proton and baryon cumulant ratios do coincide in the free hadron gas limit, where they both trivially reduce to the Skellam baseline, but this does not hold in general case.
Large differences between net proton and net baryon cumulant ratios were reported earlier in Ref. [31] for the van der Waals HRG model in the grand-canonical ensemble. Here we study these differences in the framework of the EV-HRG model constrained to lattice data and include effects of global baryon conservation and momentum smearing. Figure 6 depicts the rapidity acceptance dependence of net baryon (black squares) and net proton (blue symbols) cumulant ratios κ 2 /κ Skellam 2 , κ 4 /κ 2 , and κ 6 /κ 2 calculated using Monte Carlo sampling within the SAM. The calculations incorporate the thermal smearing and resonance decays. The results reveal large differences between net proton and net baryon cumulants ratios. Net proton cumulant ratios are considerably closer to the Skellam baseline of unity. This can be understood in the following way. By taking only a subset of baryons -the protons -one dilutes the total signal due to baryon correlations. This leads to a smaller deviation of cumulants from Poisson statistics -the limiting case of vanishing correlations.
The large difference between net proton and net baryon cumulants clearly indicates that direct comparison between the two is not justified. It is interesting that net proton cumulant ratios cross the grand-canonical value of the corresponding net baryon ratios in the grandcanonical limit (horizonal lines in Fig. 6) for a sufficiently large acceptance. This, for instance, takes place at ∆Y acc 1.4 for κ 2 /κ Skellam 2 while for κ 4 /κ 2 the crossing is at ∆Y acc 2.5. The crossings take place due to suppression of net proton cumulants from baryon number conservation. This accidental coincidence between net proton and grand-canonical net baryon cumulant ratios may be of relevance for the recent analysis of STAR data by the HotQCD collaboration in Ref. [94]. There, the net baryon lattice QCD susceptibilities were directly compared to the measured net proton cumulants and an agreement, within large error bars, was reported.
We explore also, if the method of Kitazawa and Asakawa [92,93] can be used to relate net proton and net baryon cumulants in the EV-HRG model. To do that, we calculate net proton cumulants in an alternative way, namely by registering each baryon within the rapidity acceptance with a Bernoulli probability q = p / B . For an EV-HRG model at T = 160 MeV that we use one has q 0.33. The net proton cumulants computed in this way are shown in Fig. 6 by open blue diamonds. They agree with the actual net proton cumulant ratios shown by blue circles. This confirms that cumulants of net baryon distribution can be recovered from factorial moments of net proton distribution via a binomial unfolding with probability q. The value of q in experiment can be calculated from the measured mean multiplicities of the various baryon species. The neutron yield, which is not measured, can be reconstructed from proton yields using the isospin symmetry.

D. Comparison to ALICE data
The results we have discussed so far correspond to fluctuations of baryons and protons in acceptances integrated over all transverse momenta. This has not yet been achieved experimentally. Instead, the ALICE collaboration has published measurements of the variance of net-proton distribution in Pb-Pb collisions at √ s NN = 2.76 TeV in an acceptance in a 3-momentum range 0.6 < p < 1.5 GeV/c and longitudinal pseudorapidity |η| < 0.8 [15].
The top panel of Fig. 7 depicts the comparison between the data (symbols) and the EV-HRG model with exact baryon number conservation (black line) for the ratio κ 2 / p +p of net protons. The data are described by the model within errors. However, the data are described similarly well by the ideal HRG model, where this ratio is given by the binomial baseline, (κ 2 / p +p ) bino = 1 − α p [40]. Here α p = N p acc / N B 4π where N p acc is the mean number of protons in the acceptance and N B 4π is the mean number of baryons in the full space. This implies that measurements in these acceptance windows are not very sensitive to the equation of state. The deviations from the Skellam baseline are overshadowed by the global baryon conservation. The effect of repulsive interactions in the EV-HRG model is to slightly reduce the ratio further away from the Skellam limit. This is in contrast to the baryon and proton cumulants in p T -integrated acceptances that we have shown in Figs. 5 and 6, where the effect of interactions for one unit of rapidity is already sizable. The reason is due to cuts in the transverse momentum coverage. While the presence of radial flow does induce a level of correlation between the transverse momenta and coordinates of particles, this correlation is not as strong as in the longitudinal direction given by the Bjorken flow. The p T -cuts, therefore, lead to a Poissonization of the grand-canonical fluctuations, making it challenging to extract the grand-canonical susceptibilities. This underlines the importance of expanding the acceptance for fluctuation measurements in the future runs at the LHC in order for them to be sensitive to the equation of state.
We also explore the effect of exact conservation of electric charge and strangeness of net proton fluctuations. As shown in Ref. [87], a moderate effect of these extra conservation laws on net proton cumulants is expected. To evaluate the effect, we sample the grand-canonical multiplicities of all hadrons and resonances, including mesons, in the same fashion as before, but reject, in addition to baryon number conservation, all events which do not satisfy the exact conservation of global electric charge, Q = 0, and strangeness, S = 0. These two additional rejection steps slow down the event generator procedure considerably. Therefore, we generate a smaller number of events in the BQS-canonical ensemble, equaling to about 3 · 10 7 events. For this reason we restrict the analysis within the BQS-canonical ensemble to the second and fourth order cumulants. The κ 2 / p +p ratio from the BQS-canonical EV-HRG model is depicted by a dash-dotted magenta line in the top panel of Fig. 7. The exact electric charge and strangeness conservation leads to a further reduction of κ 2 / p +p by a moderate amount. This effect is consistent with results in reported in Ref. [87] using the SAM for multiple conserved charges.
The pseudorapidity dependencies of kurtosis and hyperkurtosis of net proton fluctuations within the same ALICE acceptance are depicted in the middle and bottom panels of Fig. 7, respectively. Similar to the variance, these show a suppression with respect to the Skellam baseline, mainly due to the baryon number conservation. It is notable that the hyperkurtosis never reaches a negative value within the ALICE acceptance. Again, this is a reflection of a limited p T coverage of the acceptance as well as of measuring only a subset of all baryons.

E. Volume fluctuations
We would like to discuss another issue which may affect fluctuation measurements in heavy-ion collisions, namely fluctuations of the system volume. The volume fluctuations do not affect the behavior of the mean quantities, but they do modify the fluctuations. This effect has been studied in several works in the literature [38][39][40]. Here we follow Ref. [39] to estimate the effect of volume fluctuation on our results.
We assume that, in the absence of volume fluctuations, all the cumulants obey linear scaling with the volume, κ n ∝ V . Let us denote byκ n the cumulants which include the effect of volume fluctuations. They read [39] κ n = n l=1 V l B n,l (κ 1 /V, κ 2 /V, . . . , κ n−l+1 /V ) .
Here V l is the lth cumulant of the system volume distribution and B n,l are Bell polynomials. Let us now take into account that all odd-order cumulants of net-particle distribution at the LHC vanish, κ 2n−1 = 0. In this case the odd-order order cumulants with volume fluctuations do vanish as well,κ 2n−1 = 0. The even order cumulants up to n = 6 read Hereṽ i = V i / V i are the scaled volume cumulants. The variance of a net-particle distribution at the LHC is not influenced by volume fluctuations, as pointed out before in Refs. [39,40]. However, the volume fluctuations do influence the higher-order cumulants. The cumulant ratios read A non-zero variance of the volume distribution influences the kurtosis and hyperkurtosis. In addition, the hyperkurtosis may be affected by a non-zero skewnessṽ 3 of the volume distribution. The effect of volume fluctuations is determined by the values of the reduced cumulantsṽ i . These are mainly determined by the collision geometry and the centrality selection. To illustrate the effect of volume fluctuations we will consider net-proton fluctuations in the ALICE acceptance that we discussed in the previous subsection. For simplicity, we shall neglect the skewness of volume fluctuations,ṽ 3 = 0, which could have an influence on the hyperkurtosis, but not on the kurtosis. To fixṽ 2 we make use of the ALICE measurement of the variance of proton number distribution [15]. As the mean number of protons is non-vanishing even at the LHC energies, the varianceκ p 2 of proton number distribution is affected by the volume fluctuations, in contrast to net-proton variance which is unaffected. Following Eq. (70) the proton number scaled variance reads ALICE has measuredκ p 2 / p = 1.07 ± 0.06 and p = 18.4±0.4 in an acceptance 0.6 < p < 1.5 GeV/c and |η| < 0.8. The EV-HRG model without volume fluctuations that we use, on the other hand, predicts κ p 2 / p = 0.98. Assuming that the difference between the model and the measurements can be attributed to volume fluctuations, one can use Eq. (77) to extract the value ofṽ 2 which describes the data:ṽ 2 = 0.005 ± 0.003.
The pseudorapidity dependence of the kurtosis and hyperkurtosis of the net proton distribution in 0-5% central Pb-Pb collisions at the LHC in the EV-HRG model with baryon number conservation and volume fluctuations in depicted in Fig. 7 by the red lines with bands. The bands correspond to the error propagation of the variance of the volume distribution in Eq. (78). The volume fluctuations have a large effect on higher-order fluctuations, both the kurtosis and hyperkurtosis exceed unity in all acceptances considered, in contrast to calculations without volume fluctuations where they lie below unity. It seems, therefore, that a significant reduction of volume fluctuations will be required in the future experimental measurements to be able to reliably control this effect. As an illustration, the dotted red lines in Fig. 7 depict the cumulant ratios when the variance of volume fluctuation is decreased by an order of magnitude, i.e.ṽ 2 = 5 · 10 −4 . In this case, the results are considerably closer to the cumulant ratios without volume fluctuations, and it should be possible to reliably extract these ratios by fitting the data via Eqs. (74)- (76). This type of analysis has been performed by the HADES collaboration in Ref. [22], where the next-to-leading order volume dependence of the cumalants was additionally considered. The centrality bin width correction [95] is another possible remedy, which has been applied for net proton measurements by the STAR collaboration [11].
F. Net-Λ, net-kaon, and net-pion fluctuations Net proton cumulants are not the only fluctuation measurement performed by the ALICE collaboration. Fluctuations of net numbers of Λ's, kaons and pions are also being performed, and preliminary results were reported in Refs. [96][97][98]. Here we would like to discuss the behavior of these quantities within our approach. The main purpose here is to illustrate how the different effects like resonance decays and exact conservation of various conserved charges influence the observables semiquantitatively. Where available, we do confront our predictions with the preliminary data as well. Our analysis here is restricted to the second cumulants normalized by the Skellam baselines, which at the LHC energies are free of the influence of volume fluctuations.
To perform the analysis we sample the full EV-HRG model, including both the (anti)baryons and mesons, using the same parameters as above. The cumulants are calculated after all strong and electromagnetic decays, in the ALICE acceptance, 0.6 < p < 1.5 GeV/c and a pseudorapidity acceptance |η| < 0.5 ∆η acc , where ∆η acc is varied up to a value of 3 units. We consider three configurations for the treatment of global conservation laws: (i) global conservation laws are neglected (grandcanonical); (ii) exact conservation of baryon number is enforced (B-canonical); (iii) exact conservation of baryon number, electric charge, and strangeness is enforced (BQS-canonical). Comparing the results between the three cases allows us to distinguish the roles of different conservation laws.
Let us start with the net-Λ fluctuations. The results are depicted in the top panel of Fig. 8. The κ 2 [Λ −Λ]/ Λ +Λ ratio shows a mild suppression relative to unity as the pseudorapidity acceptance ∆η acc is increased. A small suppression exists already in the grandcanonical limit (dashed blue line), which is attributed to the presence of repulsive baryon interactions modeled by the excluded volume. A larger effect is observed when the global baryon number conservation is incorporated (solid black line). An additional suppression from exact strangeness conservation on top of baryon conservation is also observed (magenta line), although this effect is rather small. This smallness is attributed to the fact that at LHC energies the dominant part of all strange quarks is carried by kaons, with Λ's forming only a small fraction of all strange particles.
The net kaon fluctuations are interesting because they are affected by a decay φ → K + K − of the φ meson. Our calculations, as well as experimental data [99], suggest that about 6% of final state K + and K − mesons come from this decay channel. The decay generates a correlation between the numbers of K + and K − . If  both decay products fall into a measurement acceptance, this gives no contribution to the variance κ 2 [K + − K − ] as the net number of kaons is unchanged. However, the total number of charged kaons, K + + K − , increases by two. For this reason one can expect the ratio κ 2 [K + − K − ]/ K + + K − to be below unity due to resonance decays alone, even in the absence of global conservation laws. This is indeed observed in our Monte Carlo simulations depicted in the middle panel of Fig. 8: the ratio κ 2 [K + − K − ]/ K + + K − is visibly below unity in the grand-canonical calculation which we attribute to the φ → K + K − decay. The net kaon fluctuations are virtually unaffected by the exact baryon number conservation (black line). This is expected because mesons do not interact with the baryons in the EV-HRG model, hence the baryon number conservation does not have an influence on meson distribution, except for small feeddown contributions from baryonic resonances.
The kaons are affected by strangeness and, to a lesser extent, electric charge conservation. The BQScanonical calculation is depicted by the dash-dotted magenta line, showing a further suppression of the varianceover-Skellam ratio when strangeness and electric charge conservation is implemented. The resulting ∆η acc dependence of net kaon fluctuations agrees with the preliminary data of the ALICE collaboration [96,98], shown in Fig. 8 by the gray bands, although the experimental uncertainties are quite large.
The behavior of net-pion fluctuations (the bottom panel in Fig. 8) is qualitatively similar to net kaons. The pion fluctuations are affected more strongly by resonance decays than kaons. Several resonances give a notable contribution. The most notable ones are decays ρ 0 → π + π − , ω → π + π − π 0 , η → π + π − π 0 , all leading to a sizable suppression of the ratio κ 2 [π + − π − ]/ π + + π − relative to unity already in the grand-canonical limit (the blue line). Baryon conservation has a negligible influence on net pion fluctuations, similar to net kaon fluctuations. Net pion fluctuations, however, are notably suppressed by the exact conservation of electric charge, see the dash-dotted magenta line. This should not come as a big surprise, as the charged pions constitute the majority of all charged particles at the LHC, hence the sizable effect of charge conservation on pion fluctuations.
The preliminary data of the ALICE collaboration on net pion variance-over-Skellam ratio lies somewhat below our BQS-canonical model prediction, the deviations are roughly on a two-sigma level. It should be cautioned that our predictions for net pion fluctuations should be regarded as semi-quantitative, for several reasons. For instance, we use the blast-wave model parametrization from Ref. [90] which underestimates significantly the number of soft pions, p T 500 MeV/c. Also, we neglect the effect of Bose statistics, which is non-negligible for the primordial pions at the chemical freeze-out. We also neglect additional effect due to rescattering in the hadronic phase. It is known, however, that the number of ρ 0 resonances reconstructed in central Pb-Pb experimentally is about 20-25% lower than predicted by the HRG model at the chemical freeze-out [100,101]. This indicates additional dynamics in the hadronic phase involving ρ 0 resonances and their decay products, which may change the effect of ρ 0 decays on κ 2 [π + − π − ]/ π + + π − . It is true however, that both the Bose statistics as well hadronic rescattering are expected to worsen the agreement with the data rather than improve it. The Bose statistics leads to an enhancement of pion fluctuations [102], whereas the hadronic rescattering will dilute the correlations between pions from resonance decays in given acceptance, both effects thus leading to an increase of κ 2 [π + −π − ]/ π + +π − . Nevertheless, our analysis is sufficient to indicate that net pion fluctuations are affected sizably by both the resonance decays as well as exact global conservation of electric charge. Both these mechanisms should thus be taken into account in interpretations of experimental data.

G. Dynamical net-charge fluctuations
We would like to conclude our analysis of experimental data with the variance of the net-charge distribution. The corresponding measurements have been performed by the ALICE collaboration and published in Ref. [103]. There, the measurements were focused on a quantity ν (+,−,dyn) , defined as Here N +(−) is the number of positively (negatively) charged particles in the final state for a given acceptance.
In the limit N + = N − , which to a large precision holds at the LHC, ν (+,−,dyn) simplifies to Here Q ≡ N + −N − is the net charge and N ch ≡ N + + N − is the charged multiplicity. ν (+,−,dyn) is thus closely related to the so-called D-measure: The D-measure was introduced in Ref. [3] as a probe that discriminates the charge-carrier degrees of freedom in the medium. In the quark-gluon plasma (QGP), where the quarks carry fractional charges, one has D ∼ 1 − 1.5 [3] in thermal equilibrium. For a gas of hadrons and resonances, on the other hand, the baseline value is considerably larger, D ∼ 3 − 4 [35].
A direct comparison of the baselines with experimental measurements of net-charge fluctuations is complicated by several additional effects, including volume fluctuations, exact charge conservation, and acceptance cuts. in the literature. Ref. [104] advocated an additive correction: Here N 4π ch is the mean charged multiplicity in the full space. Ref. [35], on the other hand, suggested a multiplicative correction: Here 1 − α ch is the charge conservation correction factor while C µ = N + 2 / N − 2 (= 1 at the LHC) additionally corrects for the effects of finite net charge.
The ALICE measurements in Ref. [103] include charge conservation corrections and incorporate the differences between D and D as a contribution to the systematic error. Here we analyze the behavior of the D-measure within our Monte Carlo sampling of the EV-HRG model at LHC conditions. Figure 9 depicts the pseudorapidity acceptance dependence of the D-measure of the dynamical net-charge fluctuations in 0-5% central Pb-Pb collisions at the LHC calculated in the EV-HRG within various statistical ensembles. A transverse momentum cut 0.2 < p T < 5.0 GeV/c is applied. This is the same p T cut as in the ALICE measurement. The dashed blue line depicts the behavior of the D-measure in the grand-canonical version of the EV-HRG model. As the grand-canonical calculation neglects the exact charge conservation, we calculate the D-measure in this case directly using Eq. (81), without applying any of the charge conservation corrections. The resulting D-measure is a decreasing function of ∆η acc that saturates at a value of around D ∼ 2.8 in the limit ∆η acc → ∞. The suppression of D relative to the Poisson statistics baseline of D = 4 is attributed to decays of neutral resonances into a pair of charged particles, like ρ 0 → π + π − . Here the discussion of resonance decays affecting net-pion fluctuations in Sec. IV F straightforwardly applies. We note that the influence of the excluded-volume effects in the baryon sector is virtually negligible, as the majority of charged particles at the LHC are mesons. Therefore, the results shown in Fig. 9 for the EV-HRG also apply to the standard ideal HRG model.
Calculations incorporating exact conservation of various conserved charges reveal that net-charge fluctuations are affected by exact conservation of the electric charge, while the additional influence of baryon number and strangeness conservation is observed to be negligible. This observation is consistent with the results of Ref. [87], where it was shown that the variance of a conserved charge distribution is only affected by exact conservation of that charge, but not of any other conserved charge. The black lines in Fig. 9 show the results of the BQS-canonical EV-HRG model calculation where we apply the charge conservation correction in accordance with Eq. (82) [D , dash-dotted line] or (83) [D , solid line]. This is the same procedure that was performed by the ALICE collaboration in Ref. [103] to correct for global charge conservation. If these corrections were exact, one would expect to reproduce the grandcanonical result shown by the dashed blue line. Instead, we observe that both the D and D appear to overestimate the charge conservation correction, especially D at large ∆η acc . The D correction does perform better than D and stays close to the grand-canonical result for ∆η acc 4.
The experimental data of the ALICE collaboration are shown by the symbols with error bars in Fig. 9. The data points lie considerably lower than model predictions. In particular, the slope of the curve at small ∆η acc is much steeper in the data than in the model. This result is in line with the tensions of the HRG model with the preliminary data for net-pion fluctuations discussed in Sec. IV F. The visibly stronger effect obtained for the D-measure can be attributed to a significantly larger transverse momentum coverage for the net-charge fluctuations relative to those for net pions. As discussed in Sec. IV F, the effects that we neglected in our calculations, such as the Bose-Einstein statistics for pions or hadronic rescattering, would be expected to enhance the D-measure and thus even further worsen the disagreement with the data. At this point we do not see a con-ceivable mechanism to explain the ALICE data within a purely hadronic description. The measurement, therefore, points to the suppression of net-charge fluctuations in central heavy-ion collisions at the LHC relative to the hadronic scenario. One tantalizing possibility here is the QGP formation, where a suppression of the D-measure is expected [3]. We hope that future measurements and analyses will shed more light on whether the observation of a suppressed D-measure constitutes a QGP signature.

V. DISCUSSION AND SUMMARY
In this work we introduced a subsensemble sampler -a novel particlization routine for heavy-ion collisions which preserves the thermal fluctuations and correlations in an interacting hadron resonance gas on a local level. It also takes into account global conservation laws on an eventby-event basis. The key of the procedure lies in partitioning the particlization hypersurface into locally grandcanonical subvolumes. In each subvolume, the hadron numbers are sampled from the grand-canonical multiplicity distribution, while their momenta follow from a thermal distribution imposed on a collective velocity profile. The global conservation laws are enforced via a subsequent rejection sampling step. The procedure allows to evaluate event-by-event fluctuations of various particle numbers in a momentum space acceptance, as appropriate for experiment, within a fluid dynamical picture of a local thermodynamic equilibrium at particlization.
The partition into subvolumes is not unique, the choice can be optimized for the applications on hand. Certain restrictions do apply. On the one hand, each subvolume V i has to be sufficiently large such that the cumulants of hadron multiplicity distribution are in the regime where they scale linearly with V i . On the other hand, the partition should be sufficiently fine grained, both relative to the acceptance where measurements are performed as well as to capture the coordinate space inhomogeneities in the distribution of thermal parameters. In the present work we considered the partition along the space-time rapidity axis (Fig. 1), which is appropriate to study the rapidity dependence of fluctuations integrated over the transverse momenta. Other partition schemes can be considered in a more general case.
As a first application of our new particlization routine, we studied event-by-event fluctuations in Pb-Pb collisions at the LHC, with a focus on the rapidity acceptance dependence of cumulants of the net baryon distribution. To that end, we utilized a hadron resonance gas model with excluded volume interactions in the baryonic sector, which matches well the available lattice QCD data on cumulants of net baryon distribution at a particlization temperature of T = 160 MeV. We used a blast-wave flow velocity profile and neglected any dynamics in the hadronic phase except for strong and electromagnetic decays of resonances.
Our Monte Carlo simulations reveal how baryon inter-actions, global baryon conservation, thermal smearing, and resonance decays affect the behavior of cumulants in a momentum acceptance around midrapidity, as is appropriate for experimental measurements. In the absence of thermal smearing and resonance decays, net baryon cumulants follow the analytic baseline established within a subensemble acceptance method (SAM) in Ref. [77]. One can therefore use the SAM to correct experimental measurements for the effects of global baryon conservation. However, for this to work the experimental acceptance needs to have a sufficiently large rapidity coverage, roughly ∆Y acc 1, and capture the entire transverse momentum range. The reason for that is the thermal smearing, which dilutes the signal for small acceptance and causes the cumulant ratios to approach the binomial distribution baseline (see the red points in Fig. 5). We do observe that this effect is well described at LHC by a simplified analytic model which assumes the thermal smearing in kinematical rapidity to be Gaussian (see Appendix). The resulting expressions are somewhat more involved than the simple formulas of the pure SAM framework, but it is possible they can be used to subtract the effect of thermal smearing from the data in addition to global conservation. The effect of an additional rapidity smearing of baryons due to decays of resonances is found to be negligible.
We find large differences between experimentally measurable net proton cumulants and the theoretically calculated net baryon cumulants. The net proton cumulants generally lie much closer to the Skellam baseline than the net baryon cumulants. This is a reflection of the fact that protons form a subset of all baryons. Measuring a subset, as opposed to the full set, dilutes the strength of correlations, which is reflected by the difference between net proton and net baryon cumulants in Fig. 6. We do observe that net proton and net baryon cumulants are indeed related to each other by a binomial (un)folding, as advocated by Kitazawa and Asakawa [92,93]. This result, however, does not by any means imply that one can directly compare net proton cumulant ratios with the net baryon ones. Such a comparison is not only unjustified, but is likely to lead to misleading interpretations and conclusions. For meaningful comparisons one has to reconstruct the net baryon cumulants from net proton ones through the binomial unfolding procedure described in [92,93]. This has not yet been achieved in the present experiments although the procedure is doable and, in fact, straightforward, requiring the use of the factorial moments of (anti)proton distributions that are readily accessible in experiment. On the other hand, the factorial moments of baryons and antibaryons are not directly accessible in lattice QCD, therefore, applying the method of Kitazawa and Asakawa to construct net proton cumulants from the lattice results on net baryon cumulants requires model assumptions. This observation underscores the importance of measuring the factorial moments of (anti)proton distribution in addition to net proton cumulants, as only in this case one can reconstruct the net baryon cumulants and make the comparisons with various theoretical predictions meaningful.
We confronted the predictions of our event generator with the available experimental data of the ALICE collaboration on the variance of net proton, net pion, net kaon, and net charge distributions. We find good agreement of our event generator with the net-proton data. However, the data are described similarly well by the binomial distribution baseline that corresponds to an ideal hadron gas model with baryon number conservation. In other words, the currently available measurements, performed in a 3-momentum and pseudorapidity acceptance, do not allow to distinguish the subtle effects associated with the QCD chiral crossover transition. The variances of net-pion and net-kaon fluctuations are not sensitive to the interactions in the baryonic sector and global baryon conservation, but they are affected by resonance decays and exact conservation of electric charge and strangeness. Our model describes the preliminary ALICE data on net kaon fluctuations within error bars. The model also describes the trends seen in the pseudorapidity acceptance of net-pion fluctuations although the preliminary data are overestimated roughly on a two-sigma level.
The HRG model we employ does not describe the AL-ICE data on net-charge fluctuations. The experimental data on the D-measure is significantly below the model predictions (Fig. 9). It seems doubtful that the measurement can be described within a purely hadronic description. A suppression of the D-measure, on the other hand, is expected in quark-gluon plasma phase [3]. In fact, this has been the primary motivation for the corresponding measurements. It remains to be seen whether the ALICE measurement is indeed a signal of QGP.
Measurements of higher-order cumulant ratios are affected by volume fluctuations. We estimated the effect for 0-5% central 2.76 TeV Pb-Pb collisions based on the available data of the ALICE collaboration on the first two proton number cumulants and the volume fluctuations formalism of Ref. [39]. We found the effect to be quite large for the kurtosis and hyperkurtosis of net proton fluctuations in the ALICE acceptance, changing the qualitative nature of the pseudorapidity window dependence of these observables. Therefore, removing the contribution of volume fluctuations will be essential for interpreting the experimental data, and our results indicate that the centrality selection should be optimized in the future LHC measurements of the higher-order net proton fluctuations.
The formalism developed in this work has many future applications. One natural extension are the studies of fluctuations at lower collision energies probed by beam energy scan programmes at RHIC and SPS facilities [10]. There, the effects of finite (baryo)chemical potentials, nonuniform rapidity distribution of thermal parameters, and absence of the longitudinal boost invariance will play an additional role [105,106]. One can also consider a particlization hypersurface and flow velocity profile emerging from a full (3+1)-dimensional hydro simulation as op-posed to the blast-wave model that we used here, which may additionally include viscous corrections. It would also be of interest to analyze the effect of rescaterrings in the hadronic phase which would enhance the effect of momentum smearing and thus dilute the signal [43]. This can be achieved by coupling the particlization to a subsequent hadronic afterburner such as UrQMD or SMASH. Here we present a simplified analytic model to account for the effect of momentum smearing on the cumulants of net baryon distribution measured in a p T -integrated acceptance. The formalism here is applicable for interacting HRG models where correlations between numbers of baryons and antibaryons are absent in the grandcanonical ensemble. This is the case, for instance, for the EV-HRG model that we use in this study.
The model consists of two steps: (i) the effect of momentum smearing is evaluated in the grand-canonical ensemble, i.e. neglecting the exact baryon number conservation; (ii) the SAM framework is applied to the result obtained in the first step to incorporate the exact baryon number conservation.
Let us start with the first part of the procedure. Consider all baryons and antibaryons at particlization that have a longitudinal space-time coordinate η s within a narrow range [η 0 s − ∆η s /2, η 0 s + ∆η s /2]. We assume that the physical volume corresponding to this range is large enough to capture all the physics associated with the correlation length, i.e. dV /dη s ∆η s ξ 3 . This means that, in the absence of exact baryon number conservation, the distribution of (net) baryons from this space-time rapidity range is independent from all other particles outside this range and is determined by the grand-canonical susceptibilities, namely κ B,gce n (η 0 s ) = dV /dη s ∆η s T 3 χ B n , |η s − η 0 s | < ∆η s .
In the absence of correlations between numbers of baryons and antibaryons in the grand-canonical ensem-ble that we assumed, the susceptibilities and cumulants are partitioned as follows: Here χ B + n and χ B − n are the grand-canonical susceptibilities of baryon and antibaryon number, respectively.
Consider now the baryons which end up in a longitudinal rapidity acceptance |Y | < ∆Y acc /2. Since the contributions from the different ∆η s slices are independent, the resulting cumulants of the accepted particles are a sum of the contributions from the individual slices. Therefore, let us calculate the contribution from a single slice |η s − η 0 s | < ∆η s . Let us denote by p(η 0 s , ∆Y acc ) the probability that a baryon with a space-time rapidity η 0 s at particlization ends up in this acceptance. This probability is determined by thermal smearing. Assuming that all (anti)baryons at a given space-time rapidity η 0 s end up in the kinematical acceptance independently from each other and approximating the probability p(η 0 s , ∆Y acc ) to be uniform in a range ∆η s , the cumulants of distribution of (anti)baryons in acceptance |Y | < ∆Y acc /2 that came from the space-time rapidity range |η s − η 0 s | < ∆η s are obtained by applying a binomial filter with the Bernoulli probability p(η 0 s , ∆Y acc ) to the space-time rapidity cumulants κ B ± ,gce l (η 0 s ): Here k bino n [p, {κ l }] is a nth order cumulant of particle number distribution obtained by applying the binomial filter with probability p to a distribution described by a set of cumulants {κ l } with l = 1 . . . n. The cumulant generating function C k (t) for the cumulants k n after the binomial filter can be expressed in terms of the corresponding cumulant generating function C κ (t) for cumulants κ before the filter [91,107]: The explicit result for the first four cumulants reads As already mentioned, the full result for cumulants κ B ± ,gce n (∆Y acc ) of all (anti)baryons in the rapidity acceptance is obtained by summing the contributions from all η s slices. One obtaines κ B ± ,gce n (∆Y acc ) = dη s k bino n p(η s , ∆Y acc ), dκ B ± ,gce n dη s , where dκ B ± ,gce Note that dV /dη s , T , and χ B ± l can all depend on η s in general case. The net baryon cumulant is then simply κ B,gce n (∆Y acc ) = κ B + ,gce n (∆Y acc ) + (−1) n κ B − ,gce n (∆Y acc ) .
How to evaluate the binomial probability p(η s , ∆Y acc )? We shall assume that the kinematical rapidity of each baryon is smeared around the space-time rapidity coordinate η s in accordance with a Gaussian distribution with a width σ y . The width can be estimated by analyzing the flow velocity and temperature profiles at the particlization hypersurface. For the blast-wave model that we use at the LHC one has σ y ≈ 0.3. The binomial probability reads p(η s , ∆Y acc ) = Equations (95) and (96) allow to calculate the influence of thermal smearing in the grand-canonical ensemble. In order to incorporate the exact conservation of baryon number we apply the SAM framework of Ref. [77]. The canonical ensemble cumulants that include both the effect of thermal smearing and global baryon conservation read Here α is the fraction of the total volume which is covered by the acceptance and β ≡ 1 − α. For the LHC energies that we study in this paper α = ∆Y acc /(2η max s ). It is also implied κ B,ce(gce) n ≡ κ B,ce(gce) n (∆Y acc ) in the equations above.
We would like to emphasize again that the thermal smearing model here is based on the assumption that numbers of baryons and antibaryons are uncorrelated in the grand-canonical limit. While this is the case for the EV-HRG model that we use in the present paper, this is not necessarily the case for other theories. Modifying the smearing model to the general case should be possible, and will require the use of binomial filtering applied to factorial moments of baryon and antibaryon distribu-tions, as discussed in [108] in the context of acceptance corrections to net baryon and net proton cumulants. It should also be possible to generalize the model to p Tdifferential measurements and smearing based on thermal distributions superimposed on a realistic flow velocity profile and a 3-dimensional particlization hypersurface, as appropriate for the differential momentum distribution measurements in experiment.