Multi-particle Integral and Differential Correlation Functions

This paper formalizes the use of integral and differential cumulants for measurements of multi-particle event-by-event transverse momentum fluctuations, rapidity fluctuations, as well as net charge fluctuations. This enables the introduction of multi-particle balance functions, defined based on differential correlation functions (factorial cumulants), that suppress two and three prong resonance decays effects and enable measurements of underlying long range correlations obeying quantum number conservation constraints. These multi-particle balance functions satisfy simple sum rules determined by quantum number conservation. It is additionally shown that these multi-particle balance functions arise as an intrinsic component of high-order net charge cumulants. This implies that the magnitude of these cumulants, measured in a specific experimental acceptance, is strictly constrained by charge conservation and primarily determined by the rapidity and momentum width of these balance functions. The paper also presents techniques to reduce the computation time of differential correlation functions up to order $n=$10 based on the methods of moments.


I. INTRODUCTION
A variety of two-particle integral and differential correlation functions have been developed and deployed towards the analysis of particle production in heavy-ion collisions [1][2][3][4].Correlation functions formulated as functions of rapidity and azimuth angle differences have enabled, in particular, the discovery of away-side jet suppression in large collision systems [5,6].Subsequent measurements of long range particle correlations in both large and small collision systems have additionally enabled detailed studies of the collision dynamics.Of particular interest are two-particle number correlations, R 2 , based on normalized two-particle differential cumulants [3,4], transverse momentum (p T ) correlation functions, P 2 , designed to study p T fluctuations [7], and G 2 , designed for the study of viscous effects based on their longitudinal broadening in A-A collisions [8][9][10].Differential two-particle correlation functions have been studied for charge inclusive (CI) particle pairs as well as charge dependent (CD) pairs.The latter are of particular interest because they relate to charge balance functions, which provide a tool to study the correlation length, in momentum space, of charge balancing particles.Balance functions were initially developed to investigate the presence of isentropic expansion and delayed hadronization in large collision systems [11][12][13].However, balance functions were later shown to also be sensitive to quark susceptibilities [14] as well as the diffusivity of light quarks [15,16].Recent works have shown they are best measured in terms of differential two-particle cumulants and feature simple sum-rule that might be instrumental in the study of thermal hadron production [17][18][19].
Two-particle correlation functions are by construction dominated by contributions from two particle correlated production processes such as hadronic decays and jet production.They are also nominally sensitive to higher order particle correlations even though they cannot specifically discriminate such higher order processes.In many instances, it is these higher order correlation processes that are of interest to investigate particle production and properties of the matter produced in A-A collisions.As a first example, consider measurements of transverse momentum fluctuation deviates.Correlators of transverse momentum deviates were initially invoked as a proxy to study temperature fluctuations in A-A collisions [20].However, it may be argued that fluctuations observed in A-A collisions are strictly driven by initial state conditions [21].There is thus interest in establishing the magnitude of such collective fluctuations.The use of the multiple particle p T integral correlator has then been advocated to study the effects of initial fluctuations [22,23].However note that ongoing studies may suffer from short-range correlations (a.k.a.non-flow) and it might thus be desirable to develop differential analysis techniques that enable rapidity gaps designed to suppress short range correlations.As a second example of interest, consider charge, strangeness, and baryon balance functions.In this case also, one expects the correlations to receive sizable contributions from two and three prong particle decays as well as short range correlation from jets.There is thus an interest in obtaining higher order correlation functions that are sensitive to charge (strangeness, baryon) balance but suppress contributions from decays and jets and it is the primary purpose of this work to develop such multi-particle balance functions.As in studies of anisoptropic flow measurements, multi-particle cumulants, of order n=4, 6, . .., shall be used to suppress lower order correlations and obtain multi-particle balance functions.Unfortunately, measurements and calculations of higher order cumulants nominally require several nested loops to include contributions from all n-tuplets of interest on an event-by-event basis.Although such calculations are conceptually simple and remain practical for collisions and experimental acceptances featuring a modest particle multiplicity, they become prohibitively CPU intensive for large multiplicities.Fortunately, a number of techniques have been developed, particularly in the context of anisotropic flow analyses, to reduce the computational challenge of handling large multiplicities and high-order cumulants [24].One such application successfully deployed in the context of anisotropic flow studies and towards the computation of higher moment deviates relies of the method of moments [22].A second purpose of this paper is to further develop this method towards measurements of differential multi-particle correlations of inclusive as well as identified particle species.
This paper is organized as follows.First, sec.II presents a discussion detailing the need for higher order integral and differential correlators based on several examples.Higher order correlators are then introduced in sec.III in terms of integral and differential cumulants of arbitrary order as well as expectation values of the form ⟪q 1 q 2 ⋯q n ⟫ and ⟪∆q 1 ∆q 2 ⋯∆q n ⟫, where q i , i = 1, . . ., n represent particle variables of interest (e.g., transverse momentum, charge, etc) and ∆q i = q i − ⟪q⟫, i = 1, . . ., n, are their deviates relative to their event ensemble average ⟪q⟫.Techniques to compute these expectation values based on the method of moments are presented in secs.III C and III D. Equipped with these different correlators and computing tools, the notion of multi-particle balance function is then introduced in sec.IV.Finally, multi-particle balance functions are considered in sec.V in the context of measurements of net-charge cumulants and it is shown that net charge cumulants of all orders are explicitly constrained by charge conservation.The paper is summarized in sec.VI.Given much of the calculations performed in this work are somewhat tedious and lengthy, details of these calculations are presented in several appendices.Appendix A presents calculation methods and formula for moments, cumulants, factorial moments, and factorial cumulants for single and multi-variable systems, as well as for the net-charge of particles measured in a specific acceptance Ω. Appendix B extends these formula to differential correlations.General formula for the computations of deviates of the form ⟪∆q 1 ∆q 2 ⋯ ∆q m ⟫ are derived in appendix C while equations for the calculation of correlators of the form ⟪q 1 q 2 ⋯ q m p 1 ⋯ p n ⟫ are listed in appendix.D. Finally, appendix E lists definitions of multi-particle balance functions up to order n = 10.

II. MOTIVATIONS
Let us consider single particle densities of particles of type α, denoted ρ α 1 (⃗ p), and n-particle densities of mixed species α 1 , . . ., α n , denoted ρ α1⋯αn n (⃗ p 1 , . . ., ⃗ p n ).In general, mixed n-particle densities, ρ α1⋯αn n (⃗ p 1 , . . ., ⃗ p n ), correspond to the yield (per event) of n-tuplets of particles of types α 1 , α 2 , . . ., α n at momenta ⃗ p 1 , . . ., ⃗ p n .Such n-tuplets may arise from a single process yielding n correlated (mixed) particles, or combinations of processes jointly yielding n-particles.In order to focus a study on correlated particles exclusively, one commonly relies on the notions of integral and differential correlation function cumulants.Indeed recall that, by construction, integration of ρ α 1 (⃗ p) over a specific kinematic range Ω yields the average number of particles of type α in this acceptance, whereas integration of two-, three-, or n mixed particle densities yield the average number of pairs, triplets, and more generally n-tuplets of such groupings of particles.These integrals do not discriminate correlated from uncorrelated particles.
Differential cumulants and their integrals are of particular interest because they identically vanish in the absence of n or more particle correlations.They are thus an essential tool for the study of particle production.They can also be straightforwardly corrected for uncorrelated particle losses (detection efficiency).This feature is exploited in the formulation of number correlation function ratios such as with C 2 the second order cumulant and F 1 and F 2 the first and second order factorial cumulants, which is said to be robust against particle losses, i.e., independent of efficiencies provided these are approximately constant within the acceptance of a measurement [25].Ratios of factorial cumulants, similarly formulated, have also been used in recent studies [2,4].They too feature the property of robustness against particle (efficiency) losses.Particular combinations of differential and integral ratios R αβ 2 have also been used in the context of relative yield fluctuation studies and balance functions [26].
Integral and differential correlation functions have also been used to study fluctuations of specific kinematic variables.Fluctuations of event-wise total transverse momentum, in particular, have been proposed to study temperature and energy fluctuations in the initial stage of heavy-ion collisions.Voloshin et al. [20] showed fluctuations measures are best formulated in terms of transverse momentum deviates ∆p T,1 ∆p T,2 according to This integral correlator is commonly reported in terms of a dimensionless ratio ⟪∆p T,1 ∆p T,2 ⟫/⟪p T ⟫ 2 [27,28].
A differential version of this dimensionless correlator was also used [2,4,7,29] and can be written according to A generalization of ⟪∆p T,1 ∆p T,2 ⟫ to four particle correlations was first proposed by Voloshin [27] towards the study of temperature and energy fluctuations.The study of third and fourth p T moments, defined according to was proposed to probe initial stage fluctuations [22].Measurements of these correlators were recently reported by the ALICE collaboration [23].Clearly, it is trivial to also consider differential versions of these two correlators and such generalizations of P 2 (y 1 , φ 1 , y 2 , φ 2 ) might be useful to carry out higher moment analyses with finite rapidity gaps.Additionally, as we discuss in the next sections, extension to particle correlators of this form to n > 4 particles are readily accessible based on the methods of moments presented in sec.III C. The integral and differential correlation functions, Eqs.(2)(3)(4)(5), may also be applicable to the study of other types of fluctuations.For instance, replacing ∆p T,i by rapidity deviates ∆y i ≡ y i − ⟪y⟫, where y i are the rapidities of particles i = 1, . . ., N of an event, it becomes possible to study event-by-event fluctuations in the rapidity of particles.Such fluctuations might provide an alternative way to probe the longitudinal correlation length (rapidity) of produced particles.
The multi-particle correlators and the set of tools for their extraction presented in Sec.III provide a basis for the extension of former techniques used for measurements of transverse momentum correlations, rapidity correlations, as well as charge correlations (including baryon and strangeness numbers correlations) heretofore completed mostly at low orders n ≤ 4.These techniques also connect to measurements of anisotropic flow and correlations between flow and other variables discussed elsewhere [24].
The techniques, whether used with a single or several variables, are nominally very powerful because they enable joint measurements involving many particles simultaneously.However, it is also clear that the complexity of such measurements can quickly grow out of hand.For instance, assuming an interest in the rapidity, transverse momentum, and azimuth of particles, one would nominally get differential observables ⟪∆q 1 ∆q 2 ⋯∆q n ⟫ featuring 3 × n degrees of freedom.Measuring such features would then require collecting data in as many as 3 × n dimensions.Clearly, a considerable reduction of this "feature" space is required to enable the feasibility of measurements both in terms of data volumes (i.e., capturing sufficiently many n-tuples of particles to cover all partitions of the feature space with meaningful values) and in terms of its representation and interpretation.While we do not wish to preclude or dismiss possibilities of complex multidimensional analyses, we focus the discussion, as a kind of extended motivation, on some basic applications of the formalism and methods presented later in this work towards some specific physics analyses.
On general grounds, one can classify analyses of potential interest based on the number of kinematic partitions being used (i.e., partition of the 3 × n momentum space), the number of observables of interest (e.g., transverse momentum, p T , rapidity, y, charge, anisotropic flow coefficients, etc) and the number of particle types or species being considered (e.g., inclusive charged particles, positively vs. negatively charged particles, specific species such as pions, kaons, etc, and so on).We thus organize the discussion in terms of few use cases, beginning with the simplest case involving a single variable q, and next considering progressively more and more complex use cases involving two variables: q, p, as well as several variables.
Analyses based on a single variable q are already quite popular and have featured studies of fluctuations of transverse momentum, net-charge, etc.However, most prior analyses have been limited to two particles [7][8][9]29] and only few recent works have undertaken higher number of particles [23,27].Notable exceptions to this statement evidently include measurements of anisotropic flow based on multi-particle cumulants [30].
Measurements of p T (alternatively y or q, etc) correlations ⟪∆p T,1 ∆p T,2 ⋯∆p T,n ⟫ involving n ≥ 4 particles of a given type of particle in a specific acceptance can be readily undertaken based on the methods discussed in sec.III.Various types of scaling are evidently possible to obtain dimensionless observables and assess the evolution of n-order momentum correlators as a function of the A-A collision centrality or produced particle multiplicity.An obvious choice is the inclusive momentum average ⟪p T ⟫ already used in several studies [7][8][9] but other choices of scaling have also been discussed [22].Differential measurements can then be achieved, for instance, by studying the magnitude of ⟪∆p T,1 ∆p T,2 ⋯∆p T,n ⟫ vs. the width of the acceptance in rapidity.
Analyses involving two variables, q and p, are of interest, for instance, towards the study of some specific observable (e.g., anisotropic flow, transverse momentum fluctuations, or net-charge fluctuations) in two kinematic partitions separated by a finite size rapidity gap.
Considering examples of p T fluctuation studies, let q i and p i represent the transverse momentum of particles measured in two distinct rapidity acceptance ranges Ω A and Ω B of equal widths separated by a finite rapidity gap ∆η, as schematically illustrated in Fig. 1.One can then measure correlators of the form ⟪∆q 1 ∆p 1 ⟫, ⟪∆q 1 ∆q 2 ∆p 1 ∆p 2 ⟫, etc., at any order to determine the strength of n = 2, 4, etc, transverse momentum correlations as a function of the width of the rapidity gap.Measurements of transverse momentum fluctuations and correlations have been thus far mostly limited to two-particle studies and implemented with a single acceptance bin [31][32][33][34] or in a fully differential manner [8,9].The methods discused in Sec.III, however, enable differential measurements involving multiple n ≥ 4 particles.It then becomes possible to study momentum correlations arising from initial state fluctuations [22] while suppressing the influence of short range correlations (aka non-flow) associated with hadron decays and jet fragmentation.Letting q i , i = 1, . .., represent the p T of particles in partition A, and p i , i = 1, . . .represent the p T of particles in partition B, the described methods allow to obtain ⟪q 1 q 2 ⋯q n p 1 p 2 ⋯p n ⟫ and the n-order deviates ⟪∆q 1 ∆q 2 ⋯∆q n ∆p 1 ∆p 2 ⋯∆p n ⟫ corresponding to p T correlators involving n particles from partition A and n particles from partition B. As in flow studies, one expects that it becomes possible to progressively suppress resonance and jet contributions by increasing the rapidity gap ∆η between partitions A an B. The analysis can also be made more differential by also using bins in azimuth as illustrated in panels (c) and (d) of Fig. 1 thereby enabling the suppression of or focus on back-to-back jet contributions.
The methodology can readily be adopted also for measurements of charge correlations and, as it will be introduced, multi-particle balance functions.In this case, q i and p i represent the charge of particles in partitions A and B. Then generic charge correlators ⟪q 1 q 2 ⋯q n p 1 p 2 ⋯p n ⟫ and their deviates ⟪∆q 1 ∆q 2 ⋯∆q n ∆p 1 ∆p 2 ⋯∆p n ⟫ can be obtaind and, it will be seen, these can then be related to multiparticle balance functions.
Two other use cases based on two variables q i and p i are worth mentioning.One involves the study of two distinct physics observables (e.g., charge, p T , rapidity, etc) in a single kinematic partition whereas the other involves the measurement of a specific particle observable, e.g., the p T , for two types of particle species.In the first case, the variable q i and p i represent the two observables of interest whereas in the second they tag the species of interest.These latter use cases enable multi-particle correlations with specific species or between identical particles in two distinct p T ranges.
The examples discussed in the previous paragraphs are readily extended towards the computation of correlation functions involving three or more kinematic partitions and particle types.Of particular interest is the determination of multiple particle balance functions.Although it may not be practical to conduct analyses involving explicit computation of more than three or four kinematic partitions or species, it remains possible to consider balance functions involving large number of particles towards the study of long range multi-particle correlations constrained by charge conservation (or other quantum number conservation laws).Let us first consider measurements of four-particle balance functions based on two kinematic partitions A and B separated by a finite rapidity gap, as illustrated in panels (a) and (b) of Fig. 1.The partitions A and B could be azimuthally symmetric (i.e., with full azimuth coverage 0 ≤ φ < 2π), or feature partial coverage to suppress contributions from back-to-back jets, as schematically illustrated in panels (c) and (d) of the same figure.Panel (a) illustrates a measurement involving two positively charged particles in A and two negatively particles in B. A measurement of the multi-particle balance function shall then be sensitive to the strength (or probability) of processes featuring four correlated particles separated by a finite rapidity gap.Since 4-prong resonance decays are relatively few, this would reveal the likelihood of long range correlations determined by string-like fragmentation processes.In contrast, the analysis illustrated in panel (b) would focus on correlated quartets featuring two nearby pairs of unlike sign particles.These could be produced by string-like fragmentation processes yielding four or more correlated particles, but they could also result from string fragmentation producing two neutral objects, each decaying into pairs of +ve and −ve particles.An explicitly selection of the charge states to be measured in partitions A and B might thus enable a discriminating study of the relative yields of distinct processes.Indeed, an analysis of the dependence of the relative strengths of processes depicted in panel (a) and (b) could shed additional light on particle production process in elementary collisions.An analysis of the correlation strength performed as a function of the rapidity gap might then provide better sensitivity to the correlation length of string break up processes.Additionally, such analyses conducted as a function of collision centrality and beam energy in large systems (A-A), and comparisons to dependencies observed in small systems (e.g., pp and p-A), might then reveal whether this correlation length evolves with energy density, system size, collision energy, etc.Clearly, the position and size of measurement bins can be varied.Panel (c) illustrates a measurement geometry emphasizing back-to-back jet emission whereas panel (d) suppresses such processes and thus enables the study of long range non-jet and not resonance decay processes such as longitudinal string fragmentation.Obviously, a wide variety of other detection geometries can be implemented to explicitly favor or inhibit  In general, a large number of such particle species combinations could be combined to study charge (Q), strangeness (S), baryon number (B), and isospin I 3 balancing.The tools could also be applied to charmness (C) and/or bottomness (B) balancing.

III. MULTI-PARTICLE CORRELATION FUNCTIONS A. Cumulants and factorial cumulants
As already mentioned, in order to focus a study of correlated particles exclusively, one relies on the notions of integral and differential correlation function cumulants.Differential cumulants have been discussed in details elsewhere [24,30].In the context of this work, it suffices to remember they can be "reversed engineered" by listing all cluster decompositions of n-tuple densities.As such, single-, two-, and threeparticle mixed density cumulants may be obtained by writing [26] where the notation ∑ (k) stands for a sum over k ordered permutations of the labels α 1 , α 2 , and α 3 , as well as the corresponding momentum vectors ⃗ p 1 , ⃗ p 2 , and ⃗ p 3 .For n = 3, substitution of first and second order cumulants and expansion of the permutations yield Higher order cumulants are computed and listed in Appendix A.
Mixed cumulants C α1⋯αn n (⃗ p 1 , . . ., ⃗ p n ) nominally feature 3 × n degrees of freedom and are, as such, challenging to measure and visually represent.The dimensionality, and thus the number of degrees of freedom, can fortunately be reduced by integrating over several coordinates.When all coordinates are integrated within acceptances Ω k , k = 1, 2, . . ., n, one obtains mixed factorial moments.For instance, the integration of ρ α 1 (⃗ p) over a specific kinematic range Ω yields the average number of particles of type α in this acceptance, whereas integration of two-, three-, or n mixed particle densities yield the average number of pairs, triplets, and more generally n-tuplets of such groupings of particles.These averages are known as mixed factorial moments, herewith denoted f α1⋯αn n , and computed according to where ⟨N α ⟩, ⟨N α1 (N α2 − δ α1α2 )⟩, and so on, denote the ensemble average of the number of mixed tuplets of the corresponding order.Being integrals of densities, factorial moments do not discriminate correlated from uncorrelated particles and it is thus also convenient to consider mixed factorial moment cumulants, herein denoted F α1 ⋯ αn n .Mixed factorial moment cumulants, hereafter simply called factorial cumulants, are readily expressed in terms of integrals of the mixed cumulants C α1⋯αn n (⃗ p 1 , . . ., ⃗ p n ) over acceptances Ω k , k = 1, . . ., n, but can also be computed based on generating functions, as discussed in Appendix A. Lowest orders yield while formula for higher orders are also listed in Appendix A.

B. Joint Moments of Observables and their Deviates
As we discussed in Secs.II and III A, a wide variety of multi-particle correlation functions can be formulated in terms of the expectation value of products of particle observables of the form ⟪q 1 q 2 . . .q n ⟫ or their deviates ⟪∆q 1 ∆q 2 . . .∆q n ⟫.In this section, we first introduce such expectations values based on moments of the sum ∑ N i=0 q i computed for all selected particles of an event with all self-correlations removed.We then consider the expectation value of off-diagonal products of deviates of the form ⟪∆q 1 ∆q 2 ⋯ ∆q n ⟫.First note that these expressions are totally general and thus applicable for any type of particle observables, e.g., transverse momentum, rapidity, electric charge, or other quantum numbers.Together with formula introduced in sec.III C and Appendix C, these correlators enable the formulation of both integral and differential measurements of multiple particle correlations of basic particle observables.
In order to obtain expressions sought for, first consider a particle observable of interest q.This observable could be the transverse momentum of the particle, its rapidity, its charge, some other quantum numbers, or simply unity (i.e., to count the particles).We will denote the value of this observable for a specific particle q i , with the index i spanning all selected particles (i.e., satisfying specific kinematic and quality selection criteria) in a given event.We are interested in computing moments of q and its deviates ∆q i ≡ q i − ⟪q⟫, where ⟪q⟫ is the inclusive event ensemble average of q for specific collision conditions (i.e., events satisfying specific selection criteria).
Inclusive event ensemble (joint) averages of products of qs are defined according to where the sums run over all N selected particles in a given event.The notation ∑ i1≠i2 indicates the sums are computed for distinct particles, i.e., distinct values of i 1 , i 2 , etc.As such, ∑ N i=1 q i represents the sum of the q i in a given event, and ⟨∑ N i=1 q i ⟩ is the event ensemble average of this sum across all events of a selected data sample.Similarly, ⟨∑ N i≠j=1 q i q j ⟩, and higher orders, represent ensemble averages of products of the qs evaluated for n-tuplets of particles.However, the fact that sums proceeds on i ≠ j implies autocorrelations are explicitly removed.Herewith, we will use inclusive averages (i.e., averages computed over an event ensemble) but it is trivial to change the definition to event-by-event averages [22].See for instance Eqs.(D14 -D15).Additionally, note that if it is possible to pre-scan the dataset of interest to determine the inclusive average ⟪q⟫, then one can readily replace q i by ∆q i = q i − ⟪q⟫ in Eqs.(17)(18)(19) instead of carrying out the factorization discussed in details below.
We proceed to compute event ensemble averages of moments of ∆q i in terms of moments of q i .The inclusive first moment of ∆q is calculated according to and vanishes by construction.In order to compute moments of order n = 2 and n = 3, we write ∆q i ∆q j =q i q j − ⟪q⟫(q i + q j ) + ⟪q⟫ 2 , ( ∆q i ∆q j ∆q k =q i q j q k − ⟪q⟫ (q i q j + q i q k + q j q k ) + ⟪q⟫ 2 (q i + q j + q k ) − ⟪q⟫ 3 .
whereas higher order products can be written where the notation ∑ ( n k ) indicates a sum over all ( n k ) permutations of the indices i 1 , i 2 , etc. Ensemble averages of products of order n = 2, 3, 4 yield whereas higher orders can be computed according to where ⟪q i1 ⋯ q i k ⟫ = ⟪q⟫ for k = 1 and ⟪q i1 ⋯ q i k ⟫ = 1 for k = 0.
The computation of moments of order n of qs amounts to sums of the form ∑ N i≠j=1 q i q j , ∑ N i≠j≠k=1 q i q j q k , etc, which nominally require up to n nested loops for each event considered.Though conceptually trivial, such calculations involve a computation time proportional to N n that becomes quickly prohibitive for large values of the multiplicity N .Fortunately, the method of moments, which we discuss in the next two subsections, enables the computation of averages of these sums based on a single loop (per event), i.e., with a computation time proportional to N .

C. Method of Moments
The method of moments is introduced to facilitate the computation of moments ⟪q 1 q 2 . . .q n ⟫ and deviates ⟪∆q 1 ∆q 2 . . .∆q n ⟫ based on a single loop over all particles of an event rather than using nested loops.

Method of Moments for a single variable
Let Q n represent an event-wise sum of the nth power of the variable q i of the N selected particles of an event The method of moments relies on the evaluation of ensemble averages of Q n and products of the form and one obviously obtains Computation of the ensemble average of products etc, requires one properly handles the expansions of the sums.For instance, product of two and three Qs yield where the notation ∑ (3) represents a sum spanning all three ordered permutations of n 1 , n 2 , and n 3 .
Clearly, calculation of the ensemble average of This expression contains a term of the form ⟨N ⟩⟪q n1+n2 i1 ⟫ that can be readily replaced by ⟪Q n1+n2 ⟫ based on Eq. (30).Solving for ⟪q n1 1 q n2 2 ⟫, one then gets which for n 1 = n 2 = 1 evidently simplifies to Proceeding similarly for the ensemble average of The first term, ⟨N ⟩⟪q n1+n2+n3 i1 ⟫, is equal to ⟪Q n1+n2+n3 ⟫, while the next three terms are of the form of Eq. (36).Substituting these terms and solving for ⟪q n1 i1 q n2 i2 q n3 i3 ⟫, one gets The above calculation can be repeated iteratively for higher order products of qs and thus yield expressions for ⟪q 1 q 2 ⋯ q n ⟫ at arbitrarily high order n.In practice, such calculations become rather tedious for n > 4 and are best computed programmatically, as discussed in Appendix D.

Higher Order Moments of Mixed Acceptance Variates
In the previous sections, we considered calculations of higher moments ⟪q 1 q 2 ⋯q n ⟫ computed for a single acceptance or particle species (i.e., for a single kinematic bin or for a specific species or both).In order to compute differential correlation functions involving several kinematic bins or species, we now proceed to compute expectation values of the form ⟪q 1 q 2 ⋯q n p 1 ⋯p m r 1 ⋯r o ⟫ where q, p, and r represent distinct kinematic bins or species.The discussion is here limited to three bins (or species) for simplicity's sake but it is trivially extended to an arbitrary number of bins and species.The particle multiplicities in each bin are denoted N i , i = 1, . . ., 3. Deviates are denoted ∆q i , ∆p j , and ∆r k for particles in bins 1, 2, and 3, respectively.Moments for all particles detected in a single bin are given by expressions of the form (36 -39) already considered in sec.III C 1. We thus need to consider mixed moments only in this section.Lowest mixed moments are given by expressions of the form More generally, considering moments of order m 1 , m 2 , and m 3 for bins 1, 2, and 3, one gets expressions of the form where the sums and ∑ N3 k1≠⋯≠km 3 =1 span tuplets of distinct particles in bins 1, 2, and 3, respectively.Also note that the normalization corresponds to the average number of such tuplets: Event ensembles of mixed moments of deviates are defined in a similar fashion.At lowest orders, one gets expressions of the form and higher order moments are discussed in Appendix C. The computation of mixed moments and their deviates nominally requires nested loops over all particles and bins of interest.But as for the single bin case discussed in the previous section, it is advantageous to introduce event-wise sums of the variables q i , p i , and r i according to It then becomes possible, as discussed in Appendix C, to compute the event ensemble averages ⟪q 1 ⋯ q m1 p 1 ⋯ p m2 r 1 ⋯ r m3 ⟫ and their corresponding deviates ⟪∆q 1 ⋯ ∆q m1 ∆p 1 ⋯ ∆p m2 ∆r 1 ⋯ ∆r m3 ⟫ based on recursive formula of the moments of Q n , P n , and R n .

D. Differential Measurements of Multiple-Particle Correlations
The multi-particle correlators ⟪q 1 q 2 ⋯q n ⟫ and ⟪∆q 1 ∆q 2 ⋯∆q n ⟫ presented in sec.III A, equipped with the method of moments discussed in sec.III C, provide a basis for the extension of former techniques used for measurements of transverse momentum correlations, rapidity correlations, as well as charge correlations (including baryon and strangeness numbers correlations).
The method of moments, whether used with a single or several variables, is nominally very powerful because it enables joint measurements involving many particles simultaneously.However, it is also clear that the complexity of such measurements can quickly grow out of hand.
As was described in Sec.II, one can classify analyses of potential interest based on the number of kinematic bins being used (i.e., partitions of the 3× momentum space), the number of observables of interest (e.g., transverse momentum, p T , rapidity, y, charge q, anisotropic flow coefficients, etc) and the number of particle types or species being considered (e.g., inclusive charged particles, positively vs. negatively charged particles, specific species such as pions, kaons, etc, and so on).We thus organize the discussion of this section in terms of few use cases, beginning with the simplest case involving a single variable q, with event-wise variable Q n , and next considering progressively more complex use cases involving two variables: q, p with event-wise variables Q n and P n , as well as more complex analyses based on several variables.

One variable (q)
Measurements of p T (alternatively y or q, etc) correlations ⟪∆p T,1 ∆p T,2 ⋯∆p T,n ⟫ involving n ≥ 4 particles of a given type of particle in a specific acceptance can be readily undertaken based on the method of moments discussed in sec.III C. If it is possible or practical to carry out the analysis in two or more passes on the data, than one can use the first pass to determine ⟪p T ⟫.In the second pass, one can then define and compute Q n = ∑ N i (p T − ⟪p T ⟫) event by event and then use Eqs.(D4 -D10) to obtain n-oder moments ⟪∆p T,1 ∆p T,2 ⋯∆p T,n ⟫.If the determination of ⟪p T ⟫ in a first pass is not practical, then one can define and compute Q n = ∑ N i p T event by event, use Eqs.(D4 -D10) to obtain ⟪p T,1 p T,2 ⋯p T,n ⟫ and Eqs.(C1 -C5) or Eq.(C6) to obtain the n-order deviates ⟪∆p T,1 ∆p T,2 ⋯∆p T,n ⟫.

Two variables (q and p)
We first discuss examples of p T fluctuation studies.Let q i and p i represent the transverse momentum of particles measured in two distinct rapidity acceptance ranges Ω A and Ω B of equal widths separated by a finite rapidity gap ∆η, as schematically illustrated in Fig. 1.One can then measure correlators of the form ⟪∆q 1 ∆p 1 ⟫, ⟪∆q 1 ∆q 2 ∆p 1 ∆p 2 ⟫, etc., at any order to determine the strength of n = 2, 4, etc, transverse momentum correlations as a function of the width of the rapidity gap.The method of moments, however, enables differential measurements involving multiple n ≥ 4 particles.Let q i , i = 1, . .., represent the p T of particles in bin A, and p i , i = 1, . . .represent the p T of particles in bin B. One then defines eventwise variables T,i (from bin A) and P n = ∑ N i p n T,i (from bin B).Equations (D4 -D10) are then used to obtain ⟪q 1 q 2 ⋯q n p 1 p 2 ⋯p n ⟫ and Eqs.(C1 -C5) or Eq.(C6) are used to obtain the n-order deviates ⟪∆q 1 ∆q 2 ⋯∆q n ∆p 1 ∆p 2 ⋯∆p n ⟫ corresponding to p T correlators involving n particles from bin A and n particles from bin B.
The method described above can readily be adopted also for measurements of charge correlations and multi-particle balance functions.In this case, q i and p i represent the charge of particles in bins A and B. Applications of Eqs.(D4 -D10) then yield generic charge correlators ⟪q 1 q 2 ⋯q n p 1 p 2 ⋯p n ⟫ and Eqs.(C1 -C5) or Eq.(C6) can be used to compute deviates.
Two additional use cases based on two variables q i and p i (with the corresponding event-wise variables Q n and P n ) are worth mentioning.One involves the study of two distinct physics observables (e.g., charge, p T , rapidity, etc) in a single rapidity bin whereas the other involves the measurement of a specific particle observable, e.g., the p T , for two types of particle species.In the first case, the variable q i and p i represent the two observables of interest whereas in the second they tag the species of interest.The determination of correlators between two observables or two types of particle species then proceeds in the manner already described.First get ⟪q 1 q 2 ⋯q n p 1 p 2 ⋯p n ⟫ with Eqs.(D4 -D10) and next used Eqs.(C1 -C5) to compute the deviates of interest.
3. Three or more variables (q, p, r, ...) The examples discussed in the previous paragraph are readily extended towards the computation of factorial cumulants or correlation functions involving three or more kinematic bins and particle types.Of particular interest is the determination of multiple particle balance functions.Although it may not be practical to conduct analyses involving explicit computation of more than three or four kinematic bins or species, it remains possible to consider balance functions involving large number of particles towards the study of long range multi-particle correlations constrained by charge conservation (or other quantum number conservation laws).
As mentioned in Sec.II, we consider here the study of four-particle balance functions using two kinematic bins A and B separated by a finite rapidity gap, as was illustrated in panels (a) and (b) of Fig. 1.The bins A and B could be azimuthally symmetric (i.e., with full azimuth coverage 0 ≤ φ < 2π), or feature partial coverage to suppress contributions from back-to-back jets, as shown in panels (c) and (d) of the same figure.Panel (a) illustrates a measurement involving two positively charged particles in A and two negatively particles in B. A measurement of B 2+2− 4 shall then be sensitive to the strength (or probability) of processes featuring four correlated particles with two +ve and two −ve particles separated by a finite rapidity gap.By contrast, the analysis illustrated in panel (b) would focus on correlated quartets featuring two nearby pairs of +ve and −ve particles.These could be produced by string-like fragmentation processes yielding four or more correlated particles, but it could also result from string fragmentation producing two neutral objects, each decaying into pairs of +ve and −ve particles.

IV. MULTI-PARTICLE BALANCE FUNCTIONS
Yet another application of integral and differential correlation functions of the form Eqs. (2)(3)(4)(5), involves the study of (net)-charge (and other quantum numbers) fluctuations.We first show how Eq. ( 2) can be used to measure net-charge fluctuations and how it connects to measurements of balance functions [11,18,19].We also remind the reader how second moments of the charge are connected to the ν dyn observable [25] and differential charge balance functions [18,19].This then provides a convenient mechanism to introduce higher order balance functions.
In order to express net-charge fluctuations, we re-write Eq. ( 2) by replacing the transverse momentum p T by the charge q i of particles.
where deviates are defined as ∆q i ≡ q i − ⟪q⟫.The variables q i are considered implicit functions of the momentum of the particles and thus cannot be factorized out of the integral.To build this point, let us consider the expression of the average charge in the acceptance Ω of interest where ρ1 (⃗ p) represents a charge density, i.e., not the number density ρ 1 (⃗ p 1 ).Equation ( 50) may be computed based on number densities if ρ1 (⃗ p) is replaced (symbolically) by q α ρ α 1 (⃗ p) where q α is the charge of the particle species of interest, α.A similar development can be done for strangeness or baryon quantum numbers by considering strangeness and baryon densities instead of charge densities.Here for the sake of simplicity, let us restrict the discussion to three types of charged particles: positively charged, neutral, and negatively charged hadrons.The single particle density is then ρ 1 (⃗ p) = ρ (+) Substituting this expression in Eq. ( 50) and inserting the values q i = +1, 0, −1 for each of the three types, one gets where, in the second line, we applied Eq. ( 10).In the following, we assume neutral particles are not measured and consider the total multiplicity N as the sum of the multiplicities of positively and negatively charged particles, i.e.N = N (+) + N (−) .
The calculation of ⟪∆q 1 ∆q 2 ⟫, which corresponds to the covariance of the charges of two measured particles, proceeds in the same way.First expand ∆q 1 ∆q 2 and compute its ensemble average which may be expressed according to where in the second line, we used the expression of the second cumulant C 2 given by Eq. ( 7).Expanding qρ 1 as in Eq. ( 51) and q 1 q 2 ρ 2 according to the integration yields At LHC, A-A collisions produces approximately equal multiplicities (and densities) of positively and negatively charged particle: ⟨N + ⟩ ≈ ⟨N − ⟩.The above expression for ⟪∆q 1 ∆q 2 ⟫ can thus be written within which one recognizes the expression of ν +− dyn [25] ν where in the second line we used the definition of factorial cumulants F α 1 and F α1α2 2 and in the third line, we used the normalized cumulant ratios R α1α2 2 , defined in Eq. ( 1), with α 1 = + and α 2 = −.One then obtains the useful result A similar development can be carried out with a differential version of ⟪∆q 1 ∆q 2 ⟫ and one finds where ) are bound unified balance functions [18] defined according to The functions are constructed in such a way that their respective integral each converge to unity in the full acceptance limit [18] 1 .The fluctuations ⟪∆q 1 ∆q 2 ⟫(⃗ p 1 , ⃗ p 2 ) thus have an upper bound ⟨N ⟩/⟨N (N − 1)⟩.And since ⟨N (N − 1)⟩ → ⟨N ⟩ 2 in the large N (and Poisson) limit, one concludes that ⟪∆q 1 ∆q 2 ⟫(⃗ p 1 , ⃗ p 2 ) should scale in inverse proportion of the system size and the produced particle multiplicity.This expectation is verified from a number of recent measurements of charge fluctuations [35][36][37].
It is natural to seek to extend Eq. ( 60) to higher moments by considering expressions of the form ⟪∆q 1 ∆q 2 ⋯ ∆q n ⟫.We begin, in this section, with a discussion of four-particle balance functions based on an expansion of ⟪∆q 1 ∆q 2 ∆q 3 ∆q 4 ⟫.Extensions to higher orders n = 6, . . ., 10 are presented in Appendix E.
We observe the above expression nearly matches the fourth cumulant expansion, Eq. ( A33), but misses a term of the form 3⟪∆q 1 ∆q 2 ⟫⟪∆q 3 ∆q 4 ⟫.This additional term corresponds to the square of two-particle contributions and needs to be subtracted to eliminate such trivial contributions to the four particle correlator.Subtracting 3⟪∆q 1 ∆q 2 ⟫⟪∆q 3 ∆q 4 ⟫, the integral of C 4 , denoted I +− 4 , may then be written according to We thus proceed to use the integral I +− 4 and the differential cumulant C 4 (⃗ p 1 , . . ., ⃗ p 4 ) to introduce four-particle balance functions.To this end, we write the charged particle 4-tuplet decomposition of C 4 as follows where we assumed the densities are symmetric under permutations of indices, e.g., ρ , to compute each of the coefficients.The integral, Eq. ( 64), becomes where in the last line we used the definition, Eq. (A33), of the 4-particle factorial cumulant.By analogy to Eq. ( 60), one can then introduce 4-particle differential "balance functions" according to where We use the notations B +− 4 and B −+ 4 , to denote n particle balance functions involving, the balancing of negatively and positively charged particles by positively and negatively charged particles, respectively.Inclusion of the ratio of factorial coefficients (4!/2!) insures the integral of B ±∓ 4 converge to unity in the limit of full acceptance.Similar coefficients are introduced, in Appendix E to achieve proper normalization of the higher order balance functions.To indeed verify the integrals of B +− 4 and B −+ 4 integrate to unity over the full acceptance, let P (n) represent the probability of a process involving the production of n pairs of positively and negative charged particles.Mixed factorial moments are thus trivially given by Assuming these are in fact fourth order (or higher) correlations, the integral of B +− 4 (⃗ p 1 , . . ., ⃗ p 4 ) over the full momentum volume thus yields which is indeed equal to unity.By symmetry, the integral of B −+ 4 also converges to unity in a full acceptance measurement.We thus have two equivalent four particle balance functions to study charge balanced four particle correlations.
As for two-particle balance functions, the above expressions can also be generalized to cross species balance function and these shall feature simple sum rules similar to those satisfied by B 2 functions [18].
Given they are based on 4-particle cumulants, the functions ) shall suppress, by construction, contributions from 4-tuplets of uncorrelated particles.Only 4-tuplets featuring genuinely correlated particles would contribute to the strength of the correlators.Tuplets involving only two or three correlated particles would have vanishing contributions.As such the 4-cumulant components of B 4 should indeed suppress contributions from resonance decays resulting into two or three correlated particles.Contributions from hadron resonance decays would then be limited to four prong decays.Considering these typically have small probabilities, one would expect the magnitude of B 4 is then determined by other processes such as "string fragmentation", jet fragmentation, and other multi-particle production processes.Contributions from jet fragmentation processes can be singled out by using a cone-shaped acceptance with, e.g., φ ≈ 1 and y ≈ 1. Conversely, jets can be suppressed by using at least two relatively narrow kinematic bins separated by a sizable rapidity gap.Removal of two-and three-prong decays, as well as contributions from jets, enables a direct study of the underlying processes, such as string fragmentation, that lead to particle production over extended ranges of rapidity.Long range correlations have already been observed in the context of anisotropic flow measurements in pp, pA, and AA collisions.These measurements show that the long correlations are largely dominated by the geometry and fluctuations of the geometry of collisions.They however say little about the underlying nature of the correlations or the correlation length.Measurements of multi-particle charge (or baryon) balance functions would change the focus from the transverse geometry to the longitudinal structure of these correlations and might then shed light on the nature and origin of these correlations.
By construction, n-cumulants are non vanishing only if particle correlations of n-particle are present in the system considered.Consequently, should observations yield vanishing 4-cumulants, it would imply that correlation between balancing charges are only limited to second order contributions.If the 4-cumulants are non-vanishing, they would indicate more intricate production mechanisms.Either way, measurements of B 4 would provide new and valuable information on the structure of particle production dynamics.
Higher order balance functions, B ±∓ n , with 6 ≤ n ≤ 10 can be constructed in a similar way as B ±∓ 4 and are listed in Appendix E. As for B ±∓ 4 , higher order balance functions would suppress contributions from lower order correlations.A systematic study of B ±∓ n for n = 4, 6, 8, etc, would then provide sensitivity to increasingly more complicated production processes featuring a growing number of correlated particles.Such measurements should then provide additional and powerful constrains on multi-particle production models.

V. RELATION TO NET-CHARGE CUMULANTS
Cumulants of the net charge Q (as well as net baryon B and net strangeness S) of the particles measured in a specific acceptance Ω nominally provide a probe of the susceptibilities of the matter formed in nucleusnucleus collisions [38][39][40][41][42]. Several measurements of lower order cumulants (as well as mixed cumulants) have been reported in the recent literature.Ratios of lower order cumulants have been studied, in particular, in the context of the beam energy scans recently performed at the Relativistic Heavy Ion Colliders (RHIC) to identify signatures of a critical point of nuclear matter [43,44].Given the potential significance of such critical point, considerable theoretical and experimental efforts have been deployed to obtain relations between the properties of nuclear matter, net-quantum number cumulants, as well as robust techniques to measure these observables [41].In this context, note that it was recently shown that a simple relation exists between the second net charge cumulant, κ Q 2 and net-charge balance functions B [45] (also see discussions in [46]).This relation is of particular interest because it expresses the magnitude of the non trivial part (non-Poissonian) of κ Q 2 in terms of an integral of the charge balance function B across a specific experimental acceptance (i.e., a specific kinematic range).Given this integral converges to unity, by construction, in the full acceptance limit, it implies the magnitude of κ Q 2 is determined by the shape and width of the balance function relative to the width of the acceptance.This is critical for the beam energy scan because, although the acceptance can be kept fixed, the shape of the B is known to evolve with produced species, system size, nucleus-nucleus collision centrality, and beam energy [1,37,[47][48][49][50][51].As such, the magnitude of κ Q 2 thus constitutes a poorly defined reference in the search of a critical point of nuclear matter.That said, it is also of interest to consider how higher cumulants might be impacted by charge conservation, the size of the experimental acceptance of a measurement, and the dynamics of collisions.
We saw, in sec.IV, that balance functions naturally arise in the calculation of moments ⟪∆q 1 ∆q 2 ⟫ and ⟪∆q 1 ∆q 2 ∆q 3 ∆q 4 ⟫ and yield expressions proportional to factorial cumulants of the particle multiplicities.We thus seek to determine relations between net-charge cumulants, net-charge factorial cumulants, and factorial cumulants of multiplicities of positively and negatively charge particles.Details of the derivations are presented in Appendix A. In this section, we summarize results of interest which indicate that highest order contributions of net-charge cumulants are identical to integrals of multi-particle balance functions of same order.
It is well known that moments, cumulants, factorial moments, and factorial cumulants are readily computed based on their respective generating functions which, herewith, we denote by where the sub-indices m, c, f , and F indicate generating functions of moments, cumulants, factorial moments, and factorial cumulants, respectively.As discussed in Appendix A 3, moments m Q n are obtained by computing n-th derivatives of G m w.r.t.θ Q , evaluated at θ Q = 0, while cumulants κ Q n are obtained from n-th order derivatives of G c w.r.t.θ Q also evaluated at θ Q = 0. Given G c = ln G m , it is then straightforward to compute κ Q n in terms of moments m Q n ′ , with n ′ ≤ n (See Eqs.(A41 -A43).Similarly, factorial and factorial cumulants can be also obtained by n-th order derivatives of G f and G F w.r.t. to s Q .It is however more useful to express the generating functions in terms of multiplicities of positively and negatively charged particles N + and N − and their associated dummy variables θ + and θ − for moments and cumulants calculations and dummy variables s + and s − for factorial moments and factorial cumulant calculations.It is then possible to obtain relations between net-charge cumulants and (mixed) cumulants of N + and N − as well as with factorial cumulants of these multiplicities.As shown in detail in Appendix A 3, one finds even order n-cumulants are given by and so on for higher orders.Intermediate terms of order 1 < n ′ < n were omitted for the sake of clarity.
Comparing the above expressions, as well as Eqs.(A61, A62), with Eqs.(E3 -E7), we observe that the cumulants κ Q n feature a dependence on the mixed factorial moments that exactly matches the expression of the balance functions of order n ≥ 4. Indeed, as for κ Q 2 , we find that the non-trivial component of higher order κ Q n are exactly proportional to integrals of balance functions B +− n , B −+ n introduced in sec.IV.Given these balance functions are governed by sum rules, i.e., their full acceptance integrals are entirely determined by charge conservation.We conclude that as for second order cumulants κ Q 2 , the non-trivial components of higher cumulants, κ Q 2 , n ≥ 4, are determined by integral of functions whose full acceptance limit is solely driven by charge conservation.As for basic balance functions, Eqs.(61, 62), one expects that these higher order balance function integrals feature a strong dependence on the rapidity and transverse momentum coverage of the measurements, as well as the details of the particle production processes at play in the collisions being studied [41].Additionally, as for basic balance functions, it stands to reason that these higher balance function might feature some dependence on collision centrality and beam energy.Such dependencies might thus be better probed with differential balance functions.This suggests that rather than measuring cumulants κ Q n , which only feature information on the integrals of balance functions, it would be better advised to measure differential balance functions.Techniques to compute multi-particle balance functions without the drawbacks of multiple nested loops on particles of an event were discussed in sec.III C whereas kinematic configurations of measurements of potential interest were presented in sec.III D.

VI. SUMMARY
We first advocated, in sec.II, for measurements of integral and differential of multi-particle correlation functions as tools to extract characteristics of heavy ion collisions and the matter they produce heretofore somewhat neglected and susceptible of enhancing the understanding of the physics of these complex systems.We next explicitly presented detailed formula of such multiparticle correlations as well as techniques to compute them in finite time (i.e., single loop on all particles of interest) based on the methods of moments.This set the stage for the development of what we called multi-particle balance functions.These higher order balance functions were introduced based on expectation values of the form ⟪∆q 1 ∆q 2 ⋯∆q n ⟫ but are best computed in terms of combinations of nth order cumulants (or integral factorial cumulants).Much like the original balance functions B 2 introduced by Pratt et al. [11], these new balance functions are defined in such a way that they integrate to unity in full phase space (i.e., all rapidities and p T ≥ 0).As such they too provide a measure of the fraction of charge (or other quantum number) balanced when measured in a finite acceptance.This fraction is expected to be rather sensitive to the details of the (charge conserving) particle production and transport.Indeed, given they are constructed based on n-particle cumulants, they should probe the particle production rapidity and momentum correlation length scales and the details of the particle production mechanisms.
We additionally showed these higher order balance functions have integrals, formulated in terms of factorial cumulants, that are equal to the higher order contributions of net charge cumulants κ Q n .This is an important result that pertains to measurements of net charge (as well as net strangeness and net baryon number) fluctuations based on cumulants κ Q n and their evolution with beam energy in the context of the beam energy scan (BES) at the Relativistic Heavy Ion Collider.The magnitude of these cumulants cannot be corrected for charge conservation given the balance functions integrate to unity in full acceptance.Indeed the magnitude of the integral of the balance functions B ±∓ n measured in a specific acceptance (in rapidity and transverse momentum) is determined by the details of the particle production and transport (e..g., presence of radial flow) and has thus relatively little to do with the intrinsic properties of the matter they originate from (i.e., the susceptibilities of the QGP).
The formalism developed in this work for the deployment of multiple particle correlations is in many ways similar to the techniques used in the context of measurements of anisotropic flow.It is thus likely that the correlation functions discussed in this work might be calculable, with minor or no adaptation, to existing generic frameworks of anisotropic flow measurements.Of particular interest, however, are practical implementations of ⃗ p dependent acceptance and efficiency corrections at the single particle level.Also of interest are efficiency losses related to correlated detector effects that likely manifest themselves differently in the context of the correlation functions discussed in this work.The authors thus plan to follow up this work with additional studies of these practical effects.cumulants in terms of log of their respective generating functions.
It is straightforward (and convenient) to express cumulants as combinations of factorial cumulants if one notices that G c (θ) = G F (s) given s = e θ .Taking n-order derivatives ∂ n θ of the l.h.s.yields cumulants κ n , while derivatives on the r.h.s. are computed based on ∂ θ = (∂s/∂θ)∂ s = s∂ s and yield expressions in terms of F n ′ with n ′ ≤ n.The ten lowest orders are First note that cumulants of a given order k feature a linear combination of all factorial cumulants of lower order k ′ ≤ k.Second, remember that in the context of particle correlation measurements, one can conclude there are correlations of k or more particles only when a factorial cumulant F k is non vanishing.Consequently, if a factorial F k is consistent with zero, within statistical uncertainties, there is no point in measuring κ k or higher order cumulants κ k ′ , with k ′ > k since these do not carry additional experimental information about the system under study.Indeed, in such cases, the magnitude of κ k is primarily determined by factorial cumulants of the lowest orders involving few or, possibly, no particle correlations.

Multi-Variable Systems
Given a function P ( ⃗ N ) stipulating the probability of jointly observing m variables ⃗ N ≡ (N 1 , N 2 , . . ., N m ) corresponding to categories ⃗ α = (α 1 , α 2 , . . ., α m ), which in the context of this work corresponds to kinematic bins or particle species or both, mixed algebraic moments of ⃗ N , denoted m ⃗ α ⃗ n , are defined as and calculable based on a mixed moment generating functions G m ( ⃗ θ) ≡ ⟨e ∑ m i=1 θiNi ⟩ according to where ⃗ α represents all the categories for which moments are evaluated.For instance, a double derivative ∂ θ1 ∂ θ1 would yield a second moment of N 1 , whereas ∂ θ1 ∂ θ2 would yield a mix moment of N 1 and N 2 .Proceeding as for a single variable, one defines mixed cumulants according to where ⃗ θ = 0 specifies derivatives are evaluated with θ i = 0, for i = 1, . . ., m.Similarly, mixed factorial moments, f ⃗ α n , are defined based on mixed moment generating functions G f (⃗ s) ≡ ⟨∏ m i=1 s Ni ⟩ according to where ⃗ s = 1 specifies derivatives are evaluated with s i = 1, for i = 1, . . ., m. Factorial cumulants are defined as derivatives of the factorial cumulant generating functions G F (⃗ s) ≡ ln G f (⃗ s) according to Mathematica [52] enables a speedy and reliable computation of F ⃗ α n .The lowest orders are found to be where the notation ∑ (k) stands for a sum over k (ordered) permutations of the labels α 1 , α 2 , and α 3 , . ... As in the case of a single variable, one can express mixed cumulants κ ⃗ α n in terms of factorial cumulants F ⃗ α n based on the equality G c ( ⃗ θ) = G F (⃗ s) by taking derivatives on l.h.s.relative to θ i whereas derivative are taken relative to ∂ θi = (∂s i /∂θ i )∂ si = s i ∂ si on the r.h.s..

Net Charge Q
Let Q = N + −N − and N = N + +N − define the net-charge and total charged particle multiplicity, respectively, detected in a given event, with N + and N − respectively representing the number of positively and negatively charged particles in that event.The number of positively and negatively charged particles are expected to fluctuate on an event-by-event basis both owing the stochastic nature of the particle production and variations in the processes yielding particles.Moments and cumulants of Q are of interest because they nominally relate to charge susceptibility of the matter formed in heavy-ion collisions [38][39][40]42].Moments m Q n of the net-charge are defined, as in Eq. (A1), according to where P (Q, N ) represents the probability of observing a net-charge Q and total multiplicity N in a particular event.Moments m Q n can evidently be computed based on a generating function of the form G m (N, Q) = ⟨e N θ N +Qθ Q ⟩ but it is of greater interest to obtain the moments, the cumulants, and so on, in terms of moments of the multiplicities N + and N − .Clearly, a simple change of variable enables the definition of P (N + , N − ) as the joint probability of observing events with N + and N − positively and negatively charged particles.The moment generating functions of (mixed) moments of the multiplicities can then be written and successive derivatives of G m w.r.t.θ + and θ − , evaluated at θ + = θ − = 0 yield moments and mixed moments of N + and N − .Introducing the notations ∂ + = ∂/∂θ + and ∂ − = ∂/∂θ − , one computes lowest order associated to charge conservation and the widths of the measurement acceptance.A comprehensive study of the cumulants κ Q with beam energy and system size thus requires a detailed understanding of the evolution of the factorial cumulants F n with beam energy and system size.Given it is likely that multi-particle balance functions of produced particles have intricate dependencies on beam energy, and in particular the growing impact of nuclear stopping at lower energy, we advocate that differential measurements of balance functions provide better insight in the impact of effects associated to the collision dynamics that may otherwise impede the studies of the properties being sought for.There is indeed a one-to-one relation between densities ρ k (⃗ p 1 , . . ., ⃗ p k ) and factorial moments f k as well as between functional cumulants (or cumulant densities) C k (⃗ p 1 , . . ., ⃗ p k ) and factorial cumulants F k .
Once again, close inspection of the above expressions reveal straightforward patterns and one obtains a generic formula in terms of two binomial coefficients as follows q, p, and r, respectively.Proceeding as in sec.III C, the lowest order moments are found to be ⟪Q n P m ⟫ =⟨N q N p ⟩⟪q n i p m j ⟫, (D17) ⟪Q n Q m P o⟫ =⟨N q N p ⟩⟪q n+m i p o j ⟫ + ⟨N q (N q − 1)N p ⟩⟪q n i q m j p o k ⟫, (D18) ⟪Q n Q m Q o P p ⟫ =⟨N q N p ⟩⟪q n+m+o i p p j ⟫ + ⟨N q (N q − 1)N p ⟩⟪q n+m i q o j p p k ⟫ + ⟨N q (N q − 1)N p ⟩⟪q n+o i q m j p p k ⟫ (D19) ⟨N q (N q − 1)N p ⟩⟪q n i q m+o j p p k ⟫ + ⟨N q (N q − 1)(N q − 2)N p ⟩⟪q n i q m j q o k p p l ⟫, ⟪Q n Q m P o P p ⟫ =⟨N q N p ⟩⟪q n+m i p o+p j ⟫ + ⟨N q N p (N p − 1)⟩⟪q n+m i p o j p p k ⟫ + ⟨N q (N q − 1)N p ⟪q n i q m j p o+p k ⟫ (D20) + ⟨N q (N q − 1)N p (N p − 1)⟪⟪q n i q m j p o k p p l ⟫, ⟪Q n P m R o ⟫ =⟨N q N p N r ⟩⟪q n i p m j r o k ⟫, (D21) ⟪Q n Q m P o R p ⟫ =⟨N q N p N r ⟩⟪q n+m i p o j r p l ⟫ + ⟨N q (N q − 1)N p N r ⟩⟪q n i q m j p o k r p l ⟫. (D22) Based on the above expressions, one can iteratively obtain expressions for the moments ⟪q n i q m j p o k r p l ⟫ in terms of moments of products of Qs, P s, and Rs.Given the sum over q and p factorize, the moments become simple combinations of Qs and P s, and one gets at lowest orders

Figure 1 .
Figure 1.Examples of acceptance configurations and particle selections are discussed in the text.

Figure 2 .
Figure 2. Examples of multi-particles correlation for charge, baryon and strange balance studied using a rapidity gap.