{'fraction_non_alphanumeric': 0.02855839542180999, 'fraction_numerical': 0.009924940906514797, 'mean_word_length': 4.565779350669334, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Charge Conservation and Higher Moments of Charge Fluctuations

Higher moments of distributions of net charge and baryon number in heavy-ion collisions have been proposed as signals of fundamental QCD phase transitions. In order to better understand background processes for these observables, models are presented which enable one to gauge the effects of local charge conservation, decays of resonances and clusters, Bose symmetrization, and volume fluctuations. Monte Carlo methods for generating samplings of particles consistent with local charge conservation are presented, and are followed by a review of simple analytic models involving a single type of charge with a constant experimental efficiency. The main model consists of thermal emission superimposed onto a simple parameterization of collective flow, known as a blast-wave, with emission being consistent with individual canonical ensembles. The spatial extent of local charge conservation is parameterized by the size and extent over which charge is conserved. The sensitivity of third and fourth order moments, skewness and kurtosis, to these parameters, and to beam energy and baryon density is explored. Comparisons with STAR data show that a significant part of the observed non-Poissonian fluctuations in net-proton fluctuations are explained by charge and baryon-number conservation, but that measurements of the STAR collaboration for fluctuations of net electric charge significantly differ from expectations of the models presented here.


I. INTRODUCTION
The fluctuation of conserved charges is a standard means by which to investigate and classify phase transitions. At the critical point correlation lengths diverge, which results in peaks in charge fluctuations as one approaches the critical point. For systems with first-order phase transitions, fluctuations turn into phase separation and fluctuation measures are no longer intensive quantities. The growth of fluctuations becomes increasingly dramatic as one considers progressively higher-order fluctuations. In a volume V , fluctuations of a charge Q can be defined as when particles have unit charge. The measure is increasingly sensitive to the tails of the multiplicity distribution, P n , as n increases. The free energy, F ∼ a 2 (n − n) 2 + a 3 (n − n) 3 · · · , is minimized for n = n, but the quadratic part vanishes at the critical point, a 2 → 0, which allows the fluctuations to grow and be dominated by higher-order terms. The shapes of the tails of the distribution are then profoundly altered. The properties of the QCD transition, deconfinement and the restoration of chiral symmetry, are not well understood at finite baryon density. There exists the possibility that this transition is a true phase transition, with a critical point at several times nuclear density and with a critical temperature close to the pion mass. If this is the case, it begs the question as to whether the conditions for phase separation or for critical phenomena can be reproduced in the laboratory. Heavy ion collisions at high energy, measured at the Relativistic Heavy Ion Collider (RHIC) or at the LHC, can produce mesoscopic regions at temperatures of a few hundred MeV, which is well above the expectations for a critical temperature, and densities of several times nuclear matter density. These densities might or might not be sufficiently high to investigate the phase transition.
High-energy heavy ion collisions are characterized by strong explosive collective flow. Measurement is confined to the outgoing asymptotic momenta, but because of strong flow, correlations in coordinate space manifest themselves as correlations in relative momentum. Thus, measurements of correlations binned by relative momentum or charge fluctuations within some defined region of momentum space serve as surrogates for the corresponding observables in coordinate space. Indeed, measurements of charge and baryon number fluctuations have been performed at RHIC. By adjusting the beam energy of the colliding nuclei, experiments at RHIC have explored conditions at which novel phase phenomena might occur. Fluctuations of electric charge and baryon number have been especially popular. An initial scan of beam energies  was rather inconclusive, but measurements with greatly improved statistics are currently being undertaken and analyzed. The main thrust of these studies is to search for evidence of a QCD phase transition with a critical point at finite baryon density, with the hope that the phenomena one expects for an idealized equilibrated system [22,23] becomes manifest in the measured debris of heavy-ion collisions. Even if no true phase transition exists, charge fluctuations, or equivalently susceptibilities, are fundamental properties of the quark gluon plasma and can be investigated with lattice gauge theory [24].
In addition to the finite system size (event multiplicities might number in the thousands), the novel states of matter created in heavy-ion collisions persist for 10 fm/c. This severely limits the degree to which phases can separate or to which critical fluctuations can grow. This also limits the extent to which conserved charges can separate from one another. For example, if a strange and an anti-strange quanta are produced together their separation is limited by diffusion, which is difficult to calculate in lattice gauge theory [25,26], but can be roughly extracted experimentally and is gaining attention theoretically [27][28][29][30]. To justify thermal models of fluctuations of conserved charge based on equilibration, sufficient time is required for particles to enter and exit some defining volume.
The first goal of this paper is to gauge the degree to which charge-balance correlations affect higher-order correlations. Charge balance functions, which are two-particle correlations related to charge conservation, have been measured extensively [31][32][33][34][35][36][37][38][39][40][41], and modeled theoretically [42][43][44][45][46][47][48][49]. In addition to making it difficult for phases to separate or for critical correlations to emerge, local charge conservation also represents its own source of correlation, which needs to be understood as a potential source of background before making firm arguments to have observed phenomena related to phase transitions. It is well known that charge-balance correlations are readily measurable at the two particle level, N = 2 in Eq. (1), and that they explain the bulk of the N = 2 fluctuation measurement. However, their impact on N = 3, 4 fluctuations has not be investigated in great detail. For instance, the four-particle measure of correlation, which is based on cumulants, subtracts much of the contribution to M 4 coming from purely two-particle correlations. However, charge conservation can involve multiple particles, and the degree to which a cumulant-based measure, like the kurtosis, is affected by charge conservation is not fully understood. Relations based on a uniform acceptance probability and for a single type of charge were worked out in [50], and provide significant insight into how higherorder correlations are affected by charge conservation. The goal of this study is to extend such ideas to a more realistic picture, which takes into account the conservation of all three types of charge (baryon number, electric charge and strangeness) and applies a more realistic model of experimental acceptance and efficiency. The interplay of charge conservation with chemical equilibrium, decays, and Bose statistics are all considered.
To understand the role of chemical equilibrium and decays, a model is presented which creates small volumes in which the net charges B, Q and S are each fixed at some value. Even if the net charges are all zero, charged particles exist in combinations that conserve the net charge. Theoretical methods for exact calculation of the canonical ensemble and a method for Monte Carlo generation of statistically independent events are presented here. The method lends itself to including decays and accounting for experimental acceptance and efficiency. The physical picture of treating small volumes as independent patches was previously done for calculation of charge balance functions in [51][52][53], and was also applied in [54]. Following the terminology in [54], we sometimes refer to these sub-volumes as patches. The method presented here creates perfectly independent samplings, and can generate billions of such patches within a few hours. This enables highly accurate calculations of higher moments with minimal numeric cost.
After a brief review of cumulants and the definitions of skewness and kurtosis in Sec. II, the method for exact calculation of the canonical ensemble describing a multi-component, multi-charge hadron gas is presented in Sec. III. These techniques extend those used for canonical ensembles used to study isospin fluctuations of a hadron gas [55], nuclear fragmentation [56], the level density of a Fermi gas [57], and the effect of restricting a quark-gluon plasma (QGP) to having fixed charge, including being in an overall color singlet [58]. Exact methods for calculating correlations up to fourth-order are presented. Unfortunately, when including realistic acceptance effects and complex decays, the exact expressions are no longer viable. However, as shown in Sec. IV, the exact expressions show how sample events can be generated. Each event, defined by a set of particles and momenta, is generated with perfect independence from the others, and perfectly reproduces the canonical expressions from Sec. III. Section V extends the previous sections to show how Bose correlations can be included. Section VI considers the case of a single type of charge and a uniform efficiency, i.e. the probability of any particle being recorded is set to some fixed value. Much of this discussion repeats what is said in [50], and is included for completeness. This provides for a physical discussion of how charge conservation, volume fluctuations, chemical equilibrium, decays, clustering, and Bose corrections should affect higher moments.
The heart of this study is presented in Sec. VII. Here, the patches are assigned collective velocities consistent with the collective flow deduced from heavy-ion collisions. A canonical sampling of particles is generated from each patch, followed by a simulation of their decays. The particles are then overlaid onto the acceptance of the STAR detector at RHIC. Each patch is uncorrelated with any other patch, so moments can be calculated by averaging over the independent contributions from the patches. Results are displayed alongside results from the STAR Collaboration. The size and sign of the fluctuations of the net-proton distributions are consistent with observations, but the calculations of the net charge distributions differs qualitatively from STAR observations. A detailed discussion of the lessons derived from this study is presented in Sec. VIII.

II. CUMULANTS, SKEWNESS AND KURTOSIS
For the manuscript to be more self-contained, a brief review of cumulants and the definition of skewness and kurtosis is presented here. Cumulants of a charge distribution are defined by where Q is the net charge. Here, Q might refer to baryon number, to strangeness, or to the electric charge measured in units of e. Rather than showing the cumulants C n , ratios are presented to help minimize trivial dependences on system size. The skewness, S, is a measure of the third moment, This definition has the advantage in being dimensionless, but it does not become independent of volume in the limit of large volumes. Thus, it is more common to consider the ratio which becomes an intensive measure in the limit of larger volumes. However, in this study we consider the ratio which is also intensive, and approaches unity for a uncorrelated emission, i.e. a Skellam distribution. The kurtosis is a measure of four-particle correlations, but instead of K, one typically chooses where σ 2 = C 2 , to find an intensive measure of the fluctuation. For a measure to be intensive, it should be independent of volume in the large-volume limit. For small volumes, charge conservation alters the average densities of various species, which is known as canonical suppression. Canonical suppression also distorts the higher moments for smaller volumes.
The ratios C 4 /C 2 and C 3 /C 1 approach simple values in the limit that the distributions would be Poissonian. For Poissonian emission the observation of a charge in one region of momentum space is uncorrelated with the emission into any other space. Thus, particles are correlated only with themselves. If charges appear only in integral positive units, one can apply the usual expression for the Poissonian moments where the mean is η, If there exist both positive and negative charges, the distribution of the net charge can be derived by convoluting the two distributions. Convoluting two Poissonians results in a Skellam distribution. If the mean number of positives is η + and the mean number of negatives is η − , the distribution of net charge for a Skellam distribution, Q = n + − n − , yields the following cumulants Thus, if charges are produced in an uncorrelated fashion in increments of ±1, the skewness and kurtosis become Even though most of the literature focuses on Sσ = C 3 /C 2 , this study presents results for C 3 /C 1 so that one can better understand deviations from the uncorrelated baseline. Moments depend on the efficiency α with which particles are measured. In the limit of vanishing efficiency all distributions of positives or of negatives tend to become Poissonian [59], and the distribution of the net charge will thus become Skellam. Then Eq. (10) shows that as α → 0, C 4 /C 2 = C 3 /C 1 = 1. This can be understood by seeing that as α → 0, the moments are dominated by the probability of observing either zero charges or a single charge. The probability of observing a single charge is α → 0, while the probability of observing two charges is proportional to α 2 , which is negligible. This assumption would fall through if multitple charges were observed on individual particle, but for final-state hadrons the charges are only ±1.
If the net charge is fixed, a non-perfect efficiency is required to produce fluctuations. For fixed charge, the efficiency divides particles into two sets, the measured and non-measured. Each set fluctuates equally, but oppositely, relative to the mean. Thus, all the even moments of net charge will have an even reflection symmetry about an efficiency of 1/2, and the odd moments will have an odd symmetry, C n (1/2 + δα) = C n (1/2 − δα), n = 2, 4, 6 · · · −C n (1/2 − δα), n = 3, 5, 7 · · · . (12)

III. RECURSIVE TECHNIQUES FOR GENERATING CANONICAL PARTITION FUNCTIONS
For non-interacting particles the canonical partition function can be calculated exactly, or at least to the level that all partitions of A ≤ A max hadrons are taken into account, with the exact solution being reached at A max = ∞. For our case, we conserve three quantities: the electric charge Q, the baryon number B and the strangeness S. For states i with energies E i , the partition function, where Q i , B i and S i are the discrete values of the conserved quantities for the state i, can be calculated recursively. The function Z A (Q, B, S) refers to the subset of states with A hadrons, The recursive procedure begins with the canonical partition function of the vacuum. The contribution for a given A, Z A (Q, B, S), can be written as where z h is the single-particle partition function for hadron species h, which has charges q h , b h and s h . This was proved in [56], and can be understood by realizing that one can count all the ways to arrange A hadrons with a given charge by considering all the ways to arrange one hadron multiplied by all the ways to arrange the remaining hadrons.
To avoid double counting, a factor of 1/A is applied. For a fixed charge the probability to have A hadrons is, In practice, the sum over A is cut off at some A max , but in our studies here that cutoff is made large enough that contributions to Z for A > A max are negligible. Thus, once one builds the partition function from A = 0 to A max one has the partition function for all Q, B, S. Once the partition function is calculated one can also calculate the multiplicities and moments of observing specific species. For example, the multiplicity of species h in a system with charge Q, B, S is This also provides expressions for the various charges, e.g., Spectra can also be calculated. For species h with spin j h , Second-order moments can also be calculated exactly, It is straightforward to extend this expression to higher-order fluctuations. These expressions can also be extended to consider non-additive conservation laws. Net isospin conservation of a hadron gas was invoked in [55], i.e. restricting the states to being in an iso-singlet. Quark-gluon states restricted to being in both an iso-singlet and a color singlet were addressed in [58]. Bose and Fermi corrections are discussed in Sec. V.
For the case of a single kind of charge, one can see how the the recursive method above yields the same result as what one would expect by writing down the partition function for a system with (A − Q)/2 negative charges and (A + Q)/2 positive charges as shown in [50], where z is the partition function of a single charge. One can readily see that this is consistent with the recurrence relations, This result is also equivalent to expectations based on setting reaction rates equal. If one assumes that pairs are created with some rate β, and that they are destroyed with some rate αN + N − , where N + and N − are the number of positive and negative charges, N + + N − = A. Setting the rates equal, One can see that if one chooses β/α = z that Eq.s (25) and (22) are consistent. If the net charge is zero, the result is even simpler, Aside from the assumptions that Q is fixed and that there exists only one kind of charge, Eq.(25) also requires that Bose and Fermi quantum statistical corrections are negligible, and that only unit charges exist. Despite these shortcomings, this picture is useful in that it allows one to see how multiplicity fluctuations are affected by charge conservation in a simple model.

IV. GENERATION OF UNCORRELATED SAMPLE EVENTS
Complicated experimental acceptances are difficult to incorporate into expressions for the moments. It is then easiest to generate entire events via Monte Carlo, and filter the events through the acceptance. The Monte Carlo procedure involves choosing a hadron proportional to the number of ways the system might have such hadron, i.e. a product of the partition function of the individual hadron multiplied by the partition function of the remainder. The procedure becomes: 1. Calculate and store the partition function, Z A (Q, B, S), up to some size A ≤ A max for values of Q, B, S that might ultimately couple back to a given A = A max for the given total values Q, B, S.
Bose effects can be included by altering the second and third steps above. This is addressed in Sec. V. Storing the partition function can require substantial memory for large A max because the indices Q, B and S must also vary over a range of order ±A max , so memory usage roughly scales with A 4 max . Because one is usually interested in calculations with total charge near zero, one can ignore partition functions for charges that cannot couple back to the fixed overall charge at A max . Once A exceeds A max /2, the calculations here cut off values of Q, B and S that could not ultimately affect the Q = B = S = 0 partition function for A = A max . Even with this savings, partition functions with A max = 250 could require approximately 12 GB of memory, and need on the order of 10 minutes to calculate on a single processor. For A max = 125, less than a GB of memory was needed and partition functions could be calculated in less than a minute. For hadron gases at temperatures of 150 MeV, A max = 250 was sufficient for patch volumes 700 fm 3 . If multiple patch volumes are to be explored for the same temperature, computational time can also be saved by realizing that the partition functions scale as Ω A . Thus, if one performs a calculation for some initial volume Ω 0 , scaling can provide results for new volumes with minimal computation.
Once the partition function is calculated, event generation is remarkably fast. The time to generate an event scales linearly with the volume, or equivalently, linearly with the average number of particles generated. Running sufficient events to generate a million individual particles can be accomplished within a few seconds on a single CPU. Unlike Metropolis methods where events are modified by considering small changes to existing events, such as in [54], each event in this method is perfectly independent of previous events.

V. BOSE AND FERMI STATISTICS
Including Bose and Fermi statistics into the recursive relations for partition functions is straight-forward, and was shown in [55,57]. The method is related to that used for calculating the effects of multi-boson interference for pion interferometry [60]. In a fixed volume the partition function can be first treated as the usual procedure of accounting for n identical particles being in different single-particle states. This includes the 1/n! term to account for the fact that the particles are indistinguishable, i.e. the Gibbs paradox. If m indistinguishable particles are in the same single-particle state , one must correct the weight by a factor of m ! for each level, which can also be thought of as the analog of the symmetrized relative wave function with all the momenta being equal. For fermions, the weight becomes (−1) −1 m . As demonstrated in [57], the recurrence relation to the partition function then becomes where z h,n is the partition function for n particles in some level, and refers to single-particle levels of energy . The ±1 refers to bosons or fermions. For hadron gases in the high temperature environments of relativistic heavy ion collisions, only pions have a non-negligible correction from quantum degeneracy. The correction of order n for any level is of the order e −β lower than the previous term. This factor is largest for zero momentum, and for pions becomes e −βm , where m is the pion mass. For the zero-momentum level at a temperature of 150 MeV, the factor is e −m/T ≈ 0.4, and as the system cools the factor falls slightly [61]. For a more characteristic thermal momentum the factor is ≈ 0.1. For heavier particles the factor is always small in the context of relativistic heavy ion collisions. For example, for a zero-temperature ρ meson the factor is a fraction of a percent. Given that symmetrization is only being applied to pions, which are bosons, one can incorporate these corrections into the Monte Carlo procedure outlined in Sec. IV. For fermions this might be problematic because of the negative weights coming from the (−1) n−1 factors in Eq (27), but fortunately this is unneccessary because the degeneracy of fermions is negligible in the systems considered here. For pions, the algorithm is adjusted by treating each value of n as being a different species, with charges nq h and with the partition function calculated with a reduced temperature T → T /n. If one picks such a species in step 3 of the algorithm, n pions are generated, all with the same momentum. For a finite system, the n pions would be assigned small relative momenta on the order of the inverse system size.
It is well known that bosonic effects can broaden multiplicity distributions, consistent with negative binomial distributions [62,63], making them super-Poissonian. One of the goals of this study is to discern how bosonic statistics alter the skewness and kurtosis.

VI. MODELS WITH UNIFORM EFFICIENCY AND FIXED CHARGE
Even for a volume of fixed charge, finite efficiency and acceptance leads to non-zero fluctuations. The degree to which these fluctuations affect the skewness and kurtosis was worked out in [50] for emission of a fixed charge where the probability of any charge being observed was a constant α. If α were zero or unity, there would be no fluctuations, and because the charge on those particles that are not observed must fluctuate exactly opposite to the charge that is observed, the even moments must be symmetric about α = 1/2, and the odd moments must be anti-symmetric. One of the most important results of [50] is that for fixed Q and α the ratios of cumulants depend only on α and the variance and mean of the underlying multiplicity distribution. Even though Q is fixed, the net number of charged particles M can fluctuate.
From [50], the probability that M charged particles with total charge Q will result in a measured charge q = n + −n − due to a uniform efficiency α is the convolution of two binomial distributions After some tedious algebra, one can find the cumulants for fixed multiplicity, Here, the primes emphasize that the averages · · · M . denote that they consider only those events with fixed base multiplicities M . The even moments are all linear in M , while the odd moments are linear in Q. It might appear that due to the linearity one might replace M with M for a fluctuating base multiplicity, but the second term in the fourth cumulant Instead, This provides a contribution to C 4 , and the cumulants and their ratios become [50] , The relative variance of the multiplicity of the base distribution, ω M , is unity for a Poissonian distribution. In that case C 3 /C 1 and C 4 /C 2 fall below unity for non-zero α. The assumptions for this relation are that the charge and volume are fixed, that the efficiency is uniform, and that Bose and Fermi symmetrization is ignored. This relation is important as it implies that if one understands the efficiency and the second moment of the base multiplicity distribution, one can construct a baseline of cumulant ratios, and attribute any deviation from the baseline as due to fluctuations in Q, which is precisely the goal. A reasonable value of α can be taken from balance function analysis. If a charge is observed, there should be a balancing charge emitted nearby, and detected with probability α. This should correspond to the integrated strength of the charge balance function, which for electric charge is in the neighborhood of 0.35 assuming the full acceptance of the STAR TPC. Of course, α is not a constant. If a charge is observed in the center of the detector, its balancing charge has a better chance of being observed than for one observed near the periphery of the acceptance. Nonetheless, for the purposes of roughly setting expectations, this suggests that C 4 /C 2 and C 3 /C 1 can significantly differ from unity. In the next five subsections, VI A thrugh VI E, we consider how various effects, aside from fluctuating the base charge Q, might push ω M to be either super-Poissonian, ω M > 1 or sub-Possonian, ω M < 1. We discuss the effects of volume fluctuations, charge conservation with a single type of charge, decays, Bose condensation, and finally the sensitivity to considering a realistic collection of resonances accounting for all three types of conserved charge.

A. Volume Fluctuations
Heavy ion experiments measure collisions spanning a range of impact parameters. Even for a fixed impact parameter, energy deposition might significantly vary depending on how many nucleons actually collided or how many jets were produced. If more energy is deposited into a fixed volume, it might expand further before it hadronizes, resulting in larger volumes when the system hadronizes. Experimental analyses attempt to minimize these fluctuations by constraining a given fluctuation measurement to a specific centrality bin, where "centrality" might be defined by multiplicity, transverse energy, or energy deposition in a forward calorimeter. To reduce auto-correlation, centrality measurements are usually constrained to particles other than those used to construct the moments. Nonetheless, it is inevitable that a range of initial conditions is explored within any centrality bin, and one might thus expect the base multiplicity distribution to broaden. If the patch volumes are fixed at some constant value, volume fluctuations can be thought of as a fluctuation in the number of patches.
It is tempting to expect that the ratio of cumulants would be independent of the number of patches, and thus impervious to volume fluctuations. Independently, each cumulant C n scales linearly with the number of patches because there are no cross correlations between patches. Thus, if P p (N ) is the probability of having N patches, and if c n is the cumulant for a single patch, the cumulant for the overall system is and the ratio of cumulants is Thus, the ratio of cumulants is independent of the number of patches, and so independent of the overall volume, though the ratio can still depend on the volume of an individual patch. Unfortunately, although the ratios of cumulants are independent of the overall volume, they are not impervious to fluctuations of the overall volume, or equivalently, to fluctuations of the number of patches [64][65][66]. To illustrate this, one can consider a system with a distribution P p (N ) as in Eq. (34), where the average number of patches is N . For n > 1 the cumulants are expressed in terms of δQ = Q − Q. Because Q varies with the number of patches we rewrite it as where Q 1 = c 1 is the average charge in a single patch. Here δQ N is the fluctuation of the charge for a specific number N of patches, and δN = N − N . Inserting these definitions into the expressions for the cumulants, where κ n are the cumulants of P p (N ), If the number of patches is fixed, i.e. the overall volume does not fluctuate, then κ n>1 = 0 and each cumulant satisfies C n = N c n . The ratio of cumulants then cancels the factor N . Once the volume fluctuates, the ratios of cumulants of the charge distribution depend on the ratios of cumulants of P p (N ). Methods have been constructed to reduce the dependence on volume fluctuations [65][66][67][68], and have been applied to STAR data [69,70]. These methods are built on the assumption that the observable used to identify the volume is correlated to the fluctuating charge only through the fact that they both scale with the volume. If the phase space used for the centrality measure is clearly distinct of the phase space over which charge is measured, this should be a good approximation. Equations (37) and (38) are built on the assumption that the average charge scales linearly with the number of patches. This differs from the assumptions going into Eq. (33), where it was assumed that charge was fixed within the overall volume. If the total charge in the volume is absolutely fixed despite fluctuations in the number of patches, then volume fluctuations affect the answer in that the relative variance of the multiplicity distribution, ω M , is increased by volume fluctuations as described in Eq. (33).

B. Chemically Equilibrated Canonical Distribution
One can calculate the relative variance, ω M , for an equilibrated system of a single type of charge, where the measurement directly samples the equilibrated canonical distribution. This ignores decays, which in later stages proceed without the regeneration needed to maintain equilibrium. The results of the equilibrated canonical ensemble in Eq. (26) can be used to find ω M , the relative variance of the multiplicity distribution, and thus generate the moments using Eq. (33). For large average multiplicities M , the value of ω M approaches unity, the Poissonian limit. For small M it approaches two. This is expected, because for a fixed charge, you can only add particles pair-wise. In fact, as will be shown in the next section, if one generates a Poissonian number of pairs, the relative variance will be ω M = 2 for all multiplicities. Because ω M varies from from two to unity in a canonical ensemble, as the system size increases from zero to infinity, the values of C 4 /C 2 and C 3 /C 1 stay below unity as shown in Fig. 1.
The fact that ω M is above unity for small systems (due to only being able to sample even numbers of charges) is known as canonical suppression. This lowers the mean multiplicity and raises the relative variance, ω M . As seen in Fig. 1 it also raises the ratio C 4 /C 2 relative to its value for a large system.

C. Decays
A system of uncorrelated neutral particles that decay to charged particles also leads to fixed zero net charge. But unlike the results for the equilibrated canonical ensemble above, this can lead to ratios of C 4 /C 2 > 1. If each neutral The relative variance, ωM , of the charged multiplicity distribution is shown for three cases for a system carrying only a single type of ±unit charge. For a neutral equilibrated system (red dashed line) in a canonical ensemble, ωM approaches unity, the Poissonian limit for higher mean multiplicities M . In the low multiplicity limit the multiplicity will be either M = 0 or M = 2, which gives ωM = 2. This is in contrast to a system where one has a Poissonian distribution of neutral particles which all decay into pairs (green dashed line) which gives ωM = 2 for all multiplicities, or if the neutral particles all decay to four charges, which gives ωM = 4. Lower panel: According to Eq. (33), which was derived in [50], the ratio C4/C2 is determined by ωM and the acceptance probability α. It is unity for ωM = 2, below unity for ωM < 2, and above unity for ωM > 2. This shows that charge conservation in an equilibrated system pushes C4/C2 below unity, whereas if an equilibrated system of neutral particles decays, and if the decay products do not reform into the resonances, the resulting ratio of C4/C2 is unity for two-body decays, or above unity for four-body decays.
particle decays into N bodies charged particles, where the charges are ±1, the charged particle multiplicity distributions become If the emission of neutrals is Poissonian, the result is simple, ω M = N bodies . Fig. 1 illustrates how this picture affects C 4 /C 2 . In this case, because ω depends only on N bodies and does not change with multiplicity or system size, the cumulants are also independent of multiplicity. For N bodies = 2 the cumulant ratios do not even depend on the efficiency. For N bodies > 2 the ratio C 4 /C 2 exceeds unity, so if charge creation proceeded through the creation and decay of neutral clusters it would be easy to generate large values of C 4 /C 2 . This has been discussed at length in [71]. In an equilibrated system, decays and recombination have equal rates. However, as the system decouples recombination stops and decays proceed until only stable hadrons remain. During the hadronic phase the number of charged particles nearly doubles due to these decays. Thus, if the systems do equilibrate, then decay after chemical freeze-out, one would expect the ratio C 4 /C 2 to lie somewhere between the value for the canonical ensemble in Fig. 1 and the value for pure decays with N bodies = 2. STAR's experimental results for net proton fluctuations indeed satisfy this expectation, but their measured C 4 /C 2 for net charge fluctuations exceed unity, which can only be attained if N bodies > 2. Very few hadronic decays proceed via more than one charged pair, but one could have decays of clusters or hot spots. Other possible explanations for having C 4 /C 2 > 1 are volume fluctuations, which were discussed previously, or Bose symmetrization, which will be discussed ahead.

D. Bose correlations
It has long been understood that Bose correlations induce super-Poissonian fluctuations [62,63]. Here, we illustrate how Bose effects combine with charge conservation to determine C 4 /C 2 . Partition functions for a gas of positive and negative pions, m π = 139.57 MeV/c 2 , were considered to be kinetically equilibrated but with a chemical potential enforcing a fixed average density. If a gas is created at chemical equilibrium one expects µ = 0, but if it cools while maintaining a fixed number of pions, and if the pion number is fed by decays, an effective chemical potential, µ, is required. This changes the average number of pions by adjusting the phase space density as Because the system is out of equilibrium, and because the number of positives and negatives are nearly the same, the effective chemical potential has the same sign for both π + and π − . At decoupling one expects kinetic temperatures to fall near 100 MeV and the pion chemical potential to grow to perhaps as high as 75 MeV [61]. This estimate can be understood by the fact that the phase space occupancies should stay roughly constant for fixed entropy per particle for an expansion at fixed entropy, which suggests a chemical potential of approximately 50 MeV by the time the system cools to 100 MeV. Decay products further feed the phase space occupancy. Pion condensation occurs when this chemical potential reaches the pion mass. In the absence of Bose effects, this requires roughly doubling the phase space density of pions as compared to the expectations above. If such conditions were realized Bose effects could result in super-radiance [60], which should be accompanied by large multiplicity fluctuations. As a function of the effective chemical potential µ, the partition function for a pion gas in a volume of 500 fm 3 and at a temperature of 100 MeV was calculated from Eq. (27). For µ in the range of of 75 MeV the relative variance of the charged multiplicity distribution is ω M ≈ 1.2, which modestly increases C 4 /C 2 as shown in Fig. 2 for a fixed efficiency of α = 0.3. More dramatic results for C 4 /C 2 are not expected unless the chemical potential reaches within a few MeV of the pion mass. As discussed above, this is not expected. However, if the number of pion sources fluctuated wildly from one event to another, and if there were some events with twice the number of sources emitting into the same phase space, super-radiance might occur in some small fraction of the events. Such behavior would strongly contradict expectations based on chemical equilibrium.

E. Hadron Gas
Thus far, all the simple examples presented in this section considered a system with one type of conserved charge, but in actuality a hadron gas obeys the conservation of three types of charge: baryon number, strangeness and electric charge. This invalidates the use of the recurrence relation of Eq. (33) and requires the application of Eq. (16), or if Bose corrections are included, Eq. (27). Here, we calculate the canonical ensemble, Z A (Q, B, S), using the recurrence relation, then apply the Monte Carlo techniques of Sec. IV to generate sample sets of particles. The calculation is based on a large number, ∼ 300, of hadron resonances listed in [72]. After generating the particles, unstable resonances are decayed. For this section, a uniform efficiency α is assumed, independent of species, momenta, or whether the products came from a weak decay (charged pions and kaons were not decayed).
Emission is assumed to come from sub-volumes, or patches, of fixed size. Because each sub-volume is independent, there are no correlations between the various sub-volumes. Also, because cumulants scale linearly with the number of sub-volumes, ratios of cumulants can depend only on the patch volume, not on the overall volume. However, if the number of patches fluctuates, the result is modified by volume fluctuations as shown in Sec. VI A.
Here, the sensitivity to three parameters is investigated. The three varied parameters are the patch volume Ω, the baryon density ρ B , and the fixed efficiency α. When one parameter is varied, the other two are set at default values. The defaults are set at T = 150 MeV, ρ B = 0, α = 0.3 and the default patch volume is Ω = 200 fm 3 . The electric charge density is set to half of the given baryon density. After these sensitivities are studied, calculations at finite baryon density are modified so that the net baryon and electric charges responsible for the non-zero charge density are allowed to fluctuate across sub-volumes according to a Poissonian distribution. This is motivated by the fact that charges transported far away from mid-rapidity at early times cause this non-zero baryon density, so local charge conservation should play little role in the distribution of net charge across patches.
The dependence on patch size is exhibited in Fig. 3. For small patches the emission of baryons is discouraged because each baryon must be accompanied by the existence of an anti-baryon. Thus, the Boltzmann factor, e −M/T , for the mass is compounded by the necessity of having a second accompanying massive particle. Once the volumes become larger, and the mean number of baryon pairs exceeds unity, the observation of a baryon might be accompanied by seeing one fewer baryon in the remainder of the system, rather than seeing an extra anti-baryon. In fact, this is exactly what happens in a canonical ensemble in the large volume limit. In that limit, if a baryon is observed in some small amount of phase space, the mean number of baryons observed in the remainder of the system is decreased by one half, while the mean number of anti-baryons is increased by one half. It is expected that the dependence of cumulant ratios on patch size should disappear as the patch size increases beyond the threshold for the mean baryon density to be unity. Figure 3 displays this behavior by displaying the ratio C 4 /C 2 as a function of patch volume. Because the mean multiplicity for charged particles exceeds unity at a smaller volume, the ratio levels off faster for net charge than for net protons. Qualitatively, the same behavior was seen for a system with a single type of charge in Fig. 1.
The sensitivity to the acceptance probability, α, is displayed in Fig. 4. As α approaches zero, distributions tend to become Poissonian, because the number of particles measured is dominantly either zero or unity. Eq. (33) describes similar behavior for a system of a single type of charge. For net charge, the ratio is symmetric about α = 0.5. This is because when net charge is conserved, any charges not measured fluctuate exactly opposite to those that are measured. The same property is apparent in Eq. (33). However, the distribution of net protons considers only protons and anti-protons, so it does not represent the entirety of particles that carry baryon number or electric charge. Thus, that distribution is not symmetric about α = 0.5 and fluctuations are large even for α → 1. Figure 5 presents C 4 /C 2 and C 3 /C 1 as a function of baryon density. The patch volume is fixed at 200 fm 3 and the densities chosen correspond to a fixed number, 4, 8, 12, · · · , of baryons. The choice is made to display C 3 /C 1 rather than C 3 /C 2 because C 3 /C 1 would be unity for uncorrelated emission. The ratios all decrease, moving further from unity, as baryon density increases. The densities displayed cover the range of what might be reached in heavy-ion systems at the point when temperatures fall below 150 MeV.
For the calculations illustrated in Fig. 5, the net baryon number, electric charge, and strangeness were all fixed. Because of local charge conservation and the limits of diffusion, balancing charges created after the collision are constrained to stay within the same neighborhood. This constraint is adjusted by setting the size of the patch volume. However, even for measurements at mid-rapidity, non-zero charges can arise due to the transfer of baryon number and electric charge from the projectile and target rapidities. These intruder charges are not typically balanced by charges within the acceptance of the detector, and thus their fluctuation should be considered separately. Figure  6 shows how C 4 /C 2 and C 3 /C 1 behave as a function of baryon density, if the non-zero baryon number and electric As the fixed acceptance probability α approaches zero, distributions become Poissonian, and the ratios C4/C2 approach unity. Even for perfect acceptance, the fluctuations are non-zero for a multi-charge system because the conserved charges carried by a proton can be balanced by an array of other species, and charge within the proton and anti-proton sector is not conserved. Thus, the ratio is not symmetric about α = 0.5 as one might infer from Eq. (33), because protons and anti-protons represent only a fraction of the species that carry electric charge and baryon number.
charge are chosen to fluctuate within the patch according to a Poissonian distribution. As expected, the ratios rise as compared to the fixed case shown in Fig. 5, but they do not exceed unity. Discerning these fluctuations of the actual net charges is of particular interest. These results suggest that the starting point for the ratios C 4 /C 2 and C 3 /C 1 is below unity for random uncorrelated net charges accompanied by an ensemble of baryon charges.

VII. BLAST WAVE MODEL WITH A FULL HADRON GAS AND COMPARISON TO STAR RESULTS
The calculations of the previous section were based on a simple picture, where each sub-volume emitted particles whose probability of being observed was uniform, denoted by α. In practice, this probability depends on where the sub-volume is located within the overall reaction volume. A sub-volume in the region with spatial rapidity, |η s | > 1, emits particles that are only rarely recorded by detectors that specialize in mid-rapidity measurements. Even for a sub-volume centered at mid-rapidity, thermal motion can spread some of its emission to rapidities outside the acceptance. This is especially true when the experiments narrow their acceptance. For example, STAR's acceptance nominally covers ±1 units of pseudo-rapidity, but the acceptance for identified particles is confined to ±0.9 units. For real rapidity, the effective acceptance is narrowed further, due to the fact that the true rapidity y is less than 0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28 FIG. 6. The ratios C4/C2 and C3/C1 are shown for net charge and for net baryon number as a function of the baryon density. Calculations differ from those in Fig. 5 in that the net charge fluctuates according to a Poissonian distribution. The ratios rise relative to those in the fixed-charge case, but they remain below unity. the pseudo-rapidity η. This difference is magnified for more massive or slower particles. In fact, STAR analyses of identified particles often enforce cuts that only consider particles with true rapidities −0.5 < y < 0.5. Of course, even particles within the rapidity acceptance must exceed some minimum transverse momentum, and the efficiency for being detected is imperfect. For particles identified only by charge, the efficiencies are typically 80%, and for identified particles the efficiency falls by another few tens of percent.
Collective flow plays a critical role. First, longitudinal flow is what allows the measurement of rapidity to correlate to a measurement in coordinate-space rapidity. In the Bjorken model of boost-invariant flow, the mapping is simple in that a fluid element with spatial rapidity η s moves with a rapidity η s , and particles emitted from a sub-volume with rapidity η s have rapidities y ≈ η s . This simple mapping is smeared by thermal motion. Collective radial flow and cooling combine to better align the spatial rapidity η s and the measured rapidity y. The thermal spread for pions is ≈ 0.6 units of rapidity, while that for protons is ≈ 0.25.
The extent of the region over which particles are emitted in spatial rapidity, η s , should affect the result. If the extent in η s is small, there is an enhanced probability that a charge and its balancing charge will both be identified and cancel one another when assigning the net charge for an event. If the region extends over a large rapidity range, the effects of charge balance are minimized because there is a better chance that one charge will be observed while the other is outside the acceptance. The extent of the region in the transverse direction is less important. For emitting regions at the edge of the fireball, which have more collective velocity, there is a modest increase of having balancing charges both pushed into the acceptance. In addition to the size of the region over which particles from a given sub-volume are emitted, the overall size of the sub-volume at the point where chemical freeze-out occurs also plays a role as it sets the degree to which canonical suppression affects the results. As discussed in the previous section, this matters only for for smaller sub-volumes.
For the model in this section, canonical sub-volumes are overlaid onto a blast-wave parametrization of collective flow. A filter is applied to the calculations, representing the experimental STAR acceptance and efficiency. This should be sufficiently realistic to make meaningful comparisons to STAR data. The four main parameters for the blast-wave model describe the radial flow u ⊥ , the kinetic freeze-out temperature, T k , and the baryon chemical potential and temperature at chemical freeze-out, µ c and T c . The chemical freeze-out temperature µ c and chemical potential, T c , are chosen to fit relative particle yields, while T k and u ⊥ are determined by simultaneously fitting the spectra of species with varying mass, typically π, K, p. A variety of parameterizations exist, such as having the velocity increase linearly from the origin, or the transverse rapidity, or having a sharp cutoff in radius vs. assuming a Gaussian profile. Depending on the choice, the value of T k and u ⊥ vary, but are typically in the neighborhood of T k ≈ 100 MeV and u ⊥ ≈ 0.6. For increasing multiplicities, the reaction volumes can expand and cool further, which leads to modestly increased values of u ⊥ and modestly decreased values of T k for either more central or for more energetic collisions. For this study the chemical freeze-out parameters were taken from [73,74], which extracted T c and µ c for a variety of beam energies. For this section, Bose statistics were ignored. Over 300 hadron species were included in the analysis. The collective flow velocities were chosen to reproduce the mean transverse momenta of both pions and protons.
The distribution of spatial rapidities η s over which all particles are emitted is chosen to be Gaussian, with a width σ 0 that depends on beam energy as where ±y beam are the rapidities of the incoming beams. This choice reproduces the rapidity widths measured at RHIC to roughly 10% accuracy [75]. The distribution of spatial rapidities for particles of a given sub-volume was also spread according to a Gaussian with a width σ η . The parameter σ η describes how far particles created from the same sub-volume may have separated by the time of emission. For increasingly larger values of σ η , the chance that any two observed particles are correlated decreases. The rapidity of a sub-volume was then distributed according to a Gaussian with width σ 2 0 − σ 2 η , so that the emission overall was characterized by the width σ 0 . Whereas the other blast-wave parameters were chosen to describe spectra and yields, the parameter σ η is related to charge conservation. If charge is created early, and especially if the diffusion constant is large, the width might be close to one unit of spatial rapidity, whereas if all the charge were to come after hadronization, the width might more likely be a few tenths. Blast-wave models of charge balance functions, which are also driven by charge conservation and diffusion, suggest widths should be of the order of a third, but variations of a factor of two were not ruled out [51]. One goal of this section is to investigate the sensitivity to σ η .
For this blast-wave model, each sub-volume was also assigned a radial transverse velocity according to a Gaussian, and a small spread in u ⊥ . However, given that for this study, cuts are not being considered in transverse momenta or azimuthal angle, varying u ⊥ has little effect. Particles from each differential volume element were assigned momenta using the Monte Carlo algorithm described in [76]. The calculations of this section are filtered through a simplified model of the STAR detector's acceptance and efficiency. For unidentified particles, pseudo-rapidities are required to be between ±0.5 and transverse momenta are constrained to being above 200 MeV/c and below 2 GeV/c. For identified particles, rapidities were restricted to being between ±0.5, and transverse moment were required to be between 200 MeV/c and 1.6 GeV/c for pions or kaons, and between 400 MeV/c and 2 GeV/c for protons and antiprotons. Because the STAR data are corrected for efficiency, the efficiencies were set to unity. In order to compare results to measurement from STAR, the baryon densities and chemical freeze-out temperatures were mapped to beam energy according to the analysis of [74], which extracted chemical freeze-out temperatures and chemical potentials by considering ratios of particle yields. The parameter u ⊥ , which controls the transverse radial flow, was chosen along with the kinetic breakup temperature, T k , to simultaneously roughly fit the mean transverse momenta of both the protons and pions reported [77]. The kinetic freeze-out parameters were chosen according to the parameterization, where √ s 0 = 7.7 GeV and √ s f = 200 GeV. Decays were simulated, and the products of weak decays (aside from pions or charged kaons) were included in the analysis. Undoubtedly, a more realistic model of the acceptance might change the ratios, but given that these are ratios, and that the overall efficiency and acceptances are not wildly off, a more rigorous model of the acceptance is unlikely to change the result by more than a few percent.
The summation over sub-volumes was performed by summing over 400 values of u ⊥ and η s consistent with the distributions described above. For each value, emission of N sample = 2 × 10 5 sub-volumes was simulated so that cumulants could be calculated for each value of η s and u ⊥ . Because emission from different sub-volumes is uncorrelated, the cumulants for the total emission are simply the sum of the cumulants of each sub-volume. Further, because the net cumulants behave linearly in N sample , the ratios of net cumulants is independent of N sample . (color online) Ratios of cumulants from blast wave models, begining with C1/C2 (panels a -c), then skewness (panels d -f) and finally, kurtosis (panels g -i), are plotted for for different values of the sub-volume Ω. For smaller sub-volumes the ratios C1/C2, C3/C1 and C4/C2 increase. For net charge, the ratios approach an asymptotic value once Ω begins to pass ∼ 50 fm 3 , whereas for net protons the ratios appear to approach the limit at somewhat higher values of Ω. STAR measurements for net protons are not wholly dissimilar to the blast-wave calculations here, but those for net charge differ greatly. Additional physics from volume fluctuations might explain how C4/C2 might exceed unity, but it is difficult to explain how this might happen for net-charge distributions while leaving C4/C2 of the net-proton distribution unchanged. Figure 7 illustrates how results are sensitive to the size of the sub-volume, Ω. This is the effective volume at which charge conservation is enforced at chemical freeze-out. For these calculations σ η was fixed at 0.3. For a volume of 100 fm 3 the average number of hadrons is several dozen. For small sub-volumes, where in a grand canonical ensemble the typical number of charges would be zero or one, the thermodynamic cost of having a second charge to balance the first charge reduces the probability of having any charges. As mentioned earlier, this is known as canonical suppression, and it lowers both the multiplicities and the moments. For larger sub-volumes, this thermodynamic cost vanishes as the system might have had fewer balancing charges of the same sign, rather than an extra balancing charge of the opposite sign. The characteristic volume for the disappearance of canonical suppression is the volume where the mean number of pairs exceeds unity. In general, for charged particles this sets in at around 50 fm 3 , but for baryons the characteristic volume is closer to 100 fm 3 because of their being heavier and fewer. Thus, the proton moments for volumes of 50 fm 3 and 200 fm 3 differ more noticeably. The time for a fluid element to expand and cool to the point where it reaches chemical freeze-out tends to be on the order of 5 fm/c. These times are shorter for matter on the periphery, at a lower beam energy or centrality. The times are longer for fluid elements at the center, at higher beam energy, or in a more central collision. The maximum transverse distance a charge can travel before chemical freeze-out is ≈ 10 fm if they move in opposite directions, and if charge moves diffusively or if the charge is created later in the reaction, the separation should be significantly less. As discussed below, charge can separate further along the beam axis, and because that is not well understood, the size of the sub-volume carries a large uncertainty. Anywhere from 50 fm 3 to a few hundred fm 3 might be reasonable.
Balancing charges separate from one another, and depending on when they are created, they can diffuse apart from one another. This separation, represented by the parameter σ η , accounts for the separation of balancing charges both before and after hadronization, or both before and after chemical freezeout. After chemical freeze-out charges might separate and mix between the sub-volumes. Because some of the separation comes after chemical freeze-out, this distance might exceed the scales representing the size of the sub-volume Ω. Charge can spread further in the longitudinal direction due to the strong initial longitudinal collective flow at early times. This enhancement to the separation depends sensitively on when charge pairs are created. Matter thermalizes at an early time, where large collective velocity gradients along the beam axis are expected. If the motion is diffusive, the separation along the beam axis depends logarithmically on the ratio of the final time to the initial creation time [42]. Thus, if a pair is created at 0.2 fm/c, the separation for times, 0.2 < τ < 1.0 fm/c is as important as the additional separation they gain during the times 1.0 < τ < 5.0 fm/c. For this reason the size of charge spread, σ η , in spatial rapidity might be anywhere between a few tenths of a unit of rapidity to a full unit. Figure 8 shows the sensitivity of the moments to this parameter. For large σ η the observation of a charge is less likely to influence the observation of a second charge, which is similar to having a lower efficiency. For low efficiencies one expects the behavior to be more Poissonian, and C 4 /C 2 and C 3 /C 1 to be closer to unity. Indeed, this is the case, but the dependence on σ η , as shown in Fig. 8, is negligible for net charge and modest for net protons.
The beam energy dependence mainly derives from the fact that the baryon density is higher for lower beam energies. Given that C 1 falls with increasing beam energy, while C 2 increases because multiplicity increases, it is no surprise that the ratio C 1 /C 2 falls with increasing beam energy. The skewness, C 3 /C 1 , should approach unity if emission is random, i.e. a Skellam distribution. The effects of local charge conservation keep C 3 /C 1 < 1. In the limit that there are no anti-baryons, which is the limit of high net baryon density and equivalently low beam energy, the assumption that baryons are deposited amongst the sub-volumes according to a Poissonian distribution also drives C 3 /C 1 for net protons closer to unity. Combined with higher canonical suppression for heavier particles, the sensitivity to Ω, shown in Fig. 7, or to σ η , illustrated in Fig. 8, is more pronounced for blast-wave calculations of net protons than of net kaons or of net charge.
The measure of kurtosis, C 4 /C 2 , varies only modestly with beam energy. For net charge, the ratio does not vary far from 0.8 in model calculations. In contrast, C 4 /C 2 from model calculations for net protons varies from 0.6 to 0.9 depending on the values of Ω and σ η . The models also exhibit a modest dependence on beam energy for C 4 /C 2 . The ratio rises approximately 20 percent as the beam energy increases to 20 GeV, then plateaus and falls 10% until it becomes flat. Similar to the C 3 /C 1 ratio, this ratio stays below the Skellam limit. Figures 7 and 8 also display results from the STAR Collaboration for data recorded for central collisions, 0-5% centrality [10,20,21]. For net proton fluctuations, the blast-wave model results are similar to STAR measurements for all three ratios. It is difficult to conclude the meaning of the solid agreement shown in the C 1 /C 2 ratios. This agreement covers net charge, net kaons and net protons, as long as the sub-volume Ω is of the order of 50 fm 3 or greater. Interpreting the meaning of this agreement is difficult because it is mainly driven by having the model correctly match the particle yields with multiplicity. Unlike the C 3 /C 1 and C 4 /C 2 ratios there is no simple baseline for a Skellam distribution. The Skellam baseline would depend on knowing the charged particle multiplicity distributions, i.e. the moments of the multiplicity distributions for protons plus antiprotons, kaons plus antikoans, or all charged particles.
For the net-proton distribution, the behavior of C 3 /C 1 as a function of beam energy seems consistent with the experimental uncertainties, but the statistical errors in the experimental data are so large, that little can be concluded aside from the fact that the models predict that C 3 /C 1 should be in the neighborhood of 0.75. These ratios from the models can show a modest sensitivity to beam energy, but any such trends are overwhelmed by the statistical errors of the experimental results at the moment.
For both net kaons and for net charge, C 3 /C 1 and C 4 /C 2 lie above the range of model predictions, but the large experimental error bars forbid one from stating this with great confidence. Should the experimental results with improved statistics confirm this discrepancy, it will be difficult to explain unless a significant number of charges are emitted from large clusters. In contrast, the model should always give ratios below unity for C 4 /C 2 , regardless of the choice of parameters. These failures of the blast-wave model for net kaons and for net charge are discussed in the upcoming summary, Sec. VIII.

VIII. SUMMARY
The principal goals of this study were to clarify background contributions for higher moments of charge distributions measured by the STAR Collaboration at RHIC, and to state the degree to which current experimental results are either consistent or inconsistent with these contributions. By background, this refers to sources of fluctuations besides those that arise from baryon number or charge clustering due to processes such as phase separation. The list of such sources includes charge conservation, Bose corrections, volume fluctuations, and the decays of resonances. In order to gain better insight both simple semi-analytic models with a single type of conserved charge, similar to the work performed in [50], and a more realistic blast-wave model which includes a more realistic accounting of the STAR acceptance, similar to what was applied in [54], were investigated. By using a highly efficient algorithm for Monte Carlo generation of particles according to the canonical ensemble, results were produced with small statistical uncertainties. This enabled the exploration of sensitivities to critical parameters of the model.
Of the various background correlations, charge conservation provided the strongest non-Poissonian contributions. Because fourth-order cumulants were defined to subtract contributions from second-order correlations, one might have expected a small contribution. Consistent with the result from [50] for a single type of charge, it was found that generating sets of particles from a sub-volume equilibrated according to the canonical ensemble produced values of C 3 /C 1 and C 4 /C 2 which were significantly lower than the Skellam value of unity, which is what one would expect for uncorrelated emission. In contrast, the contribution from two-particle decays does not change either C 3 /C 1 or C 4 /C 2 . This conclusion persisted for the more realistic blast-wave model which incorporated the conservation of all three charges and included a filter of the STAR acceptance. The correlation varied modestly according to the size of the canonical sub-volume for small sub-volumes due to canonical suppression. A modest sensitivity was also found to the spatial extent of this volume along the beam axis, as the overlay of collective longitudinal flow onto the finite acceptance in rapidity effectively lowers the probability for two balancing charges to both be observed. Due to the fact that baryon density falls with increasing beam energy, the strength of such correlations do depend on beam energy. But this sensitivity was not dramatic. Thus, even though this background contribution is rather large, it is quite smooth with respect to beam energy, so if sharp non-monotonic structures are observed experimentally with respect to beam energy, such structures are not likely to be driven by charge conservation.
A second source of background arises from multi-particle symmetrization of the outgoing pions. By extending the recursive techniques applied for the canonical ensemble to include symmetrization, it was found that such effects should not affect the skewness or kurtosis unless the pion phase space density were to become surprisingly large. In order for the symmetrization effects to become large, the pion phase space density in the absence of symmetrization would have to double. If this were the case, the pion spectra would be more dramatically altered by Bose effects and the measured HBT radii would have to be significantly altered.
Volume fluctuations are somewhat of a wildcard for background processes. As shown in Sec. VI A, such fluctuations can significantly increase the C 4 /C 2 ratio. The STAR Collaboration invested great effort in minimizing their impact, but it is difficult to gauge the degree to which such effects might have persisted. Volume fluctuations also increase similar moments of the multiplicity distribution, which is constructed by counting charged particles rather than net charge. It is critical for the experiments to simultaneously present moments of the multiplicity distribution alongside those of the net-charge distributions. This sensitivity was also illustrated with the simple results from a system with one charge and uniform acceptance, restated in Eq. (33) from the previous work of [50]. Because volume fluctuations should similarly increase the C 4 /C 2 ratio for net charge, net strangeness and net baryon number, behavior of such ratios for one type of charge that are not seen in other types of charge can be considered as originating from some other effect.
Due to the way in which cumulants are constructed, the third and fourth-order cumulants should be impervious to the effects of two-particle decays. However, decays of clusters that produce four or more charged particles do contribute to C 4 . Here, it was found that if a significant fraction of charged particles come from such decays, the higher moments can be profoundly altered, as was already known from the work in [71]. However, the decays of such clusters would also affect the multiplicity distribution, which again underscores the importance of the experiments to simultaneously analyze fluctuation of net charge and of the the multiplicity distribution. If such effects were important, it would suggest novel contributions to the dynamics of charge production, outside of the usual paradigm of creating equilibrated distributions of hadrons.
The results of the blast-wave calculations were displayed alongside STAR results in Sec. VII. The experimental results had much larger statistical errors than the calculations, which limits the conclusions that can be drawn. The fluctuations of net protons were not far from the range of those calculated here. This is consistent with charge conservation being the dominant source of non-Poissonian behavior, i.e. C 4 /C 2 = 1 and C 3 /C 1 = 1. By no means does this suggest that this is evidence for a lack of more novel sources of correlation in the baryons, such as that arising from phase separation. The experimental uncertainties may be currently too large to unmask such phenomena.
The observed moments of the net-charge distribution are perplexing. Although there are large statistical uncertainties in the experimental data, it appears that C 4 /C 2 exceeds expectations of the blast-wave calculations, and even lies above unity. Given that the net-proton distributions are roughly in line with the model predictions, this discrepancy cannot be explained by volume fluctuations, or equivalently by the systematics of event binning. Because phase transition phenomena are expected to manifest themselves mainly in the net-proton distributions, this would suggest that the decays of larger clusters into four or more charged particles might be present. If such processes were present, it would motivate a rethinking of models of chemical evolution and charge production in heavy-ion collisions. In order to confirm or dismiss this hypothesis, it is imperative that a simultaneous analysis of charge (not net charge) multiplicity distributions be undertaken from the same data sets with the same cuts on centrality. Given the greatly improved data sets currently being analyzed by STAR in the Beam Energy Scan II program at RHIC [78], this puzzle should be clarified in the next year or two.
Finally, the studies here help point the way to future improvements in modeling. The picture of independent canonical sub-volumes is crude. It does incorporate the truth that charge conservation is enforced locally, over some length scale, and is sufficient to provide the understanding of how large such effects might be. However, in reality baryon number, electric charge and strangeness are created and evolve in different ways. Strangeness tends to be created early in the collision, thus allowing the balancing strange and anti-strange quarks to separate before the emission of the hadrons to which they are asymptotically assigned. In contrast, up and down quarks are more likely to be produced later in the reaction. Thus, the characteristic canonical volumes should have different sizes and different longitudinal extents depending on whether one is considering up, down or strange quarks. This is also true for off-diagonal correlations, e.g. correlations between strange and up. These only appear in the hadronic phase. A diagrammatic formalism has been developed for evolving such two-, three-, and four-body correlations as a function of the positions of each charge. These equations are based on knowing the chemical evolution of the system and the diffusion constants for each type of charge. Such calculations have been performed for two-body correlations and compared to experimentally measured charge-balance functions [79,80] or fluctuations [27]. Unfortunately, the method for three-and four-body correlations is challenging to implement numerically [81]. It is tractable, but would require significant effort. The study presented here suggests that such calculations would be warranted if one wanted to reproduce the moments of these distributions for each type of charge, and especially if one wants to consider cross terms [82], such as moments involving powers of both charge and net baryon number. Otherwise, given that charge conservation effects are expected to evolve smoothly with beam energy, one could simply see whether measurements of the ongoing Beam Energy Scan at RHIC unveil any sharp features, and assign such features to more novel types of physics.