Correcting event-by-event fluctuations in heavy-ion collisions for exact global conservation laws with the generalized subensemble acceptance method

We introduce the subensemble acceptance method 2.0 (SAM-2.0) -- a procedure to correct cumulants of a random number distribution inside a subsystem for the effect of exact global conservation of a conserved quantity to which this number is correlated, with applications to measurements of event-by-event fluctuations in heavy-ion collisions. The method expresses the corrected cumulants in terms of the cumulants inside and outside the subsystem that are not subject to the exact conservation. The derivation assumes that all probability distributions associated with the cumulants are peaked at the mean values but are otherwise of arbitrary shape. The formalism reduces to the original SAM [V. Vovchenko et al., Phys.Lett.B 811 (2020) 135868 [arXiv:2003.13905]] when applied to a coordinate space subvolume of a uniform thermal system. As the new method is restricted neither to the uniform systems nor to the coordinate space, it is applicable to fluctuations measured in heavy-ion collisions at various collision energies in different momentum space acceptances. The SAM-2.0 thus brings the experimental measurements and theoretical calculations of event-by-event fluctuations closer together, as the latter are typically performed without the account of exact global conservation.


I. INTRODUCTION
Event-by-event fluctuations are among the key observables for probing the QCD phase structure experimentally [2], with measurements performed at experiments in a broad collision energy range from HADES-GSI to ALICE-LHC [3][4][5][6][7]. Proper theoretical modeling of fluctuations in heavy-ion collisions is crucial for correct interpretation of experimental data. Exact conservation of quantities like energy-momentum or the conserved QCD charges has a significant influence on many fluctuation observables [8][9][10][11]. On one hand, the effect of conservation laws prevents meaningful direct comparisons of experimental data to theoretical calculations based on the grand-canonical ensemble, like lattice QCD [12,13] or effective theories [14,15]. On the other hand, a proper accounting for exact conservation laws in dynamical models of heavy-ion collisions such as hydrodynamics is challenging, and most of the available implementations are restricted specifically to the model of ideal hadron resonance gas [16][17][18][19][20] (see, however, a recent Ref. [21] for a Monte Carlo sampling of an interacting hadron resonance gas).
Recently, a subensemble acceptance method (SAM) has been developed [1,22] that allows one to correct cumulants of conserved charges for the effect of exact global conservation laws. More specifically, the SAM expresses the cumulants measured in a coordinate space subsystem of a thermal system in terms of the grand-canonical susceptibilities and a fraction α of the whole volume covered by the subsystem. The method has been applied in Ref. [1] to describe baryon number fluctuations at the LHC energies, where cuts in momentum rapidity can be identified with the subsystem in the longitudinal coordinate space due to Bjorken flow. The applicability of the SAM to measurements at lower collision energies, like e.g. beam energy scan at RHIC, however, is less clear: the created system is not uniform, with matter properties differing strongly between central and forward-backward rapidity regions [23] while the absence of the longitudinal boost invariance does not allow one to properly correlate the experimental momentum acceptance with a coordinate space subvolume.
This work introduces the SAM-2.0, which addresses the limitations of the original SAM. Namely, the method allows one to calculate the effect of exact global conservation on cumulants measured in an arbitrary subsystem, for instance in momentum space as appropriate for the experiment. Also, the system need not be uniform or equilibrated, the method is applicable in any scenario where all the relevant distributions are highly peaked around the means. In addition, the SAM-2.0 allows one to evaluate the effect of a global conservation law, such as that of net baryon number, on the cumulants of a nonconserved quantity like the net proton number. This is particularly relevant for the experiment, where it is usually not possible to measure all particles contributing to a given conserved charge. A typical use case of the SAM-2.0 would be correcting the theoretical calculations of the "grand-canonical" event-by-event fluctuations within dynamical models of heavy-ion collisions, like hydrodynamics, for global conservation laws. In the limit of uniform thermal system and fluctuations in the coordinate space subvolume, the SAM-2.0 reduces to the results of the original SAM [1].
The paper is structured as follows. The SAM-2.0 is presented in Sec. II, covering fluctuations of both the conserved and non-conserved quantities. The applications of SAM-2.0 are illustrated in Sec. III on the example of net proton and net baryon fluctuations at the LHC. Limitations of the method are discussed in Sec. IV. Discussion and summary in Sec. V close the article. arXiv:2106.13775v2 [hep-ph] 6 Jan 2022

II. METHOD
Consider fluctuations of a random number, to which one can refer to as a conserved quantity (charge) B, that characterizes the system. The whole system is partitioned into two parts: the in-system (the acceptance) and the out-system (the complement). The charge B is additive between the subsystems, i.e. the total charge is the sum First, consider the fluctuations of B in the absence of constraints on the total value of the charge. The superscript "gce" will be used to denote this scenario, in analogy to the grand-canonical ensemble in statistical mechanics. It should be emphasized, however, that these fluctuations do not necessarily have to come from the grand-canonical statistical mechanics. Instead, the method is solely based on the following two assumptions regarding the unconstrained distributions: 1. The probability distributions P gce in (B in ) and P gce out (B out ) of the charges in the two subsystems are independent, i.e. the joint probability distribution factorizes: 2. P gce in (B in ) and P gce out (B out ) are unimodal distributions that are highly peaked at the mean, such that the maximum term method is applicable. To be more precise, it is assumed that, for any t in the neighborhood of t = 0, the following property holds Here B max in/out (t) corresponds to the maximum term of a generalized (unnormalized) probability distributionP gce in/out (B in/out ; t) = e tB in/out P gce in/out (B in/out ), thus B max in/out (t = 0) corresponds to the maximum of P gce in/out (B in/out ). One can show that the distributions P gce in/out are peaked at the corresponding mean values, i.e. B max in/out (t = 0) = B in/out (Fig. 1). Indeed, The schematic view of the "grandcanonical" (dashed lines) and the "canonical" (solid line) distributions of the conserved quantity B inside (blue) and outside (red) the acceptance.
where Eq. (3) was used, as well as the fact that dP gce in/out (B in/out )/dB in/out = 0 at B in/out = B max in/out . Without the loss of generality, one can represent the probability distributions as Here A gce in/out is the normalization factor while the shape of the whole distribution is fully encoded by f in/out (B in/out ). The fact that the distributions are peaked at the mean can be represented as dP gce or, more generally, for t = 0 one can define B in/out (t) through dP gce in/out /dB in/out [ B in/out (t)] = 0: The relation (7) can be used to express cumulants κ gce n,in/out of the P gce in/out distribution in terms of f in/out and its derivatives. This can be achieved by differentiating Eq. (7) and observing that, by definition, Taking a single derivative with respect to t allows one to calculate κ gce 2,in/out : The higher-order cumulants can be obtained by continuously differentiating Eq. (9) with respect to t where it is implied B in/out → B in/out (t). Appendix B presents a recursive procedure to express the derivatives of f in terms of κ gce n,in/out , which will be useful for the application of SAM.
The exact global conservation of charge B is achieved by rejecting all configurations that do not satisfy this conservation law where B = B in + B out . Due to the Kroenecker symbol δ Bin+Bout,B , the distributions of B in and B out that are subject to global conservation are no longer independent of one another. Here the goal is to evaluate the cumulants κ ce n,in of the constrained distribution of B in in terms of the cumulants κ gce n,in/out of the unconstranined B in/out distributions. The distributions P gce in/out (B in/out ) and P ce in (B in ) are depicted schematically in Fig. 1. A particular physical example which satisfies the assumptions of the method is a statistical mechanics system with a U (1) conserved charge B in the thermodynamic limit, where the implementation of global conservation laws corresponds to the transition from the grandcanonical to the canonical ensemble characterizing the full system. In this case the two subsystems may be identified with a partition in the coordinate space whereas the probability functions are proportional to the canonical partition functions, P ∝ Z ce ∝ e − V T f , where f is the free energy density. Such a setup forms the basis of the original SAM framework [1]. Note, however, that the SAM-2.0 presented here is more general and applies to any scenario that satisfies the two assumptions listed above. In particular, the system does not have to be uniform and in perfect equilibrium whereas the partition into subsystems need not be performed solely in the coordinate space.

A. Cumulants of a conserved quantity
The goal is to evaluate the cumulants κ ce n,in of the constrained distribution of B in in terms of the cumulants κ gce n,in/out of the unconstranined B in/out distributions. Such a mapping shall be denoted by It follows from Eq. (10) that the probability function P ce in (B in ) can be written as The distribution P ce in (B in ) is peaked at the "grandcanonical" mean B in given by Eq. (4). Indeed, by maximizing P ce in (B in ) with respect to B in one obtains an equation The first order t-dependent cumulantκ ce 1,in (t) ≡ B ce in (t) = ∂G Bin (t)/∂t corresponds to the maximum of a generalized probability functionP ce in (B in ; t) = e tBin P ce in (B in ). MaximizingP ce in (B in ; t) yields the equation where B ce out (t) = B − B ce in (t). The higher order cumulants κ ce n,in correspond to the derivatives ofκ ce 1,in (t) = B ce in/out (t) evaluated at t = 0: These cumulants can be evaluated recursively, by successively applying the t derivatives to Eq. (15) and evaluating them at t = 0.
The nth order derivative of Eq. (15) with respect to t can be obtained by applying Faá di Bruno (FdB) formula [24] to the right-hand-side (see Appendix A for details): where Eq. (16) was used and f (k+1) in/out [ B in/out ] is the (k + 1)th derivative of the function f in/out [Eq. (5)] evaluated at the mean B in/out . Here B n,k are Bell polynomials. One can now separate the term k = 1 in Eq. (17) from the rest of the sum. Using the property B n,1 (x 1 , . . . , x n ) = x n one obtains thus, κ ce n+1,in reads Equation (19) is a recursive relation for κ ce n+1,in , expressing this quantity in terms of the lower order cumulants up to κ ce n,in and the "grand-canonical" functions f in/out and their derivatives, which themselves can be expressed in terms of the "grand-canonical" cumulants κ gce n,in (see Appendix B). This recursive procedure thus defines the mapping function S in Eq. (11).
It is instructive to evaluate the leading cumulants explicitly. The first-order cumulant corresponds to the mean value B in which is unaffected by the conservation laws, thus κ ce 1,in = κ gce 1,in = B in . Setting n = 1 one obtains the variance where Eq. (9) was used for f in/out . For n = 2 one has Here f in/out was taken from Eq. (B6). The fourth-order cumulant (n = 3) is Any higher-order cumulant can be obtained from Eq. (19) in the same fashion. A Mathematica notebook to express cumulant κ ce n,in of arbitrary order n in terms of {κ gce n,in/out } is available via [25].
It is instructive to consider the following specific case: a uniform thermal system of volume V and temperature T which is partitioned into two subsystems of volumes V 1 = αV and V 2 = (1 − α)V in the coordinate space. This was studied in the original SAM paper [1]. In this case the unconstrained cumulants are determined by the grand-canonical susceptibilities χ B n as follows: κ gce n,in = αV T 3 χ B n and κ gce n,out = (1 − α)V T 3 χ B n . Substituting these relations into Eqs. (20) and (21) one obtains Here β = 1 − α. These relations reproduce the results of Ref. [1]. It should be emphasized however, that the SAM-2.0 is more general than the original SAM, and allows one to perform the correction for exact conservation for any two distributions that satisfy the assumptions (2) and (3). In particular, the method is applicable to momentum space acceptances.

B. Cumulants of a non-conserved quantity
The global conservation affects not only the cumulants of a quantity which is conserved, but also that of nonconserved numbers that are correlated to the conserved one. A common example encountered in heavy-ion collisions are cumulants of net proton number. Even though this quantity is not exactly conserved, the fact that (anti)protons form a subset of all (anti)baryons make the net proton number correlated to the net baryon number and thus being affected by baryon conservation [10]. This issue is of high practical relevance for heavy-ion collisions because the protons, and not the baryons, can be measured directly in the experiment. The argument can be extended to other measurements and conserved charges, such as the net kaon number affected by the strangeness conservation. Fluctuations of non-conserved quantities that are correlated to globally conserved charges have earlier been studied in Ref. [22] in uniform thermal systems on the level of second order cumulants. Here these considerations are extended to arbitrary highly peaked distributions and high-order cumulants.
The non-conserved number will be denoted by N . This number is correlated to a conserved quantity (charge) B. In the absence of exact global conservation of B, the joint distribution of numbers N in/out and B in/out inside/outside acceptance is characterized by a joint probability function W gce in/out (B in/out , N in/out ). From the definition it follows that As before, it is assumed that the joint distributions W gce in (B in , N in ) and W gce out (B out , N out ) are uncorrelated to one another, and that both distributions are highly peaked at the means, such that the maximum term method is applicable. The distributions W gce in/out (B in/out , N in/out ) are rewritten in the exponential form as follows The exact global conservation of B is achieved via the Kronecker delta: The joint probability function W ce in (B in , N in ) of the numbers B in and N in inside the subsystem reads The distribution W ce in (B in , N in ) depends on the "grand-canonical" distribution of the conserved charge B out outside of the subsystem but not on the N out distribution of the non-conserved number.
The cumulants of the joint distribution W ce in (B in , N in ) shall be denoted by κ ce n,m;in . The goal is to express these cumulants in terms of the joint "grand-canonical" cumulants κ gce n,m;in/out of the distributions W gce in/out (B in/out , N in/out ), i.e. to find the mappingS such that κ ce n,m;in =S κ gce n,m;in , κ gce n,m;out .
The cumulant generating function for the distribution W ce in (B in , N in ) reads Maximizing a generalized probability functioñ with respect to B in and N in gives the equations defining B ce in (t B , t N ) and N ce in (t B , t N ): where the following short-hands were introduced , .
The cumulants κ ce n,m;in can be evaluated via a recursive procedure. To do that a derivative ∂ n+m /∂(t B ) n ∂(t N ) m is applied to Eqs. (34) and (35). These derivatives can be evaluated with the help of the multivariate Faá di Bruno's formula in the combinatorial form, which presents them as a sum over partitions of a set (see Appendix A for the details on the notation). Applying the FdB formula to Eq. (34) yields The first sum runs over all partitions Π 2 n+m of a set {1, . . . , n + m} into blocks of two colors while the second sum corresponds to partitions Π n+m into blocks of a single color. One can separate the terms from both sums that correspond to the partitions consisting of single blocks: Note that the sums in the r.h.s of Eq. (39) contain cumulants κ ce in only up to order n + m. (35) and performing the same manipulations one obtains Equations (39) and (40) can be cast in the form of a system of linear equations for κ ce n+1,m;in and κ ce n,m+1;in : [g (2,0) in Solving the system of equations allows one to express the "canonical" cumulants κ ce n+1,m;in and κ ce n,m+1;in in terms of lower-order cumulants as well the derivatives of the "grand-canonical" distribution functions g in (B in , N in ) and f out (B out ). The derivatives of f out and g in can be expressed in terms of the "grand-canonical" cumulants κ gce in/out (see Appendices B & C, respectively), i.e. this defines the mapping functionS in Eq. (31). This completes the recursive procedure to evaluate the joint cumulants κ ce in of the (B in , N in ) distribution constrained by exact global conservation of charge B.
The leading order cumulants can be written down as follows. The first order cumulants -the means -are unaffected by the exact global conservations, thus To illustrate the formalism it is applied here to correct the cumulants of the net-proton and net-baryon number in different rapidity windows in heavy-ion collisions at the LHC for baryon number conservation. More specifically Pb-Pb collisions at √ s NN = 2.76 TeV are considered here. In Ref. [21] the fluctuations were analyzed using Monte Carlo sampling of hadrons at a blast-wave particlization hypersurface, with account for effects of baryon excluded volume, thermal smearing, and exact global conservation of baryon number. Here the fluctuations are first calculated analytically in the grand-canonical limit, and then the SAM-2.0 is applied to perform the correction for baryon conservation. Comparisons with the Monte Carlo data of Ref. [21] are then performed to verify the validity of the method. The fluctuations of net proton and net baryon number studied here correspond to measurements around midrapidity as a function of rapidity cut |y| < ∆Y acc /2. To account for the thermal smearing, it is assumed that the longitudinal rapidity y of each baryon is smeared at particlization around its space-time rapidity coordinate η s by a Gaussian. Under this assumption it is possible to calculate analytically the cumulants κ B ± ,gce n,in (∆Y acc ) of baryons and antibaryons in the given rapidity acceptance ∆Y acc in the grand-canonical limit. The cumulants κ B ± ,gce n,out (∆Y acc ) of (anti)baryons outside the acceptance ∆Y acc are evaluated in the same fashion. This procedure is described in detail in the Appendix of Ref. [21] [see Eq. (93) there] and the results of that paper are used here.
In order to apply the SAM-2.0 to correct net proton and net baryon cumulants for global baryon conservation one first needs to evaluate the corresponding grand-canonical joint cumulants κ B,p,gce n,m;in (∆Y acc ) and κ B,p,gce n,m;out (∆Y acc ) inside and outside the acceptance ∆Y acc . The procedure to evaluate κ B,p,gce n,m;in/out (∆Y acc ) is the same both for inside and outside the acceptance, thus the focus here is on the former case.
Let P (N B ) be the probability to have N B baryons in the acceptance. This distribution is characterized by the cumulants κ B + ,gce n,in (∆Y acc ). Each baryon ends up as a proton in the final state with a probability q ≈ 0.33 [21]. As the model does not contain isospin-dependent correlations at particlization, the selection of final state protons from all baryons is described by the binomial distribution [26,27]: where The cumulant generating function for the joint (N B ,N p ) distribution reads Here G B is the cumulant generating function of N B distribution and To obtain the joint cumulants κ B,p,gce n,m;in/out (∆Y acc ) of the accepted net-proton/net-baryon distribution one makes use of the following property of the model: there are no grand-canonical correlations between baryons and antibaryons. In this case one obtains κ B,p,gce n,m;in/out (∆Y acc ) = κ B + ,p + ,gce n,m;in/out (∆Y acc ) + (−1) n+m κ B − ,p − ,gce n,m;in/out (∆Y acc ) .
Note that in cases where the grand-canonical correlations between baryons and antibaryons are present, the procedure to evaluate κ B,p,gce n,m;in/out (∆Y acc ) will be more involved and will require the use factorial moments of baryon and antibaryon distributions [28].
After evaluating κ B,p,gce n,m;in/out (∆Y acc ), these cumulants are inserted into the recursive procedure defined by Eqs. (41) and (42) to evaluate the cumulants κ B,p,ce n,m;in (∆Y acc ) of accepted net-proton/net-baryon distribution that are constrained by exact baryon conservation. The following three ratios of cumulants are analyzed, for both the net protons,  [21], shown by the symbols. It is seen that the SAM-2.0 results are consistent with the Monte Carlo sampling procedure, for all values of ∆Y acc . This validates the method for conditions realized at the LHC energies and thus allows one to forgo the computationally expensive Monte Carlo sampling for calculating the highorder cumulants of (net-)proton number distribution.
The results also illustrate some of the advantages of the SAM-2.0 over the original SAM. In particular, the SAM-2.0 made it possible to compute analytically the net proton cumulants, which cannot be done using the original SAM since net proton is not a conserved quantity. For the net baryons the advantages are not as evident at the LHC, as the created system is largely uniform due to the longitudinal boost invariance, making also the original SAM applicable for net baryons. The differences between SAM-2.0 and SAM become more important at lower collision energies where the created fireball is no longer uniform. Applications of the SAM-2.0 framework for Au-Au collisions at RHIC beam energy scan can be found in Ref. [29]. , (b) κ ce 4,in /κ ce 2,in , and (c) κ ce 6,in /κ ce 2,in in 0-5% central Pb-Pb collisions at the LHC evaluated using within the model of Ref. [21]. The symbols depict the Monte Carlo results of Ref. [21]. The solid lines correspond to the analytical calculations where the correction for baryon conservation has been applied within the SAM-2.0 framework.

IV. LIMITATIONS
The SAM-2.0 is based on the two assumptions introduced at the beginning of Sec. II, namely that the unconstrained conserved charge distributions inside and outside the acceptance are independent of one another [Eq. (2)], and that these distributions are highly peaked at the mean values, such that the maximum term method is applicable [Eq. (3)]. In practice, the conditions (2) and (3) may not always be satisfied fully, and if that is the case, the true results may deviate from the SAM-2.0 expressions. These possible deviations are discussed in this section.

A. The peakedness of the distributions
The relation (3) is accurate in large statistical systems and it becomes exact in the infinite volume limit. In particular, this is the case when the total conserved charge B is large compared to unity. Note that the method is also accurate when B is small (or even vanishing like the net strangeness in heavy-ion collisions), but the system volume is sufficiently large. To see why this is the case one can consider a system of particles and antiparticles with vanishing net conserved charge, B = N B − NB = 0 [9]. Even though in this case the net charge is zero, the numbers of particles and antiparticles are large, N B , NB 1, and the B distribution is highly peaked around this small value of B.
As far as heavy-ion applications are concerned, the finite volume effects are essentially negligible in central collisions of heavy ions across a broad energy range [10,21]. The only possible exceptions to this are the vicinity of a critical point, where the multiplicity distributions become broad, or a first-order phase transition, where the distributions may be bimodal. In such cases, as well as in the case of small systems or net strangeness at low energies [30], the accuracy of the relation (3) should be verified explicitly.

B. The independence of the in/out distributions
The assumption (2) of independence of P gce in (B in ) and P gce out (B in ) is more subtle. As follows from Fig. 2, this assumption appears to hold true for baryon number distributions in rapidity acceptances at the LHC energies, where there is strong correlation the longitudinal momenta and coordinates due to Bjorken flow. One can now consider the extreme opposite: No correlations between momenta and coordinates of particles. To be more specific, consider a static thermal system of particles where B plays the role of the conserved particle number in the system. 1 It is assumed that correlations between particles may only depend on their coordinates, but not on their momenta. As shown below, this represents an example of the worst-case scenario for the accuracy of SAM-2.0. Consider the fluctuations of particle number in a particular momentum acceptance ∆p. As the absence of momentum-dependent correlations between particles has been assumed, the probability α that each particle ends 1 Compared to the baryon number in heavy-ion collisions, here the production of antiparticles is not considered. Such as a setup thus resembles nucleon number fluctuations in heavy-ion collisions at moderate energies, where antibaryons can be neglected.
up in the acceptance ∆p is independent of all other particles and corresponds to the ratio of mean number of accepted particles B in to the total number B: The B in cumulants are then given by the binomial distribution and the first three cumulants read: These expressions are independent of the equation of state, i.e. they are independent of the grand-canonical cumulants κ gce n,in/out that would characterize the system in the absence of exact global conservation.
One can now calculate these cumulants using the SAM-2.0. In the absence of global conservation, the total particle number B fluctuates and this is characterized by the grand-canonical cumulants κ gce n . As the probability that a given particle ends up inside or outside the momentum space acceptance ∆p is independent from all other particles, the grand-canonical cumulants κ gce n,in/out are obtained by convoluting the underlying distributions described by κ gce n with the binomial distribution [27,28,31]. For example, the leading three κ gce n,in read The expressions for κ gce n,out are obtained by a substitution α → 1 − α.
If the distributions of B in and B out were independent, the sum of the in/out cumulants would give the total cumulant. One obtains instead κ gce 2,in + κ gce 2,out = κ gce 2 (1 − 2αβ) + 2 αβ κ gce 1 = κ gce 2 , (67) κ gce 3,in + κ gce 3,out = κ gce 3 (1 − 3αβ) + 3 αβ κ gce 2 = κ gce 3 . (68) Here β = 1 − α. The in/out cumulants are not additive for n ≥ 2 in the general case, thus the distributions of B in and B out are not independent. The only partial case where the cumulants are additive is when the distribution of particle number B follows the Poisson distribution, i.e. κ gce n = B for all n ≥ 1.
Because the assumption of the independence of the B in and B out distributions does not hold in the presence of grand-canonical correlations between particles, the SAM-2.0 will not yield the expected exact result for the canonical cumulants of B in given by Eqs. (60)-(62). The difference between the SAM-2.0 and the exact result for the second order cumulant can be evaluated. One has κ gce 1 = B and κ gce 2 = B + δ, where δ = 0 corresponds to the presence of correlations between particles. The difference between the SAM-2.0 result, evaluated using (20), and the exact result (61) is It follows from Eq. (69) that, for |δ| < B, the sign of κ ce 2,in SAM−2.0 − κ ce 2,in exact is determined by the sign of δ, thus the SAM-2.0 overestimates the effect of correlations (which in this example should not contribute at all) to κ ce 2,in . Therefore, when the independence of the distributions P gce in (B in ) and P gce out (B in ) does not hold, the results of the SAM-2.0 may be taken as an estimate for the maximum possible effect of grand-canonical correlations in the cumulants of accepted particles constrained by global conservation. Note that the largest deviations occur α = 1/2 and disappear as α → 0 or α → 1.
In the context of heavy-ion applications, it should be noted that the correlation between the momenta and coordinates of particles due to Bjorken flow diminishes as collision energy is decreased. The applications of SAM-2.0 at moderate collision energies, for instance below those of RHIC-BES ( √ s NN < 7.7 GeV), should thus be performed with care.

V. DISCUSSION AND SUMMARY
This work introduced a generalized subensemble acceptance method -the SAM-2.0 -to correct cumulants of both conserved and non-conserved quantities in a subsystem (acceptance) for exact global conservation of a single conserved quantity. The method expresses the corrected cumulants in terms of the uncorrected cumulants inside and outside the acceptance. In contrast to the original SAM [1], the new method is applicable to non-uniform systems and momentum space acceptances, making it suitable for applications in heavy-ion collisions at various collision energies.
The practical use of the SAM-2.0 is in applying a correction for global charge conservation to theoretical calculations of various event-by-event fluctuations. One can take hydrodynamical simulations of heavy-ion collisions as an example: 1. First one calculates the cumulants of interest, both inside and outside the acceptance, without the effect of global conservation. In this case the standard grand-canonical particlization is performed, where emission of particles from every hypersurface element is calculated independently.
2. Then, one applies the SAM-2.0 to correct the cumulants for exact global conservation.
In such a way one calculates the cumulants affected by global conservation laws without the need to apply a time-consuming, if not intractable, procedure to account for exact global conservation laws at particlization. The method requires that the unconstrained probability distributions inside and outside the acceptance are independent [Eq. (2)] and highly peaked [Eq. (3)], but have otherwise arbitrary properties. It should be noted that only the total baryon number (as well as electric charge) is conserved exactly in heavyion collisions. This number contains contributions from participants and spectators, which both fluctuate. In the context of hydrodynamic simulations, the correction for baryon conservation should thus be applied to calculations for a fixed number of participants and then folded with participant fluctuations, as discussed in Ref. [32]. Note that the effect of participant fluctuations is less relevant for central collisions, where the number of spectators is small, as well as in those cases where the experimental data have been corrected for participant fluctuations.
It has been shown here that the SAM-2.0 accurately reproduces the effect of global baryon conservation on cumulants of net proton number fluctuations in Pb-Pb collisions at the LHC. This has been achieved by comparing the SAM-2.0 results with the direct Monte Carlo sampling. Applications of the SAM-2.0 to proton number cumulants in Au-Au collisions at RHIC-BES can be found in Ref. [29].
Here only a single globally conserved quantity was considered. This is sufficient in cases where a single conservation law has a dominant effect on the observable of interest. The method can be extended to multiple conserved quantities by utilizing cumulants of their joint distributions, similarly as it was done in Ref. [22] for the original SAM. Faà di Bruno's (FdB) formula is an identity generalizing the chain rule to higher derivatives. In case of multivariate derivatives, the "combinatorial" form of the FdB formula is particularly useful, which involves summation over partitions of a set.

Partitions of a set into uncolored blocks
A partition of a set is a grouping of its elements into non-empty subsets -the blocks. One can denote Π n as the set of all partitions of a set {1, . . . , n} [33]. For example, the set {1} has only a single partition, thus Π 1 = {{{1}}}. On the other hand, there are two different partitions possible for the set {1, 2}, either both elements are in the same block, or they are split into two different blocks: All partitions of a set Π n+1 can be constructed recursively, by utilizing the previously calculated partitions Π n . To construct Π n+1 one iterates over all partitions π in Π n and either adds the element n + 1 to one of the blocks of π or adds a new block to π which contains the element n + 1. In this way one iterates over all valid partitions Π n+1 of a set {1, . . . , n + 1}. The partitions of a set {1, . . . , n} can be used to evaluate the following nth order derivative of a composite function f [g(x)]: Here π runs over all partititions Π n of a set {1, . . . , n}, |π| is the number of blocks in the partition π, B runs over all blocks in the partition π, and |B| is the number of elements in the block B. Equation (A2) can be derived by successive applications of the chain, product, and sum rules, that are necessary to compute the higher derivatives. This expression is equivalent to the more common form of the FdB formula in terms of the Bell polynomials: , . . . , g (n−k+1) (x) . (A3)

Partitions of a set into colored blocks
Partitions of a set into colored (or labeled ) blocks generalize the standard partitions Π n into a single type of blocks [34]. This extension allows one to generalize the FdB formula to derivatives of multivariate composite functions. Assume that each block can be labeled by one of l distinct colors and denote Π l n as the set of all partitions of a set {1, . . . , n} into such colored blocks. E.g. for l = 2 one has Here the subscripts in the right-hand-sides encode the label (l = 1 or 2) of each block. Similarly to the partitions into unlabeled blocks, Π l n+1 can be computed recursively, by iterating over all partitions Π l n and either adding the element n + 1 to an existing block or adding a new block containing the single element n + 1. The only difference is that when creating a new block, all the possible l colors for this block have to be considered.
The partitions of a set into l labeled blocks allow one to evaluate the following nth order derivative of a multivariate composite function f [g 1 (x 1 , . . . , x n ), . . . , g l (x 1 , . . . , x n )] [35,36]: 2 Here π j is the subset of π which contains all blocks of color j, |π j | is the dimension of this set and |π| = l j=1 |π j |. It is instructive to consider a special case of l = 2 and f = f [g 1 (x, y), g 2 (x, y)], which is relevant for the applications studied in the present work. The derivative ∂ n+m f /∂x n ∂y m reads Here B 1 is the number of elements in block B that are smaller or equal than n and B 2 is the number of all other elements of B.

(B7)
Appendix C: Derivatives of g in/out in terms of grand-canonical cumulants The highly peaked joint distribution W gce in/out (B in/out , N in/out ) of B in/out and N in/out from Sec. II B reads W gce in/out (B in/out , N in/out ) = C gce in/out e −g in/out (B in/out ,N in/out ) .
One can derive the relations between the function g in/out and the "grand-canonical" cumulants κ gce in/out . Introducing the cumulant generating function G gce B in/out ,N in/out (t B , t N ) and applying the maximum term method one obtains the following equations defining B in/out (t B , t N ) and N in/out (t B , t N ): The "grand-canonical" cumulants are defined as κ gce n+1,m;in/out = κ gce n,m+1;in/out = Using Eqs. (C2) and (C3) one can rewrite the identities B in/out (t B , t N ) = B in/out and N in/out (t B , t N ) = N in/out as follows: B in/out g (1,0) in/out ( B in/out , N in/out ), g (0,1) in/out ( B in/out , N in/out ) = B in/out , N in/out g (1,0) in/out ( B in/out , N in/out ), g (0,1) in/out ( B in/out , N in/out ) = N in/out .
The recursive procedure to evaluate g The sum {π1,π2}∈Π 2 n+m contains two terms which correspond to partitions of the set {1, . . . , n + m} into a single block, of colors 1 and 2, respectively. Separating these terms from the sums in Eqs. (C8) and (C9) allows one to rewrite these expressions as a linear system of equations for g (1+n,m) in/out and g HereΠ 2 n+m = Π 2 n+m \ {|π 1 | + |π 2 | = 1} corresponds to all partitions that contain more than a single block. The sums in the r.h.s of Eqs. (C10) and (C11) contain the derivatives of g of order up to n + m. Thus, solving this linear system of equations provide a recursive relation to calculate g (1+n,m) in/out and g (n,1+m) in/out in terms of the cumulants κ gce in/out . The recursive procedure starts from n + m = 1.