# Internal Cumulants for Femtoscopy with Fixed Charged Multiplicity

Authors: H. C. Eggers; B. Buschbeck
Journal: Advances in High Energy Physics
 (2013)

---

## Abstract

A detailed understanding of all effects and influences on higher-order correlations is essential. At low charged multiplicity, the effect of a non-Poissonian multiplicity distribution can significantly distort correlations. Evidently, the reference samples with respect to which correlations are measured should yield a null result in the absence of correlations. We show how the careful specification of desired properties necessarily leads to an average-of-multinomials reference sample. The resulting internal cumulants and their averaging over several multiplicities fulfill all requirements of correctly taking into account non-Poissonian multiplicity distributions as well as yielding a null result for uncorrelated fixed-N samples. Various correction factors are shown to be approximations at best. Careful rederivation of statistical variances and covariances within the frequentist approach yields errors for cumulants that differ from those used so far. We finally briefly discuss the implementation of the analysis through a multiple event buffer algorithm.


---

## Body

## 1. Introduction and Motivation

The understanding of hadronic collisions is now considered an essential baseline for ultrarelativistic heavy-ion collisions. Given the correspondingly low final-state multiplicities, there are significant deviations, even for inclusive samples, from assumptions commonly made both in the general theory and in the definition of experimentally measured quantities such as a non-Gaussian shape of the correlation function and non-Poissonian multiplicity distributions. Constraints such as energy-momentum conservation [1, 2] would also play a role in at least some regions of phase space. Multiplicity-class and fixed-multiplicity analysis differ increasingly from Poissonian and inclusive distributions, and with the good statistics now available, measurements have become accurate enough to require proper understanding and treatment of these assumptions and deviations, which play an ever larger role with increasing order of correlation.


### 1.1. Correlations as a Function of Charged Multiplicity

There are a number of reasons to study correlations at fixed charged multiplicityN or, if necessary, charged-multiplicity classes. Firstly, the physics of multiparticle correlations will evidently change with N, and indeed the multiplicity dependence of various quantities such as the intercept parameter and radii associated with Gaussian parametrisations is under constant scrutiny [3–12]. Measurement of many observables as a function of multiplicity class, regarded a proxy for centrality dependence, has been routine for years. Corresponding theoretical considerations, for example, in the quantum optical approach go back a long time [13]. Secondly, correlations for fixed-N are the building blocks which are combined into multiplicity-class and inclusive correlations [14].However, such fixed-N correlations have been beset by an inconsistency in that they are nonzero even when the underlying sample is uncorrelated and do not integrate to zero either. This has been recognised from the start [15], and various attempts have been made to fix the problem.Combining events from several fixed-N subsamples into multiplicity classes does not solve these problems. To quote an early reference [16], “Averaging over multiplicities inextricably mixes the properties of the correlation mechanism with those of the multiplicity distribution. Instead, the study of correlations at fixed multiplicities allows one to separate both effects and to investigate the behaviour of correlation functions as a function of multiplicity.” Under the somewhat inappropriate name of “Long-Range Short-Range correlations” [15, 17], an attempt was made to separate these multiplicity-mixing correlations from the fixed-N correlations, but the inconsistencies inherent in the underlying fixed-N correlations were not addressed. Building on [18], we propose doing so now.


### 1.2. Cumulants in Multiparticle Physics

Multiparticle cumulants have entered the mainstream of analysis, as shown by the following incomplete list of topics. In principle, the considerations presented in this paper would apply to any, and all such cumulants to the degree that their reference distribution deviates from a Poisson process or that the type of particle kept fixed differs from the particle being analysed.Integrated cumulants of multiplicity distributions have a long history in multiparticle physics [19]. Second-order differential cumulants, normally termed “correlation functions,” have likewise been ubiquitous for decades [7] both in charged-particles correlations [15] and in femtoscopy since they provide information on spacetime characteristics of the emitting sources, most recently at the LHC [10, 11, 20]. Differential three-particle cumulants generically measure asymmetries in source geometry and exchange amplitude phases [21]. They also provide consistency checks [22] and a tool to disentangle the coherence parameter from other effects [23, 24]. Three-particle cumulants are also sensitive to differences between longitudinal and transverse correlation lengths in the Lund model  [25]. Inclusive three-particle cumulants have been measured, albeit with different methodologies, in, for example, hadronic [26–29], leptonic [30–33], and nuclear collisions [34–37]. They play a central role in direct QCD-based calculations [38–40], in some recent theory, and in experiment of azimuthal and jet-like correlations [41–45]. Net-charge and other charge combinations are considered probes of the QCD phase diagram [46–48]. Cumulants of order 4 or higher are, of course, increasingly difficult to measure, and so early investigations were largely confined to their scale dependence [49–52]. The large event samples now available have, however, made feasible measurements of fourth- and higher-order cumulants in other variables as proposed in [13, 53–56] as, for example, recently measured by ALICE [57]. Reviews of femtoscopy theory range from [58–60] to more recent ones such as [8].


### 1.3. Outline of This Paper

It has long been obvious that the root cause of the problems and inconsistencies set out in Section1.1 was the reference sample [61]. Insofar as cumulants are concerned, the solution was outlined in [18] as a subtraction of the reference sample cumulant from the measured one; important pieces of the puzzle were, however, still missing at that stage. In this paper, we clarify and extend the basic concept of internal cumulants and consider in detail the case of second- and third-order differential cumulants in the invariant Q=-(p1-p2)2 for fixed charged multiplicity N. The method may be implemented for other variables without much fuss.A second cornerstone of the present paper is the recognition that then particles which enter a correlation analysis are usually only a subset of the N charged pions. While in the case of charged-particle correlations, all N particles are used in the analysis, Bose-Einstein correlations, for example, would use only the n≡n+ positive pions (and, in a separate analysis, only the n-=N-n+ negatives). In addition, there may be reasons to restrict the analysis itself to subregions of the total acceptance Ω in which N was measured, as exemplified in this paper by restriction to a “good azimuthal region” subinterval around the beam axis, 𝒜⊂[0,2π], in which detection efficiency is high. 𝒜 can, however, be reinterpreted generically as any restriction in momentum space compared to Ω and/or as a selection such as charge or particle species. Even when setting 𝒜=Ω, that is, doing the femtoscopy analysis in the full acceptance, n still does not equal N but fluctuates around N/2. The trivial observation that n≠N fundamentally changes the analysis: identical-particle correlations at fixedNand charged-particle correlations at fixedNrequire different definitions.As we will show, ad hoc prescriptions such as simply inserting prefactors or implementing event mixing using only events of the sameN do alleviate the effect of the overall non-Poissonian multiplicity distribution in part but fail to remove them completely. The same issues will, of course, arise in any other correlation type of, for example, nonidentical particles or net charge correlations. The formalism set out here can be easily extended to such cases. A refined version of the abovementioned Long-Range-Short-Range method, which we term “Averaged-Internal” cumulants, will be presented in Section 5. Along the way, we document in Section 2 extended versions of the particle counters [62, 63] which we need as the basis for correlation studies and in Section 4 demonstrate from first principles that statistical errors for cumulants used so far have captured only some of the terms and with partly incorrect prefactors. Section 6 outlines the implementation of event mixing for fixed-N analysis. While experimental results will be published elsewhere, preliminary results in Figures 2 and 3 show that, in third and even in second order, corrections due to proper treatment of fixed-N reference samples can be large.


## 1.1. Correlations as a Function of Charged Multiplicity

There are a number of reasons to study correlations at fixed charged multiplicityN or, if necessary, charged-multiplicity classes. Firstly, the physics of multiparticle correlations will evidently change with N, and indeed the multiplicity dependence of various quantities such as the intercept parameter and radii associated with Gaussian parametrisations is under constant scrutiny [3–12]. Measurement of many observables as a function of multiplicity class, regarded a proxy for centrality dependence, has been routine for years. Corresponding theoretical considerations, for example, in the quantum optical approach go back a long time [13]. Secondly, correlations for fixed-N are the building blocks which are combined into multiplicity-class and inclusive correlations [14].However, such fixed-N correlations have been beset by an inconsistency in that they are nonzero even when the underlying sample is uncorrelated and do not integrate to zero either. This has been recognised from the start [15], and various attempts have been made to fix the problem.Combining events from several fixed-N subsamples into multiplicity classes does not solve these problems. To quote an early reference [16], “Averaging over multiplicities inextricably mixes the properties of the correlation mechanism with those of the multiplicity distribution. Instead, the study of correlations at fixed multiplicities allows one to separate both effects and to investigate the behaviour of correlation functions as a function of multiplicity.” Under the somewhat inappropriate name of “Long-Range Short-Range correlations” [15, 17], an attempt was made to separate these multiplicity-mixing correlations from the fixed-N correlations, but the inconsistencies inherent in the underlying fixed-N correlations were not addressed. Building on [18], we propose doing so now.


## 1.2. Cumulants in Multiparticle Physics

Multiparticle cumulants have entered the mainstream of analysis, as shown by the following incomplete list of topics. In principle, the considerations presented in this paper would apply to any, and all such cumulants to the degree that their reference distribution deviates from a Poisson process or that the type of particle kept fixed differs from the particle being analysed.Integrated cumulants of multiplicity distributions have a long history in multiparticle physics [19]. Second-order differential cumulants, normally termed “correlation functions,” have likewise been ubiquitous for decades [7] both in charged-particles correlations [15] and in femtoscopy since they provide information on spacetime characteristics of the emitting sources, most recently at the LHC [10, 11, 20]. Differential three-particle cumulants generically measure asymmetries in source geometry and exchange amplitude phases [21]. They also provide consistency checks [22] and a tool to disentangle the coherence parameter from other effects [23, 24]. Three-particle cumulants are also sensitive to differences between longitudinal and transverse correlation lengths in the Lund model  [25]. Inclusive three-particle cumulants have been measured, albeit with different methodologies, in, for example, hadronic [26–29], leptonic [30–33], and nuclear collisions [34–37]. They play a central role in direct QCD-based calculations [38–40], in some recent theory, and in experiment of azimuthal and jet-like correlations [41–45]. Net-charge and other charge combinations are considered probes of the QCD phase diagram [46–48]. Cumulants of order 4 or higher are, of course, increasingly difficult to measure, and so early investigations were largely confined to their scale dependence [49–52]. The large event samples now available have, however, made feasible measurements of fourth- and higher-order cumulants in other variables as proposed in [13, 53–56] as, for example, recently measured by ALICE [57]. Reviews of femtoscopy theory range from [58–60] to more recent ones such as [8].


## 1.3. Outline of This Paper

It has long been obvious that the root cause of the problems and inconsistencies set out in Section1.1 was the reference sample [61]. Insofar as cumulants are concerned, the solution was outlined in [18] as a subtraction of the reference sample cumulant from the measured one; important pieces of the puzzle were, however, still missing at that stage. In this paper, we clarify and extend the basic concept of internal cumulants and consider in detail the case of second- and third-order differential cumulants in the invariant Q=-(p1-p2)2 for fixed charged multiplicity N. The method may be implemented for other variables without much fuss.A second cornerstone of the present paper is the recognition that then particles which enter a correlation analysis are usually only a subset of the N charged pions. While in the case of charged-particle correlations, all N particles are used in the analysis, Bose-Einstein correlations, for example, would use only the n≡n+ positive pions (and, in a separate analysis, only the n-=N-n+ negatives). In addition, there may be reasons to restrict the analysis itself to subregions of the total acceptance Ω in which N was measured, as exemplified in this paper by restriction to a “good azimuthal region” subinterval around the beam axis, 𝒜⊂[0,2π], in which detection efficiency is high. 𝒜 can, however, be reinterpreted generically as any restriction in momentum space compared to Ω and/or as a selection such as charge or particle species. Even when setting 𝒜=Ω, that is, doing the femtoscopy analysis in the full acceptance, n still does not equal N but fluctuates around N/2. The trivial observation that n≠N fundamentally changes the analysis: identical-particle correlations at fixedNand charged-particle correlations at fixedNrequire different definitions.As we will show, ad hoc prescriptions such as simply inserting prefactors or implementing event mixing using only events of the sameN do alleviate the effect of the overall non-Poissonian multiplicity distribution in part but fail to remove them completely. The same issues will, of course, arise in any other correlation type of, for example, nonidentical particles or net charge correlations. The formalism set out here can be easily extended to such cases. A refined version of the abovementioned Long-Range-Short-Range method, which we term “Averaged-Internal” cumulants, will be presented in Section 5. Along the way, we document in Section 2 extended versions of the particle counters [62, 63] which we need as the basis for correlation studies and in Section 4 demonstrate from first principles that statistical errors for cumulants used so far have captured only some of the terms and with partly incorrect prefactors. Section 6 outlines the implementation of event mixing for fixed-N analysis. While experimental results will be published elsewhere, preliminary results in Figures 2 and 3 show that, in third and even in second order, corrections due to proper treatment of fixed-N reference samples can be large.


## 2. Raw Data, Counters, and Densities

### 2.1. Raw Data

The starting point for experimental correlation analysis is the inclusive sample𝒮, made up of ℰ events a=1,…,ℰ. Each event consists of a varying number of final-state elementary particles and photons; for our purposes, we consider only the N(a) charged pions of event a in Ω, the maximal acceptance region used. Each pion i=1,…,N(a) is characterised by a data vector (Pia,sia,eia) containing its measured information, including the three components of its momentum, Pia=(pixa,piya,piza)or(yia,ϕia,ptia), while its discrete attributes such as mass and charge are captured in a data vector sia of discrete values; for the moment, we will consider only the charge, sia→cia. From the sample's raw data, we can immediately find derived quantities such as the total charged multiplicity N(a) and the total transverse energy, and such derived quantities are hence considered part of the raw data. The list of particle attributes should be augmented by an error vector eia containing the measurement errors for each track, but we will not consider detector resolution errors here. In summary, the inclusive data sample is fully described in terms of 𝒮={Pa,sa}a=1ℰ consisting of lists of vectors in continuous and discrete spaces
(1)Pa={P1a,P2a,…,PN(a)a},sa={s1a,s2a,…,sN(a)a}.


### 2.2. Data after Conditioning and Cuts

For a particular analysis, the inclusive sample is invariably subdivided and modified through “conditioning,” the statistics terminology for semiinclusive or triggered analysis. From the total sample of events, a subsample is selected according to some restriction or precondition. In our case, this conditioning proceeds in the following steps.(i)Conditioning into Fixed-N Subsamples. For the fixed-multiplicity analyses that form the subject of this paper, 𝒮 is subdivided into a set of fixed-N subsamples 𝒮N, each of which contains only events a whose measured multiplicity N(a) is equal to the constant N characterising 𝒮N={Pa,sa∣δ(N,N(a))}, N=0,1,2,… We use the vertical bar ∣ here and everywhere in the usual sense of “conditioning” whereby the events in sample 𝒮N must satisfy the condition that their charged multiplicity must equal the specified constant N, denoted in this case by the Kronecker delta δ(N,N(a)). Quantities to the right of the vertical bar are generally considered known and fixed, while quantities to the left of the bar are variable or unknown. The number of events in 𝒮N equals the δ-restricted sum over the inclusive sample,
(2)ℰN=∑a=1ℰ‍δ(N,N(a)),∑N=0∞‍ℰN=ℰ.
The usual multiplicity distribution is the list of relative frequencies,1(3)ℛN=ℰNℰ,∑N=0∞‍ℛN=1.
While desirable, it is not easy to measure the total multiplicity of final-state charged pions, a quantity which approximately tracks the variation in the physics. Choosing charged pions measured within the maximal detector acceptance N=N(Ω) as marker is in any case only an approximation because it excludes charged particles outside the primary cuts and also ignores final-state particles other than charged pions. Nevertheless, we expect N to be a reasonable measure of the multiplicity dependence of the physics. Alternatively, the multiplicity density in pseudorapidity at central rapidities dN/dη can be used as a model-dependent proxy for N.(ii)Azimuthal Cut. While N(a) is the charged multiplicity measured in Ω, there is no a priori reason why the correlation analysis itself may not be conducted within a restricted part 𝒜⊂Ω of momentum space within which the actual analysis is done. In the case of the UA1 detector from which the data used in the examples below was drawn, 𝒜 refers to azimuthal regions within which measurement efficiency was high, and pions found in the low-efficiency azimuthal regions were excluded. Correspondingly, the multiplicity n(a) which enters the correlation analysis itself differs from N(a) and will, for a given fixed-N, fluctuate with relative frequency
(4)ℛnN=ℰnNℰNforeachfixedN=1,2,…,
where ℰnN is the number of events for which N(a)=N and n(a)=n. The outcome space for n(a) will depend on its definition; in the present case where only positive (or only negative) pions within 𝒜 are used in the analysis, it will be [0,1,…,N] so that the relative frequency is normalised by
(5)∑n=0N‍ℛnN=1∀N=1,2,….
With approximate charge conservation n+≃n-, we expect the fixed-N average for positive (or negative) pions in 𝒜 to hover around
(6)〈n〉N≃N2volumeof𝒜volumeofΩ.
An example of the resulting relative frequencies (conditional normalised multiplicity distributions) is shown in Figure 1. Since n≤N, these conditional multiplicity distributions are almost always sub-Poissonian, that is, narrower than a Poisson distribution with the same 〈n〉N would be.Figure 1
Conditional normalised multiplicity distribution (relative frequency) of the number of positive pionsn in restricted azimuthal region 𝒜={20∘≤|ϕ|≤160∘} for, respectively, fixed charged multiplicity N= 3, 5, 10, 20, 30, 40, and 50, for the UA1 dataset used in [65]. In accordance with (6), 〈n〉N≃0.39N.(a) Correction factorsF2N=〈n(n-1)〉N/〈n〉N2 (red squares) and F3N=〈n(n-1)(n-2)〉N/〈n〉N3 (blue circles) as defined in (57) as well as G3N (green triangles) of (70), for the UA1 dataset used in [65]. (b) Inverse factors as used in the normalisation of internal cumulants (74)-(75) and (83)-(84).
(a)
(b)Figure 3
Comparison, for the same UA1 data and cuts, of(1/F2N)-1 (red filled squares), (1/F3N)-1 (blue filled circles), and both averages over n for fixed N, with prefactors [N2/N(N-1)]-1 (red open squares) and [N3/N(N-1)(N-2)]-1 (blue open circles). The Poissonian limit is represented by zero on the y-axis; that is, 100 times the y-scale represents the percentage deviation from the Poissonian limit.(iii)Generalisation. While in this paper the analysis will be carried out for the n(a) positive pions of event a falling into 𝒜, the same formalism obviously applies to negative pions and may equally refer to any other particles such as kaons, baryons, and photons, in any combination which depends on N. There is no a priori connection between the definitions of N and n.(iv)Identical-Particle versus Multispecies Analysis. While we do not develop the formalism for correlations between two or three particles of different species or charge, the methodology developed here can be easily modified to deal with such cases. For example, positive-negative pion combinations and “charge balance correlations” [64] can be handled by inserting delta functions δ(c,cia), where c is the desired charge and cia the measured charge of track i in event a, into the definitions of the counters in Section 2.3.


### 2.3. Counters and Densities for FixedN

This section is based on an old formalism [53, 62, 63] which must, however, be updated to accommodate the issues being considered here. The basic building block of correlation analysis is the counter; it is a particular projection of the raw data particularly suited to the construction of histograms. Eventwise counters ρ^ for a given event a are averaged to give sample counters ρ.We take the simple case where eventa contains N(a) tracks with three-momenta Pa={P1a,…,PN(a)a}, no discrete attributes s and no further cuts or selection. For each point p1 in momentum space, only that particle i (if any) whose momentum Pia happens to coincide with p is to be counted,
(7)ρ^(p1∣Pa)=∑i=1N(a)δ(p1-Pia).
Such counters always appear under an integral over some region of the P space, so that the delta functions fulfill the purpose of counting those particles falling within that region. Alternatively, one can consider the delta functions here and below to represent small nonoverlapping intervals around the specified momenta. The integral over the full momentum space Ω yields
(8)∫Ω‍dp1ρ^(p1∣Pa)=N(a),
while an integral over some subspace or bin Ωb⊂Ω will yield the number of particles of event a in bin Ωb. The second-order eventwise counter for event a is2(9)ρ^(p1,p2∣Pa)=∑i≠j=1N(a)‍δ(p1-Pia)δ(p2-Pja),
with the inequality i≠j ensuring that a single particle is not counted as a “pair.” The counter integrates to
(10)∫Ω‍dp1dp2ρ^(p1,p2∣Pa)=N(a)(N(a)-1)=N(a)2_,
using the falling factorial notation
(11)Nr_=N(N-1)⋯(N-r+1)r=1,2,…,
as contrasted to the rising factorial (Pochhammer symbol) Nr¯=N(N+1)⋯(N+r-1). The single-particle counter is a projection of ρ^(p1,p2∣Pa) because
(12)∫Ω‍dp2ρ^(p1,p2∣Pa)=[N(a)-1]ρ^(p1∣Pa).
The most general eventwise counter which enters the exclusive cross-section for events with charged multiplicity N(a)(13)ρ^(p1,p2,…,pN(a)∣Pa)=∑i1≠i2≠⋯≠iN(a)=1N(a)‍∏d=1N(a)‍δ(pd-Pida)
fully describes the event, including any and all correlations between its particles. It integrates to the factorial of the event multiplicity
(14)∫Ω‍dp1dp2⋯dpN(a)ρ^(p1,p2,…,pN(a)∣Pa)=N(a)!
and contains all counters of lower order by projection. An rth-order counter ρ^(p1,p2,…,pr|Pa) is zero whenever there are more observation points than particles being observed, r>N(a).3To distinguish eventwise counters for nonfixed-N from eventwise counters for fixed N, we define the separate eventwise counter for fixed N by specifying an additional Kronecker delta,
(15)ρ^(p1,…,pr∣Pa,N)=δ(N,N(a))∑i1≠⋯≠irN‍∏d‍δ(pd-Pida)r=1,2,…,N.While the counters and densities defined above and below are clearly frame dependent, it is easy to define corresponding Lorentz-invariant versions by supplementing each delta function in 3 momenta with the corresponding energy; thus (7) would become, for example,
(16)ρ^(p1∣Pa)=∑i=1N(a)‍E(p1)δ(p1-Pia),
with E(p1)=E1=p12+m2 the on-shell energy, and in general
(17)ρ^(p1,p2,…,pr∣Pa)=∑i1≠i2≠⋯≠ir=1N(a)‍∏d=1N(a)‍Edδ(pd-Pida)r=1,…,N(a),ρ^(p1,…,pr∣Pa,N)=δ(N,N(a))∑i1≠⋯≠irN∏d‍Edδ(pd-Pida)r=1,…,N,
which are manifestly invariant. Because such counters and densities are, however, always integrated over some Ωb by ∏d‍(dpd/Ed), the additional factors Ed always cancel and play no role on this level of analysis and will be ignored for the time being. The bin boundaries of Ωb do, however, remain frame dependent.Charge-, spin-, or species-specific counters are defined in the same way, that is, by supplying appropriate Kronecker deltas to the counters; for example, the particle counter for pions with chargec1 at p1 for fixed N is
(18)ρ^(p1∣c1,Pa,N)=δ(N,N(a))∑i=1N(a)‍δ(c1,cia)δ(p1-Pia),
while the two-particle counter for charge combination (c1,c2) at momenta (p1,p2) for any N is, for example,
(19)ρ^(p1,p2∣c1,c2,Pa)=∑i,j=1N(a)‍δ(c1,cia)δ(p1-Pia)δ(c2,cja)δ(p2-Pja).
In contrast to (18), charge counters rather than particle counters would be
(20)ρ^c(p1∣c1,Pa)=∑i=1N(a)‍c1δ(c1,cia)δ(p1-Pia),
so that ρ^c(p1∣+1,Pa)+ρ^c(p1∣-1,Pa) represents the net charge of event a at p1. The two-particle counter for charges (c1,c2) at momenta (p1,p2) is
(21)ρ^c1c2≡ρ^c(p1,p2∣c1,c2,Pa)=∑i,j=1N(a)‍c1δ(c1,cia)δ(p1-Pia)c2δ(c2,cja)δ(p2-Pja),
and “charge flow” correlations can be constructed from this (for rapidities (y,y′) in the case of [66]) such as
(22)Φ(y,y′)=-〈∑i≠j‍ciacjaδ(y-Yia)δ(y′-Yja)〉,
which can be expressed as Φ=〈ρ^+-+ρ^+--ρ^++-ρ^--〉 and the related “charge balance functions” described in for example [46].Returning to the fixed-N case, eventwise counters will usually be combined with similar events to form event averages. The simplest average is the fixed-N density for the subsample of fixed 𝒮N,
(23)ρ(p1∣𝒮N)=1ℰN∑a‍ρ^(p1∣Pa,N),
with the Kronecker delta in (15) ensuring that only events in 𝒮N are considered, so we need not further specify the individual terms or limits of the a-sum. Using (2), it is immediately clear that
(24)∫Ω‍dp1ρ(p1∣𝒮N)=N
compared to the integral over the corresponding eventwise counter
(25)∫Ω‍dpρ^(p∣Pa,N)=Nδ(N,N(a))
and to the integral (8); similarly
(26)∫Ω‍dp1⋯dprρ(p1,…,pr∣𝒮N)=N(N-1)⋯(N-r+1)=Nr_.
The inclusive averaged density ρ(p1∣𝒮) is the weighted average over all N of the fixed-N averages,
(27)ρ(p1∣𝒮)=∑N=1∞‍ℛNρ(p1∣𝒮N).
Using (3), (15), and (23), this can be written as
(28)ρ(p1∣𝒮)=1ℰ∑a‍∑i‍δ(p1-Pia),
and for general r=1,2,…(29)ρ(p1,…,pr∣𝒮)=∑N=r∞‍ℛNρ(p1,…,pr∣𝒮N)=1ℰ∑a‍∑i1≠⋯≠ir‍∏d=1r‍δ(pd-Pida),
keeping in mind that ρ^(p1,…,pr∣𝒮N) will be zero whenever N(a)<r. The integral of any rth-order inclusive averaged density is the rth-order factorial moment of the multiplicity distribution,
(30)∫Ω‍dp1⋯dprρ(p1,…,pr∣𝒮)=∑N=r∞‍ℛNNr_=〈Nr_〉,
with simple angle brackets denoting inclusive averaging.The averaged counters are of course directly related to the traditional definitions in terms of cross-sections. Ifℒ is the integrated luminosity of incoming particles, the topological cross-section is σN=ℰN/ℒ, the inelastic cross-section is σI=ℰ/ℒ, and the inclusive cross-section is σincl=∑N‍NσN=〈N〉σI while the relative frequency (multiplicity distribution) can be written as usual as ℛN=σN/σI. The relation between the differential cross sections and our counters is
(31)ρincl(p1,…,pN)=ρ(p1,…,pN∣𝒮)=1σId3Nσincldp1⋯dpN,ρ(p1,…,pN∣𝒮N)=1σNd3Nσexcldp1⋯dpN,
and so as usual inclusive and exclusive densities are related by [14]
(32)ρ(p1,…,pN∣𝒮)=∑N=r∞‍ℛN(N-r)!∫‍ρ(p1,…,pN∣𝒮N)dpr+1⋯dpN,
while the semi-inclusive cross sections and counters follow by the usual projections.


### 2.4. Counters and Densities for Fixed(N,n)

Our choice of a basic counter is motivated by the experimental situation set out in Section1: we wish to work in event subsamples of fixed total charged multiplicity N(a) in the entire momentum space Ω but do the differential correlation analysis using only those pions n which fall into the restricted space 𝒜 and of a particular charge +1 or -1. This requires the use of “sub-subsamples” for which both N and n are kept fixed,
(33)𝒮nN={Pa,withaconstrainedbyδ(n,n+(a))δ(N,N(a))},
with n+(a) being the number of positive pions of event a in 𝒜, and eventwise sub-subsample counters
(34)ρ^(p1,…,pr∣n,N,Pa)=δ(n,n+(a))δ(N,N(a))×∑i1≠⋯≠ir=1n‍δ(p1-Pi1a)⋯δ(pr-Pira).
As in (2), the number of events in a sub-subsample ℰnN=∑a‍δ(n,n+(a))δ(N,N(a)) enters the relevant event averages
(35)ρ(p1,…,pr∣𝒮nN)=1ℰnN∑a‍ρ^(p1,…,pr∣n,N,Pa),
where once again the double Kronecker deltas in (34) ensure selection of events in 𝒮nN only. Integrals of the counters over the good-azimuth region 𝒜 yield, for the eventwise and 𝒮nN-averaged counters,
(36)∫𝒜‍dp1⋯dprρ^(p1,…,pr∣n,N,Pa)=nr_δ(n,n+(a))δ(N,N(a)),∫𝒜‍dp1⋯dprρ(p1,…,pr∣𝒮nN)=nr_.
Bearing in mind that observation points p1,p2,… refer to positive pions in 𝒜 only, the event-averaged counters for fixed N but any n are given by the average weighted in terms of the relative frequency ℛnN=ℰnN/ℰN,
(37)ρ(p1,…,pr∣𝒮N)=∑n=rN‍ℛnNρ(p1,…,pr∣𝒮nN)=〈ρ(p1,…,pr∣𝒮nN)〉N,
for r=1,2,3,…, which integrate to
(38)∫𝒜‍dp1⋯dprρ(p1,…,pr∣𝒮N)=〈nr_〉N.


## 2.1. Raw Data

The starting point for experimental correlation analysis is the inclusive sample𝒮, made up of ℰ events a=1,…,ℰ. Each event consists of a varying number of final-state elementary particles and photons; for our purposes, we consider only the N(a) charged pions of event a in Ω, the maximal acceptance region used. Each pion i=1,…,N(a) is characterised by a data vector (Pia,sia,eia) containing its measured information, including the three components of its momentum, Pia=(pixa,piya,piza)or(yia,ϕia,ptia), while its discrete attributes such as mass and charge are captured in a data vector sia of discrete values; for the moment, we will consider only the charge, sia→cia. From the sample's raw data, we can immediately find derived quantities such as the total charged multiplicity N(a) and the total transverse energy, and such derived quantities are hence considered part of the raw data. The list of particle attributes should be augmented by an error vector eia containing the measurement errors for each track, but we will not consider detector resolution errors here. In summary, the inclusive data sample is fully described in terms of 𝒮={Pa,sa}a=1ℰ consisting of lists of vectors in continuous and discrete spaces
(1)Pa={P1a,P2a,…,PN(a)a},sa={s1a,s2a,…,sN(a)a}.


## 2.2. Data after Conditioning and Cuts

For a particular analysis, the inclusive sample is invariably subdivided and modified through “conditioning,” the statistics terminology for semiinclusive or triggered analysis. From the total sample of events, a subsample is selected according to some restriction or precondition. In our case, this conditioning proceeds in the following steps.(i)Conditioning into Fixed-N Subsamples. For the fixed-multiplicity analyses that form the subject of this paper, 𝒮 is subdivided into a set of fixed-N subsamples 𝒮N, each of which contains only events a whose measured multiplicity N(a) is equal to the constant N characterising 𝒮N={Pa,sa∣δ(N,N(a))}, N=0,1,2,… We use the vertical bar ∣ here and everywhere in the usual sense of “conditioning” whereby the events in sample 𝒮N must satisfy the condition that their charged multiplicity must equal the specified constant N, denoted in this case by the Kronecker delta δ(N,N(a)). Quantities to the right of the vertical bar are generally considered known and fixed, while quantities to the left of the bar are variable or unknown. The number of events in 𝒮N equals the δ-restricted sum over the inclusive sample,
(2)ℰN=∑a=1ℰ‍δ(N,N(a)),∑N=0∞‍ℰN=ℰ.
The usual multiplicity distribution is the list of relative frequencies,1(3)ℛN=ℰNℰ,∑N=0∞‍ℛN=1.
While desirable, it is not easy to measure the total multiplicity of final-state charged pions, a quantity which approximately tracks the variation in the physics. Choosing charged pions measured within the maximal detector acceptance N=N(Ω) as marker is in any case only an approximation because it excludes charged particles outside the primary cuts and also ignores final-state particles other than charged pions. Nevertheless, we expect N to be a reasonable measure of the multiplicity dependence of the physics. Alternatively, the multiplicity density in pseudorapidity at central rapidities dN/dη can be used as a model-dependent proxy for N.(ii)Azimuthal Cut. While N(a) is the charged multiplicity measured in Ω, there is no a priori reason why the correlation analysis itself may not be conducted within a restricted part 𝒜⊂Ω of momentum space within which the actual analysis is done. In the case of the UA1 detector from which the data used in the examples below was drawn, 𝒜 refers to azimuthal regions within which measurement efficiency was high, and pions found in the low-efficiency azimuthal regions were excluded. Correspondingly, the multiplicity n(a) which enters the correlation analysis itself differs from N(a) and will, for a given fixed-N, fluctuate with relative frequency
(4)ℛnN=ℰnNℰNforeachfixedN=1,2,…,
where ℰnN is the number of events for which N(a)=N and n(a)=n. The outcome space for n(a) will depend on its definition; in the present case where only positive (or only negative) pions within 𝒜 are used in the analysis, it will be [0,1,…,N] so that the relative frequency is normalised by
(5)∑n=0N‍ℛnN=1∀N=1,2,….
With approximate charge conservation n+≃n-, we expect the fixed-N average for positive (or negative) pions in 𝒜 to hover around
(6)〈n〉N≃N2volumeof𝒜volumeofΩ.
An example of the resulting relative frequencies (conditional normalised multiplicity distributions) is shown in Figure 1. Since n≤N, these conditional multiplicity distributions are almost always sub-Poissonian, that is, narrower than a Poisson distribution with the same 〈n〉N would be.Figure 1
Conditional normalised multiplicity distribution (relative frequency) of the number of positive pionsn in restricted azimuthal region 𝒜={20∘≤|ϕ|≤160∘} for, respectively, fixed charged multiplicity N= 3, 5, 10, 20, 30, 40, and 50, for the UA1 dataset used in [65]. In accordance with (6), 〈n〉N≃0.39N.(a) Correction factorsF2N=〈n(n-1)〉N/〈n〉N2 (red squares) and F3N=〈n(n-1)(n-2)〉N/〈n〉N3 (blue circles) as defined in (57) as well as G3N (green triangles) of (70), for the UA1 dataset used in [65]. (b) Inverse factors as used in the normalisation of internal cumulants (74)-(75) and (83)-(84).
(a)
(b)Figure 3
Comparison, for the same UA1 data and cuts, of(1/F2N)-1 (red filled squares), (1/F3N)-1 (blue filled circles), and both averages over n for fixed N, with prefactors [N2/N(N-1)]-1 (red open squares) and [N3/N(N-1)(N-2)]-1 (blue open circles). The Poissonian limit is represented by zero on the y-axis; that is, 100 times the y-scale represents the percentage deviation from the Poissonian limit.(iii)Generalisation. While in this paper the analysis will be carried out for the n(a) positive pions of event a falling into 𝒜, the same formalism obviously applies to negative pions and may equally refer to any other particles such as kaons, baryons, and photons, in any combination which depends on N. There is no a priori connection between the definitions of N and n.(iv)Identical-Particle versus Multispecies Analysis. While we do not develop the formalism for correlations between two or three particles of different species or charge, the methodology developed here can be easily modified to deal with such cases. For example, positive-negative pion combinations and “charge balance correlations” [64] can be handled by inserting delta functions δ(c,cia), where c is the desired charge and cia the measured charge of track i in event a, into the definitions of the counters in Section 2.3.


## 2.3. Counters and Densities for FixedN

This section is based on an old formalism [53, 62, 63] which must, however, be updated to accommodate the issues being considered here. The basic building block of correlation analysis is the counter; it is a particular projection of the raw data particularly suited to the construction of histograms. Eventwise counters ρ^ for a given event a are averaged to give sample counters ρ.We take the simple case where eventa contains N(a) tracks with three-momenta Pa={P1a,…,PN(a)a}, no discrete attributes s and no further cuts or selection. For each point p1 in momentum space, only that particle i (if any) whose momentum Pia happens to coincide with p is to be counted,
(7)ρ^(p1∣Pa)=∑i=1N(a)δ(p1-Pia).
Such counters always appear under an integral over some region of the P space, so that the delta functions fulfill the purpose of counting those particles falling within that region. Alternatively, one can consider the delta functions here and below to represent small nonoverlapping intervals around the specified momenta. The integral over the full momentum space Ω yields
(8)∫Ω‍dp1ρ^(p1∣Pa)=N(a),
while an integral over some subspace or bin Ωb⊂Ω will yield the number of particles of event a in bin Ωb. The second-order eventwise counter for event a is2(9)ρ^(p1,p2∣Pa)=∑i≠j=1N(a)‍δ(p1-Pia)δ(p2-Pja),
with the inequality i≠j ensuring that a single particle is not counted as a “pair.” The counter integrates to
(10)∫Ω‍dp1dp2ρ^(p1,p2∣Pa)=N(a)(N(a)-1)=N(a)2_,
using the falling factorial notation
(11)Nr_=N(N-1)⋯(N-r+1)r=1,2,…,
as contrasted to the rising factorial (Pochhammer symbol) Nr¯=N(N+1)⋯(N+r-1). The single-particle counter is a projection of ρ^(p1,p2∣Pa) because
(12)∫Ω‍dp2ρ^(p1,p2∣Pa)=[N(a)-1]ρ^(p1∣Pa).
The most general eventwise counter which enters the exclusive cross-section for events with charged multiplicity N(a)(13)ρ^(p1,p2,…,pN(a)∣Pa)=∑i1≠i2≠⋯≠iN(a)=1N(a)‍∏d=1N(a)‍δ(pd-Pida)
fully describes the event, including any and all correlations between its particles. It integrates to the factorial of the event multiplicity
(14)∫Ω‍dp1dp2⋯dpN(a)ρ^(p1,p2,…,pN(a)∣Pa)=N(a)!
and contains all counters of lower order by projection. An rth-order counter ρ^(p1,p2,…,pr|Pa) is zero whenever there are more observation points than particles being observed, r>N(a).3To distinguish eventwise counters for nonfixed-N from eventwise counters for fixed N, we define the separate eventwise counter for fixed N by specifying an additional Kronecker delta,
(15)ρ^(p1,…,pr∣Pa,N)=δ(N,N(a))∑i1≠⋯≠irN‍∏d‍δ(pd-Pida)r=1,2,…,N.While the counters and densities defined above and below are clearly frame dependent, it is easy to define corresponding Lorentz-invariant versions by supplementing each delta function in 3 momenta with the corresponding energy; thus (7) would become, for example,
(16)ρ^(p1∣Pa)=∑i=1N(a)‍E(p1)δ(p1-Pia),
with E(p1)=E1=p12+m2 the on-shell energy, and in general
(17)ρ^(p1,p2,…,pr∣Pa)=∑i1≠i2≠⋯≠ir=1N(a)‍∏d=1N(a)‍Edδ(pd-Pida)r=1,…,N(a),ρ^(p1,…,pr∣Pa,N)=δ(N,N(a))∑i1≠⋯≠irN∏d‍Edδ(pd-Pida)r=1,…,N,
which are manifestly invariant. Because such counters and densities are, however, always integrated over some Ωb by ∏d‍(dpd/Ed), the additional factors Ed always cancel and play no role on this level of analysis and will be ignored for the time being. The bin boundaries of Ωb do, however, remain frame dependent.Charge-, spin-, or species-specific counters are defined in the same way, that is, by supplying appropriate Kronecker deltas to the counters; for example, the particle counter for pions with chargec1 at p1 for fixed N is
(18)ρ^(p1∣c1,Pa,N)=δ(N,N(a))∑i=1N(a)‍δ(c1,cia)δ(p1-Pia),
while the two-particle counter for charge combination (c1,c2) at momenta (p1,p2) for any N is, for example,
(19)ρ^(p1,p2∣c1,c2,Pa)=∑i,j=1N(a)‍δ(c1,cia)δ(p1-Pia)δ(c2,cja)δ(p2-Pja).
In contrast to (18), charge counters rather than particle counters would be
(20)ρ^c(p1∣c1,Pa)=∑i=1N(a)‍c1δ(c1,cia)δ(p1-Pia),
so that ρ^c(p1∣+1,Pa)+ρ^c(p1∣-1,Pa) represents the net charge of event a at p1. The two-particle counter for charges (c1,c2) at momenta (p1,p2) is
(21)ρ^c1c2≡ρ^c(p1,p2∣c1,c2,Pa)=∑i,j=1N(a)‍c1δ(c1,cia)δ(p1-Pia)c2δ(c2,cja)δ(p2-Pja),
and “charge flow” correlations can be constructed from this (for rapidities (y,y′) in the case of [66]) such as
(22)Φ(y,y′)=-〈∑i≠j‍ciacjaδ(y-Yia)δ(y′-Yja)〉,
which can be expressed as Φ=〈ρ^+-+ρ^+--ρ^++-ρ^--〉 and the related “charge balance functions” described in for example [46].Returning to the fixed-N case, eventwise counters will usually be combined with similar events to form event averages. The simplest average is the fixed-N density for the subsample of fixed 𝒮N,
(23)ρ(p1∣𝒮N)=1ℰN∑a‍ρ^(p1∣Pa,N),
with the Kronecker delta in (15) ensuring that only events in 𝒮N are considered, so we need not further specify the individual terms or limits of the a-sum. Using (2), it is immediately clear that
(24)∫Ω‍dp1ρ(p1∣𝒮N)=N
compared to the integral over the corresponding eventwise counter
(25)∫Ω‍dpρ^(p∣Pa,N)=Nδ(N,N(a))
and to the integral (8); similarly
(26)∫Ω‍dp1⋯dprρ(p1,…,pr∣𝒮N)=N(N-1)⋯(N-r+1)=Nr_.
The inclusive averaged density ρ(p1∣𝒮) is the weighted average over all N of the fixed-N averages,
(27)ρ(p1∣𝒮)=∑N=1∞‍ℛNρ(p1∣𝒮N).
Using (3), (15), and (23), this can be written as
(28)ρ(p1∣𝒮)=1ℰ∑a‍∑i‍δ(p1-Pia),
and for general r=1,2,…(29)ρ(p1,…,pr∣𝒮)=∑N=r∞‍ℛNρ(p1,…,pr∣𝒮N)=1ℰ∑a‍∑i1≠⋯≠ir‍∏d=1r‍δ(pd-Pida),
keeping in mind that ρ^(p1,…,pr∣𝒮N) will be zero whenever N(a)<r. The integral of any rth-order inclusive averaged density is the rth-order factorial moment of the multiplicity distribution,
(30)∫Ω‍dp1⋯dprρ(p1,…,pr∣𝒮)=∑N=r∞‍ℛNNr_=〈Nr_〉,
with simple angle brackets denoting inclusive averaging.The averaged counters are of course directly related to the traditional definitions in terms of cross-sections. Ifℒ is the integrated luminosity of incoming particles, the topological cross-section is σN=ℰN/ℒ, the inelastic cross-section is σI=ℰ/ℒ, and the inclusive cross-section is σincl=∑N‍NσN=〈N〉σI while the relative frequency (multiplicity distribution) can be written as usual as ℛN=σN/σI. The relation between the differential cross sections and our counters is
(31)ρincl(p1,…,pN)=ρ(p1,…,pN∣𝒮)=1σId3Nσincldp1⋯dpN,ρ(p1,…,pN∣𝒮N)=1σNd3Nσexcldp1⋯dpN,
and so as usual inclusive and exclusive densities are related by [14]
(32)ρ(p1,…,pN∣𝒮)=∑N=r∞‍ℛN(N-r)!∫‍ρ(p1,…,pN∣𝒮N)dpr+1⋯dpN,
while the semi-inclusive cross sections and counters follow by the usual projections.


## 2.4. Counters and Densities for Fixed(N,n)

Our choice of a basic counter is motivated by the experimental situation set out in Section1: we wish to work in event subsamples of fixed total charged multiplicity N(a) in the entire momentum space Ω but do the differential correlation analysis using only those pions n which fall into the restricted space 𝒜 and of a particular charge +1 or -1. This requires the use of “sub-subsamples” for which both N and n are kept fixed,
(33)𝒮nN={Pa,withaconstrainedbyδ(n,n+(a))δ(N,N(a))},
with n+(a) being the number of positive pions of event a in 𝒜, and eventwise sub-subsample counters
(34)ρ^(p1,…,pr∣n,N,Pa)=δ(n,n+(a))δ(N,N(a))×∑i1≠⋯≠ir=1n‍δ(p1-Pi1a)⋯δ(pr-Pira).
As in (2), the number of events in a sub-subsample ℰnN=∑a‍δ(n,n+(a))δ(N,N(a)) enters the relevant event averages
(35)ρ(p1,…,pr∣𝒮nN)=1ℰnN∑a‍ρ^(p1,…,pr∣n,N,Pa),
where once again the double Kronecker deltas in (34) ensure selection of events in 𝒮nN only. Integrals of the counters over the good-azimuth region 𝒜 yield, for the eventwise and 𝒮nN-averaged counters,
(36)∫𝒜‍dp1⋯dprρ^(p1,…,pr∣n,N,Pa)=nr_δ(n,n+(a))δ(N,N(a)),∫𝒜‍dp1⋯dprρ(p1,…,pr∣𝒮nN)=nr_.
Bearing in mind that observation points p1,p2,… refer to positive pions in 𝒜 only, the event-averaged counters for fixed N but any n are given by the average weighted in terms of the relative frequency ℛnN=ℰnN/ℰN,
(37)ρ(p1,…,pr∣𝒮N)=∑n=rN‍ℛnNρ(p1,…,pr∣𝒮nN)=〈ρ(p1,…,pr∣𝒮nN)〉N,
for r=1,2,3,…, which integrate to
(38)∫𝒜‍dp1⋯dprρ(p1,…,pr∣𝒮N)=〈nr_〉N.


## 3. Construction of Correlation Quantities

### 3.1. Criteria

Correlation measurements of any sort are only meaningful if a reference baseline signifying “independence” or “lack of correlation” is defined quantitatively; indeed, many different kinds of correlations may be defined and measured on the same data, depending on which particular physical and mathematical scenario is considered to be known or trivial and taken to be the baseline [61]. In our case, we require the reference distribution to have the following properties.(1)
The number of charged pions in all phase spaceNis an important parameter as a measure of possibly different physics, but only thenpositive pions in𝒜are to be considered in the differential analysis.(2)
For a given(N,n), the momenta of the reference densityρref(p1,…,pr∣𝒮nN)should be mutually independent for any order1≤r≤n. This and the previous requirement imply that the reference should be a n multinomial distributed over continuous momentum space; see Section 3.2.1.(3)
Given fixedN, the reference densityρref(p1,…,pr∣𝒮N)must reproduce then-multiplicity structure of the sub-subsamples𝒮nNas embodied inℛnN. As set out further in Section3.2.2, this translates into an average of multinomials,(39)ρref(p1,…,pr∣𝒮N)=∑n=rN‍ℛnNρmult(p1,…,pr∣α,𝒮nN)=〈ρmult(p1,…,pr∣α,𝒮nN)〉N.(4)
The reference density should reproduce the measured one-particle density in momentum space. This can in principle be satisfied by three different expressions for the multinomial's parameters α: see Section 3.2.3.(5)
Measures of correlation must reduce to zero even on a differential basis whenever the data is, in fact, uncorrelated. While this may seem self-evident, this requirement is often ignored or not satisfied in the literature. We address the resulting proper baseline through the use of internal cumulants in Section 3.3.(6)
The measure of correlation should be insensitive to the one-particle distribution. This is addressed as usual by normalisation; see Section 3.3.


### 3.2. The Reference Distribution

#### 3.2.1. Multinomials in Discrete and Continuous Spaces

Before (39) can be developed further, it is necessary to take a detour into discrete outcome spaces before tackling the continuous outcome space defined by p and Pa. The reason is that multinomial distributions for continuous arguments p can be written only as a limit of the discrete precursor.Let there be binsΩb,b=1,…,B with the corresponding set of Bernoulli probabilitiesα={α(b)}b=1B of a single particle falling into bin Ωb, normalised by ∑b‍α(b)=1. Independent tossing of n particles into these bins results in the multinomial for the bin counts n={nb}b=1B,
(40)p(n∣α,n)=n!∏b=1B‍α(b)nbnb!,
with normalisation
(41)∑U(n)‍p(n∣α,n)=1,
where the sum must be taken over the “universal set”
(42)U(n)={n∣nb≥0;∑b‍nb=n}.
The multivariate factorial moment generating function (FMGF) for this multinomial for the set of source parameters λ={λ(b)}b=1B can be solved in closed form,
(43)Qmult(λ∣α,n)=∑U‍p(n∣α,n)∏b‍(1-λ(b))nb=[1-∑b‍λ(b)α(b)]n.
The FMGF Q(λ) can generally be used to find multivariate factorial moments ρ(bi1,bi2,…,bir) and factorial cumulants κ(bi1,bi2,…,bir) for any selection of bins (bi1,bi2,…,bir)∈(1,…,B), including repeated indices, by differentiation
(44)ρ(bi1,bi2,…,bir)=(-1)r∂rQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0,κ(bi1,bi2,…,bir)=(-1)r∂rlnQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0.
For the multinomial case (43), the factorial moments and cumulants are therefore
(45)ρmult(bi1,bi2,…,bir∣α,n)=nr_∏d=1r‍α(id)∀r≤n,κmult(bi1,bi2,…,bir∣α,n)=(-1)r-1(r-1)!·n∏d=1r‍α(id).
The multinomial for variable p in continuous outcome space ℝ is derived by keeping n constant while taking the limit B→∞ with bin sizes tending to zero and changing to a Bernoulli probability density α(b)→dpα(p) normalised by ∫𝒜‍dpα(p)=1. The result is the point process where the probability for the count n(p) in the infinitesimal “bin” around any p to be larger than 1 becomes negligible; that is, we have at most one particle at a given p. While the multinomial probability itself can be written only as a limit, the FMGF can be written analytically as the functional [67]
(46)Qmult[λ(p)∣α(p),n]=[1-∫𝒜‍dpλ(p)α(p)]n.
Factorial moments and factorial cumulants are found generically from functional derivatives [14]
(47)ρ(pi1,pi2,…,pir)=(-1)rδrQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,κ(pi1,pi2,…,pir)=(-1)rδrlnQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,
which for the multinomial Q[λ(p)]=Qmult[λ(p)∣α(p),n] of (46) yield
(48)ρmult(pi1,pi2,…,pir∣α(p),n)=nr_∏k=1r‍α(pik),1≤r≤n,(49)=κmult(pi1,pi2,…,pir∣α(p),n)=(-1)r-1(r-1)!·n∏k=1r‍α(pik).


#### 3.2.2. Multinomial Reference for FixedN

Applying the above general case to our reference distribution (39), we must rewrite (46) to make provision for the fact that α may in general depend not only on N but also on n,
(50)Qmult[λ(p)∣α(p∣𝒮nN)]=[1-∫‍dpλ(p)α(p∣𝒮nN)]n.
Inserting (50) into (39), we find the FMGF for the reference distribution of subsample 𝒮N to be
(51)Qref[λ(p)∣α(p),𝒮N]=∑n‍ℛnNQmult[λ(p)∣α(p∣𝒮nN)]=〈[1-∫‍dpλ(p)α(p∣𝒮nN)]n〉N.
Using (48), the reference factorial moments are therefore
(52)ρref(p1,p2,…,pr∣𝒮N)=〈nr_∏k=1r‍α(pk∣𝒮nN)〉N,∀r≤n,
with corresponding expressions for the reference factorial cumulants.


#### 3.2.3. Reproducing the One-Particle Distribution

The set of functionsα(p∣𝒮nN) are as yet undetermined, apart from the general constraints α(p∣𝒮nN)≥0 and ∫𝒜‍dpα(p∣𝒮nN)=1. In multinomials of all kinds, the Bernoulli probabilities α are fixed parameters and therefore are the conveyers of whatever remains constant in the outcomes while the detailed outcomes fluctuate as statistical outcomes do. The “field” α(p∣𝒮nN) can and must therefore be seen as the quantity encompassing the “physics” of the one-particle distributions, which, in the absence of additional external information, is embodied by our experimental data sample: the experimental densities ρ(p1,…,pN∣𝒮N) “are” the physics, including all correlations, and their first-order projections ρ(p1∣𝒮N) “are” the one-particle physics. The question immediately arises whether α(p∣𝒮nN) should be fixed by ρ(p∣𝒮nN) or the n-average ρ(p∣𝒮N)=〈ρ(p∣𝒮nN)〉N. Three possible choices come to mind.(1)
It is tempting to define it in terms of the density for each sub-subsample𝒮nN,
(53)α(p∣𝒮nN)=ρ(p∣𝒮nN)n∀(N,n),which is correctly normalised since∫‍dpρ(p∣𝒮nN)=n. As this choice would attribute physical significance to n, it would be appropriate whenever n is associated with additionally measured experimental information. If, however, n fluctuates randomly from event to event based in part on unmeasured or unmeasurable properties such as an event's azimuthal orientation, use of (53) makes no sense.Ifn is deemed physically relevant, correlations in terms of ρ(p1,…,pr∣𝒮nN) of (35) may be feasible, conditional on the availability of a sufficient number of events ℰnN. Where sample sizes do not permit this, one could nevertheless attempt to measure what have historically been termed “short-range correlations” but in this case not in the traditional sense of fixed-N correlations versus inclusive ones but rather of fixed-n-fixed-N correlations versus fluctuating-n-fixed-N correlations. See Section 3.6.(2)
A second choice(54)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)n〉Nwould be properly normalised but fails to satisfy the crucial relations (71)–(73) below and is hence discarded.(3)
While remaining open-minded towards Choice 1, we therefore choose the third possibility, the ratio of the average density divided by the average, all for fixedN,
(55)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)〉N〈n〉N=ρ(p∣𝒮N)〈n〉N,which would be appropriate for samples whereℰnN is too small or physical significance can be attributed only to N but not to n. According to (38), it is also correctly normalised and ensures that the Bernoulli parameters are the same for all events in 𝒮N, independent of n. Substituting this into (52), the differential reference factorial moments orders become
(56)ρref(p1,p2,…,pr∣𝒮N)=〈nr_〉N〈n〉Nr∏d=1r‍ρ(pd∣𝒮N)=FrN∏d=1r‍ρ(pd∣𝒮N),where we identify the prefactor as the normalised factorial moments of then-multiplicity distribution for given N,
(57)FrN=〈nr_〉N〈n〉Nr,while the generating functional (51) becomes (see also [13])
(58)Qref[λ(p)∣α(p),𝒮N]=〈[1-∫‍dpλ(p)ρ(p∣𝒮N)〈n〉N]n〉N.Taking functional derivatives of the logarithm of (58), the first-, second-, and third-order cumulants of the reference density are
(59)κref(p1∣𝒮N)=ρref(p1∣𝒮N)=ρ(p1∣𝒮N),κref(p1,p2∣𝒮N)=ρref(p1,p2∣𝒮N)-ρ(p1∣𝒮N)ρ(p1∣𝒮N)(60)=(〈n2_〉N〈n〉N2-1)ρ(p1∣𝒮N)×ρ(p2∣𝒮N),(61)κref(p1,p2,p3∣𝒮N)=ρref(p1,p2,p3∣𝒮N)-[3]ρref(p1,p2∣𝒮N)ρ(p3∣𝒮N)+2ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N)(62)=(〈n3_〉N〈n〉N3-3〈n2_〉N〈n〉N2+2)×ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
where the square bracket [3] indicates the number of distinct permutations which must be taken into account. The terms in the rounded brackets are readily recognised as the normalised factorial cumulants of the n distribution for a given fixed N(63)KrN=1〈n〉Nr(-∂)r∂Λrln(∑n=0N‍ℛnN(1-Λ)n)Λ=0,
and so generalisation to arbitrary orders is immediate,
(64)κref(p1,…,pr∣𝒮N)=KrN∏k=1r‍ρ(pk∣𝒮N).
This can be proven generally by defining the functional Λ[λ(p)]=∫‍dpλ(p)ρ(p∣𝒮N)/〈n〉N which has only a first nonzero functional derivative δΛ/δλ(p1)=ρ(p∣𝒮N)/〈n〉N and the multiplicity generating function 𝒵(Λ∣𝒮N)=∑n‍ℛnN(1-Λ)n, in terms of which Qref[λ]=𝒵[Λ[λ]].


### 3.3. Internal Cumulants for Fixed𝒮N

Equation (64) shows that the differential cumulants of the reference distribution are directly proportional to the integrated cumulants KrN of n, which are zero only if ℛnN is Poissonian. For fixed N, neither the integrated cumulants KrN nor the differential ones are zero. While this has long been recognised in the literature [15], the inevitable consequence was not drawn; namely, that “Poissonian” cumulants for fixed N(65)κ(p1,p2∣𝒮N)=ρ(p1,p2∣𝒮N)-ρ(p1∣𝒮N)ρ(p2∣𝒮N),(66)κ(p1,p2,p3∣𝒮N)=ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)+2ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
and so forth cannot possibly represent true correlations because they are nonzero even when the momenta are fully independent. It is known that the theory of cumulants needs improvement on a fundamental level which reaches well beyond the scope of this paper  [68, 69], but those difficulties are irrelevant here. A first step which does address the above concerns was taken in [18], where it was shown very generally on the basis of generating functionals that correlations for samples of fixed N are best measured using the internal cumulantsκI, which are defined as the differences between the measured and the reference cumulants of the same order
(67)κI(p1,…,pr∣𝒮N)=κ(p1,…,pr∣𝒮N)-κref(p1,…,pr∣𝒮N).
For our averaged-multinomial reference case, the internal cumulants of second and third-orders are given by the differences between (65) and (60) and between (66) and (62), resulting in
(68)κI(p1,p2∣𝒮N)=ρ(p1,p2∣𝒮N)-F2Nρ(p1∣𝒮N)ρ(p1∣𝒮N),(69)κI(p1,p2,p3∣𝒮N)=ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)+G3Nρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
with
(70)G3N=3F2N-F3N=3〈n(n-1)〉N〈n〉N2-〈n(n-1)(n-2)〉N〈n〉N3,
and so on for higher orders. These internal cumulants are identically zero if and when the measured densities for fixed 𝒮nN are multinomials since then from (56)
(71)ρ(p1,p2∣𝒮N)⟶ρref(p1,p2∣𝒮N)=F2Nρ(p1∣𝒮N)ρ(p2∣𝒮N),
so that κI(p1,p2∣𝒮N)→0 whenever the data is multinomial, while
(72)ρ(p1,p2,p3∣𝒮N)⟶F3Nρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N)
ensures that κI(p1,p2,p3∣𝒮N)→0 in the same case. On another level, the internal cumulants always integrate to zero over the full good-azimuth space 𝒜, irrespective of the presence of correlations,
(73)∫𝒜‍dp1dp2κI(p1,p2∣𝒮N)=∫𝒜‍dp1dp2dp3κI(p1,p2,p3∣𝒮N)=0,
and so on for all orders. Both properties will remain valid after transformation from three-momentum to invariant four-momentum differences in Section 3.4. In the case of Poissonian statistics, FrN=1 for all r, so that the above internal cumulants revert to their usual definitions.As stated in Section3.1, the measured correlations may in addition be made insensitive to the one-particle distribution through normalisation. As set out in [18], such normalisation is achieved for fixed N by dividing the internal cumulants by the corresponding reference distribution density, which for the case at hand is given by (56). This leads to the second-order normalised internal cumulant
(74)KI(p1,p2∣𝒮N)=1F2Nρ(p1,p2∣𝒮N)ρ(p1∣𝒮N)ρ(p2∣𝒮N)-1,
while in third-order we get
(75)KI(p1,p2,p3∣𝒮N)=1F3N×(ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N))+3F2NF3N-1.


### 3.4. Correlation Integrals for Momentum Differences

In femtoscopy, correlations are mostly expressed in terms of pair variablesK=(1/2)(p1+p2) and difference q=p1-p2 or the invariant four-momenta [70] Q=-(p1-p2)2=(p1-p2)2-(E1-E2)2 where the energies are on-shell, Er=pr2+m2. As shown in [53], the formulation of eventwise counters as sums and products of Dirac delta functions makes it easy to change variables. Writing ρrN(p1,…,pr) as shorthand for ρ(p1,…,pr∣𝒮N) and so forth, the second-order unnormalised internal cumulant in terms of Q is, for example, found from the identity ∫‍dQκ2NI(Q)=∫‍dQ∫𝒜‍dp1dp2κ2NI(p1,p2)δ(Q-(p1-p2)2-(E1-E2)2) to be
(76)κ2NI(Q)=〈∑i≠j‍δ(Q-Qijaa)〉Na-F2N〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na=ρ2N(Q)-F2Nρ1N⊗ρ1N(Q),
where the counters in the second line are defined by the terms in the first, while Qijaa=[-(Pia-Pja)2]1/2 and Qijab=[-(Pia-Pjb)2]1/2 are four-momentum differences between sibling pairs aa and event mixing pairs ab, respectively. It is easy to show that ∫0∞‍dQρ2N(Q)=〈n2_〉N and ∫0∞‍dQρ1N⊗ρ1N(Q)=〈n〉N2 and hence, as before, ∫0∞‍dQκ2NI(Q)=0, which will be true for any correlation whatsoever. The double event averages in the product term
(77)ρ1N⊗ρ1N(Q)=〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na
are the theoretical foundations of event mixing [53]; the inner b-average is usually shortened to a smaller “moving average tail” subsample of 𝒮N.In third-order, the “GHP average” invariant is defined as the average of three two-momentum differences over all pairs (with or without the3),
(78)Qa=-(p1-p2)2-(p2-p3)2-(p3-p1)23,
and it is related to the invariant mass of three pions M3=(p1+p2+p3)2 by Qa2=(1/3)M32-m2. Other “topologies” such as the “GHP max” Qm=max[-(p1-p2)2,-(p2-p3)2,-(p3-p1)2] can also be employed. For large multiplicities, the “Star” topology may be preferred [71], but we will not pursue it here. For the GHP average, the third internal cumulant is given by
(79)κ3NI(Qa)=〈∑i≠j≠k‍δ(Qa-Qikjaaa)〉Na-3〈〈∑i≠j‍∑k‍δ(Qa-Qijkaab)〉Nb〉Na+G3N〈〈〈∑i,j,k‍δ(Qa-Qijkabc)〉Nc〉Nb〉Na,
with Qijkabc=(1/3)[-(Pia-Pjb)2-(Pjb-Pkc)2-(Pkc-Pia)2] and similarly for Qikjaaa and Qikjaab. Second- and third-order cumulants are normalised by, respectively,
(80)F2Nρ1N⊗ρ1N(Q)=F2N〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na,F3Nρ1N⊗ρ1N⊗ρ1N(Qa)=F3N〈〈〈∑i,j,k‍δ(Qa-Qijkabc)〉Nc〉Nb〉Na.
After transforming from momenta to Q, the formulae of Section 3.3 become
(81)κ2NI(Q)=ρ2N(Q)-F2Nρ1N⊗ρ1N(Q),(82)κ3NI(Qa)=ρ3N(Qa)-[3]ρ2N⊗ρ1N(Qa)+G3Nρ1N⊗ρ1N⊗ρ1N(Qa),
while the normalised cumulants of Section 3.3 become
(83)K2I(Q∣𝒮N)=ρ2N(Q)F2Nρ1N⊗ρ1N(Q)-1,(84)K3I(Qa∣𝒮N)=ρ3N(Qa)-[3]ρ2N⊗ρ1N(Qa)F3Nρ1N⊗ρ1N⊗ρ1N(Qa)+3F2NF3N-1.


### 3.5. Effect of Fixed-N Correction Factors

To get a feeling for the size of the corrections involved, we measured the correction factorsFrN and G3N with the same UA1 dataset and the same cuts as in Figure 1. As shown in Figure 2, the consequence of the clearly sub-Poissonian multiplicity distributions shown in Figure 1 is that these factors are significantly less than 1, in contrast to the usual factorial moments of the charged multiplicity distribution which are super-Poissonian with factors exceeding 1. For very low multiplicities N<10, normalised internal cumulants are hence larger than their Poissonian counterparts but converge to them with increasing N. Nevertheless, up to N≃30 corrections of more than 5% for K2I and more than 20% for K3I compared to their uncorrected counterparts can be expected. By contrast, the additive correction G3N does not deviate much from the Poissonian limit of 2 except for very small N. By contrast, unnormalised internal cumulants (81) and (82) are far less sensitive to the multinomial correction.It is of interest to zoom in on the approach to the Poissonian limit of 1 and to compare these corrections to the equivalent charged-multiplicity-based ones, which for the case of fixedN, would be just N(N-1)/N2 and N(N-1)(N-2)/N3. In Figure 3, the Poissonian limit corresponds to zero on the y-axis. It is clear that the fixed-N factors go some way to correct for the fixed-N conditioning; the gap between them is approximately determined by 〈nr_〉N/Nr_, that is, by the exact definition and outcome space for n.


### 3.6. Eliminating Fluctuations inn

We return briefly to the first choice in Section3.2.3, that is, α(p∣𝒮nN)=ρ(p∣𝒮nN)/n which would permit different physics for each (N,n) combination. If we were willing and able to do analyses for each 𝒮nN, we would use the fixed-(N,n) equivalent of (68) and (69) derived from κref(p1,…,pr∣𝒮nN)=(-1)r-1(r-1)!n∏k=1r‍ρ(pk∣𝒮nN),
(85)κI(p1,p2∣𝒮nN)=ρ(p1,p2∣𝒮nN)-(1-1n)ρ(p1∣𝒮nN)ρ(p2∣𝒮nN),κI(p1,p2,p3∣𝒮nN)=ρ(p1,p2,p3∣𝒮nN)-[3]ρ(p1,p2∣𝒮nN)ρ(p3∣𝒮nN)+2(1-1n2)ρ(p1∣𝒮nN)×ρ(p2∣𝒮nN)ρ(p3∣𝒮nN),
and normalise by ρref(p1,…,pr∣𝒮nN)=(nr_/nr)∏k=1r‍ρ(pk|𝒮nN). Where that is not possible, we can still average over the above to form “Averaged Internal” (AI) correlations (see Section 5), but in this case averaging over n for fixed N,
(86)κSRC(p1,p2∣𝒮N)=〈κI(p1,p2∣𝒮nN)〉N,κSRC(p1,p2,p3∣𝒮N)=〈κI(p1,p2,p3∣𝒮N)〉N,
and normalise if necessary by the moment ρSRC(p1,…,pr∣𝒮N)=〈(nr_/nr)∏k=1r‍ρ(pk∣𝒮nN)〉N. Given that this involves products of moments in the sub-subsample 𝒮nN, event mixing would have to be restricted to the same sub-subsamples also, for example
(87)〈(1-1n)ρ(p1∣𝒮nN)ρ(p2∣𝒮nN)〉N=∑n=2N‍ℰnNℰN(1-1n)1(ℰnN)2×∑a,b∈𝒮nN‍ρ^(p1∣n,N,Pa)ρ^(p2∣n,N,Pb).
The transformation to pair variables works in the same way as in previous sections.


## 3.1. Criteria

Correlation measurements of any sort are only meaningful if a reference baseline signifying “independence” or “lack of correlation” is defined quantitatively; indeed, many different kinds of correlations may be defined and measured on the same data, depending on which particular physical and mathematical scenario is considered to be known or trivial and taken to be the baseline [61]. In our case, we require the reference distribution to have the following properties.(1)
The number of charged pions in all phase spaceNis an important parameter as a measure of possibly different physics, but only thenpositive pions in𝒜are to be considered in the differential analysis.(2)
For a given(N,n), the momenta of the reference densityρref(p1,…,pr∣𝒮nN)should be mutually independent for any order1≤r≤n. This and the previous requirement imply that the reference should be a n multinomial distributed over continuous momentum space; see Section 3.2.1.(3)
Given fixedN, the reference densityρref(p1,…,pr∣𝒮N)must reproduce then-multiplicity structure of the sub-subsamples𝒮nNas embodied inℛnN. As set out further in Section3.2.2, this translates into an average of multinomials,(39)ρref(p1,…,pr∣𝒮N)=∑n=rN‍ℛnNρmult(p1,…,pr∣α,𝒮nN)=〈ρmult(p1,…,pr∣α,𝒮nN)〉N.(4)
The reference density should reproduce the measured one-particle density in momentum space. This can in principle be satisfied by three different expressions for the multinomial's parameters α: see Section 3.2.3.(5)
Measures of correlation must reduce to zero even on a differential basis whenever the data is, in fact, uncorrelated. While this may seem self-evident, this requirement is often ignored or not satisfied in the literature. We address the resulting proper baseline through the use of internal cumulants in Section 3.3.(6)
The measure of correlation should be insensitive to the one-particle distribution. This is addressed as usual by normalisation; see Section 3.3.


## 3.2. The Reference Distribution

### 3.2.1. Multinomials in Discrete and Continuous Spaces

Before (39) can be developed further, it is necessary to take a detour into discrete outcome spaces before tackling the continuous outcome space defined by p and Pa. The reason is that multinomial distributions for continuous arguments p can be written only as a limit of the discrete precursor.Let there be binsΩb,b=1,…,B with the corresponding set of Bernoulli probabilitiesα={α(b)}b=1B of a single particle falling into bin Ωb, normalised by ∑b‍α(b)=1. Independent tossing of n particles into these bins results in the multinomial for the bin counts n={nb}b=1B,
(40)p(n∣α,n)=n!∏b=1B‍α(b)nbnb!,
with normalisation
(41)∑U(n)‍p(n∣α,n)=1,
where the sum must be taken over the “universal set”
(42)U(n)={n∣nb≥0;∑b‍nb=n}.
The multivariate factorial moment generating function (FMGF) for this multinomial for the set of source parameters λ={λ(b)}b=1B can be solved in closed form,
(43)Qmult(λ∣α,n)=∑U‍p(n∣α,n)∏b‍(1-λ(b))nb=[1-∑b‍λ(b)α(b)]n.
The FMGF Q(λ) can generally be used to find multivariate factorial moments ρ(bi1,bi2,…,bir) and factorial cumulants κ(bi1,bi2,…,bir) for any selection of bins (bi1,bi2,…,bir)∈(1,…,B), including repeated indices, by differentiation
(44)ρ(bi1,bi2,…,bir)=(-1)r∂rQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0,κ(bi1,bi2,…,bir)=(-1)r∂rlnQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0.
For the multinomial case (43), the factorial moments and cumulants are therefore
(45)ρmult(bi1,bi2,…,bir∣α,n)=nr_∏d=1r‍α(id)∀r≤n,κmult(bi1,bi2,…,bir∣α,n)=(-1)r-1(r-1)!·n∏d=1r‍α(id).
The multinomial for variable p in continuous outcome space ℝ is derived by keeping n constant while taking the limit B→∞ with bin sizes tending to zero and changing to a Bernoulli probability density α(b)→dpα(p) normalised by ∫𝒜‍dpα(p)=1. The result is the point process where the probability for the count n(p) in the infinitesimal “bin” around any p to be larger than 1 becomes negligible; that is, we have at most one particle at a given p. While the multinomial probability itself can be written only as a limit, the FMGF can be written analytically as the functional [67]
(46)Qmult[λ(p)∣α(p),n]=[1-∫𝒜‍dpλ(p)α(p)]n.
Factorial moments and factorial cumulants are found generically from functional derivatives [14]
(47)ρ(pi1,pi2,…,pir)=(-1)rδrQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,κ(pi1,pi2,…,pir)=(-1)rδrlnQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,
which for the multinomial Q[λ(p)]=Qmult[λ(p)∣α(p),n] of (46) yield
(48)ρmult(pi1,pi2,…,pir∣α(p),n)=nr_∏k=1r‍α(pik),1≤r≤n,(49)=κmult(pi1,pi2,…,pir∣α(p),n)=(-1)r-1(r-1)!·n∏k=1r‍α(pik).


### 3.2.2. Multinomial Reference for FixedN

Applying the above general case to our reference distribution (39), we must rewrite (46) to make provision for the fact that α may in general depend not only on N but also on n,
(50)Qmult[λ(p)∣α(p∣𝒮nN)]=[1-∫‍dpλ(p)α(p∣𝒮nN)]n.
Inserting (50) into (39), we find the FMGF for the reference distribution of subsample 𝒮N to be
(51)Qref[λ(p)∣α(p),𝒮N]=∑n‍ℛnNQmult[λ(p)∣α(p∣𝒮nN)]=〈[1-∫‍dpλ(p)α(p∣𝒮nN)]n〉N.
Using (48), the reference factorial moments are therefore
(52)ρref(p1,p2,…,pr∣𝒮N)=〈nr_∏k=1r‍α(pk∣𝒮nN)〉N,∀r≤n,
with corresponding expressions for the reference factorial cumulants.


### 3.2.3. Reproducing the One-Particle Distribution

The set of functionsα(p∣𝒮nN) are as yet undetermined, apart from the general constraints α(p∣𝒮nN)≥0 and ∫𝒜‍dpα(p∣𝒮nN)=1. In multinomials of all kinds, the Bernoulli probabilities α are fixed parameters and therefore are the conveyers of whatever remains constant in the outcomes while the detailed outcomes fluctuate as statistical outcomes do. The “field” α(p∣𝒮nN) can and must therefore be seen as the quantity encompassing the “physics” of the one-particle distributions, which, in the absence of additional external information, is embodied by our experimental data sample: the experimental densities ρ(p1,…,pN∣𝒮N) “are” the physics, including all correlations, and their first-order projections ρ(p1∣𝒮N) “are” the one-particle physics. The question immediately arises whether α(p∣𝒮nN) should be fixed by ρ(p∣𝒮nN) or the n-average ρ(p∣𝒮N)=〈ρ(p∣𝒮nN)〉N. Three possible choices come to mind.(1)
It is tempting to define it in terms of the density for each sub-subsample𝒮nN,
(53)α(p∣𝒮nN)=ρ(p∣𝒮nN)n∀(N,n),which is correctly normalised since∫‍dpρ(p∣𝒮nN)=n. As this choice would attribute physical significance to n, it would be appropriate whenever n is associated with additionally measured experimental information. If, however, n fluctuates randomly from event to event based in part on unmeasured or unmeasurable properties such as an event's azimuthal orientation, use of (53) makes no sense.Ifn is deemed physically relevant, correlations in terms of ρ(p1,…,pr∣𝒮nN) of (35) may be feasible, conditional on the availability of a sufficient number of events ℰnN. Where sample sizes do not permit this, one could nevertheless attempt to measure what have historically been termed “short-range correlations” but in this case not in the traditional sense of fixed-N correlations versus inclusive ones but rather of fixed-n-fixed-N correlations versus fluctuating-n-fixed-N correlations. See Section 3.6.(2)
A second choice(54)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)n〉Nwould be properly normalised but fails to satisfy the crucial relations (71)–(73) below and is hence discarded.(3)
While remaining open-minded towards Choice 1, we therefore choose the third possibility, the ratio of the average density divided by the average, all for fixedN,
(55)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)〉N〈n〉N=ρ(p∣𝒮N)〈n〉N,which would be appropriate for samples whereℰnN is too small or physical significance can be attributed only to N but not to n. According to (38), it is also correctly normalised and ensures that the Bernoulli parameters are the same for all events in 𝒮N, independent of n. Substituting this into (52), the differential reference factorial moments orders become
(56)ρref(p1,p2,…,pr∣𝒮N)=〈nr_〉N〈n〉Nr∏d=1r‍ρ(pd∣𝒮N)=FrN∏d=1r‍ρ(pd∣𝒮N),where we identify the prefactor as the normalised factorial moments of then-multiplicity distribution for given N,
(57)FrN=〈nr_〉N〈n〉Nr,while the generating functional (51) becomes (see also [13])
(58)Qref[λ(p)∣α(p),𝒮N]=〈[1-∫‍dpλ(p)ρ(p∣𝒮N)〈n〉N]n〉N.Taking functional derivatives of the logarithm of (58), the first-, second-, and third-order cumulants of the reference density are
(59)κref(p1∣𝒮N)=ρref(p1∣𝒮N)=ρ(p1∣𝒮N),κref(p1,p2∣𝒮N)=ρref(p1,p2∣𝒮N)-ρ(p1∣𝒮N)ρ(p1∣𝒮N)(60)=(〈n2_〉N〈n〉N2-1)ρ(p1∣𝒮N)×ρ(p2∣𝒮N),(61)κref(p1,p2,p3∣𝒮N)=ρref(p1,p2,p3∣𝒮N)-[3]ρref(p1,p2∣𝒮N)ρ(p3∣𝒮N)+2ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N)(62)=(〈n3_〉N〈n〉N3-3〈n2_〉N〈n〉N2+2)×ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
where the square bracket [3] indicates the number of distinct permutations which must be taken into account. The terms in the rounded brackets are readily recognised as the normalised factorial cumulants of the n distribution for a given fixed N(63)KrN=1〈n〉Nr(-∂)r∂Λrln(∑n=0N‍ℛnN(1-Λ)n)Λ=0,
and so generalisation to arbitrary orders is immediate,
(64)κref(p1,…,pr∣𝒮N)=KrN∏k=1r‍ρ(pk∣𝒮N).
This can be proven generally by defining the functional Λ[λ(p)]=∫‍dpλ(p)ρ(p∣𝒮N)/〈n〉N which has only a first nonzero functional derivative δΛ/δλ(p1)=ρ(p∣𝒮N)/〈n〉N and the multiplicity generating function 𝒵(Λ∣𝒮N)=∑n‍ℛnN(1-Λ)n, in terms of which Qref[λ]=𝒵[Λ[λ]].


## 3.2.1. Multinomials in Discrete and Continuous Spaces

Before (39) can be developed further, it is necessary to take a detour into discrete outcome spaces before tackling the continuous outcome space defined by p and Pa. The reason is that multinomial distributions for continuous arguments p can be written only as a limit of the discrete precursor.Let there be binsΩb,b=1,…,B with the corresponding set of Bernoulli probabilitiesα={α(b)}b=1B of a single particle falling into bin Ωb, normalised by ∑b‍α(b)=1. Independent tossing of n particles into these bins results in the multinomial for the bin counts n={nb}b=1B,
(40)p(n∣α,n)=n!∏b=1B‍α(b)nbnb!,
with normalisation
(41)∑U(n)‍p(n∣α,n)=1,
where the sum must be taken over the “universal set”
(42)U(n)={n∣nb≥0;∑b‍nb=n}.
The multivariate factorial moment generating function (FMGF) for this multinomial for the set of source parameters λ={λ(b)}b=1B can be solved in closed form,
(43)Qmult(λ∣α,n)=∑U‍p(n∣α,n)∏b‍(1-λ(b))nb=[1-∑b‍λ(b)α(b)]n.
The FMGF Q(λ) can generally be used to find multivariate factorial moments ρ(bi1,bi2,…,bir) and factorial cumulants κ(bi1,bi2,…,bir) for any selection of bins (bi1,bi2,…,bir)∈(1,…,B), including repeated indices, by differentiation
(44)ρ(bi1,bi2,…,bir)=(-1)r∂rQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0,κ(bi1,bi2,…,bir)=(-1)r∂rlnQ(λ)∂λ(bi1)∂λ(bi2)⋯∂λ(bir)|λ=0.
For the multinomial case (43), the factorial moments and cumulants are therefore
(45)ρmult(bi1,bi2,…,bir∣α,n)=nr_∏d=1r‍α(id)∀r≤n,κmult(bi1,bi2,…,bir∣α,n)=(-1)r-1(r-1)!·n∏d=1r‍α(id).
The multinomial for variable p in continuous outcome space ℝ is derived by keeping n constant while taking the limit B→∞ with bin sizes tending to zero and changing to a Bernoulli probability density α(b)→dpα(p) normalised by ∫𝒜‍dpα(p)=1. The result is the point process where the probability for the count n(p) in the infinitesimal “bin” around any p to be larger than 1 becomes negligible; that is, we have at most one particle at a given p. While the multinomial probability itself can be written only as a limit, the FMGF can be written analytically as the functional [67]
(46)Qmult[λ(p)∣α(p),n]=[1-∫𝒜‍dpλ(p)α(p)]n.
Factorial moments and factorial cumulants are found generically from functional derivatives [14]
(47)ρ(pi1,pi2,…,pir)=(-1)rδrQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,κ(pi1,pi2,…,pir)=(-1)rδrlnQ[λ(p)]δλ(pi1)δλ(pi2)⋯δλ(pir)|λ(p)=0,
which for the multinomial Q[λ(p)]=Qmult[λ(p)∣α(p),n] of (46) yield
(48)ρmult(pi1,pi2,…,pir∣α(p),n)=nr_∏k=1r‍α(pik),1≤r≤n,(49)=κmult(pi1,pi2,…,pir∣α(p),n)=(-1)r-1(r-1)!·n∏k=1r‍α(pik).


## 3.2.2. Multinomial Reference for FixedN

Applying the above general case to our reference distribution (39), we must rewrite (46) to make provision for the fact that α may in general depend not only on N but also on n,
(50)Qmult[λ(p)∣α(p∣𝒮nN)]=[1-∫‍dpλ(p)α(p∣𝒮nN)]n.
Inserting (50) into (39), we find the FMGF for the reference distribution of subsample 𝒮N to be
(51)Qref[λ(p)∣α(p),𝒮N]=∑n‍ℛnNQmult[λ(p)∣α(p∣𝒮nN)]=〈[1-∫‍dpλ(p)α(p∣𝒮nN)]n〉N.
Using (48), the reference factorial moments are therefore
(52)ρref(p1,p2,…,pr∣𝒮N)=〈nr_∏k=1r‍α(pk∣𝒮nN)〉N,∀r≤n,
with corresponding expressions for the reference factorial cumulants.


## 3.2.3. Reproducing the One-Particle Distribution

The set of functionsα(p∣𝒮nN) are as yet undetermined, apart from the general constraints α(p∣𝒮nN)≥0 and ∫𝒜‍dpα(p∣𝒮nN)=1. In multinomials of all kinds, the Bernoulli probabilities α are fixed parameters and therefore are the conveyers of whatever remains constant in the outcomes while the detailed outcomes fluctuate as statistical outcomes do. The “field” α(p∣𝒮nN) can and must therefore be seen as the quantity encompassing the “physics” of the one-particle distributions, which, in the absence of additional external information, is embodied by our experimental data sample: the experimental densities ρ(p1,…,pN∣𝒮N) “are” the physics, including all correlations, and their first-order projections ρ(p1∣𝒮N) “are” the one-particle physics. The question immediately arises whether α(p∣𝒮nN) should be fixed by ρ(p∣𝒮nN) or the n-average ρ(p∣𝒮N)=〈ρ(p∣𝒮nN)〉N. Three possible choices come to mind.(1)
It is tempting to define it in terms of the density for each sub-subsample𝒮nN,
(53)α(p∣𝒮nN)=ρ(p∣𝒮nN)n∀(N,n),which is correctly normalised since∫‍dpρ(p∣𝒮nN)=n. As this choice would attribute physical significance to n, it would be appropriate whenever n is associated with additionally measured experimental information. If, however, n fluctuates randomly from event to event based in part on unmeasured or unmeasurable properties such as an event's azimuthal orientation, use of (53) makes no sense.Ifn is deemed physically relevant, correlations in terms of ρ(p1,…,pr∣𝒮nN) of (35) may be feasible, conditional on the availability of a sufficient number of events ℰnN. Where sample sizes do not permit this, one could nevertheless attempt to measure what have historically been termed “short-range correlations” but in this case not in the traditional sense of fixed-N correlations versus inclusive ones but rather of fixed-n-fixed-N correlations versus fluctuating-n-fixed-N correlations. See Section 3.6.(2)
A second choice(54)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)n〉Nwould be properly normalised but fails to satisfy the crucial relations (71)–(73) below and is hence discarded.(3)
While remaining open-minded towards Choice 1, we therefore choose the third possibility, the ratio of the average density divided by the average, all for fixedN,
(55)α(p∣𝒮nN)=〈ρ(p∣𝒮nN)〉N〈n〉N=ρ(p∣𝒮N)〈n〉N,which would be appropriate for samples whereℰnN is too small or physical significance can be attributed only to N but not to n. According to (38), it is also correctly normalised and ensures that the Bernoulli parameters are the same for all events in 𝒮N, independent of n. Substituting this into (52), the differential reference factorial moments orders become
(56)ρref(p1,p2,…,pr∣𝒮N)=〈nr_〉N〈n〉Nr∏d=1r‍ρ(pd∣𝒮N)=FrN∏d=1r‍ρ(pd∣𝒮N),where we identify the prefactor as the normalised factorial moments of then-multiplicity distribution for given N,
(57)FrN=〈nr_〉N〈n〉Nr,while the generating functional (51) becomes (see also [13])
(58)Qref[λ(p)∣α(p),𝒮N]=〈[1-∫‍dpλ(p)ρ(p∣𝒮N)〈n〉N]n〉N.Taking functional derivatives of the logarithm of (58), the first-, second-, and third-order cumulants of the reference density are
(59)κref(p1∣𝒮N)=ρref(p1∣𝒮N)=ρ(p1∣𝒮N),κref(p1,p2∣𝒮N)=ρref(p1,p2∣𝒮N)-ρ(p1∣𝒮N)ρ(p1∣𝒮N)(60)=(〈n2_〉N〈n〉N2-1)ρ(p1∣𝒮N)×ρ(p2∣𝒮N),(61)κref(p1,p2,p3∣𝒮N)=ρref(p1,p2,p3∣𝒮N)-[3]ρref(p1,p2∣𝒮N)ρ(p3∣𝒮N)+2ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N)(62)=(〈n3_〉N〈n〉N3-3〈n2_〉N〈n〉N2+2)×ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
where the square bracket [3] indicates the number of distinct permutations which must be taken into account. The terms in the rounded brackets are readily recognised as the normalised factorial cumulants of the n distribution for a given fixed N(63)KrN=1〈n〉Nr(-∂)r∂Λrln(∑n=0N‍ℛnN(1-Λ)n)Λ=0,
and so generalisation to arbitrary orders is immediate,
(64)κref(p1,…,pr∣𝒮N)=KrN∏k=1r‍ρ(pk∣𝒮N).
This can be proven generally by defining the functional Λ[λ(p)]=∫‍dpλ(p)ρ(p∣𝒮N)/〈n〉N which has only a first nonzero functional derivative δΛ/δλ(p1)=ρ(p∣𝒮N)/〈n〉N and the multiplicity generating function 𝒵(Λ∣𝒮N)=∑n‍ℛnN(1-Λ)n, in terms of which Qref[λ]=𝒵[Λ[λ]].


## 3.3. Internal Cumulants for Fixed𝒮N

Equation (64) shows that the differential cumulants of the reference distribution are directly proportional to the integrated cumulants KrN of n, which are zero only if ℛnN is Poissonian. For fixed N, neither the integrated cumulants KrN nor the differential ones are zero. While this has long been recognised in the literature [15], the inevitable consequence was not drawn; namely, that “Poissonian” cumulants for fixed N(65)κ(p1,p2∣𝒮N)=ρ(p1,p2∣𝒮N)-ρ(p1∣𝒮N)ρ(p2∣𝒮N),(66)κ(p1,p2,p3∣𝒮N)=ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)+2ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
and so forth cannot possibly represent true correlations because they are nonzero even when the momenta are fully independent. It is known that the theory of cumulants needs improvement on a fundamental level which reaches well beyond the scope of this paper  [68, 69], but those difficulties are irrelevant here. A first step which does address the above concerns was taken in [18], where it was shown very generally on the basis of generating functionals that correlations for samples of fixed N are best measured using the internal cumulantsκI, which are defined as the differences between the measured and the reference cumulants of the same order
(67)κI(p1,…,pr∣𝒮N)=κ(p1,…,pr∣𝒮N)-κref(p1,…,pr∣𝒮N).
For our averaged-multinomial reference case, the internal cumulants of second and third-orders are given by the differences between (65) and (60) and between (66) and (62), resulting in
(68)κI(p1,p2∣𝒮N)=ρ(p1,p2∣𝒮N)-F2Nρ(p1∣𝒮N)ρ(p1∣𝒮N),(69)κI(p1,p2,p3∣𝒮N)=ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)+G3Nρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N),
with
(70)G3N=3F2N-F3N=3〈n(n-1)〉N〈n〉N2-〈n(n-1)(n-2)〉N〈n〉N3,
and so on for higher orders. These internal cumulants are identically zero if and when the measured densities for fixed 𝒮nN are multinomials since then from (56)
(71)ρ(p1,p2∣𝒮N)⟶ρref(p1,p2∣𝒮N)=F2Nρ(p1∣𝒮N)ρ(p2∣𝒮N),
so that κI(p1,p2∣𝒮N)→0 whenever the data is multinomial, while
(72)ρ(p1,p2,p3∣𝒮N)⟶F3Nρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N)
ensures that κI(p1,p2,p3∣𝒮N)→0 in the same case. On another level, the internal cumulants always integrate to zero over the full good-azimuth space 𝒜, irrespective of the presence of correlations,
(73)∫𝒜‍dp1dp2κI(p1,p2∣𝒮N)=∫𝒜‍dp1dp2dp3κI(p1,p2,p3∣𝒮N)=0,
and so on for all orders. Both properties will remain valid after transformation from three-momentum to invariant four-momentum differences in Section 3.4. In the case of Poissonian statistics, FrN=1 for all r, so that the above internal cumulants revert to their usual definitions.As stated in Section3.1, the measured correlations may in addition be made insensitive to the one-particle distribution through normalisation. As set out in [18], such normalisation is achieved for fixed N by dividing the internal cumulants by the corresponding reference distribution density, which for the case at hand is given by (56). This leads to the second-order normalised internal cumulant
(74)KI(p1,p2∣𝒮N)=1F2Nρ(p1,p2∣𝒮N)ρ(p1∣𝒮N)ρ(p2∣𝒮N)-1,
while in third-order we get
(75)KI(p1,p2,p3∣𝒮N)=1F3N×(ρ(p1,p2,p3∣𝒮N)-[3]ρ(p1,p2∣𝒮N)ρ(p3∣𝒮N)ρ(p1∣𝒮N)ρ(p2∣𝒮N)ρ(p3∣𝒮N))+3F2NF3N-1.


## 3.4. Correlation Integrals for Momentum Differences

In femtoscopy, correlations are mostly expressed in terms of pair variablesK=(1/2)(p1+p2) and difference q=p1-p2 or the invariant four-momenta [70] Q=-(p1-p2)2=(p1-p2)2-(E1-E2)2 where the energies are on-shell, Er=pr2+m2. As shown in [53], the formulation of eventwise counters as sums and products of Dirac delta functions makes it easy to change variables. Writing ρrN(p1,…,pr) as shorthand for ρ(p1,…,pr∣𝒮N) and so forth, the second-order unnormalised internal cumulant in terms of Q is, for example, found from the identity ∫‍dQκ2NI(Q)=∫‍dQ∫𝒜‍dp1dp2κ2NI(p1,p2)δ(Q-(p1-p2)2-(E1-E2)2) to be
(76)κ2NI(Q)=〈∑i≠j‍δ(Q-Qijaa)〉Na-F2N〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na=ρ2N(Q)-F2Nρ1N⊗ρ1N(Q),
where the counters in the second line are defined by the terms in the first, while Qijaa=[-(Pia-Pja)2]1/2 and Qijab=[-(Pia-Pjb)2]1/2 are four-momentum differences between sibling pairs aa and event mixing pairs ab, respectively. It is easy to show that ∫0∞‍dQρ2N(Q)=〈n2_〉N and ∫0∞‍dQρ1N⊗ρ1N(Q)=〈n〉N2 and hence, as before, ∫0∞‍dQκ2NI(Q)=0, which will be true for any correlation whatsoever. The double event averages in the product term
(77)ρ1N⊗ρ1N(Q)=〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na
are the theoretical foundations of event mixing [53]; the inner b-average is usually shortened to a smaller “moving average tail” subsample of 𝒮N.In third-order, the “GHP average” invariant is defined as the average of three two-momentum differences over all pairs (with or without the3),
(78)Qa=-(p1-p2)2-(p2-p3)2-(p3-p1)23,
and it is related to the invariant mass of three pions M3=(p1+p2+p3)2 by Qa2=(1/3)M32-m2. Other “topologies” such as the “GHP max” Qm=max[-(p1-p2)2,-(p2-p3)2,-(p3-p1)2] can also be employed. For large multiplicities, the “Star” topology may be preferred [71], but we will not pursue it here. For the GHP average, the third internal cumulant is given by
(79)κ3NI(Qa)=〈∑i≠j≠k‍δ(Qa-Qikjaaa)〉Na-3〈〈∑i≠j‍∑k‍δ(Qa-Qijkaab)〉Nb〉Na+G3N〈〈〈∑i,j,k‍δ(Qa-Qijkabc)〉Nc〉Nb〉Na,
with Qijkabc=(1/3)[-(Pia-Pjb)2-(Pjb-Pkc)2-(Pkc-Pia)2] and similarly for Qikjaaa and Qikjaab. Second- and third-order cumulants are normalised by, respectively,
(80)F2Nρ1N⊗ρ1N(Q)=F2N〈〈∑i,j‍δ(Q-Qijab)〉Nb〉Na,F3Nρ1N⊗ρ1N⊗ρ1N(Qa)=F3N〈〈〈∑i,j,k‍δ(Qa-Qijkabc)〉Nc〉Nb〉Na.
After transforming from momenta to Q, the formulae of Section 3.3 become
(81)κ2NI(Q)=ρ2N(Q)-F2Nρ1N⊗ρ1N(Q),(82)κ3NI(Qa)=ρ3N(Qa)-[3]ρ2N⊗ρ1N(Qa)+G3Nρ1N⊗ρ1N⊗ρ1N(Qa),
while the normalised cumulants of Section 3.3 become
(83)K2I(Q∣𝒮N)=ρ2N(Q)F2Nρ1N⊗ρ1N(Q)-1,(84)K3I(Qa∣𝒮N)=ρ3N(Qa)-[3]ρ2N⊗ρ1N(Qa)F3Nρ1N⊗ρ1N⊗ρ1N(Qa)+3F2NF3N-1.


## 3.5. Effect of Fixed-N Correction Factors

To get a feeling for the size of the corrections involved, we measured the correction factorsFrN and G3N with the same UA1 dataset and the same cuts as in Figure 1. As shown in Figure 2, the consequence of the clearly sub-Poissonian multiplicity distributions shown in Figure 1 is that these factors are significantly less than 1, in contrast to the usual factorial moments of the charged multiplicity distribution which are super-Poissonian with factors exceeding 1. For very low multiplicities N<10, normalised internal cumulants are hence larger than their Poissonian counterparts but converge to them with increasing N. Nevertheless, up to N≃30 corrections of more than 5% for K2I and more than 20% for K3I compared to their uncorrected counterparts can be expected. By contrast, the additive correction G3N does not deviate much from the Poissonian limit of 2 except for very small N. By contrast, unnormalised internal cumulants (81) and (82) are far less sensitive to the multinomial correction.It is of interest to zoom in on the approach to the Poissonian limit of 1 and to compare these corrections to the equivalent charged-multiplicity-based ones, which for the case of fixedN, would be just N(N-1)/N2 and N(N-1)(N-2)/N3. In Figure 3, the Poissonian limit corresponds to zero on the y-axis. It is clear that the fixed-N factors go some way to correct for the fixed-N conditioning; the gap between them is approximately determined by 〈nr_〉N/Nr_, that is, by the exact definition and outcome space for n.


## 3.6. Eliminating Fluctuations inn

We return briefly to the first choice in Section3.2.3, that is, α(p∣𝒮nN)=ρ(p∣𝒮nN)/n which would permit different physics for each (N,n) combination. If we were willing and able to do analyses for each 𝒮nN, we would use the fixed-(N,n) equivalent of (68) and (69) derived from κref(p1,…,pr∣𝒮nN)=(-1)r-1(r-1)!n∏k=1r‍ρ(pk∣𝒮nN),
(85)κI(p1,p2∣𝒮nN)=ρ(p1,p2∣𝒮nN)-(1-1n)ρ(p1∣𝒮nN)ρ(p2∣𝒮nN),κI(p1,p2,p3∣𝒮nN)=ρ(p1,p2,p3∣𝒮nN)-[3]ρ(p1,p2∣𝒮nN)ρ(p3∣𝒮nN)+2(1-1n2)ρ(p1∣𝒮nN)×ρ(p2∣𝒮nN)ρ(p3∣𝒮nN),
and normalise by ρref(p1,…,pr∣𝒮nN)=(nr_/nr)∏k=1r‍ρ(pk|𝒮nN). Where that is not possible, we can still average over the above to form “Averaged Internal” (AI) correlations (see Section 5), but in this case averaging over n for fixed N,
(86)κSRC(p1,p2∣𝒮N)=〈κI(p1,p2∣𝒮nN)〉N,κSRC(p1,p2,p3∣𝒮N)=〈κI(p1,p2,p3∣𝒮N)〉N,
and normalise if necessary by the moment ρSRC(p1,…,pr∣𝒮N)=〈(nr_/nr)∏k=1r‍ρ(pk∣𝒮nN)〉N. Given that this involves products of moments in the sub-subsample 𝒮nN, event mixing would have to be restricted to the same sub-subsamples also, for example
(87)〈(1-1n)ρ(p1∣𝒮nN)ρ(p2∣𝒮nN)〉N=∑n=2N‍ℰnNℰN(1-1n)1(ℰnN)2×∑a,b∈𝒮nN‍ρ^(p1∣n,N,Pa)ρ^(p2∣n,N,Pb).
The transformation to pair variables works in the same way as in previous sections.


## 4. Statistical Errors

While the various versions of internal cumulants constructed above may all be relevant at some point, we concentrate on finding expressions for experimental standard errors for the unnormalised and normalised internal cumulants of (81)–(84). This turns out to be more subtle than merely applying a generic root-mean-square prescription. We will show in this section that standard errors implemented thus far may have been underestimated even in the standard two-particle case.The calculations performed in this section belong to the “frequentist” view of probability; a proper Bayesian analysis, which can be expected to rest on more solid foundations, is beyond the scope of this paper. The two viewpoints can reasonably be expected to yield similar results in the limit of large bin contents and sample sizes.In this section, we often simplify notation by writingρr(Q∣𝒮N)→ρrN and ρ1⊗ρ1(Q∣𝒮N)→ρ1N2 and so forth, since the formulae apply to samples and variables of any kind.


### 4.1. Propagation of Statistical Errors

Because cumulants can be measured only through the moments that enter their definitions, the first task is to identify which moment variances and covariances are needed. By means of standard error propagation [68], we find the sample variances for second-order internal cumulants of (81) and (83) to be
(88)var(κ2I(Q∣𝒮N))=var(ρ2N-F2Nρ1N2)=var(ρ2N)+F2N2var(ρ1N2)-2F2Ncov(ρ2N,ρ1N2),(89)var(K2I(Q∣𝒮N))=var(1F2Nρ2Nρ1N2-1)(90)=[1F2Nρ2Nρ1N2]2×[var(ρ2N)ρ2N2+var(ρ1N2)ρ1N4-2cov(ρ2N,ρ1N2)ρ2N·ρ1N2],
under the assumption that var(F2N) is much smaller than the other variances, so that F2N can be treated as a constant; this is the case if there are many bins for Q. Similarly, from (69), the variance of the unnormalised internal cumulant is, assuming G3 of (70) to be constant,
(91)var(κ3I(Qa∣𝒮N))=var[ρ3N-3ρ2Nρ1N+G3ρ1N3](92)=var(ρ3N)+9var(ρ2Nρ1N)+G32var(ρ1N3)-6cov(ρ3N,ρ2Nρ1N)+2G3cov(ρ3N,ρ1N3)-6G3cov(ρ2Nρ1N,ρ1N3),
while the normalised version has variance (again assuming var(G3)≪varρr)
(93)var(K3I(Qa∣𝒮N))=var[ρ3N-3ρ2Nρ1N+G3ρ1N3F3Nρ1N3]=1F3N2var[ρ3N-3ρ2Nρ1Nρ1N3+G3](94)=[ρ3N-3ρ2Nρ1NF3Nρ1N3]2×[var(ρ3N-3ρ2Nρ1N)(ρ3N-3ρ2Nρ1N)2+var(ρ1N3)ρ1N6-2cov(ρ3N-3ρ2Nρ1N,ρ1N3)(ρ3N-3ρ2Nρ1N)ρ1N3]=[ρ3N-3ρ2Nρ1NF3Nρ1N3]2×[var(ρ3N)+9var(ρ2Nρ1N)-6cov(ρ3N,ρ2Nρ1N)(ρ3N-3ρ2Nρ1N)2+var(ρ1N3)ρ1N6+6cov(ρ2Nρ1N,ρ1N3)-2cov(ρ3N,ρ1N3)(ρ3N-3ρ2Nρ1N)ρ1N3].
The unnormalised cumulants (81) and (82) and their variances (88) and (92) require knowledge of the multiplicity factorial moments F2N,F3N, so that the individual terms must be accumulated until the entire sample has been analysed. By contrast, the normalised cumulants (83) and (84) and their variances (90) and (94) contain the multiplicity moments only as prefactors.


### 4.2. Expectation Values of Counters

While we will not make direct use of the results in this section, it is nevertheless useful briefly to consider what we might mean by an “expectation value of experimental counters and densities.” For any scalar functionf(p) of the momenta, the theoretical expectation value E[f] is defined as the integral over the entire outcome space Ω of f weighted by a “parent distribution” P(p), an abstract entity supposedly containing everything there is to know on this level,
(95)E[f(p)]=∫Ω‍dpP(p)f(p).
Purely theoretical concepts such as P(p) and E[f] should be given little or no room in a strongly experimentally-oriented study. In calculating standard errors on counters below, we will, however, make use of the exact factorisation that expectation values provide whenever two variables x,y are statistically independent, E[xy]=∫‍dxdyP(x,y)xy=E[x]E[y].Expectation values for pairwise variables such as the four-momentum differenceQ we are considering here must be based on the underlying physics. We can deduce some properties of the parent distribution based on the usual definition of the femtoscopic correlation function:
(96)C2(Q)≡K2(Q)+1=ρ2sibling(Q)ρ2reference(Q)=〈ρ^(Qaa)〉〈〈ρ^(Qab)〉〉,
which is a function of two entirely different quantities: the four-momentum differences Qaa of “sibling” tracks taken from the same event a and one constructed from the mixed-event sample using tracks from different events, written as Qab, Qbc, and so forth. For second-order correlations, the parent distribution is therefore necessarily a two-variable probability4  P(Qaa,Qbc) which, depending on whether the cases b=a and c=a occur, may or may not factorise into “sibling” and “mixed” marginal probabilities
(97)P(Qaa,Qbc)=Ps(Qaa)Pm(Qbc)iffa≠b≠c,
but (unless a=b=c), the marginals will always be
(98)Ps(Qaa)=∫‍dQbcP(Qaa,Qbc),Pm(Qbc)=∫‍dQaaP(Qaa,Qbc).
The shapes of Ps(Q) and Pm(Q) must necessarily be different since it is precisely this difference that leads to a nontrivial signal in (96). In terms of this joint probability, we can write expectation values of eventwise counters (separately for inclusive, fixed-N, or fixed-n cases) as
(99)E[ρ^(Qaa)]=∫Ω‍dQaadQbcP(Qaa,Qbc)ρ^aa=∑i≠j‍∫Ω‍dQaaPs(Qaa)δ(Qaa-Qijaa)=∑i≠j‍Ps(Qijaa),E[ρ^(Qbc)]=∫Ω‍dQaadQbcP(Qaa,Qbc)ρ^bc=∑i,j‍∫Ω‍dQbcPm(Qbc)δ(Qbc-Qijbc)=∑i,j‍Pm(Qijbc).
Later, we will meet expectation values for cases such asa=c,
(100)E[ρ^(Qaa)ρ^(Qab)]=∑i≠j‍∑k,ℓ‍∫Ω‍dQaadQabP(Qaa,Qab)×δ(Qaa-Qijaa)δ(Qab-Qkℓab)=∑i≠j‍∑k,ℓ‍P(Qijaa,Qkℓab),
which definitely do not factorise. The above expressions can be simplified because we know that the parent distribution is not a function of the individual track indices i,j,k, and ℓ(101)P(Qijaa,Qkℓbc)=P(Qaa,Qbc)∀i,j,k,ℓ
and similarly Ps(Qijaa)=Ps(Qaa) and Pm(Qkℓbc)=Pm(Qbc). For the event-averaged counters, this results in
(102)E[ρaa(Q)]=〈na2_〉NPs(Qaa),E[ρbc(Q)]=〈nb〉N〈nc〉NPm(Qbc),
or in terms of the notation of Section 3.4,
(103)E[ρ2N(Q)]=〈n2_〉NPs(Qaa),E[ρ1N⊗ρ1N(Q)]=〈n〉N2Pm(Qbc).
As mentioned, we do not need the factorisation (97) of P(Qaa,Qbc) as long as we keep careful track of the equal-event-indices cases. Whenever a≠b or a≠c, independence of the events ensures that expectation values of products of any functions f(Qaa) and g(Qab) of the pair variables do factorise,
(104)E[f(Qaa)g(Qab)]=E[f(Qaa)]E[g(Qab)]a≠b.
For third-order correlations, the parent distribution is a function of three different variables Qaaa, Qbbc, and Qdef containing, respectively, three, two, or one track from the same event, and corresponding considerations regarding equal and unequal event indices apply there, too.


### 4.3. Statistical Error Calculation from First Principles

It was shown in Section4.2 that expectation values would have well-defined meanings in terms of underlying parent distributions and their marginals if  their parent distributions were known, which, however, they are not. We are therefore forced to revert from expectation values E[·] to sample averages 〈·〉 after completing a calculation. The real use of such expectation values in frequentist statistics has been in the form of a gedankenexperiment which we now reproduce from Kendall [68]. Let x be any generic eventwise counter or any other eventwise statistic. Since the formulae in this section remain true for inclusive and fixed-N samples, we omit any notation related to N in this derivation. In this simplified notation, the well-known standard error of the sample mean 〈x〉 is given by (simplifying ℰ-1→ℰ)
(105)σ(〈x〉)=var(〈x〉)=1ℰ[〈x2〉-〈x〉2],
which follows from the combinatorics of equal and unequal event indices by the above artificial use of expectation values, reverting from expectation values E[·] to sample means 〈·〉 in the last step:
(106)var(〈x〉)=E[〈x〉2]-E[〈x〉]2=1ℰ2∑a,b‍[E[xaxb]-E[xa]E[xb]](107)=1ℰ2∑a=b‍[E[xaxb]-E[xa]E[xb]]+1ℰ2∑a≠b‍[E[xa]E[xb]-E[xa]E[xb]]=1ℰ2∑a‍[E[x2]-E[x]2]+0=1ℰ[〈x2〉-〈x〉2],
where we have used the fact that E[xaxb]=E[xa]E[xb] for all a≠b and assumed that all x are identically distributed, E[xa]=E[x] for all a. Equality or inequality of event indices is thus crucial. We will follow the same approach below, keeping careful track of equal and unequal event indices, factorising expectation values for unequal event indices, and reverting to sample means in the last step.


### 4.4. Variances and Covariances for Multiple Event Averages

#### 4.4.1. Statistical Errors for Second-Order Cumulants

According to (88)–(94), we must handle variances and covariances of products of several event averages. To derive these, we will use the following shortened notation: letting δijab≡δ(Q-Qijab) and so forth, then ρ^aa=∑i≠j‍δ(Q-Qijaa)=∑i≠j‍δijaa is the eventwise pair counter for event a, while ρ^bc=∑i,j‍δ(Q-Qijbc)=∑i,j‍δijbc is the mixed-event counter of events b and c (with b≠c≠a assumed), so that ρ2(Q)=〈ρ^aa〉a while ρ1⊗ρ1(Q)=〈ρ^ab〉ab is a double event average. We reserve the event index a for the “sibling” event whose correlations are currently being analysed and use indices b,c,…,t,u,v,w,… for events entering the event-mixing parts.5 All quantities are assumed to be measured within a particular subsample 𝒮N, but we omit the N-subscript and the argument. The event-index subscripts such as 〈·〉bc above are included or omitted depending on whether they convey relevant information on the specific averaging.In this notation, the method that led to (107) reads
(108)var(ρ2)=var(〈ρ^aa〉a)=1ℰ2∑a,b‍[E(ρ^aaρ^bb)-E(ρ^aa)E(ρ^bb)]=1ℰ[〈(ρ^aa)2〉-〈ρ^aa〉2].
The same method of disentangling the combinatorics of equal and unequal event indices is applied consistently to all variances and covariances below. The ρ12 term in second-order cumulants has variance
(109)var(ρ1⊗ρ1)=var(〈ρ^bc〉bc)=1(ℰ2_)2∑b≠c‍∑d≠e‍[E(ρ^bcρ^de)-E(ρ^bc)E(ρ^de)].
The case b≠c≠d≠e yields zero, but the cases b=d≠c≠e and three other equivalent combinations yield
(110)var(ρ1⊗ρ1)=4ℰ3_(ℰ2_)2[E[ρ^bcρ^be]-E[ρ^bc]E[ρ^be]]⟶4ℰ3_(ℰ2_)2[〈ρ^bcρ^be〉bce-〈ρ^bc〉〈ρ^be〉]⟶4ℰ[〈ρ^bcρ^be〉-(ρ1⊗ρ1)2],
where in the second step we reverted E[·]→〈·〉 and in the third6 assumed ℰ≫1. Note that the requirement b≠c≠e implies that 〈ρ^bcρ^be〉bce=〈〈∑i≠j‍δijbc〉c〈∑k≠l‍δkℓbe〉e〉b cannot be simplified to the square of a single counter 〈〈∑i≠j‍δijbc〉c2〉b: the event mixing involves three different events, not two. Secondly, the combinations of two equalities b=d≠c=e and b=e≠c=d in (109) yield another term of order ℰ-2,
(111)2ℰ2_[〈ρ^bdρ^bd〉bd-(ρ1⊗ρ1)2],
which we can safely neglect when 〈n2〉/ℰ≪1 except when there are few bins or large multiplicities even in small bins. It is worth emphasising that the extra factor 4 which appears in (110) arises from the same method that has been used for decades to justify use of (105). We find, by the same method, that the covariance between ρ2 and ρ1⊗ρ1 is given by
(112)cov(ρ2,ρ1⊗ρ1)=cov(〈ρ^dd〉,〈ρ^bc〉)=1ℰℰ2_∑d‍∑b≠c‍[E(ρ^ddρ^bc)-E(ρ^dd)E(ρ^bc)]=2ℰ[〈ρ^ddρ^db〉-〈ρ^dd〉〈ρ^dc〉]=2ℰ[〈ρ^ddρ^dc〉-(ρ2)(ρ1⊗ρ1)],
so that we must, in addition, accumulate, for every event d, the product of the counters
(113)ρ^ddρ^dc=∑i≠j‍δijdd∑k‍〈∑ℓ‍δkℓdc〉c.
Note that there is no restriction on track indices k≠i or k≠j in the d-event, meaning that events with n(a)=2 contribute to this counter which would otherwise not be the case. Combining these, we find, to leading order in ℰ-1 and renaming mixed-event indices,
(114)var(κ2I)=1ℰ{(ρ2-2F2Nρ1⊗ρ1)2〈(ρ^aa-2F2Nρ^ac)(ρ^aa-2F2Nρ^ad)〉a-(ρ2-2F2Nρ1⊗ρ1)2},
with all event indices strictly unequal and c- and d event averages understood where appropriate.7 Contrasting this with the traditional way to calculate the same variance,
(115)var(κ2I)=1ℰ{(ρ2-F2Nρ1⊗ρ1)2〈(ρ^aa-F2Nρ^ac)(ρ^aa-F2Nρ^ad)〉-(ρ2-F2Nρ1⊗ρ1)2},
it is clear that in previous analyses the two possible ways to set a equal to b or c were overlooked, while normal (non-internal) cumulants also omit the F2N.


#### 4.4.2. Statistical Errors for Third-Order Cumulants

In third-order, we will needδijkabc≡δ(Qa-Qijkabc) and similar quantities and the notation for counters ρ^aaa, ρ^aab, and ρ^abc corresponding to the event averages ρ3(Qa)=〈ρ^aaa〉, ρ2⊗ρ1(Qa)=〈ρ^aab〉, and ρ1⊗ρ1⊗ρ1(Qa)=〈ρ^abc〉, respectively. Clearly, a≠b≠c must hold in the third-order case. We obtain for the necessary third-order quantities (shuffling and/or renaming indices if necessary)
(116)var(ρ3)=1ℰ[〈(ρ^aaa)2〉a-ρ32]cov(ρ3,ρ2⊗ρ1)=1ℰℰ2_∑r‍∑s≠t‍[E(ρ^rrrρ^sst)-E(ρ^rrr)E(ρ^sst)]=1ℰ[〈ρ^rrrρ^rrs〉+〈ρ^rrrρ^rss〉-2ρ3(ρ2⊗ρ1)],
with ρ^rrrρ^rrs=∑i≠j≠k‍δijkrrr∑ℓ≠m‍〈∑n‍δℓmnrrs〉s and ρ^rrrρ^rss=∑i≠j≠k‍δijkrrr∑ℓ‍〈∑m≠n‍δℓmnrss〉s. The remaining variances and covariances needed for third-order correlations with GHP topology are, after renaming of indices,
(117)var(ρ2⊗ρ1)=ℰ3_(ℰ2_)2∑g≠e‍∑c≠d‍[E(ρ^ggeρ^ccd)-E(ρ^gge)E(ρ^ccd)]=ℰ3_(ℰ2_)2[〈ρ^ggd[ρ^ggc+ρ^gcc+ρ^ddc+ρ^dcc]〉-4(ρ2⊗ρ1)2],
while we neglect
(118)1ℰ2_[ρ^ggcρ^ggc+ρ^ggcρ^gcc-2(ρ2⊗ρ1)2].
The next term is simpler,
(119)cov(ρ3,ρ13)=1ℰℰ3_∑t‍∑u≠v≠w‍[E(ρ^tttρ^uvw)-E(ρ^ttt)E(ρ^uvw)]=3ℰ[〈ρ^tttρ^tuv〉-ρ3ρ13],
but the following is not,
(120)cov(ρ2⊗ρ1,ρ13)=1ℰ2_ℰ3_∑u≠v‍∑x≠y≠z‍[E(ρ^uuvρ^xyz)-E(ρ^uuv)E(ρ^xyz)]=3ℰ4_ℰ2_ℰ3_[〈ρ^wwxρ^wyz〉+〈ρ^wxxρ^wyz〉-2(ρ2ρ1)(ρ13)]+6ℰ2_[〈ρ^wwxρ^wxy〉-(ρ2ρ1)(ρ13)],
and the large number of combinations makes the variance of ρ13 particularly complicated,
(121)var(ρ1⊗ρ1⊗ρ1)=1(ℰ3_)2∑d≠e≠h‍∑b≠p≠q‍[E(ρ^dehρ^bpq)-E(ρ^deh)E(ρ^bpq)]=9ℰ5_(ℰ3_)2[〈ρ^behρ^bpq〉-(ρ13)(ρ13)]+18ℰ4_(ℰ3_)2[〈ρ^behρ^beq〉-(ρ13)(ρ13)]+6ℰ3_[〈ρ^behρ^beh〉-(ρ13)(ρ13)].
For large ℰ, the leading order terms will usually dominate, so that we can neglect the subleading terms.8 To leading order, we therefore obtain after substitution in (91) and again omitting brackets for non-a event averages
(122)var(κ3I)=1ℰ{〈ρ^aaa2+9ρ^aab(ρ^aac+ρ^acc+ρ^bbc+ρ^bcc)+9G32ρ^aaa(ρ^abcρ^ade)-6ρ^aaa(ρ^aab+ρ^abb)+6G3ρ^aaaρ^abc-18G3ρ^aab(ρ^acd+ρ^bcd)ρ^aaa2〉a-[ρ32+36(ρ2⊗ρ1)2+9G32(ρ13)2-12ρ3(ρ2⊗ρ1)+6G3ρ3(ρ1)3-36G32(ρ2⊗ρ1)(ρ1)3]},=1ℰ{〈(ρ^aaa-3ρ^aab-3ρ^abb+3G3ρ^abc)×(ρ^aaa-3ρ^aad-3ρ^add+3G3ρ^ade)〉a-[ρ3-6ρ2⊗ρ1+3G3(ρ1)3]2}.
While the factorised form is again instructive, it cannot be calculated in this form within the a-loop in the analysis since the G3N constants are known only on completion of the entire sample analysis. Rather, the full palette of product counters ρ^ρ^ has to be accumulated and averaged and combined only in the final phase of the analysis.


## 4.1. Propagation of Statistical Errors

Because cumulants can be measured only through the moments that enter their definitions, the first task is to identify which moment variances and covariances are needed. By means of standard error propagation [68], we find the sample variances for second-order internal cumulants of (81) and (83) to be
(88)var(κ2I(Q∣𝒮N))=var(ρ2N-F2Nρ1N2)=var(ρ2N)+F2N2var(ρ1N2)-2F2Ncov(ρ2N,ρ1N2),(89)var(K2I(Q∣𝒮N))=var(1F2Nρ2Nρ1N2-1)(90)=[1F2Nρ2Nρ1N2]2×[var(ρ2N)ρ2N2+var(ρ1N2)ρ1N4-2cov(ρ2N,ρ1N2)ρ2N·ρ1N2],
under the assumption that var(F2N) is much smaller than the other variances, so that F2N can be treated as a constant; this is the case if there are many bins for Q. Similarly, from (69), the variance of the unnormalised internal cumulant is, assuming G3 of (70) to be constant,
(91)var(κ3I(Qa∣𝒮N))=var[ρ3N-3ρ2Nρ1N+G3ρ1N3](92)=var(ρ3N)+9var(ρ2Nρ1N)+G32var(ρ1N3)-6cov(ρ3N,ρ2Nρ1N)+2G3cov(ρ3N,ρ1N3)-6G3cov(ρ2Nρ1N,ρ1N3),
while the normalised version has variance (again assuming var(G3)≪varρr)
(93)var(K3I(Qa∣𝒮N))=var[ρ3N-3ρ2Nρ1N+G3ρ1N3F3Nρ1N3]=1F3N2var[ρ3N-3ρ2Nρ1Nρ1N3+G3](94)=[ρ3N-3ρ2Nρ1NF3Nρ1N3]2×[var(ρ3N-3ρ2Nρ1N)(ρ3N-3ρ2Nρ1N)2+var(ρ1N3)ρ1N6-2cov(ρ3N-3ρ2Nρ1N,ρ1N3)(ρ3N-3ρ2Nρ1N)ρ1N3]=[ρ3N-3ρ2Nρ1NF3Nρ1N3]2×[var(ρ3N)+9var(ρ2Nρ1N)-6cov(ρ3N,ρ2Nρ1N)(ρ3N-3ρ2Nρ1N)2+var(ρ1N3)ρ1N6+6cov(ρ2Nρ1N,ρ1N3)-2cov(ρ3N,ρ1N3)(ρ3N-3ρ2Nρ1N)ρ1N3].
The unnormalised cumulants (81) and (82) and their variances (88) and (92) require knowledge of the multiplicity factorial moments F2N,F3N, so that the individual terms must be accumulated until the entire sample has been analysed. By contrast, the normalised cumulants (83) and (84) and their variances (90) and (94) contain the multiplicity moments only as prefactors.


## 4.2. Expectation Values of Counters

While we will not make direct use of the results in this section, it is nevertheless useful briefly to consider what we might mean by an “expectation value of experimental counters and densities.” For any scalar functionf(p) of the momenta, the theoretical expectation value E[f] is defined as the integral over the entire outcome space Ω of f weighted by a “parent distribution” P(p), an abstract entity supposedly containing everything there is to know on this level,
(95)E[f(p)]=∫Ω‍dpP(p)f(p).
Purely theoretical concepts such as P(p) and E[f] should be given little or no room in a strongly experimentally-oriented study. In calculating standard errors on counters below, we will, however, make use of the exact factorisation that expectation values provide whenever two variables x,y are statistically independent, E[xy]=∫‍dxdyP(x,y)xy=E[x]E[y].Expectation values for pairwise variables such as the four-momentum differenceQ we are considering here must be based on the underlying physics. We can deduce some properties of the parent distribution based on the usual definition of the femtoscopic correlation function:
(96)C2(Q)≡K2(Q)+1=ρ2sibling(Q)ρ2reference(Q)=〈ρ^(Qaa)〉〈〈ρ^(Qab)〉〉,
which is a function of two entirely different quantities: the four-momentum differences Qaa of “sibling” tracks taken from the same event a and one constructed from the mixed-event sample using tracks from different events, written as Qab, Qbc, and so forth. For second-order correlations, the parent distribution is therefore necessarily a two-variable probability4  P(Qaa,Qbc) which, depending on whether the cases b=a and c=a occur, may or may not factorise into “sibling” and “mixed” marginal probabilities
(97)P(Qaa,Qbc)=Ps(Qaa)Pm(Qbc)iffa≠b≠c,
but (unless a=b=c), the marginals will always be
(98)Ps(Qaa)=∫‍dQbcP(Qaa,Qbc),Pm(Qbc)=∫‍dQaaP(Qaa,Qbc).
The shapes of Ps(Q) and Pm(Q) must necessarily be different since it is precisely this difference that leads to a nontrivial signal in (96). In terms of this joint probability, we can write expectation values of eventwise counters (separately for inclusive, fixed-N, or fixed-n cases) as
(99)E[ρ^(Qaa)]=∫Ω‍dQaadQbcP(Qaa,Qbc)ρ^aa=∑i≠j‍∫Ω‍dQaaPs(Qaa)δ(Qaa-Qijaa)=∑i≠j‍Ps(Qijaa),E[ρ^(Qbc)]=∫Ω‍dQaadQbcP(Qaa,Qbc)ρ^bc=∑i,j‍∫Ω‍dQbcPm(Qbc)δ(Qbc-Qijbc)=∑i,j‍Pm(Qijbc).
Later, we will meet expectation values for cases such asa=c,
(100)E[ρ^(Qaa)ρ^(Qab)]=∑i≠j‍∑k,ℓ‍∫Ω‍dQaadQabP(Qaa,Qab)×δ(Qaa-Qijaa)δ(Qab-Qkℓab)=∑i≠j‍∑k,ℓ‍P(Qijaa,Qkℓab),
which definitely do not factorise. The above expressions can be simplified because we know that the parent distribution is not a function of the individual track indices i,j,k, and ℓ(101)P(Qijaa,Qkℓbc)=P(Qaa,Qbc)∀i,j,k,ℓ
and similarly Ps(Qijaa)=Ps(Qaa) and Pm(Qkℓbc)=Pm(Qbc). For the event-averaged counters, this results in
(102)E[ρaa(Q)]=〈na2_〉NPs(Qaa),E[ρbc(Q)]=〈nb〉N〈nc〉NPm(Qbc),
or in terms of the notation of Section 3.4,
(103)E[ρ2N(Q)]=〈n2_〉NPs(Qaa),E[ρ1N⊗ρ1N(Q)]=〈n〉N2Pm(Qbc).
As mentioned, we do not need the factorisation (97) of P(Qaa,Qbc) as long as we keep careful track of the equal-event-indices cases. Whenever a≠b or a≠c, independence of the events ensures that expectation values of products of any functions f(Qaa) and g(Qab) of the pair variables do factorise,
(104)E[f(Qaa)g(Qab)]=E[f(Qaa)]E[g(Qab)]a≠b.
For third-order correlations, the parent distribution is a function of three different variables Qaaa, Qbbc, and Qdef containing, respectively, three, two, or one track from the same event, and corresponding considerations regarding equal and unequal event indices apply there, too.


## 4.3. Statistical Error Calculation from First Principles

It was shown in Section4.2 that expectation values would have well-defined meanings in terms of underlying parent distributions and their marginals if  their parent distributions were known, which, however, they are not. We are therefore forced to revert from expectation values E[·] to sample averages 〈·〉 after completing a calculation. The real use of such expectation values in frequentist statistics has been in the form of a gedankenexperiment which we now reproduce from Kendall [68]. Let x be any generic eventwise counter or any other eventwise statistic. Since the formulae in this section remain true for inclusive and fixed-N samples, we omit any notation related to N in this derivation. In this simplified notation, the well-known standard error of the sample mean 〈x〉 is given by (simplifying ℰ-1→ℰ)
(105)σ(〈x〉)=var(〈x〉)=1ℰ[〈x2〉-〈x〉2],
which follows from the combinatorics of equal and unequal event indices by the above artificial use of expectation values, reverting from expectation values E[·] to sample means 〈·〉 in the last step:
(106)var(〈x〉)=E[〈x〉2]-E[〈x〉]2=1ℰ2∑a,b‍[E[xaxb]-E[xa]E[xb]](107)=1ℰ2∑a=b‍[E[xaxb]-E[xa]E[xb]]+1ℰ2∑a≠b‍[E[xa]E[xb]-E[xa]E[xb]]=1ℰ2∑a‍[E[x2]-E[x]2]+0=1ℰ[〈x2〉-〈x〉2],
where we have used the fact that E[xaxb]=E[xa]E[xb] for all a≠b and assumed that all x are identically distributed, E[xa]=E[x] for all a. Equality or inequality of event indices is thus crucial. We will follow the same approach below, keeping careful track of equal and unequal event indices, factorising expectation values for unequal event indices, and reverting to sample means in the last step.


## 4.4. Variances and Covariances for Multiple Event Averages

### 4.4.1. Statistical Errors for Second-Order Cumulants

According to (88)–(94), we must handle variances and covariances of products of several event averages. To derive these, we will use the following shortened notation: letting δijab≡δ(Q-Qijab) and so forth, then ρ^aa=∑i≠j‍δ(Q-Qijaa)=∑i≠j‍δijaa is the eventwise pair counter for event a, while ρ^bc=∑i,j‍δ(Q-Qijbc)=∑i,j‍δijbc is the mixed-event counter of events b and c (with b≠c≠a assumed), so that ρ2(Q)=〈ρ^aa〉a while ρ1⊗ρ1(Q)=〈ρ^ab〉ab is a double event average. We reserve the event index a for the “sibling” event whose correlations are currently being analysed and use indices b,c,…,t,u,v,w,… for events entering the event-mixing parts.5 All quantities are assumed to be measured within a particular subsample 𝒮N, but we omit the N-subscript and the argument. The event-index subscripts such as 〈·〉bc above are included or omitted depending on whether they convey relevant information on the specific averaging.In this notation, the method that led to (107) reads
(108)var(ρ2)=var(〈ρ^aa〉a)=1ℰ2∑a,b‍[E(ρ^aaρ^bb)-E(ρ^aa)E(ρ^bb)]=1ℰ[〈(ρ^aa)2〉-〈ρ^aa〉2].
The same method of disentangling the combinatorics of equal and unequal event indices is applied consistently to all variances and covariances below. The ρ12 term in second-order cumulants has variance
(109)var(ρ1⊗ρ1)=var(〈ρ^bc〉bc)=1(ℰ2_)2∑b≠c‍∑d≠e‍[E(ρ^bcρ^de)-E(ρ^bc)E(ρ^de)].
The case b≠c≠d≠e yields zero, but the cases b=d≠c≠e and three other equivalent combinations yield
(110)var(ρ1⊗ρ1)=4ℰ3_(ℰ2_)2[E[ρ^bcρ^be]-E[ρ^bc]E[ρ^be]]⟶4ℰ3_(ℰ2_)2[〈ρ^bcρ^be〉bce-〈ρ^bc〉〈ρ^be〉]⟶4ℰ[〈ρ^bcρ^be〉-(ρ1⊗ρ1)2],
where in the second step we reverted E[·]→〈·〉 and in the third6 assumed ℰ≫1. Note that the requirement b≠c≠e implies that 〈ρ^bcρ^be〉bce=〈〈∑i≠j‍δijbc〉c〈∑k≠l‍δkℓbe〉e〉b cannot be simplified to the square of a single counter 〈〈∑i≠j‍δijbc〉c2〉b: the event mixing involves three different events, not two. Secondly, the combinations of two equalities b=d≠c=e and b=e≠c=d in (109) yield another term of order ℰ-2,
(111)2ℰ2_[〈ρ^bdρ^bd〉bd-(ρ1⊗ρ1)2],
which we can safely neglect when 〈n2〉/ℰ≪1 except when there are few bins or large multiplicities even in small bins. It is worth emphasising that the extra factor 4 which appears in (110) arises from the same method that has been used for decades to justify use of (105). We find, by the same method, that the covariance between ρ2 and ρ1⊗ρ1 is given by
(112)cov(ρ2,ρ1⊗ρ1)=cov(〈ρ^dd〉,〈ρ^bc〉)=1ℰℰ2_∑d‍∑b≠c‍[E(ρ^ddρ^bc)-E(ρ^dd)E(ρ^bc)]=2ℰ[〈ρ^ddρ^db〉-〈ρ^dd〉〈ρ^dc〉]=2ℰ[〈ρ^ddρ^dc〉-(ρ2)(ρ1⊗ρ1)],
so that we must, in addition, accumulate, for every event d, the product of the counters
(113)ρ^ddρ^dc=∑i≠j‍δijdd∑k‍〈∑ℓ‍δkℓdc〉c.
Note that there is no restriction on track indices k≠i or k≠j in the d-event, meaning that events with n(a)=2 contribute to this counter which would otherwise not be the case. Combining these, we find, to leading order in ℰ-1 and renaming mixed-event indices,
(114)var(κ2I)=1ℰ{(ρ2-2F2Nρ1⊗ρ1)2〈(ρ^aa-2F2Nρ^ac)(ρ^aa-2F2Nρ^ad)〉a-(ρ2-2F2Nρ1⊗ρ1)2},
with all event indices strictly unequal and c- and d event averages understood where appropriate.7 Contrasting this with the traditional way to calculate the same variance,
(115)var(κ2I)=1ℰ{(ρ2-F2Nρ1⊗ρ1)2〈(ρ^aa-F2Nρ^ac)(ρ^aa-F2Nρ^ad)〉-(ρ2-F2Nρ1⊗ρ1)2},
it is clear that in previous analyses the two possible ways to set a equal to b or c were overlooked, while normal (non-internal) cumulants also omit the F2N.


### 4.4.2. Statistical Errors for Third-Order Cumulants

In third-order, we will needδijkabc≡δ(Qa-Qijkabc) and similar quantities and the notation for counters ρ^aaa, ρ^aab, and ρ^abc corresponding to the event averages ρ3(Qa)=〈ρ^aaa〉, ρ2⊗ρ1(Qa)=〈ρ^aab〉, and ρ1⊗ρ1⊗ρ1(Qa)=〈ρ^abc〉, respectively. Clearly, a≠b≠c must hold in the third-order case. We obtain for the necessary third-order quantities (shuffling and/or renaming indices if necessary)
(116)var(ρ3)=1ℰ[〈(ρ^aaa)2〉a-ρ32]cov(ρ3,ρ2⊗ρ1)=1ℰℰ2_∑r‍∑s≠t‍[E(ρ^rrrρ^sst)-E(ρ^rrr)E(ρ^sst)]=1ℰ[〈ρ^rrrρ^rrs〉+〈ρ^rrrρ^rss〉-2ρ3(ρ2⊗ρ1)],
with ρ^rrrρ^rrs=∑i≠j≠k‍δijkrrr∑ℓ≠m‍〈∑n‍δℓmnrrs〉s and ρ^rrrρ^rss=∑i≠j≠k‍δijkrrr∑ℓ‍〈∑m≠n‍δℓmnrss〉s. The remaining variances and covariances needed for third-order correlations with GHP topology are, after renaming of indices,
(117)var(ρ2⊗ρ1)=ℰ3_(ℰ2_)2∑g≠e‍∑c≠d‍[E(ρ^ggeρ^ccd)-E(ρ^gge)E(ρ^ccd)]=ℰ3_(ℰ2_)2[〈ρ^ggd[ρ^ggc+ρ^gcc+ρ^ddc+ρ^dcc]〉-4(ρ2⊗ρ1)2],
while we neglect
(118)1ℰ2_[ρ^ggcρ^ggc+ρ^ggcρ^gcc-2(ρ2⊗ρ1)2].
The next term is simpler,
(119)cov(ρ3,ρ13)=1ℰℰ3_∑t‍∑u≠v≠w‍[E(ρ^tttρ^uvw)-E(ρ^ttt)E(ρ^uvw)]=3ℰ[〈ρ^tttρ^tuv〉-ρ3ρ13],
but the following is not,
(120)cov(ρ2⊗ρ1,ρ13)=1ℰ2_ℰ3_∑u≠v‍∑x≠y≠z‍[E(ρ^uuvρ^xyz)-E(ρ^uuv)E(ρ^xyz)]=3ℰ4_ℰ2_ℰ3_[〈ρ^wwxρ^wyz〉+〈ρ^wxxρ^wyz〉-2(ρ2ρ1)(ρ13)]+6ℰ2_[〈ρ^wwxρ^wxy〉-(ρ2ρ1)(ρ13)],
and the large number of combinations makes the variance of ρ13 particularly complicated,
(121)var(ρ1⊗ρ1⊗ρ1)=1(ℰ3_)2∑d≠e≠h‍∑b≠p≠q‍[E(ρ^dehρ^bpq)-E(ρ^deh)E(ρ^bpq)]=9ℰ5_(ℰ3_)2[〈ρ^behρ^bpq〉-(ρ13)(ρ13)]+18ℰ4_(ℰ3_)2[〈ρ^behρ^beq〉-(ρ13)(ρ13)]+6ℰ3_[〈ρ^behρ^beh〉-(ρ13)(ρ13)].
For large ℰ, the leading order terms will usually dominate, so that we can neglect the subleading terms.8 To leading order, we therefore obtain after substitution in (91) and again omitting brackets for non-a event averages
(122)var(κ3I)=1ℰ{〈ρ^aaa2+9ρ^aab(ρ^aac+ρ^acc+ρ^bbc+ρ^bcc)+9G32ρ^aaa(ρ^abcρ^ade)-6ρ^aaa(ρ^aab+ρ^abb)+6G3ρ^aaaρ^abc-18G3ρ^aab(ρ^acd+ρ^bcd)ρ^aaa2〉a-[ρ32+36(ρ2⊗ρ1)2+9G32(ρ13)2-12ρ3(ρ2⊗ρ1)+6G3ρ3(ρ1)3-36G32(ρ2⊗ρ1)(ρ1)3]},=1ℰ{〈(ρ^aaa-3ρ^aab-3ρ^abb+3G3ρ^abc)×(ρ^aaa-3ρ^aad-3ρ^add+3G3ρ^ade)〉a-[ρ3-6ρ2⊗ρ1+3G3(ρ1)3]2}.
While the factorised form is again instructive, it cannot be calculated in this form within the a-loop in the analysis since the G3N constants are known only on completion of the entire sample analysis. Rather, the full palette of product counters ρ^ρ^ has to be accumulated and averaged and combined only in the final phase of the analysis.


## 4.4.1. Statistical Errors for Second-Order Cumulants

According to (88)–(94), we must handle variances and covariances of products of several event averages. To derive these, we will use the following shortened notation: letting δijab≡δ(Q-Qijab) and so forth, then ρ^aa=∑i≠j‍δ(Q-Qijaa)=∑i≠j‍δijaa is the eventwise pair counter for event a, while ρ^bc=∑i,j‍δ(Q-Qijbc)=∑i,j‍δijbc is the mixed-event counter of events b and c (with b≠c≠a assumed), so that ρ2(Q)=〈ρ^aa〉a while ρ1⊗ρ1(Q)=〈ρ^ab〉ab is a double event average. We reserve the event index a for the “sibling” event whose correlations are currently being analysed and use indices b,c,…,t,u,v,w,… for events entering the event-mixing parts.5 All quantities are assumed to be measured within a particular subsample 𝒮N, but we omit the N-subscript and the argument. The event-index subscripts such as 〈·〉bc above are included or omitted depending on whether they convey relevant information on the specific averaging.In this notation, the method that led to (107) reads
(108)var(ρ2)=var(〈ρ^aa〉a)=1ℰ2∑a,b‍[E(ρ^aaρ^bb)-E(ρ^aa)E(ρ^bb)]=1ℰ[〈(ρ^aa)2〉-〈ρ^aa〉2].
The same method of disentangling the combinatorics of equal and unequal event indices is applied consistently to all variances and covariances below. The ρ12 term in second-order cumulants has variance
(109)var(ρ1⊗ρ1)=var(〈ρ^bc〉bc)=1(ℰ2_)2∑b≠c‍∑d≠e‍[E(ρ^bcρ^de)-E(ρ^bc)E(ρ^de)].
The case b≠c≠d≠e yields zero, but the cases b=d≠c≠e and three other equivalent combinations yield
(110)var(ρ1⊗ρ1)=4ℰ3_(ℰ2_)2[E[ρ^bcρ^be]-E[ρ^bc]E[ρ^be]]⟶4ℰ3_(ℰ2_)2[〈ρ^bcρ^be〉bce-〈ρ^bc〉〈ρ^be〉]⟶4ℰ[〈ρ^bcρ^be〉-(ρ1⊗ρ1)2],
where in the second step we reverted E[·]→〈·〉 and in the third6 assumed ℰ≫1. Note that the requirement b≠c≠e implies that 〈ρ^bcρ^be〉bce=〈〈∑i≠j‍δijbc〉c〈∑k≠l‍δkℓbe〉e〉b cannot be simplified to the square of a single counter 〈〈∑i≠j‍δijbc〉c2〉b: the event mixing involves three different events, not two. Secondly, the combinations of two equalities b=d≠c=e and b=e≠c=d in (109) yield another term of order ℰ-2,
(111)2ℰ2_[〈ρ^bdρ^bd〉bd-(ρ1⊗ρ1)2],
which we can safely neglect when 〈n2〉/ℰ≪1 except when there are few bins or large multiplicities even in small bins. It is worth emphasising that the extra factor 4 which appears in (110) arises from the same method that has been used for decades to justify use of (105). We find, by the same method, that the covariance between ρ2 and ρ1⊗ρ1 is given by
(112)cov(ρ2,ρ1⊗ρ1)=cov(〈ρ^dd〉,〈ρ^bc〉)=1ℰℰ2_∑d‍∑b≠c‍[E(ρ^ddρ^bc)-E(ρ^dd)E(ρ^bc)]=2ℰ[〈ρ^ddρ^db〉-〈ρ^dd〉〈ρ^dc〉]=2ℰ[〈ρ^ddρ^dc〉-(ρ2)(ρ1⊗ρ1)],
so that we must, in addition, accumulate, for every event d, the product of the counters
(113)ρ^ddρ^dc=∑i≠j‍δijdd∑k‍〈∑ℓ‍δkℓdc〉c.
Note that there is no restriction on track indices k≠i or k≠j in the d-event, meaning that events with n(a)=2 contribute to this counter which would otherwise not be the case. Combining these, we find, to leading order in ℰ-1 and renaming mixed-event indices,
(114)var(κ2I)=1ℰ{(ρ2-2F2Nρ1⊗ρ1)2〈(ρ^aa-2F2Nρ^ac)(ρ^aa-2F2Nρ^ad)〉a-(ρ2-2F2Nρ1⊗ρ1)2},
with all event indices strictly unequal and c- and d event averages understood where appropriate.7 Contrasting this with the traditional way to calculate the same variance,
(115)var(κ2I)=1ℰ{(ρ2-F2Nρ1⊗ρ1)2〈(ρ^aa-F2Nρ^ac)(ρ^aa-F2Nρ^ad)〉-(ρ2-F2Nρ1⊗ρ1)2},
it is clear that in previous analyses the two possible ways to set a equal to b or c were overlooked, while normal (non-internal) cumulants also omit the F2N.


## 4.4.2. Statistical Errors for Third-Order Cumulants

In third-order, we will needδijkabc≡δ(Qa-Qijkabc) and similar quantities and the notation for counters ρ^aaa, ρ^aab, and ρ^abc corresponding to the event averages ρ3(Qa)=〈ρ^aaa〉, ρ2⊗ρ1(Qa)=〈ρ^aab〉, and ρ1⊗ρ1⊗ρ1(Qa)=〈ρ^abc〉, respectively. Clearly, a≠b≠c must hold in the third-order case. We obtain for the necessary third-order quantities (shuffling and/or renaming indices if necessary)
(116)var(ρ3)=1ℰ[〈(ρ^aaa)2〉a-ρ32]cov(ρ3,ρ2⊗ρ1)=1ℰℰ2_∑r‍∑s≠t‍[E(ρ^rrrρ^sst)-E(ρ^rrr)E(ρ^sst)]=1ℰ[〈ρ^rrrρ^rrs〉+〈ρ^rrrρ^rss〉-2ρ3(ρ2⊗ρ1)],
with ρ^rrrρ^rrs=∑i≠j≠k‍δijkrrr∑ℓ≠m‍〈∑n‍δℓmnrrs〉s and ρ^rrrρ^rss=∑i≠j≠k‍δijkrrr∑ℓ‍〈∑m≠n‍δℓmnrss〉s. The remaining variances and covariances needed for third-order correlations with GHP topology are, after renaming of indices,
(117)var(ρ2⊗ρ1)=ℰ3_(ℰ2_)2∑g≠e‍∑c≠d‍[E(ρ^ggeρ^ccd)-E(ρ^gge)E(ρ^ccd)]=ℰ3_(ℰ2_)2[〈ρ^ggd[ρ^ggc+ρ^gcc+ρ^ddc+ρ^dcc]〉-4(ρ2⊗ρ1)2],
while we neglect
(118)1ℰ2_[ρ^ggcρ^ggc+ρ^ggcρ^gcc-2(ρ2⊗ρ1)2].
The next term is simpler,
(119)cov(ρ3,ρ13)=1ℰℰ3_∑t‍∑u≠v≠w‍[E(ρ^tttρ^uvw)-E(ρ^ttt)E(ρ^uvw)]=3ℰ[〈ρ^tttρ^tuv〉-ρ3ρ13],
but the following is not,
(120)cov(ρ2⊗ρ1,ρ13)=1ℰ2_ℰ3_∑u≠v‍∑x≠y≠z‍[E(ρ^uuvρ^xyz)-E(ρ^uuv)E(ρ^xyz)]=3ℰ4_ℰ2_ℰ3_[〈ρ^wwxρ^wyz〉+〈ρ^wxxρ^wyz〉-2(ρ2ρ1)(ρ13)]+6ℰ2_[〈ρ^wwxρ^wxy〉-(ρ2ρ1)(ρ13)],
and the large number of combinations makes the variance of ρ13 particularly complicated,
(121)var(ρ1⊗ρ1⊗ρ1)=1(ℰ3_)2∑d≠e≠h‍∑b≠p≠q‍[E(ρ^dehρ^bpq)-E(ρ^deh)E(ρ^bpq)]=9ℰ5_(ℰ3_)2[〈ρ^behρ^bpq〉-(ρ13)(ρ13)]+18ℰ4_(ℰ3_)2[〈ρ^behρ^beq〉-(ρ13)(ρ13)]+6ℰ3_[〈ρ^behρ^beh〉-(ρ13)(ρ13)].
For large ℰ, the leading order terms will usually dominate, so that we can neglect the subleading terms.8 To leading order, we therefore obtain after substitution in (91) and again omitting brackets for non-a event averages
(122)var(κ3I)=1ℰ{〈ρ^aaa2+9ρ^aab(ρ^aac+ρ^acc+ρ^bbc+ρ^bcc)+9G32ρ^aaa(ρ^abcρ^ade)-6ρ^aaa(ρ^aab+ρ^abb)+6G3ρ^aaaρ^abc-18G3ρ^aab(ρ^acd+ρ^bcd)ρ^aaa2〉a-[ρ32+36(ρ2⊗ρ1)2+9G32(ρ13)2-12ρ3(ρ2⊗ρ1)+6G3ρ3(ρ1)3-36G32(ρ2⊗ρ1)(ρ1)3]},=1ℰ{〈(ρ^aaa-3ρ^aab-3ρ^abb+3G3ρ^abc)×(ρ^aaa-3ρ^aad-3ρ^add+3G3ρ^ade)〉a-[ρ3-6ρ2⊗ρ1+3G3(ρ1)3]2}.
While the factorised form is again instructive, it cannot be calculated in this form within the a-loop in the analysis since the G3N constants are known only on completion of the entire sample analysis. Rather, the full palette of product counters ρ^ρ^ has to be accumulated and averaged and combined only in the final phase of the analysis.


## 5.  Averaged Internal Cumulants

AsN is only an approximation for the true total event multiplicity anyway, and for cases of small sample statistics, it may be necessary or desirable to group subsamples of fixed N into multiplicity classes N∈[A,B]. It is important, however, not to simply lump all events within this multiplicity class into a single “half-inclusive” subsample, because, as has long been known [15], that results in terms entering the cumulants which arise solely to “multiplicity mixing” (MM) of events of different N. Given the arbitrary choice of [A,B], such MM correlations are spurious and avoided in favour of “Averaged-Internal” (AI) correlations9 defined as follows. Using the renormalised multiplicity distribution
(123)ℛN′=ℛN∑N=ABℛN=ℰN∑N=ABℰN,
the AI unnormalised cumulants, reference distributions, and normalised AI cumulants are
(124)κ2AI(Q∣𝒮AB)=∑N=AB‍ℛN′κ2I(Q∣𝒮N),(125)ρ1⊗ρ1(Q∣𝒮AB)=∑N=AB‍ℛN′F2Nρ1⊗ρ1(Q∣𝒮N),(126)K2AI(Q∣𝒮AB)=κ2AI(Q∣𝒮AB)ρ1⊗ρ1(Q∣𝒮AB)=∑N=ABℛN′ρ2(Q∣𝒮N)∑N=ABℛN′F2Nρ1⊗ρ1(Q∣𝒮N)-1.
Note that the correction factors F2N, which are normalised factorial moments of n for fixed N, are part of the summed normalisations.10 In third-order, we have correspondingly
(127)κ3AI(Qa∣𝒮AB)=∑N=AB‍ℛN′κ3I(Qa∣𝒮N),(128)ρ1⊗ρ1⊗ρ1(Qa∣𝒮AB)=∑N=AB‍ℛN′F3Nρ1⊗ρ1⊗ρ1(Qa∣𝒮N),(129)K3AI(Qa∣𝒮AB)=κ3AI(Qa∣𝒮AB)ρ1⊗ρ1⊗ρ1(Qa∣𝒮AB).
Note also that (128) holds for the normalisation only and not for the last term in κ3I, which is (3F2N-F3N)ρ1⊗ρ1⊗ρ1— but that is already taken care of in the formula (69) for κ3I itself. Expressions for an inclusive (all-N) multiplicity summation of internal cumulants are obtained from the above by setting A=0 and B=∞. The AI (Averaged Internal) correlations (124) and (127) represent refined versions of what has traditionally been termed “Short-Range Correlations,” differing from the original formulae [15, 17] by the F2N and G3N factors, respectively. This was originally pointed out in [18] but only for multinomials in N.Regarding variances and standard errors for AI correlations, we first note that, since subsamples𝒮N are strictly mutually independent, a variance over the [A,B] range is simply the weighted sum of the corresponding fixed-N variances. From (124) and (129), we have for all orders r(130)κrAI(Q∣𝒮AB)=∑N=AB‍ℛN′κrI(Q∣𝒮N)r=2,3,…,
and given the independence of any functions f and g of different multiplicity subsamples, E[f(𝒮N)·g(𝒮N′)]=E[f(𝒮N)]·E[g(𝒮N′)] for all N≠N′, we conclude that
(131)var[κrAI(Q∣𝒮AB)]=∑N‍(ℛN′)2var[κrI(Q∣𝒮N)],(132)var[ρ1r(Q∣𝒮AB)]=∑N‍(ℛN′)2(FrN)2var[ρ1r(Q∣𝒮N)],(133)cov[κrAI(Q∣𝒮AB),ρ1r(Q∣𝒮AB)]=∑N‍(ℛN′)2FrNcov[κrI(Q∣𝒮N),ρ1r(Q∣𝒮N)],
which are known functions in terms of Sections 4.4.1 and 4.4.2, while for the normalised cumulants in [A,B], covariances between numerator and denominator are (omitting the Q)
(134)cov[κ2I(𝒮N),ρ12(𝒮N)]=cov[ρ2(𝒮N),ρ1⊗ρ1(𝒮N)]-F2Nvar(ρ1⊗ρ1(𝒮N)),cov[κ3I(𝒮N),ρ13(𝒮N)]=cov[ρ3(𝒮N),ρ13(𝒮N)]-3cov[ρ2⊗ρ1(𝒮N),ρ13(𝒮N)]+G3var[ρ13(𝒮N)],
so that the normalised range cumulants have variances
(135)var[KrAI(Q∣𝒮AB)]=(KrAI)2·[var[κr(𝒮AB)](κr(𝒮AB))2+var[ρ1r(𝒮AB)](ρ1r(𝒮AB))2-2cov[κr(𝒮AB),ρ1r(𝒮AB)]κr(𝒮AB)·ρ1r(𝒮AB)],
and standard errors are given by11(136)σ(KrAI(𝒮AB))=var(KrAI(Q∣𝒮AB)).


## 6. Event Mixing Algorithms

“Event mixing” [72] is widely used to simulate uncorrelated or semicorrelated quantities such as ρ13 and ρ2⊗ρ1. The idea has always been to use the experimental sample at hand to simulate the baseline of independence referred to in Section 3.1 in such a way that criteria 2 (independence of momenta), 4 (reproducing the one-particle momentum space distribution), and 6 (normalisation) are all addressed simultaneously. Ideally, all effects bar the desired correlation are elegantly removed in this way.For the internal cumulants and their variances and covariances derived above, event mixing requires keeping track of counters of all orders in each of the subsamples𝒮N. A count of event indices in Section 4 shows that in a brute-force calculation we would need, for each subsample, a minimum of five independent event averages or O(ℰ5) event combinations; furthermore, caution would advise not to use the same event in calculating related counters, so that selection and use of more than five events in mixing are advisable. The resulting excessive number of full event averages, mixing every (sub) sample event with every other one, is therefore not feasible.If the order of events in the sample is random, the multiple event averages can be simplified by the use of the following multiple event buffer algorithm.(1)
A single overall event loop equivalent to the event indexa runs over the entire inclusive sample 𝒮. A given event a will have a multiplicity N=N(a), so for that particular a correlation analysis for subsample 𝒮N is advanced by one event while the others remain dormant.12 In this way, a, which always refers to the sibling event, runs over all ℰN events of every subsample 𝒮N. There is no need to either explicitly sort the inclusive sample into subsamples or to run multiple jobs for fixed N.(2)
The firstℰB events13 of a given multiplicity N are used solely to fill up the buffer without doing any analysis. Once a given buffer has been filled, event mixing analysis proceeds for that subsample as follows.
(a)
A newly reada-event is assigned to the N=N(a) buffer, the earliest event in that buffer is discarded, and sums for averages entering F2N and F3N as well as the sibling counters ρ^aa, ρ^aa, (ρ^aa)2, and (ρ^aaa)2 are updated.(b)
Event combinations for mixed-event counters are built up by picking any one of theℰB-1 other events in that buffer and calling it b, thereafter picking any one of the remaining ℰB-2 events in the buffer, calling it c and so on. While for third-order, only five events (including the sibling event) are needed to construct all the counters required; in practice it is better to use different mixing events for different counters to root out even traces of unwanted correlations between different mixing counters. The random selection of events rather than tracks for mixing is necessary to ensure that more than one track per event can be used as required for counters of Sections 3 and 4 such as ρ^bbc.(c)
For a given event setb,c,d,…, the mixing counters are incremented using all possible combinations of the na tracks in event a together with all the nb,nc,… tracks in the selected events b,c,… mixing events. For example, ρ^bbc would use all possible nb(nb-1) pairs of b-tracks14 together with all possible nc single c-tracks. The mixing of all tracks of a given event rather than just selected ones ensures that the fluctuations in n for given fixed N are automatically contained in the counters.(d)
For constanta, the process of selecting events b,c,… is repeated Cmix= 10–100 times to reduce the statistical errors, avoiding events that have been used in previous selections. Efficiency can be maximised by tuning of both the number of events ℰB stored in each buffer and by the number of resamples Cmix.(3)
Once the entire sample has been processed via thea-event loop, the b- c- d-event averages are normalised by Cmix(ℰN-ℰB), while the primary a-average is normalised by (ℰN-ℰB).(4)
Unnormalised and normalised correlation quantities, their standard errors, and correction factors are constructed by appropriate combinations of averaged counters.(5)
Results from fixed-N subsamples can then be combined into AI correlations over partial ranges of N or the entire inclusive sample at the end of the event loop using the methods outlined in Section 5.


## 7. Discussion and Conclusions

(
1
)  Correlations are only defined properly if the null case or reference distribution is defined on the same level of sophistication as the correlation itself. Translated into statistics, the six criteria set out in Section 3.1 for a reference sample for correlations at fixed charged multiplicity N lead straight to the definition of the reference sample as the average of multinomials given in (51), weighted by the conditional multiplicity distribution ℛnN. Assigning Bernoulli probabilities α(p∣𝒮nN)=〈ρ(p∣𝒮nN)〉N/〈n〉N yields the reference density (56) and generating functional (58).(
2
) From the theorem that internal cumulants are given by the difference between measured and reference cumulants, we obtain normalised and unnormalised internal cumulants which satisfy every stated criterion for proper correlations for fixed-N samples.(
3
)  We have highlighted the distinction between n, the particles entering the correlation analysis itself, and N, the particles determining the event selection criterion for a particular semi-inclusive subsample. Various correction factors are shown to be fair to good approximations of these exact results in some cases but far off the mark in others. Surprisingly, normalised cumulants are far more sensitive to these corrections, through the normalisation prefactor, than their unnormalised counterparts. To belabour the point, for any variables (x1,x2), different definitions for correction factor 1/F for fixed-N correlations of positive pions,
(137)K2(x1,x2∣𝒮N)=1Fρ(x1,x2∣𝒮N)ρ(x1∣𝒮N)ρ(x2∣𝒮N)-1,
can be very important at low multiplicities, with F=1 (Poisson) being the worst approximation, F=N(N-1)/N2 being a fair one, and 〈n(n-1)〉N/〈n〉N2 being the best.For inclusive correlations, correction factors such as〈n〉2/〈n(n-1)〉incl were proposed early on in [73, 74] in an approach based on probabilities rather than densities. Reference [75] specifically calls the inclusion of these correction factors meaningless because the theory then requires that the emission function be identically zero. We note that the argument in all those references relates to inclusive samples, while for the samples of fixed N considered in this paper the prefactor, which is an average ofnat fixedN, is a necessity. Either way, the arbitrariness of the use or nonuse of the prefactor has been eliminated here based solely on considerations related to the reference distribution.(
4
)  The problem posed in this paper, namely, the relation between charged multiplicity on the one hand and correlations based on the conditional n-distribution ℛnN has attracted little attention in the literature. Indeed, almost all theoretical work on multipion correlations, as for example summarised in [76], starts from the projection of final-state events, with all their different particle species, onto the single-species subspace of either n+ (positive pions) or n- (negative pions) correlations, to the exclusion of the other charge. It would be interesting to see a combined theory for both n+andn-, which would encompass all the work done so far plus correlations between unlike-sign pions and, of course, the issue raised by us here.(
5
)  As shown in Figures 2 and 3, correction factors for third-order are larger than the second-order ones. For higher rth-order correlations, the effect of using a fixed-N subsample is suppressed by approximately 1/〈nr-1〉N for unnormalised cumulants but actually worsens for normalised cumulants due to the normalisation prefactors 〈n〉Nr/〈nr_〉N for small n. The importance of accurate correction of normalised quantities therefore rises with order of correlation.(
6
)  The difference between Poissonian and internal cumulants is largest at small multiplicities n. The mixed-multinomial prescription will therefore be required for any correlation analysis involving small n, independently of the magnitude of N. Apart from the usual suspects of leptonic, hadronic, and low-energy collisions, the low-n case occurs both for very restricted phase space (such as in spectrometer experiments) and for correlations of rare particles such as kaons and baryons, even for large N.(
7
)  Since n fluctuates according to the conditional multiplicity distribution ℛnN, the degree to which fixed-N correlations differ from inclusive ones is strongly coupled to the character of ℛnN. In general, ℛnN is sub-Poissonian and so FrN falls below the Poisson limit of 1. The correction from fixed-N Poissonian to internal normalised cumulants is hence upward, not downward as in the case of multiplicity-mixing corrections.(
8
) While not the main subject of the present paper, some light is cast on the relationship between three levels of correlation, namely, correlations inherent in the overall multiplicity distribution, multiplicity-mixing correlations, and the true internal correlations for fixed N. Each of these can and should be treated separately. The averaged-internal correlations of Section 5 are a compromise solution which may be useful both for physics reasons and for small datasets.(
9
) The fixed-N corrections discussed here are separate and complementary to other important effects at low multiplicity. References [76, 77] highlight, for example, possible effects of “residual correlations” resulting from projecting from multipion to two-pion correlations.Energy-momentum conservation would also play a role. Borghini [78, 79] has, for example, calculated the effect of momentum conservation for normalised two- and three-particle cumulants in momenta and pt. However, the saddlepoint method used applies to the large-N limit, and the results cannot be directly applied to the low-N (and hence low-n) samples under discussion here. Indeed, momentum conservation will be near-irrelevant for cases of large N and small n as discussed above, but the small-n corrections of this paper will remain important. For the specific case of like-sign pion femtoscopy, the fact that only n~𝒜N/2Ω out of the N charged pions are used and that momentum conservation constraints include all other final-state particles both imply that momentum conservation constraints may be less important than the internal-cumulant correction introduced here.For the specific choice of correlation variablesQ and Qa for two- and three-particle cumulants, the contribution of momentum conservation to cumulants at small Q will be small since the counts will be dominated by pairs at small (Δϕ,Δy) and intermediate (pt1,pt2). As pointed out by [1], momentum conservation exerts the greatest influence at large pair or triplet momenta and hence mostly at large Q, where it may lead to moments and cumulants which do not converge to a constant as presupposed in most fits. The ad hoc method of multiplying fit parametrisations by a prefactor 1+cQ with c a free parameter does not adequately address the problem.Regarding the multiplicity dependence of the influence of energy-momentum conservation, [2] calculates the effect of energy-momentum conservation on single-particle differential observables and finds significant systematic effects. No doubt this must also be the case for multiparticle observables, although as we have pointed out above, the effect of conservation laws will be diluted by the fact that fewer than half of the final-state particles of any given event are actually used in the present analysis. Detailed investigations are beyond the scope of this paper.(
10
) We have also recalculated statistical errors for products of event averages starting from the original prescription which forms the basis of frequentist statistical error calculations. Compared to conventional statistical error calculations, new prefactors appear in our calculations — see for example (110) and the results in Section 4.4.2—which have surprisingly been missed so far.(
11
) We note that the present formalism is still in the frequentist statistics mindset, which may be inaccurate for low multiplicities and should be supplanted by a proper Bayesian analysis. The final word has certainly not been spoken about correlation analysis of small-n datasets.



---

*Source: 230515-2013-07-29.xml*