Cumulants of multiple conserved charges and global conservation laws

We analyze the behavior of cumulants of conserved charges in a subvolume of a thermal system with exact global conservation laws by extending a recently developed subensemble acceptance method (SAM) [V. Vovchenko et al., arXiv:2003.13905] to multiple conserved charges. Explicit expressions for all diagonal and off-diagonal cumulants up to sixth order that relate them to the grand canonical susceptibilities are obtained. The derivation is presented for an arbitrary equation of state with an arbitrary number of different conserved charges. The global conservation effects cancel out in any ratio of two second order cumulants, in any ratio of two third order cumulants, as well as in a ratio of strongly intensive measures $\Sigma$ and $\Delta$ involving any two conserved charges, making all these quantities particularly suitable for theory-to-experiment comparisons in heavy-ion collisions. We also show that the same cancellation occurs in correlators of a conserved charge, like the electric charge, with any non-conserved quantity such as net proton or net kaon number. The main results of the SAM are illustrated in the framework of the hadron resonance gas model. We also elucidate how net-proton and net-$\Lambda$ fluctuations are affected by conservation of electric charge and strangeness in addition to baryon number.


I. INTRODUCTION
Fluctuations and correlations of conserved charges in statistical systems carry rich information on intrinsic properties of matter. These quantities play a central role in studies of the QCD phase diagram, both in first-principle lattice QCD simulations [2,3] and in heavy-ion collision experiments [4]. Event-by-event fluctuations of different quantities are used in the search of the QCD critical point [5][6][7]. Various correlators of conserved charges, on the other hand, carry information on the relevant QCD degrees of freedom, such as the baryon-strangeness correlator [8].
Fluctuations and correlations of many different quantities, that include both the conserved charges and various hadron number distributions, have been measured in a number of experiments. These include measurements of second order cumulants, both diagonal [9][10][11][12] and off-diagonal [13][14][15], as well as higher-order fluctuation measures [16][17][18][19]. An important question is how to relate the experimental measurements to theoretical predictions. For instance, cumulants of the net-proton number cannot be computed in many of the theories, lattice gauge theory in particular, where only the conserved baryon number is accessible. In such a case one either has to reconstruct net-baryon fluctuations from net-proton measurements [20,21], or directly compare net-proton and net-baryon cumulants, accepting an inevitable systematic error stemming from such an approximation. Another problem is participant (or volume) fluctuations, which is a source of non-dynamical fluctuations affecting comparisons between theory and experiment [22,23].
Perhaps the most important issue is the choice of statistical ensemble. The vast majority of theories operate in the grand canonical ensemble, where the system can freely exchange conserved charges with a reservoir. Direct comparison of grand canonical susceptibilities with heavy-ion data is commonplace in the literature [24][25][26][27][28][29][30][31][32][33]. However, all charges are globally conserved in heavy-ion collisions. This would imply that the canonical ensemble is more appropriate than the grand canonical ensemble. The difference between ensembles does not play a major role if only mean hadron yields are considered in central collisions of heavy ions -due to the thermodynamic equivalence of statistical ensembles for the averages, the difference between hadron abundances evaluated in different statistical ensembles disappears in large systems. However, the thermodynamic equivalence of statistical ensembles does not extend to fluctuations, meaning that values of second and higher order cumulants will depend on the choice of the ensemble, no matter how large the system is.
The experimental measurements typically have a limited momentum acceptance, covering only a fraction of the total momentum space. In Ref. [34] the necessary conditions to emulate the grand canonical ensemble in heavy-ion collisions have been outlined: measurements should be performed in a rapidity acceptance ∆Y acc which is, on one hand, large enough to capture all the relevant physics, ∆Y acc ∆Y cor , where ∆Y cor characterizes the correlation range in rapidity, while on the other hand, it covers only a small fraction of the whole momentum space such that global conservation laws can be neglected, ∆Y acc ∆Y 4π . Furthermore, the measurements should cover the entire transverse momentum range.
Global conservation effects are non-negligible whenever ∆Y acc is comparable to ∆Y 4π . The magnitude of these arXiv:2007.03850v1 [hep-ph] 8 Jul 2020 effects, as well as ways to deal with them, have been studied in the past using a picture of an uncorrelated hadron gas with a single globally conserved charge in a number of papers [35][36][37][38][39][40][41][42]. The analysis in Ref. [37] indicated that the effects of global conservation are sizable already for moderate values of the acceptance fraction α ≡ ∆Y acc /∆Y 4π 0.2, especially for higher-order cumulants. In our recent work [1], we introduced a subensemble acceptance method (SAM) -a procedure to calculate the cumulants in a presence of a single conserved charge for an arbitrary equation of state. In Ref. [43] this formalism was applied to fluctuations in vicinity of a critical point. In the present work, we extend the SAM to equations of state with multiple globally conserved charges, as is appropriate e.g. for QCD with baryon number B, electric charge Q, and strangeness S. In addition to conserved charges, we also explore how cumulants of non-conserved quantities, such as e.g. net-proton number, are affected by multiple global conservation laws. The paper is organized as follows. Sec. II presents the SAM for multiple conserved charges. In Sec. III we illustrate the formalism on an example of a hadron resonance gas model. Discussion and conclusions in Sec. IV close the article.

A. Notation
We shall use a tensor notation throughout this section. Each tensor is denoted by a hat. Where applicable, the number of indices shall determine the tensor rank. We also adopt the Einstein notation, where a repetition of each index implies summation over that index.
Let us have a vectorQ = (Q 1 , . . . , Q N ) of N independent conserved charges in the system. Each conserved charge is associated with a chemical potential. The vector of chemical potentials is denotedμ = (µ 1 , . . . , µ N ). In the grand canonical ensemble, GCE, the relation between cumulantsκ gce and susceptibilitiesχ is straightforward: Here p is the pressure, T the temperature and V the volume of the system. The relation (1) applies for an arbitrary cumulant (susceptibility) of order M . Both the susceptibilitiesχ i1...i M and the cumulantsκ i1...i M are symmetric with respect to any permutation of their indices. The notation (1) for the susceptibilities is different from the one commonly used in the QCD literature [44,45]. There, the susceptibilities read The quantities in Eqs. (1) and (2) are equivalent, i.e.χ i1...i M ≡ χ Q1...Q N l1...l N , when the set of indices (i 1 . . . i M ) contains exactly l 1 elements equal to unity, exactly l 2 elements equal to two, and so on. For QCD with three conserved charges, baryon number, electric charge, and strangeness, one hasQ = (Q 1 , Q 2 , Q 3 ) = (B, Q, S). As an example, we write here a diagonal and an off-diagonal fourth order QCD susceptibilities using the two notations:

B. Subensemble acceptance method
In the present work we generalize the relation (1) to account for the presence of global conservation laws for all of the conserved charges. Let us consider a subvolume V 1 = αV of a uniform thermal system where all the charges are globally conserved. Following our earlier work [1], we assume that the subvolume V 1 as well as the remaining volume V 2 = (1 − α)V are both of a macroscopic size, i.e. they are large compared to correlation length ξ, V 1 ξ 3 and V 2 ξ 3 . Here ξ refers to any correlation length of relevance to the cumulants of multiple conserved charges under consideration. As a consequence, one can neglect all interactions at the surface separating the two subsystems, meaning that the total Hamiltonian can be expressed as a sum the subsystem Hamiltonians, i.e.
In this case the canonical ensemble partition function of the total system with total conserved charge vectorQ reads [46] Z(T, V,Q) = Q1 Z(T, αV,Q 1 ) Z(T, βV,Q −Q 1 ) .
Here β ≡ 1 − α, and the sum goes over all possible values of conserved charges in the first subsystemQ 1 . The probability P (Q 1 ) to simultaneously find all conserved charges in the subsystem with volume V 1 equal toQ 1 is proportional to the product of the canonical partition functions of the two subsystems: In the thermodynamic limit, V → ∞, we have whereρ =Q/V is the vector of densities of all the conserved charges and f (T,ρ) is the free energy density.
To evaluate the cumulantsκ i1...i M [Q 1 ] of the distribution of conserved chargesQ 1 inside the subvolume V 1 we introduce the cumulant generating function GQ 1 (t): The cumulants,κ i1,...,i M [Q 1 ], correspond to the Taylor coefficients of GQ 1 (t): Here we introduced a shorthandκ i1,...,im (t) for generalizedt-dependent cumulants. All second and higher-order cumulants can be obtained by differentiating the first order cumulants,κ i1 (t), with respect to the components oft. The first order cumulants read with the (un-normalized)t-dependentQ 1 probability function, In the thermodynamic limit, the probabilityP (Q 1 ;t) is highly peaked at the mean values ofQ 1 . The vector Q 1 i (t) of the mean values ofQ 1 is, therefore, determined by requiring that all partial derivatives ofP (Q 1 ;t) vanish, Hereμ is the reduced chemical potential as a function of temperature and densities of conserved charges, andρ 1 (t) = Q 1 (t) /(αV ) andρ 2 (t) = Q 2 (t) /(βV ) = (Q − Q 1 (t) )/(βV ). Fort = 0 the solution is: i.e. the charges are distributed uniformly between the subsystems, as should be the case by construction. Given that the first order cumulant of the i-th charge isκ i [Q 1 ] ≡ Q 1 i = αVρ i and the first order susceptibility equalŝ χ i ≡ρ i /T 3 by definition, we getκ C. Second order cumulants Given the generalized first order cumulantκ i , the second order cumulantsκ ij (t) arễ To evaluateκ ij (t) we differentiate Eq. (13) with respect tot j . We obtain where we used the chain rule multiple times. Here one observes that and, by definition, Equation (18), therefore, reduces toδ The derivatives ∂μ 1,2 i /∂ρ 1,2 j1 can be expressed in terms of second order susceptibilities,χ ij . Indeed, since the second order susceptibilities are defined by the derivatives of the conserved number densities, the inverse derivatives, ∂μ i /∂ρ j , are defined by the inverse matrix of second order susceptibilities, The matrix equation (21) for the second order cumulants, therefore, readŝ Hereχ ij1 andχ ij1 correspond to the matrix of second order conserved charge susceptibilities evaluated at an arbitrary finitet in the first and second subsystems, respectively.
The solution to Eq. (24) isκ We note, that the correction factors due to global charge conservation are the same for all second order susceptibilities. An important consequence of this result is that a ratio of any two second order cumulants of conserved charges coincides with the corresponding ratio of grand canonical susceptibilities, i.e. effects of global conservation are canceled out.

D. Third order cumulants
The third order cumulants are calculated by differentiating Eq. (25) with respect tot k : To evaluate this derivative we will make use of the following identity from the matrix calculus: Now, in Eq. (27), we expressκ ij asκ where, following Eq. (25),κ The third order cumulant defined by Eq. (27), therefore, readŝ Let us now evaluate ∂χ −1 j1j2 /∂t k . Using the identity (28) gives which, after applying the chain rule to the middle term, turns into Using all the identities we derived above [Eqs. (19), (20), and (23)] for each of the new terms, as well as the definitioñ χ ijk = ∂χ ij /∂μ k of the third order GCE susceptibilities, we get The derivation for ∂χ −1 j1j2 /∂t k that appears in the second term of the expression forκ ijk [Eq. (31)] is analogous. The only difference is the prefactor: The final result forκ ijk (t) is obtained by substituting Eqs. (34) and (35) into (31): Fort = 0 we haveχ =χ =χ for all ranks of the susceptibility tensorχ. Alsoκ ij = αV T 3 βχ ij [Eq. (26)]. This implies that every convolution ofκ ij1 withχ −1 j1j provides a Kronecker symbol times a factor, namelyκ ij1χ This simplifies the evaluation ofκ ijk considerably. In the end we obtain the following: Similar to the second order cumulants [Eq. (26)], the global conservation corrections are identical for all third order cumulants. As a consequence, effects of global conservation cancel in any ratio of any two third order cumulants of conserved charges.

E. Results up to sixth order
The fourth order and higher order cumulants are calculated by iteratively differentiating Eq. (36). For example, for the fourth order cumulants one has:κ The details of the iterative procedure for calculating higher-order cumulants are described in Appendix A. Here we present the resulting SAM expressions up to the sixth order, evaluated att = 0. For completeness we list here all cumulants starting from the first order: Here i 1 , . . . , i 6 = 1, . . . , N in all of the above equations. The notation σ∈S M corresponds to a sum over all M! permutations of a set (1, . . . , M ). σ i is the ith element of the permutation σ. In Appendix B we provide an example of a detailed calculation of a fourth order cumulant using Eq. (42). This example is useful for understanding the notation entering Eqs. (39)- (44). A Mathematica notebook to express cumulants in Eqs. (39)- (44) in terms of the susceptibilities using the QCD notation (2) is available via [47]. It follows from Eq. (42) that global conservation effects generally affect different fourth order cumulants in a different way. In contrast to the cumulants of lower order, the global conservation effects do not cancel in ratios of fourth order cumulants. An exception to this are vanishing chemical potentials,μ = 0, whereχ ijk = 0 for all i, j, k and where only the first term in the r.h.s of (42) is non-zero. The fifth order cumulants (43) have a structure similar to the fourth order cumulants. Thus, atμ = 0 the effects of global conservation cancel in any ratio of any two four order and any two five order cumulants. The sixth order cumulants (44) have a considerably more involved structure.
These expressions reproduce the results of Ref. [1], where the SAM was originally formulated for the case of a single conserved charge.

Two conserved charges B and Q
In a case of two conserved charges, say baryon number B and electric charge Q, we haveQ = (B, Q) and N = 2. Here we would like to illustrate how the cumulants of baryon number are affected by the presence of exact conservation of other conserved charges. First, we note that, following Eqs.
It is evident from Eq. (51) that the presence of a conserved electric charge influences the fourth order baryon number cumulant if there are baryon-electric charge correlations in the underlying equation of state. This implies corrections to the relation (48) derived in Ref. [1] for a system with a single conserved charge. To elucidate these corrections we rewrite Eq. (51) in the following form: Effects of BQ-correlations can be sizable when the grand canonical off-diagonal susceptibilities are comparable to the diagonal ones. On the other hand, corrections to Eq. (48) are small when χ BQ 11 , χ BQ 21 χ Q 2 or when the matterantimatter symmetry is small, i.e. for zero (small)μ. In the latter case, realized in heavy-ion collisions at LHC energies, all odd-order susceptibilities are (close to) zero, thus the second term in the r.h.s. of Eq. (52) (nearly) vanishes.
The fifth order cumulant, κ 5 [B 1 ], reads The discussion above on the behavior of κ 4 [B 1 ] largely applies to κ 5 [B 1 ] as well.
Explicit expression for κ 6 [B 1 ] can be obtained by expanding Eq. (44). The expression is very lengthy, containing many different contributions from various BQ correlators. Here we only present the explicit expression for κ 6 [B 1 ] at µ = 0, where all odd-order susceptibilities vanish and thus simplify the formula considerably: 3. Three conserved charges B, Q, and S Let us consider now three conserved charges: baryon number B, electric charge Q, and strangeness S. In this casê Q = (B, Q, S) and N = 3. The first three cumulants of baryon number B 1 are the same as in the case of a single baryon charge. The fourth cumulant reads Here D[χ 2 ] is the determinant of the matrix of second order susceptibilities: The fourth-order cumulant κ 4 [B 1 ] is affected by both the baryon-electric charge and baryon-strangeness correlations. Even the correlation between electric charge and strangeness does contribute, through a correlator χ QS 11 . It is notable that the entire second term in the r.h.s. of Eq. (55) vanishes at LHC energies (μ = 0), i.e.
We do not write here the lengthy expressions for κ 5 [B 1 ] and κ 6 [B 1 ]. These can be worked out from Eqs. (43) and (44), if desired. We will only write, for completeness, the expression for κ 6 [B 1 ] forμ = 0 (LHC energies), where it is considerably simplified:

G. Strongly intensive quantities
Our considerations in the present paper are focused on effects of global conservation of multiple conserved charges, and how they distort the grand canonical baseline in the measured cumulants. Another non-dynamical source that affects the measurements are fluctuations of the system volume, that cannot be completely avoided in heavy-ion collisions. Different methods exist to address volume fluctuations [22,23,38]. Here we explore briefly the possibility to construct quantities that are insensitive to both the global charge conservation and the volume fluctuations.
More specifically, we consider strongly intensive quantities -fluctuation measures that were developed in Ref. [22] and designed to be insensitive to volume fluctuations. These comprise of two combinations of first and second moments of two extensive quantities. Here we take two conserved charges, say Q a and Q b , both measured in a subvolume V 1 : The normalization factors C ∆ and C Σ correspond to an arbitrary extensive measure not sensitive to volume fluctuations. As an example, one can take any linear combination of mean values of Q a and Q b . Possible specific choices of C ∆ and C Σ have been discussed in Ref. [48]. The strongly intensive quantities can be written in terms of cumulants: The quantities ∆[Q a , Q b ] and Σ[Q a , Q b ] (and any combination of the two) are insensitive to fluctuations of the total system volume V . To show their sensitivity to global charge conservation we use Eqs. (39) and (40): We note that the above two equations are obtained assuming that the subvolume fraction α is constant, i.e. it is unaffected by volume fluctuations. If α does fluctuate, e.g. due to possible fluctuations in baryon stopping, the formalism will require a generalization to account for that. Given that C ∆ and C Σ are proportional to any extensive (i.e. proportional to the subvolume V 1 ) measure that is not sensitive to volume fluctuations, we can, without the loss of generality, write these factors as C ∆ = αV T 3 χ ∆ 1 and C Σ = αV T 3 χ Σ 1 where χ ∆ 1 and χ Σ 1 do not depend on V and α, e.g. they may be chosen as linear combinations of χ Qa 1 and χ Q b 1 . Then: Both strongly intensive measures ∆[Q a , Q b ] and Σ[Q a , Q b ] are affected by global charge conservation, through a common factor 1 − α. This implies that a ratio of these two quantities is neither affected by volume fluctuations nor the global charge conservation: H. Non-conserved quantities So far we focused on the cumulants of multiple conserved charges. It is, however, notoriously difficult to measure neutral particles event-by-event in heavy-ion experiments, preventing direct measurements of baryon number and strangeness fluctuations. For this reason one usually considers non-conserved quantities such as net-proton or net-kaon number as proxies for net-baryon and net-strangeness fluctuations. The measurements of electric charge fluctuations do not suffer from this problem and can be done directly. For instance, the STAR collaboration has recently reported measurements of net-proton, net-kaon, and net-charge second order cumulants [15].
Here we consider the behavior of a non-conserved quantity, such as net-proton or net-kaon number, in the presence of exact conservation of charges. For clarity, we will refer to this quantity as a net proton number, N p . The considerations below, however, are general and apply to any non-conserved quantity, not only net-proton number.
As the net-proton number is not a conserved quantity, it has a distribution for fixed values of the conserved chargeŝ Q. This implies that the canonical partition function can be written as sum over all possible values of N p Here W (T, V,Q; N p ) counts the number of configurations that yield a particular net-proton number N p in the final state. We note that W (T, V,Q; N p ) is not easily accessible in theoretical calculations, such as lattice QCD, as it requires a careful projection on asymptotic states that count the net number of protons emerging from a thermalized QCD matter created in heavy-ion collision. However, the number of (anti)protons in a given heavy-ion event is a well defined, physical observable accessible to experiment. In the framework of the hadron resonance gas (HRG), which is successfully applied to extract the chemical freeze-out temperature, the number of protons is given by the sum of primordial protons and those arising from all strong and electromagnetic decays of resonances, such as for example ∆(1232). Or in other words, the projection to asymptotic states in the HRG is equivalent to taking all resonance decays into account. For the purposes of this paper is not important how and if W (T, V,Q; N p ) can be calculated in a given theory but rather that it is a well defined quantity, which it is. Since N p is not conserved, the system can freely fluctuate from a configuration with a particular value of N p to a configuration with another value. Thus, N p behaves just as a conserved charge in a grand canonical ensemble, so that W (T, V,Q; N p ) can be regarded as a generalized canonical partition function with fixed values of conserved chargesQ and net-proton number N p . In the thermodynamic limit, V → ∞, the following representation of W (T, V,Q; N p ) holds: Heref (T,ρ, ρ p ) is a generalized free energy density that depends on net-proton density ρ p = N p /V , in addition to the conserved charge densities. In the thermodynamic limit W (T, V,Q; N p ) is highly peaked around N p , the sum (68) is determined by the maximum term at N p = N p , i.e.
where ρ 0 p = N p /V . By considering the canonical partition function Z(T, V,Q) in the form given by Eqs. (68) and (69) we can now introduce a generalized cumulant generating function for the joint distribution ofQ 1 , N 1 p , and N 2 p . Here N 1 p and N 2 p are the net-proton numbers in the first and second subsystems, respectively. The generating function reads , therefore all cumulants that involve the conserved chargesQ, but not N 1,2 p , will coincide with the results obtained in previous section.

Off-diagonal cumulants involving a single conserved charge
Let us introduce a matrix of the grand canonical second order susceptibilities for the joint (Q, N p ) distribution: Theχ ij matrix is (N + 1) × (N + 1) dimensional, and χQ iQj is the N × N matrix of second order susceptibilities of conserved chargesQ, defined in Eq. (22), χQ ip = χ pQi = ∂ρ p /∂μ i is a grand canonical correlator between N p and conserved chargeQ i , and χ pp = ∂ρ p /∂μ p is the grand canonical susceptibility for the net-proton number N p .
Differentiating Eq. (73) with respect tot j , where j ∈ 1 . . . N , we get The tilde inχ means that the susceptibilities are calculated in the first subsystem at arbitrary values oft and t 1 p , whereas the same quantity without a tilde corresponds to the susceptibilities evaluated att = 0 and t 1 p = 0. The same notation applies to the use of tilde in the notation for the cumulants κ of the ( 75) correspond to the elements of the inverseχ matrix [Eq. (74)].κ N +1,j corresponds to a mixed cumulant involving net-proton number N 1 p and conserved chargeQ 1 j , both evaluated in the first subsystem. Equation (75) can be solved forκ N +1,j , yielding: Now we sett = 0 and t 1 p = 0, so thatχ −1 →χ −1 . Furthermore, for κ j1j with j 1 , j = 1 . . . N we can use the result (40) leading to We observe that N +1 j1=1χ −1 N +1,j1χ j1j = 0 as this expression corresponds to an off-diagonal element of the (N + 1) × (N + 1) identity matrixÎ =χ −1χ . This implies N j1=1χ Given the structure ofχ [Eq. (74)] we can identifyχ N +1,j = χ pQj as a grand canonical correlator of net-proton number with a conserved chargeQ j . The corresponding off-diagonal cumulant κ pQj reads The cumulant κ pQj is affected by the global conservation of charges in the same way as all second order cumulants of conserved charges [Eq. (40)]. This means effects of conservation laws cancel out in following ratios: Furthermore, as our derivation has been obtained for an arbitrary non-conserved quantity, we can also consider correlators of the electric charge with different non-conserved quantities, such as e.g. net proton and net kaon numbers. The global conservation factors cancel κ pQj

Variance of a non-conserved quantity
Let us now differentiate Eq. (73) with respect to t 1 p : This can be solved forκ N +1,N +1 . Fort = 0, t 1 p = 0, the solution reads where we used the result (78) for κ j1,N +1 . The sum N +1,j1χ j1,N +1 corresponds to the (N +1), (N +1) element of an identity matrixÎ =χ −1χ , therefore, By definition, in the limit α → 1 κ pp reduces to the variance of net-proton distribution in the canonical ensemble κ ce pp ≡ V T 3 χ ce pp . Therefore, by setting α = 1 in Eq. (85) we obtain the canonical net-proton susceptibility: From an analysis of Eqs. (85) and (86) one observes that the net-proton susceptibility κ pp /(αV T 3 ) represents a linear combination of the grand canonical net-proton susceptibility χ pp and the canonical net-proton susceptibility χ ce pp = detχ det χ at all values of α.

A. HRG model setup
To illustrate the main features of the SAM in the presence of multiple conserved charges we shall consider the behavior of various fluctuation measures in a hadron resonance gas (HRG) model [49]. The HRG model describes the hadronic phase as a multi-component gas of known hadrons and resonances and has broad applications for describing hadrochemistry in relativistic heavy-ion collisions [50][51][52].
In the present work we take the simplest variant of the HRG model where we neglect quantum statistics, finite resonance widths, and excluded volume corrections. Even though the model in this case reduces simply to a multicomponent ideal gas of Maxwell-Boltzmann particles, the different hadron species do carry different values of the three QCD conserved charges, baryon number B, electric charge Q, and strangeness S. This induces non-trivial cross-correlations between conserved charges, making the model suitable for studying the effects of multiple global conservation laws.
For an arbitrary grand canonical HRG model the conserved charge susceptibilities are Here we used the conventional notation [Eq. (2)] for the susceptibilities. The index i sums over all hadron species in the HRG, including both particles and antiparticles. The degeneracy and mass of a hadron specie i are denoted by d i and m i , respectively, while b i , q i , and s i are its baryon number, electric charge, and strangeness. We also introduced a shorthand χ hrg i (88) for the susceptibility of particle number of hadron species i. We take into account all established light-flavor and strange hadrons and resonances from the 2014 edition of Particle Data Tables [53], as incorporated in the default particle list of the open source thermal-statistical package FIST [54].
For our calculation we use the following values of the thermal parameters: T = 160 MeV and µ B = 100 MeV. The temperature value is typical for chemical freeze-out in heavy-ion collisions [52,55], while a non-zero value of µ B allows to study effects that are absent at µ B = 0, such as mixing of cumulants of various order discussed in In addition to analytic calculations of cumulants within the SAM, we also perform Monte Carlo calculations of the same cumulants by sampling the HRG in the canonical ensemble. We employ the canonical ensemble sampler of the FIST package to generate hadron multiplicities in the full acceptance. The canonical ensemble sampler follows the efficient multi-step procedure introduced by Becattini and Ferroni [58], and its implementation in FIST is detailed in Ref. [59]. Given that hadron coordinates are uncorrelated in the HRG model, we then independently apply a Bernoulli trial with probability α to each hadron in each event to establish whether it belongs to a subvolume V 1 = αV . In our sampling procedure we take values B = 20, Q = 8, S = 0 of the globally conserved charges. The value of the total volume of V is calculated using the grand canonical HRG model such that the mean baryon number, electric charge, and strangness are equal to the canonical ensemble values listed above. The resulting value of V = 522.8 fm 3 ensures consistency between the grand canonical and canonical formulations regarding all average quantities. This value of V is, on one hand, large enough such that deviations from the thermodynamic limit are small, and on the other hand, it is small enough to obtain sufficient statistics for an accurate calculation of cumulants up to the fourth order. In total we generate 10 8 events for the present analysis.
The Monte Carlo calculations present a powerful validation of the SAM, as they make no use of the assumptions which go into that formalism.

B. Second order cumulants of conserved charges
The SAM predicts that appropriately scaled second order cumulants κ X 2 /(αβV T 3 ), where X = B, Q, S, are independent of the value of α and coincide with the corresponding grand canonical susceptibilities χ X 2 [see Eq. (40)]. The Monte Carlo calculations do indeed confirm this behavior, as shown in the left panel of Fig. 1. Measurements of the variance of conserved charges can thus be used to extract the grand canonical susceptibility if the values of the acceptance fraction α, the volume V , and the temperature T can be reliably estimated.
In practice, a reliable estimation of α, V , and T is challenging. For this reason it is useful to consider observables where these quantities do not appear. As already discussed, Eq. (40), a ratio of any two second order cumulants of conserved charges is insensitive to α and V , and coincides with the corresponding ratio of grand canonical susceptibilities. The right panel of Fig. 1 depicts the α-dependence of cumulants ratios κ BQ 11 /κ B 2 , κ QS 11 /κ S 2 , and κ BS 11 /κ S 2 , as calculated within the Monte Carlo event generator (symbols) and within the SAM using the grand canonical susceptiblities. The Monte Carlo calculations confirm the expected α-independence of these ratios, their constant values being consistent with the grand canonical baseline.

C. Third order cumulants of conserved charges
Let us consider now third order cumulants. First we analyze the so-called skewness ratio κ X 3 /κ X 2 for X = B, Q, S. The skewness is a non-Gaussian fluctuation measure that characterizes the asymmetry of a distribution around the mean value. The signs of the skewness of QCD conserved charges are thought to be sensitive probes of the QCD phase structure [60]. The SAM predicts that the skewness κ X 3 /κ X 2 scaled by (1 − 2α) is independent of acceptance and coincides with the skewness χ X 3 /χ X 2 evaluated in the grand canonical ensemble. The left panel of Fig. 2 depicts the ratios (κ X 3 /κ X 2 )/(1 − 2α) for baryon number, electric charge, and strangeness evaluated using the Monte Carlo event generator. The Monte Carlo results are consistent with (κ X 3 /κ X 2 )/(1 − 2α) being independent on α and coincide with the grand canonical (χ X 3 /χ X 2 ) susceptibility ratios. In addition, the SAM predicts that any ratio of two third order cumulants of conserved charges is insensitive to global conservation laws. This follows from Eq. (41). The right panel of Fig. 2 depicts the α-dependence of cumulant ratios κ B 3 /κ Q 3 , κ QS 21 /κ S 3 , and κ BS 12 /κ S 3 , as calculated within the Monte Carlo event generator (symbols) and within the SAM using the grand canonical susceptiblities. The Monte Carlo calculations are consistent with the α-independence of all these ratios, and in agreement with the grand canonical baseline, as predicted by the SAM. The statistical errors in the Monte Carlo calculations become large in the vicinity of α = 1/2, as clearly seen in Fig. 2. This is a consequence of the fact that third order cumulants vanish at α = 1/2, as follows from Eq. (41). A ratio of third order cumulants in the vicinity of α = 1/2 corresponds to a ratio of two small numbers, hence the large statistical uncertainties. As a consequence, it would be advisable to perform experimental analysis of third order cumulants in acceptances away from α = 1/2.

D. Fourth order cumulants of conserved charges
We turn now to fourth order cumulants. An interesting new aspect here is that fourth order cumulants are determined not only by the corresponding fourth order GCE susceptibilities but also by second and third order mixed susceptibilities, as seen from Eqs. (42) [see also Eq. (55)]. The HRG model analysis allows to estimate the importance of these mixed susceptibilities with regard to the behavior of fourth order cumulants in a finite acceptance.
We shall analyze here the fourth-to-second order ratios κ X 4 /κ X 2 of diagonal cumulants for X = B, Q, S -the socalled kurtosis of a conserved charge distribution. The kurtosis of the baryon number in a subvolume V 1 = αV is evaluated using Eqs. (55) and (46): The kurtosis of electric charge Q and strangeness S distributions are given by similar expressions, which can be explicitly derived from the general formula (42). In fact, these expressions can also be obtained from Eq. (89) by cyclic permutations in (B, Q, S). Figure 3 depicts the α-dependence of the kurtosis of baryon number (black symbols), electric charge (blue symbols), and net strangeness (red symbols), as obtained from Monte Carlo simulations. The solid lines in Fig. 3 correspond to the SAM predictions [Eq. (89)], these agree with the Monte Carlo results. To estimate the relevance of crosscorrelations between multiple conserved charges on κ X 4 /κ X 2 it is instructive to compare the results to predictions of the SAM for a single conserved charge. The ratio κ X 4 /κ X 2 in this case is obtained by dividing Eq. (48) by Eq. (46), which is depicted in Fig. 3 by dashed lines. Equation (90) deviates from the general result (89) by no more than a few percent, indicating the behavior of kurtosis of a conserved charge is primarily driven by the exact conservation of that charge, whereas the influence of exact conservation of other conserved charges is subleading. We also remind the reader that this influence of other conserved charges vanishes completely at µ B = 0 (LHC energies), where all odd-order susceptibilities is zero. These observations lead us to conclude that the simplified Eq. (90) is rather accurate for practical applications, at least for µ B 100 MeV ( √ s N N 40 GeV).

E. Off-diagonal cumulants involving non-conserved quantities
Next, we switch from cumulants of globally conserved quantities to cumulants involving quantities that are not globally conserved. This does better reflect the current experimental reality. To be more specific, we shall consider second order cumulants involving net-proton, net-kaon, and net-charge numbers. This is in part motivated by recent experimental efforts of the STAR collaboration in measuring these quantities [15]. While the net electric charge is globally conserved, the net proton and net kaon numbers do fluctuate even in the full acceptance. As follows from the results of Sec. II H, a correlation of a non-conserved quantity with a conserved charge is affected by the global conservation laws by the same factor (1 − α) as all second order cumulants of conserved charges. The implication is that measurable ratios such as κ pQ 11 /κ Q 2 and κ kQ 11 /κ Q 2 are expected to be unaffected by the global conservation laws and coincide with the corresponding ratios χ pQ 11 /χ Q 2 and χ kQ 11 /χ Q 2 of the grand canonical susceptibilities.
The grand canonical susceptibilities can be evaluated in the HRG model. Extra care should be taken to account for the large feeddown from decays of resonances to final yields of protons and kaons. The grand canonical fluctuations and correlations of final particle numbers after resonance decays in HRG have been worked out in Refs. [36,61]. The grand canonical susceptibility describing correlations between two final-state hadron species i and j reads Here the index R sums over all resonances. The quantity n i n j R is an average product of the number of hadron species i and j which result from decays of resonance R. This quantity takes into account the multinomial nature of resonance decays. χ hrg i is defined in Eq. (87). The correlator χ kQ 11 of net-kaon number with net-charge is evaluated by incorporating contributions from protons and antiprotons as well as of all charged hadrons in the final state: Here k + and k − corresponds to positively and negatively charged final state kaons, respectively, while the index j runs over all hadron species with charge q j in the final state, including both particles and antiparticles. The expression for the net-proton-net-charge correlator χ pQ 11 is analogous to Eq. (92). Evaluation of the grand canonical susceptibilities in accordance with Eqs. (91) and (92) is readily implemented in the FIST package, as documented in Ref. [54].
The Monte Carlo procedure requires some extensions to incorporate the resonance decays as well. This is done in the following way. First, we mark all the primordial hadrons and resonances as either belonging to subvolume V 1 or subvolume V 2 by doing the Bernoulli trials, as before. Then we let all unstable resonances decay until only the stable decay products are left. All decay products stemming from the same primordial resonance are assigned the same  Figure 4. Dependence of cumulant ratios κ pQ 11 /κ Q 2 and κ kQ 11 /κ Q 2 (left panel ) and κ πp 11 /κ p 2 and κ pk 11 /κ k 2 (right panel ) on the acceptance α, as calculated in the hadron resonance gas model using canonical ensemble Monte Carlo sampler (symbols) and analytically in the framework of the subensemble acceptance method (lines). Here p, k, and π stand, respectively, for netproton, net-kaon and net-pion numbers evaluated after resonance decays. The Monte Carlo sample contains 10 7 events and uses the same thermal parameters as in Fig. 1. subvolume as that primordial resonance. We generate 10 7 events where we perform resonance decays in accordance with the above description. All the cumulants of interest are then computed in the standard way, as a statistical average from final-state hadron distributions.
The left panel of Fig. 4 depicts the Monte Carlo calculation results for cumulant ratios κ pQ 11 /κ Q 2 and κ kQ 11 /κ Q 2 as functions of acceptance parameter α. The ratios do not exhibit any sensitivity to the value of α. As predicted by the SAM, these quantities coincide with the corresponding ratios χ pQ 11 /χ Q 2 and χ kQ 11 /χ Q 2 of the grand canonical susceptibilities, evaluated through Eqs. (91) and (92). Measurements of such quantities can thus directly reflect intrinsic properties of matter that are characterized by the grand canonical susceptibilities.
Second order cumulant ratios involving non-conserved quantities only, on the other hand, do depend on the size of the subvolume. Effects of global conservation laws no longer cancel out in such a case. To illustrate this aspect we show in the right panel of Fig. 4 the α-dependence of cumulant ratios κ πp 11 /κ p 2 and κ pk 11 /κ k 2 involving correlations of net proton number with net pion and net kaon numbers, respectively. These ratios clearly do exhibit α dependence, interpolating between ratios of grand canonical (α → 0) and canonical (α → 1) susceptibilities.
The SAM predicts a linear α-dependence of net-proton and net-kaon cumulants κ p 2 and κ k 2 [Eq. (85)]. Assuming the same linear α-dependence holds also for the correlators κ pπ 11 and κ pk 11 one can expresses the ratios κ πp 11 /κ p 2 and κ pk 11 /κ k 2 as follows: Here the quantities denoted as (. . .) ce are the susceptibilities evaluated in the canonical ensemble. The procedure to evaluate such quantities analytically is detailed in Refs. [36,61], and is readily available in the FIST package. The α-dependence of κ pπ 11 /κ p 2 and κ pk 11 /κ k 2 evaluated according to Eqs. (93) and (94) is shown in the right panel of Fig. 4 by dashed lines. These analytic expectations agree with the Monte Carlo results, suggesting that all second order cumulants of non-conserved quantities can be generally described by a linear function in α which interpolates between the grand canonical (α → 0) and canonical (α → 1) limits, as written in the numerator and denominator of Eqs. (93) and (94).

F. Net-proton and net-Λ fluctuations
We shall conclude our HRG model analysis by exploring fluctuations of net-proton and net-Λ number. These quantities are more accessible to experimental measurements than net-baryon number fluctuations. The grand canonical HRG model provides a natural baseline for the second order cumulants of net-proton (net-Λ) fluctuations: the net-proton (net-Λ) number fluctuations are described by the Skellam distribution, meaning that the second order cumulant simply counts the mean number of protons (Λ) and antiprotons (Λ): This baseline is not modified by resonance decays since no decays are known to generate proton-proton (ΛΛ) or proton-antiproton (ΛΛ) correlations. Naturally, both the net proton and net Λ fluctuations are affected by the exact global conservation of charges. In particular, the effect of baryon number conservation on net-proton fluctuations has extensively been studied in the literature within the HRG model with a canonical treatment of baryon number [37,38]. If the acceptance for particles and antiparticles is uniform, the effect of global baryon number conservation on net-proton fluctuations is the following: Here N acc p(p) is the mean number of (anti)protons within the acceptance where net-proton fluctuations are measured and N 4π B(B) is the mean number of (anti)baryons in the full space. The expression for net-Λ fluctuations is analogous. Note that α p B can also be defined as α p 37,38] if the acceptance factor is the same for protons and antiprotons. This is the case at the LHC due to particle-antiparticle symmetry at µ B = 0, as well as at low collision energies where the production of antibaryons can be neglected. In a general case, however, differences between the two definitions do appear because of non-zero values of µ Q and µ S . We checked that these differences do not exceed 2% at all collision energies in the HRG model, therefore, in practice either of the two definitions may be adopted.
Expression (97) has been used by the ALICE collaboration to interpret their measurements of the behavior of net-proton fluctuations in Pb-Pb collisions at the LHC as primarily driven by an effect of global baryon number conservation [11]. At the same time, fluctuations can also be influenced by exact conservation of electric charge and strangeness, in particular since protons do carry electric charge and Λ's do carry strangeness. This point has been made in a recent STAR paper on net-Λ fluctuations [12], where an ad hoc relation has been used to estimate the simultaneous effect of net baryon and net strangeness conservation. In the STAR paper [12] the baseline κ NBD 2 [Λ −Λ] corresponds to a negative binomial distribution, which is similar although slightly different than the Skellam baseline κ Sk 2 [Λ −Λ] discussed here. Here we shall employ the SAM to rigorously study the effect of multiple conserved charges on net-proton and net-Λ fluctuations. Our starting point is Eq. (85) which describes the second order cumulant of an arbitrary non-conserved quantity within a subvolume for an arbitrary equation of state. We rewrite Eq. (85) here for net-protons in a HRG model: Normalizing the net proton cumulant by the Skellam baseline, κ Sk 2 [p 1 −p 1 ] = αV T 3 χ p 2 , we get Let us first consider a single conserved charge -the baryon number B. In this case det χ hrg = χ B 2 and detχ hrg = χ B 2 χ p 2 − (χ pB 11 ) 2 . Furthermore, in the HRG model one has χ pB 11 = χ p 2 . Inserting these relations into Eq. (99) one obtains In the HRG model χ p 2 and χ B 2 are proportional to the total number of protons plus antiprotons and to the total number of baryons plus antibaryons, respectively, with the same proportionality factor, so that χ p 2 /χ B 2 = ( N 4π . The subvolume fraction α is proportional to the ratio of mean number of protons plus antiprotons in acceptance to the mean number of protons plus antiprotons in the full space, α = ( N acc in agreement with Eq. (97). One gets an analogous relation for net-Λ fluctuations using the same logic.
Let us now consider the exact conservation of electric charge Q in addition to baryon number B. In this case the susceptibilities involving Q contribute to det χ hrg and detχ hrg in the r.h.s. of Eq. (99). To simplify the resulting expression we shall make use of the fact that resonance decays generate negligibly small proton-charge correlations in addition to proton self-correlation, implying that χ pQ 11 ≈ χ p 2 holds to a large precision. 1 After some algebra, one obtains Here and N 4π ch,prim corresponds to the charged particle multiplicity at the chemical freeze-out stage, i.e. before resonance decays. If contributions of multi-charged hadrons are small, the approximate relation in Eq. (103) becomes exact. Note that the final charged multiplicity N 4π ch,fin measured in the experiment can be considerably larger than N 4π ch,prim , due to decays of neutral resonances into charged particles, e.g. the ρ 0 → π + + π − decay. This effect is significant, and our HRG calculations suggest that N 4π ch,fin can be up to a factor two larger than N 4π ch,prim at RHIC and LHC energies. For this reason a reliable estimation of α p Q can be challenging. Taking N 4π ch,fin in place of N 4π ch,prim in Eq. (103) will underestimate the value of α p Q . For the electrically neutral Λ-hyperon it makes sense to consider global conservation of strangeness instead of electric charge. The simultaneous effect of baryon number and strangeness conservation on net-Λ fluctuations is the following: Here N 4π S,prim counts the total number of strange hadrons. As in the case of net-charge in Eq. (103), N 4π S,prim corresponds to hadrons before resonance decays. The decays, however, lead only to a small distortion, primarily through a φ → K + K − decay, thus N 4π S,prim ≈ N 4π S,fin to a very good approximation. Equations (102) and (104) allow to establish the theoretical basis for the ad hoc relation proposed in Ref. [12] to take into account simultaneous global conservation of two conserved charges. That relation works for net-proton and net-Λ fluctuations if, respectively, baryon-electric and baryon-strangeness correlations in the equation of state can be neglected, i.e.
(107) 1 We note that ∆(1232) ++ → p + π + decays generate an excess of χ pQ 11 over χ p 2 . On the other hand, decays like ∆(1232) 0 → p + π − lead to a reduction of χ pQ 11 relative to χ p 2 . Our HRG model calculations reveal only percent level deviations of χ pQ 11 from χ p 2 after all resonance decays are accounted for. Below we present explicit calculations to establish the accuracy of these relations. We conclude that a quantitative HRG model analysis of net-proton and net-Λ fluctuations requires, in addition to baryon number, a canonical treatment of, respectively, electric charge and strangeness. This may be even more relevant for the higher-order fluctuations, which are expected to be more sensitive to exact conservation of multiple charges. From a practical point of view, the approximate relations (101) and (106),(107) can be used to estimate the magnitude of the global charge conservation effects, as they are found to bracket the true values of net-proton and net-Λ cumulants in an HRG model calculation.

IV. DISCUSSION AND CONCLUSIONS
The main findings of the paper can be described as follows: • The subensemble acceptance method, originally formulated in Ref. [1] for a single conserved charge, has been extended for the case of multiple conserved charges. This allowed to express cumulants of conserved charge distributions measured in a subvolume of a thermal system with a globally conserved charge in terms of the grand canonical susceptibilities for any equation of state. Explicit expressions have been provided for all diagonal and off-diagonal cumulants up to the sixth order, the formalism permits iterative computation of higher-order cumulants as well, if desired.
• A conserved charge cumulant up to third order depends only on a single grand canonical susceptibility with the same indices, but not on any other susceptibility, as follows from Eqs. (39), (40), and (41). Starting from the fourth order, however, the cumulants depend on multiple grand canonical susceptibilities [Eqs. (42)- (44)].
• The global conservation effects cancel out in any ratio of two second order cumulants and in any ratio of two third order cumulants. As follows from Eqs. (40) and (41), such ratios reduce simply to the corresponding ratios of the grand canonical susceptibilities. We verified this statement explicitly using Monte Carlo sampling of the canonical hadron resonance gas model shown in Figs. 1 and 2. We thus argue that such quantities are particularly suitable for experimental measurements, since they allow to eliminate the dependence of results on a relatively difficult-to-constrain value of the acceptance parameter α.
• The kurtosis of a conserved charge distribution is, in general, affected by conservation laws involving other conserved charges, see e.g. Eq. (55). Our HRG model analysis, however, suggests that this effect is small in heavy-ion collisions at √ s NN 40 GeV (see Fig. 3), and likely at lower energies as well. This implies that the kurtosis of a conserved charge is mainly affected by the exact conservation of that charge, while the conservation laws involving other conserved charges have a small effect.
• From an experimental point of view, it is challenging to measure fluctuations of conserved charges other than electric charge Q. In Sec. II H we have extended the SAM to incorporate fluctuations and correlations involving non-conserved quantities, such as net-proton or net-kaon number. Our main result here is Eq. (79): a correlator of non-conserved quantity, such as net-proton number, with a conserved charge, such as Q, is affected by global conservation laws by the same factor as any second order cumulant of conserved charges. As a consequence, the readily measurable cumulant ratios such as κ pQ 11 /κ Q 2 , κ kQ 11 /κ Q 2 , or κ pQ 11 /κ kQ 11 are unaffected by global conservation laws. However, similar ratios involving cumulants of two non-conserved quantities, such as κ πp 11 /κ p 2 do depend on the acceptance. This is a useful observation for present [15] and future measurements of the off-diagonal cumulants of net-particle distributions. We note that experimental measurements should be performed such that the acceptance parameter α is the same for all hadron species that go into the measurement. Ideally, this entails p T -integrated measurements in a finite rapidity Y acceptance, as opposed to the currently available measurements in a finite p T range and/or pseudorapidity η acceptance [15].
• The second order cumulant κ p 2 of a non-conserved quantity, such as e.g. net-proton number p, in acceptance α represents a linear combination between the grand canonical (α → 0) and canonical (α → 1) limits [Eq. (85)]. Furthermore, the canonical susceptibility (cumulant) is expressed solely in terms of the matrix of second order grand canonical susceptibilities involving the conserved charges and a non-conserved quantity [Eq. (74)]. In Sec. III F we used this expression to analyze the influence of the various conservation laws on net-proton and net-Λ fluctuations in a HRG model. We found that, in addition to baryon number conservation, the variances of net-proton and net-Λ fluctuations are markedly influenced by net-charge and net-strangeness conservation, respectively. This is a new element compared to prior HRG model studies [37,38,42] that considered only the effect of baryon number conservation on these quantities.
• Our studies in the present paper have been focused on a scenario where the total system volume is fixed in all events. Event-by-event fluctuations of the system volume, on the other hand, cannot be avoided completely in heavy-ion collisions and they have their own influence on fluctuations of conserved charges. Different methods exist to address these [22,23,38]. In this paper, we have shown in Sec. II G that a ratio of strongly intensive measures Σ and ∆ involving any two conserved charges is insensitive to both the global charge conservation and volume fluctuations. While the physical interpretation of strongly intensive measures is somewhat less straightforward than that of the traditional cumulants, their insensitivity to volume fluctuations makes them useful observables when volume fluctuations are difficult to control. The concept of strongly intensive cumulants, introduced in Ref. [62], can be used to extend these considerations to higher-order fluctuation measures.
• The SAM is formulated to study the effects of global charge conservation on cumulants measured in the coordinate space. Experimental measurements in heavy-ion collisions, on the other hand, are performed in the momentum space. Nevertheless, strong space-momentum correlations at the highest collision energies due to longitudinal Bjorken flow allow to associate measurements in finite rapidity space with spatial subvolumes at the freeze-out stage. Furthermore, the fluctuation measures where global conservation factors were found to cancel out can be expected to be robust probes of the grand canonical susceptibilities even in the absence of strong space-momentum correlations.
• We focused the discussion on the effects of global conservation laws. It is not unfeasible, however, that the exact conservation of charges takes place not only globally, but also in localized spatial regions, as discussed in a number of recent papers [63][64][65][66][67][68][69][70]. The SAM can be applied in such a scenario if these localized spatial regions -called patches according to the terminology developed in Refs. [65,69] -are regarded as total volumes V where the conserved charges are conserved exactly. The only requirement is that patches are sufficiently large to contain all physics associated with the correlation length.
To summarize, we developed a formalism to quantify the effect of global conservation of multiple conserved charges on fluctuation measurements in heavy-ion collisions. In particular, this has allowed us, for the first time, to construct fluctuation measures that, to a leading order, are insensitive to effects of global charge conservation. In the future we plan to apply the concepts developed in this paper to construct a sampler of an interacting hadron resonance gas that preserves local correlations and fluctuations encoded in the equation of state, and which can be used in state-of-the-art hydrodynamic simulations of heavy-ion collisions.
• The derivatives ofχ −1 2 andχ −1 2 . These have already been computed in Sec. II D. The corresponding expressions are given by Eqs. (34) and (35) which we rewrite here for completeness: (A5) • The derivatives ofχ i1...i M andχ i1...i M . These are evaluated by applying the chain rule. The result is: As follows from the rules (A3)-(A7), at-derivative of a term which comprises a convolution of an arbitrary number ofκ,χ −1 2 , andχ elements yields a sum of terms, each again being a certain convolution ofκ,χ −1 2 , andχ elements. Therefore, the rules (A3)-(A7) are sufficient to iteratively computet-derivatives ofκ ijk (t) up to arbitrary high order. Further simplifications can be achieved by observing that all theκ,χ −1 2 , andχ tensors are symmetric with respect to any permutation of their indices.
We have implemented the above rules within the Mathematica package to compute the higher order cumulantsκ. Equations (42), (43), and (44) in the main text depict the final results for the fourth, fifth, and sixth order cumulants, respectively, all evaluated att = 0. i σ3 = i 3 , and i σ4 = i 4 . For the second permutation one has σ = (1, 2, 4, 3), so that i σ1 = i 1 , i σ2 = i 2 , i σ3 = i 4 , and i σ4 = i 3 . And so on for all remaining permutations.
A Mathematica notebook to express any cumulant up to sixth order within the QCD notation is available via [47].