Multiparticle Cumulant Method in Flow Analysis

• The multiparticle cumulant method is a connected-correlation framework that subtracts lower-order contributions to isolate genuine collective behavior in collision data.
• It utilizes generating functions, Q-vectors, and subevent techniques to extract azimuthal anisotropies and suppress non-flow effects in small systems.
• The approach extends to high-order, mixed-harmonic, differential, and multiplicity analyses, providing insights into fluctuation strengths and collective dynamics.

The multiparticle cumulant method is a connected-correlation framework in which lower-order disconnected contributions are subtracted from multiparticle observables so that the remaining quantity isolates genuine collective or correlated structure. In relativistic collision studies, it is most often formulated for azimuthal anisotropy through cumulants such as $c_n\{2\}$, $c_n\{4\}$, symmetric cumulants, and asymmetric cumulants, implemented with Q-vectors and, in small systems, with pseudorapidity-separated subevents to reduce non-flow. Closely related cumulant constructions are also used for identified-particle multiplicity fluctuations, multibin long-range correlations, multi-particle balance functions, and femtoscopic many-body correlation functions (Nie et al., 2018, Bialas et al., 2011, Pruneau et al., 2023, Grande et al., 2021).

1. Statistical basis and connected correlations

At its most general level, the method is organized by generating functions and their logarithms. For multiplicity fluctuations of a species with event-by-event multiplicity $N$, the ordinary cumulants $C_n$ are defined from the moment generating function $M(t)=\langle e^{tN}\rangle$ and the cumulant generating function $K(t)=\ln M(t)$ through

$$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$

The same formalism admits factorial moments and factorial cumulants through

$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$

with

$G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$

The ordinary and factorial cumulants are related by Stirling-number expansions such as

$$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$

(Nabroth, 17 Oct 2025).

For long-range multiplicity correlations in several bins, the basic observables are the joint factorial moments

$c_n\{4\}$0

where $c_n\{4\}$1, together with the multivariate factorial moment generating function

$c_n\{4\}$2

The joint factorial cumulants are obtained from $c_n\{4\}$3 by differentiation at the origin. This construction makes explicit that mixed cumulants vanish if the bins are statistically independent, whereas nonzero mixed cumulants signal connected long-range structure (Bialas et al., 2011).

The same connected-correlation logic appears in differential correlation functions. Ordinary cumulant densities $c_n\{4\}$4 are obtained from inclusive densities $c_n\{4\}$5 by cluster expansion, and integral factorial cumulants are their phase-space integrals. In the balance-function framework, this distinction between differential cumulants and their acceptance-integrated counterparts is essential because differential objects retain the full $c_n\{4\}$6 dependence while integral cumulants compress that information into acceptance-dependent fluctuation measures (Pruneau et al., 2023).

2. Azimuthal cumulants and the Q-cumulant formalism

For azimuthal anisotropy, the starting point is the Fourier expansion

$c_n\{4\}$7

Multiparticle azimuthal correlators are built from particle angles and averaged first over distinct tuples within an event and then over events. For harmonic $c_n\{4\}$8, the fundamental two- and four-particle correlators are

$c_n\{4\}$9

$N$0

The corresponding standard cumulants are

$N$1

$N$2

In the pure-flow limit,

$N$3

so that

$N$4

(Nie et al., 2018).

The standard implementation uses event-wise flow vectors

$N$5

or, in weighted analyses,

$N$6

For unit weights, the two-particle correlator is

$N$7

For weighted reference particles, one convenient form is

$N$8

Higher-order correlators are expressed through combinations of $N$9, $C_n$0, and higher harmonics, and can be generated recursively in a single loop over particles and events (Collaboration, 2021, Nadderd et al., 2022).

This formalism extends beyond the basic two- and four-particle case. The recursive cumulant relations

$C_n$1

$C_n$2

$C_n$3

generalize systematically to higher order (Nadderd et al., 2022, Moravcova et al., 2020).

3. Subevent cumulants and the suppression of non-flow

In small systems, low multiplicity makes non-flow especially problematic. Jets, dijets, resonance decays, Coulomb/HBT correlations, and global momentum conservation can dominate standard cumulants. The subevent cumulant method addresses this by splitting the event into pseudorapidity-separated regions and requiring particles in a correlator to come from different subevents (Nie et al., 2018).

For two subevents $C_n$4 and $C_n$5,

$C_n$6

For three subevents, a particularly clean four-particle construction is

$C_n$7

$C_n$8

The same logic extends to symmetric cumulants,

$C_n$9

and to their two- and three-subevent variants (Nie et al., 2018).

The physical rationale is explicit. Requiring particles to come from different $M(t)=\langle e^{tN}\rangle$0-separated subevents removes same-jet and resonance pairs, three-subevent constructions break typical dijet topologies, and $M(t)=\langle e^{tN}\rangle$1 gaps reduce Coulomb/HBT contributions. In AMPT for $M(t)=\langle e^{tN}\rangle$2+Pb at $M(t)=\langle e^{tN}\rangle$3 TeV, the standard method gave $M(t)=\langle e^{tN}\rangle$4 that changed sign, becoming positive at smaller $M(t)=\langle e^{tN}\rangle$5 where non-flow dominates, whereas both two- and three-subevent $M(t)=\langle e^{tN}\rangle$6 remained negative across the full $M(t)=\langle e^{tN}\rangle$7 range. In the same study, standard-method $M(t)=\langle e^{tN}\rangle$8 and $M(t)=\langle e^{tN}\rangle$9 were comparable to data, while the subevent results behaved differently, indicating contamination of the standard method, especially when the number of produced particles is small (Nie et al., 2018).

Small-system measurements at RHIC reinforce the same point in a different way. In $K(t)=\ln M(t)$0+Au collisions, PHENIX observed real-valued $K(t)=\ln M(t)$1, that is $K(t)=\ln M(t)$2, at $K(t)=\ln M(t)$3, $K(t)=\ln M(t)$4, $K(t)=\ln M(t)$5, and $K(t)=\ln M(t)$6 GeV, while in $K(t)=\ln M(t)$7+Au at $K(t)=\ln M(t)$8 GeV the four-particle cumulant had the opposite sign. At $K(t)=\ln M(t)$9 GeV in $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$0+Au, $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$1 agreed with $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$2, so that $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$3, which PHENIX interpreted as evidence that nonflow effects are subdominant in the four-particle signal (Aidala et al., 2017).

4. Mixed-harmonic, differential, and high-order generalizations

The method extends naturally from single-harmonic cumulants to correlations among different harmonics. The symmetric cumulant

$$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$4

measures whether events with larger $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$5 tend to have larger or smaller $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$6 (Nie et al., 2018). The asymmetric cumulant

$$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$7

projects the nonlinear coupling between $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$8 and $$C_n=\left.\frac{\partial^n K(t)}{\partial t^n}\right|_{t=0}.$$9, while the four-particle symmetric cumulant

$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$0

probes correlations between $$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$1 and $$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$2. Their normalized forms,

$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$3

were calculated from transverse momentum conservation and flow and found to be in good agreement with recent ATLAS subevent measurements in $$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$4 collisions (Pei et al., 2024).

For rare probes, the cumulant construction becomes differential. In PbPb collisions at $$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$5 TeV, prompt $$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$6 mesons were analyzed with a two-subevent reference built from HF$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$7 and HF$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$8 towers. The differential four-particle cumulant is

$$F(s)=\langle(1+s)^N\rangle,\qquad G(s)=\ln F(s),\qquad \kappa_n=\left.\frac{\partial^n G(s)}{\partial s^n}\right|_{s=0},$$9

and the corresponding four-particle flow coefficient is

$G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$0

The reported $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$1 rose with $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$2 to a maximum near $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$3 GeV and remained systematically below $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$4, providing direct access to event-by-event fluctuations of charm-quark azimuthal anisotropy (Collaboration, 2021).

The same logic generalizes to arbitrary mixtures of particles of interest and inclusive reference particles. Mixed correlators of the form

$G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$5

lead to mixed cumulants $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$6, $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$7, $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$8, and $G(s)=K(\ln(1+s)),\qquad K(t)=G(e^t-1).$9, together with rare-probe estimators such as

$$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$0

This construction was proposed specifically to study fluctuations and correlations in the azimuthal anisotropies of rare probes (Holtermann et al., 2023).

At higher order, several general frameworks now exist. A generic recursive algorithm produced explicit new 10-, 12-, 14-, and 16-particle single-harmonic cumulants and the corresponding $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$1 coefficients (Moravcova et al., 2020). A decomposition based on mathematical induction generated exact event-level expressions for high-order Q-cumulants and tabulated the needed integer coefficients and basis elements up to 14-particle order (Nadderd et al., 2022). A later partition-based method made the fast calculation of arbitrary-order $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$2 practical and was validated up to 40th order (Nađđerđ et al., 31 May 2025). In parallel, a multidimensional generating-function analysis classified 33 distinct cumulants with orders $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$3 and harmonics $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$4, including the genuine three-particle correlation function that contains information of all third-order cumulants (Taghavi, 2020).

5. Multiplicity cumulants, balance functions, and femtoscopic connected correlations

Although azimuthal flow dominates present usage, the same cumulant logic is also central in fluctuation analyses of identified-particle multiplicities. In Ag+Ag collisions at $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$5 GeV, HADES used ordinary cumulants $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$6, factorial cumulants $$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$7, and normalized factorial cumulants together with a probabilistic particle-identification method in which each track carries a species weight

$$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$8

Event-wise proxies such as

$$C_n=\sum_{k=1}^n S(n,k)\kappa_k,\qquad \kappa_n=\sum_{k=1}^n s(n,k)C_k$$9

permit the reconstruction of factorial moments and cumulants without hard PID cuts. Because MC closure tests showed that binomial efficiency corrections fail in the studied phase space, HADES used unfolding and full detector-response corrections, and removed centrality fluctuations with a data-driven event-mixing procedure (Nabroth, 17 Oct 2025).

The differential-cumulant framework also yields multi-particle balance functions. For the two-particle case,

$c_n\{4\}$00

while higher-order balance functions are built from connected $c_n\{4\}$01 combinations so that two- and three-prong resonance decays are suppressed. In full acceptance, these balance functions satisfy simple sum rules determined by quantum number conservation, and the paper argues that they arise as an intrinsic component of high-order net-charge cumulants (Pruneau et al., 2023).

A distinct but structurally related use of cumulants appears in femtoscopy. There, the measured three-body correlation function is decomposed as

$c_n\{4\}$02

with

$c_n\{4\}$03

The lower-order contribution is obtained by projecting measured pair correlations into the triplet variable $c_n\{4\}$04 through an analytic kernel

$c_n\{4\}$05

which isolates the genuine three-body connected correlation $c_n\{4\}$06 after subtraction (Grande et al., 2021).

6. Interpretation, model discrimination, and open issues

A central virtue of the method is its sensitivity to fluctuations. Under widely used assumptions for the $c_n\{4\}$07 distribution,

$c_n\{4\}$08

so ratios such as $c_n\{4\}$09 directly probe fluctuation strength (Collaboration, 2021). In $c_n\{4\}$10+Pb, AMPT reproduced $c_n\{4\}$11 but yielded a $c_n\{4\}$12 magnitude smaller than ATLAS data, which the study interpreted as either underestimated collectivity or missing non-Gaussian $c_n\{4\}$13 fluctuations (Nie et al., 2018).

The sign structure of cumulants is physically informative but not unique to a single mechanism. A parton-model initial-state calculation in $c_n\{4\}$14 collisions reproduced negative $c_n\{4\}$15, the ordering

$c_n\{4\}$16

and the qualitative behavior of symmetric cumulants, and argued that coherent multiple scattering is essential to obtain a positive definite $c_n\{4\}$17 (Dusling et al., 2017). In a related CGC analysis, the 3-dipole correlator up to order $c_n\{4\}$18 led to the conjecture

$c_n\{4\}$19

with the dilute-limit calculation of $c_n\{4\}$20 supporting that scaling (Zhang et al., 2019). These results imply that negative $c_n\{4\}$21 and higher-order hierarchy do not by themselves uniquely identify hydrodynamic flow.

Global constraints can also bias cumulants. For transverse momentum conservation, the leading contribution to the $c_n\{4\}$22-particle cumulant was derived as

$c_n\{4\}$23

so the sign is $c_n\{4\}$24 and the magnitude falls as $c_n\{4\}$25. The interplay of TMC and flow can change the signs of $c_n\{4\}$26, $c_n\{4\}$27, and $c_n\{4\}$28 at some multiplicities (Xie et al., 2022). The same interplay was shown to describe $c_n\{4\}$29 and $c_n\{4\}$30 measured with the subevent cumulant method, indicating that subevent cumulants suppress short-range non-flow but do not remove global TMC (Pei et al., 2024).

Taken together, these developments indicate that the multiparticle cumulant method is not a single observable but a hierarchy of connected-correlation constructions. In large systems it constrains the full joint probability density of flow magnitudes and symmetry planes; in small systems it separates collectivity from jets, resonances, and global constraints as far as the chosen subevent topology allows; in multiplicity, balance-function, and femtoscopic applications it isolates connected many-body structure from acceptance, PID overlap, and lower-order backgrounds. A plausible implication is that future progress will depend less on introducing entirely new observables than on combining high-order, mixed-harmonic, differential, and subevent cumulants with response-based modeling and detector-accurate inference frameworks (Moravcova et al., 2020, Nadderd et al., 2022, Nađđerđ et al., 31 May 2025).