Revealing long-range multi-particle collectivity in small collision systems via subevent cumulants

TL;DR

This paper introduces a subevent cumulant method to better identify genuine collective flow signals in small collision systems by suppressing non-flow effects like dijets, enhancing the study of long-range azimuthal correlations.

Contribution

The paper proposes a novel subevent cumulant technique that reduces non-flow contamination, improving the detection of true collective flow in small collision systems.

Findings

Standard cumulants are dominated by non-flow effects such as dijets.

The new subevent cumulant method effectively suppresses non-flow contributions.

The method can detect flow signals as low as 4% in PYTHIA simulations.

Abstract

Multi-particle azimuthal cumulants, often used to study collective flow in high-energy heavy-ion collisions, have recently been applied in small collision systems such as $pp$ and $p$+A to extract the second-order azimuthal harmonic flow $v_2$. Recent observation of four-, six- and eight-particle cumulants with "correct sign" $c_2\{4\}<0, c_2\{6\}>0, c_2\{8\}<0$ and approximate equality of the inferred single-particle harmonic flow, $v_2\{4\}\approx v_2\{6\}\approx v_2\{8\}$, have been used as strong evidence for a collective emission of all soft particles produced in the collisions. We show that these relations in principle could be violated due to the non-Gaussianity in the event-by-event fluctuation of flow and/or non-flow. Furthermore, we show, using $pp$ events generated with the PYTHIA model, that $c_2\{2k\}$ obtained with standard cumulant method are dominated by non-flow from dijets. An alternative cumulant method based on two or more $\eta$-separated subevents is proposed to suppress the dijet contribution. The new method is shown to be able to recover a flow signal as low as 4\% imposed on the PYTHIA events, independently of how the event activity class is defined. Therefore the subevent cumulant method offers a more robust way of studying collectivity based on the existence of long-range azimuthal correlations between multiple distinct $\eta$ ranges. The prospect of using the subevent cumulants to study collective flow in A+A collisions, in particular its longitudinal dynamics, is discussed.

Tables

Table 1. Table 1: The setup used in this analysis. The first column shows the three cumulant methods and corresponding η 𝜂 \eta ranges, the other three columns shows the p T subscript 𝑝 T p_{\mathrm{T}} selection used for calculating azimuthal correlation ⟨ 2 ​ k ⟩ delimited-⟨⟩ 2 𝑘 \left\langle 2k\right\rangle , for selecting events used for c n ​ { 2 ​ k } subscript 𝑐 𝑛 2 𝑘 c_{n}\{2k\} , and for defining event activity used to represent the x 𝑥 x -axis, respectively. See the text for more detail.

Cumulant types	$X_1$: criteria for $⟨2k⟩$	$X_2$: criterial for $N_{\text{ch}}^{\text{sel}}$	$X_3$: criteria for $N_{\mathrm{ch}}$
		used in event averaging	to represent the $x$-axis
standard method	GeV      GeV	same as for $X_1$     $p_T>0.2$ GeV     $p_T>0.4$ GeV     $p_T>0.6$ GeV
two-subevent method			$p_T>0.4$ GeV
$\eta \in [-2.5,0],[0,2.5]$
three-subevent method
$\eta \in [-2.5,-0.83],[-0.83,0.83],[0.83,2.5]$

Equations

$$\text{P}(\varphi )=\frac{1}{2\pi }\sum _{n=-\infty }^{\infty }{\mathbf{v}}_n{e}^{-in\varphi },{\mathbf{v}}_n=v_n{e}^{in{\Phi }_n},$$

$$\displaystyle p({\bm{v}}_{n})=\frac{1}{2\pi\delta^{2}_{n}}e^{-\left|{\bm{v}}_{n}-v_{n}^{\;{}_{0}}\right|^{2}\big{/}\left(2\delta^{2}_{n}\right)},\;p(v_{n})=\frac{v_{n}}{\delta_{n}^{2}}e^{-\frac{(v_{n})^{2}+(v_{n}^{0})^{2}}{2\delta_{n}^{2}}}I_{0}\left(\frac{v_{n}^{0}v_{n}}{\delta_{n}^{2}}\right).$$

$${\mathbf{Q}}_n\equiv \sum _i{e}^{in{\varphi }_i}=Q_n{e}^{in{\Psi }_n},{\mathbf{q}}_n\equiv \frac{{\sum }_i{e}^{in{\varphi }_i}}{M}=q_n{e}^{in{\Psi }_n}$$

$${\mathbf{q}}_n={\mathbf{v}}_n+{\mathbf{s}}_{\mathrm{n}}+{\mathbf{s}}_{\mathrm{n}}^{\mathrm{stat}}$$

$$p({\mathbf{q}}_n)=p({\mathbf{v}}_n)\otimes p({\mathbf{s}}_n)\otimes p({\mathbf{s}}_{\mathrm{n}}^{\mathrm{stat}})$$

$$\displaystyle\left\llangle 2k\right\rrangle=\left\llangle e^{in\sum_{j=1}^{k}(\phi_{2j-1}-\phi_{2j})}\right\rrangle=\left\langle v_{n}^{2k}\right\rangle.$$

$$\displaystyle\ln\left\langle e^{z({\bm{v}}_{n}+{\bm{v}}_{n}^{*})}\right\rangle=\ln\left(\sum_{k=1}^{\infty}\frac{z^{2k}}{k!^{2}}\left\langle({\bm{v}}_{n}{\bm{v}}_{n}^{*})^{k}\right\rangle\right)=\ln\left(\sum_{k=1}^{\infty}\frac{z^{2k}}{k!^{2}}\left\llangle 2k\right\rrangle\right)\equiv\sum_{k=1}^{\infty}\frac{z^{2k}}{k!^{2}}c_{n}\{2k\}$$

$$c_n\{2\}$$

$$c_n\{4\}$$

$$c_n\{6\}$$

$$v_n\{2\}=\sqrt{c_n\{2\}},v_n\{4\}=\sqrt[4]{-c_n\{4\}},v_n\{6\}=\sqrt[6]{c_n\{6\}/4}.$$

$$⟨2⟩$$

$$⟨4⟩$$

$$c_n\{2\}$$

$$c_n\{4\}$$

$$c_n\{6\}$$

$$v_n\{2k\}$$

$$p(v_n)=2\frac{\alpha }{\kappa }\frac{v_n}{\kappa }(1-(\frac{v_n}{\kappa }{)}^2{)}^{\alpha -1}$$

$$\ln ⟨e^{z({\mathbf{q}}_n+{\mathbf{q}}_n^{\ast })}⟩=\ln ⟨e^{z({\mathbf{v}}_n+{\mathbf{v}}_n^{\ast })}⟩+\ln ⟨e^{z({\mathbf{s}}_n+{\mathbf{s}}_n^{\ast })}⟩=\sum _{k=1}^{\infty }\frac{z^{2k}}{k{!}^2}\left(c_n\{2k,v\}+c_n\{2k,s\}\right),$$

$$c_n\{2k\}=c_n\{2k,v\}+c_n\{2k,s\}.$$

$$p(x;a)$$

$${⟨2⟩}_{a|b}$$

$${⟨4⟩}_{a,a|b,b}$$

$$W_{{⟨2⟩}_{a|b}}$$

$${\mathbf{Q}}_{n,o}\equiv \sum _i{e}^{in{\varphi }_i^o},{\mathbf{q}}_{n,o}\equiv \frac{{\sum }_i{e}^{in{\varphi }_i^o}}{M_o},{\tau }_o=\frac{1}{M_o},o=a\text{mbox}orb$$

$$\displaystyle\left\llangle 2\right\rrangle_{a|b}$$

$$\displaystyle\left\llangle 4\right\rrangle_{a,a|b,b}$$

$$\displaystyle c_{n}^{a,a|b,b}\{4\}\equiv\left\llangle 4\right\rrangle_{a,a|b,b}-2\left\llangle 2\right\rrangle_{a|b}^{2}=\left\langle e^{in(\phi_{1}^{a}+\phi_{2}^{a}-\phi_{3}^{b}-\phi_{4}^{b})}\right\rangle_{c}+\left\langle{\bm{v}}^{2}_{n,a}{\bm{v}}^{*2}_{n,b}\right\rangle-2\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,b}^{*}\right\rangle^{2}\;.$$

$${⟨2⟩}_{a|a}\equiv ⟨e^{in({\varphi }_1^a-{\varphi }_2^a)}⟩=\frac{Q_{n,a}^2-M_a}{M_a(M_a-1)}=\frac{q_{n,a}^2-{\tau }_a}{1-{\tau }_a},{⟨2⟩}_{b|b}\equiv ⟨e^{in({\varphi }_1^b-{\varphi }_2^b)}⟩=\frac{Q_{n,b}^2-M_b}{M_b(M_b-1)}=\frac{q_{n,b}^2-{\tau }_b}{1-{\tau }_b}$$

$${⟨4⟩}_{a,b|a,b}\equiv ⟨e^{in({\varphi }_1^a+{\varphi }_2^b-{\varphi }_3^a-{\varphi }_4^b)}⟩=\frac{(Q_{n,a}^2-M_a)(Q_{n,b}^2-M_b)}{M_a(M_a-1)M_b(M_b-1)}=\frac{(q_{n,a}^2-{\tau }_a)(q_{n,b}^2-{\tau }_b)}{(1-{\tau }_a)(1-{\tau }_b)}$$

$$\displaystyle c_{n}^{a,b|a,b}\{4\}\equiv\left\llangle 4\right\rrangle_{a,b|a,b}-\left\llangle 2\right\rrangle_{a|a}\left\llangle 2\right\rrangle_{b|b}-\left\llangle 2\right\rrangle_{a|b}^{2}=\left\langle e^{in(\phi_{1}^{a}+\phi_{2}^{b}-\phi_{3}^{a}-\phi_{4}^{b})}\right\rangle_{c}+\left\langle{\bm{v}}^{2}_{n,a}{\bm{v}}^{*2}_{n,b}\right\rangle-\left\langle v_{n,a}^{2}\right\rangle\left\langle v_{n,b}^{2}\right\rangle-\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,b}^{*}\right\rangle^{2}\;.$$

$${⟨4⟩}_{a,a|b,c}$$

$$\displaystyle\left\llangle 4\right\rrangle_{a,a|b,c}$$

$$\displaystyle c_{n}^{a,a|b,c}\{4\}\equiv\left\llangle 4\right\rrangle_{a,a|b,c}-2\left\llangle 2\right\rrangle_{a|b}\left\llangle 2\right\rrangle_{a|c}=\left\langle e^{in(\phi_{1}^{a}+\phi_{2}^{a}-\phi_{3}^{b}-\phi_{4}^{b})}\right\rangle_{c}+\left\langle{\bm{v}}^{2}_{n,a}{\bm{v}}^{*}_{n,b}{\bm{v}}^{*}_{n,c}\right\rangle-2\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,b}^{*}\right\rangle\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,c}^{*}\right\rangle$$

$${⟨4⟩}_{a,b|a,c}$$

Full text

I Introduction

High energy heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) create a strongly-interacting nuclear matter that exhibits many interesting characteristics. One such characteristic is the collimated emission of particle pairs with small azimuthal-angle separation, $\Delta\phi$, that extends over a large range of pseudorapidity differences, $\Delta\eta$. This so called “ridge” correlation was first observed in A+A collisions Adare et al. (2008); Abelev et al. (2009); Alver et al. (2010); Aamodt et al. (2011); Aad et al. (2012); Chatrchyan et al. (2014), and later was also observed in $pp$ Khachatryan et al. (2010); Aad et al. (2016); Aaboud et al. (2017a) and $p$+A collisions Chatrchyan et al. (2013a); Abelev et al. (2013); Aad et al. (2013a); Adare et al. (2013); Aad et al. (2014a); Khachatryan et al. (2015a). In A+A collisions, the ridge is believed to be the consequence of collective emission of particles in the azimuthal direction, and the collectivity is generated in the final state after local thermalization, described by relativistic viscous hydrodynamic models Heinz and Snellings (2013); Gale et al. (2013); Luzum and Petersen (2014); Jia (2014). For small systems such as $pp$ and $p$+A collisions, the origin of the ridge is less clear Dusling et al. (2016). Current efforts are focused on understanding whether the ridge in small systems reflects global collectivity of the event Loizides (2016), and if so whether it is of hydrodynamic origin similar to A+A collisions Bozek and Broniowski (2013) or it is created in the initial state from gluon saturation Dusling and Venugopalan (2013).

The ridge signal from two-particle correlation (2PC) is characterized by a Fourier decomposition \sim 1+2v_{n}^{2}\cos(n\mbox{\Delta\phi}), where the $v_{n}$ denotes the single-particle anisotropy harmonics. The second-order coefficient $v_{2}$ is by far the largest, followed by $v_{3}$. In small collision systems, the extraction of the ridge signal requires a careful removal of a large contribution from dijets, which is estimated from 2PC in very low multiplicity events and then subtracted from higher multiplicity events Abelev et al. (2013); Aad et al. (2013a, 2014a, 2016); Aaboud et al. (2017a); Khachatryan et al. (2017). On the other hand, since collectivity is intrinsically a multi-particle phenomenon, it can be probed more directly using the multi-particle correlations (or cumulants) technique Borghini et al. (2001a). One of the perceived hallmark features of collectivity is the observation of positively defined signal from $2k$-particle correlation $v_{2}\{2k\},k\geq 2$, which has been measured for $v_{2}\{4\}$, $v_{2}\{6\}$ and $v_{2}\{8\}$ in high-multiplicity $pp$ and $p$+Pb collisions Aad et al. (2013b); Khachatryan et al. (2015a, 2017). However, this perception could be wrong in small collision systems, where the non-flow correlations can be as large as or bigger than the genuine long-range ridge correlations. As we show in this paper, the sign of $v_{n}\{2k\}$ is sensitive to the event-by-event (ebye) fluctuations of flow and non-flow, and positive definiteness of $v_{n}\{2k\}$ in general is not required for flow correlations (also discussed in Ref. Jia and Radhakrishnan (2015)), and it could also be the result of non-flow correlations.

Multi-particle cumulants suppress short-range correlations, but does not completely remove them. In fact $v_{2}\{4\}$ is observed to change sign at a smaller number of charged particle, $N_{\mathrm{ch}}$, in $pp$ and $p$+Pb collisions Aad et al. (2013b); Khachatryan et al. (2015a, 2017). Recently ATLAS observed Aaboud et al. (2017b) that the $N_{\mathrm{ch}}$ value where sign-change happens and the magnitude of $v_{2}\{4\}$ depend on how the event classes are chosen for the calculation of cumulants. In this paper, we show that the choice of event class influences the probability distribution of non-flow, and consequently the non-flow contribution to the $v_{2}\{4\}$. An improved cumulant method is proposed to further suppress non-flow, and therefore reduce the sensitivity of $v_{2}\{4\}$ to non-flow fluctuations. In this method, cumulants are constructed from particles in several subevents separated in $\eta$. A subevent cumulant idea was proposed in Ref. Adler et al. (2002), where the particles in the event are divided “randomly” into four subevents. However, the purpose there was not to suppress non-flow, but instead to circumvent a brute-force nested-loop calculation.

The structure of the paper is as follows. In sections II- IV, we discuss the contributions of flow and non-flow to multi-particle cumulants and how such contributions depend on the probability distributions (or event-by-event fluctuations) of flow and non-flow. In sections V and VI, we introduce the subevent cumulant method, and demonstrate the robustness of the method using simulations based on PYTHIA8. The performance is compared to standard cumulant method focusing on four-particle correlations. In section VII, we discuss the implication of our findings for the interpretation of collectivity in small systems, as well as possible measurements enabled by our method, in particular for probing the longitudinal dynamics of collectivity. The formula for higher-order cumulants, and for symmetric and asymmetric cumulants are given in the Appendix.

II Event-by-event distribution of flow and non-flow

The azimuthal anisotropy of the particle production in an event can be characterized by Fourier expansion of the underlying probability distribution $\emph{P}(\phi)$ in azimuthal angle $\phi$:

[TABLE]

where $v_{n}$ and $\Phi_{n}$ are magnitude and phase, respectively. In heavy ion collisions, flow harmonics ${\bm{v}}_{n}$ vary event to event, and can be described by a probability distribution $p({\bm{v}}_{n})$ Luzum and Petersen (2014); Jia (2014). In A+A collisions, the flow fluctuations are close to Gaussian or equivalently the distribution of $v_{n}$ is Bessel-Gaussian after integrating out the $\phi$ angle Voloshin and Zhang (1996):

[TABLE]

The parameter $v_{n}^{\;{}_{0}}$ reflects the component driven by the average geometry of the overlap region; it is expected to be sizable only for $n=2$.

Due to the finite number of particles $M$ produced in each event, harmonic flow can only be estimated from the observed flow vector ${\bm{Q}}_{n}$ or per-particle normalized flow vector ${\bm{q}}_{n}$:

[TABLE]

where the sum runs over the particles in the event, $\phi_{i}$ are their azimuthal angles and $\Psi_{n}$ is the event plane. The magnitude and direction of ${\bm{q}}_{n}$ differ from those for the true flow, due to azimuthal fluctuations associated with finite particle multiplicity, denoted by ${\bm{s}}_{\mathrm{n}}^{\rm{stat}}$, as well as non-flow from various short-range correlations, denoted by ${\bm{s}}_{\mathrm{n}}$,

[TABLE]

In heavy ion collisions, all three components, flow, non-flow and finite number effects fluctuate event to event. Therefore the probability distributions for ${\bm{q}}_{n}$ and ${\bm{v}}_{n}$ can be related to each other by random smearing functions that reflect the non-flow and statistical fluctuation.

[TABLE]

The statistical fluctuation component usually cancels out after averaging over many events in two- or multi-particle correlation analyses (see the discussion in the next section). The non-flow fluctuations are more complicated, and they depend on the number of short-range sources, the particle multiplicity in each source and possible correlations between different sources (for example dijets) 111In this definition, the special case where the cluster sizes are fixed and only only the number of clusters fluctuates event by event, is considered as fluctuation of non-flow..

III Contribution of flow and non-flow to multi-particle cumulants

Let us first consider the case where there are only flow correlations for events with finite multiplicity. The moment of $p(v_{n})$ distribution can be extracted from multi-particle correlations: A $2k$-particle azimuthal correlator is obtained by averaging over all unique combinations in one event then over all events Borghini et al. (2001b); Bilandzic et al. (2011):

[TABLE]

where $\left\langle x^{2k}\right\rangle\equiv\int x^{2k}p(x)dx$ is the $2k$-th moment of the probability distribution for $x$, and we have used the fact that statistical fluctuations (${\bm{s}}_{\mathrm{n}}^{\rm{stat}}$) drop out after averaging over many events. The $2k$-particle cumulant is then obtained by proper combination of correlations involving $\leq 2k$ number of particles, whose expression can be obtained with the following generating function Borghini et al. (2001b):

[TABLE]

The formula for the first three are Borghini et al. (2001b):

[TABLE]

which leads to the following cumulant-based ${\it definition}$ of harmonic flow $v_{n}$ 222This definition of $v_{n}\{2k\}$ assumes the sign of $c_{n}\{2k\}$ is negative for even $k$, and positive for odd $k$, which is not true for an arbitrary $p(v_{n})$ distribution Jia and Radhakrishnan (2015).:

[TABLE]

Eqs. 6–10 play a crucial role for understanding how the flow and non-flow contribute to multi-particle correlations, so we shall discuss them in more detail below.

Let’s consider the usual expressions for two- and four-particle correlators in terms of ${\bm{q}}_{n}$ and ${\bm{Q}}_{n}$ for one event with $M$ particles Bilandzic et al. (2011); Bilandzic (2012),

[TABLE]

The advantage of using per-particle normalized flow vector ${\bm{q}}_{n}$ is that all quantities in the equation are smaller than one, and the terms can be sorted in powers of $\tau=1/M\ll 1$ and $q_{kn}\sim q_{n}^{k}$ Jia and Radhakrishnan (2015). The event-by-event weights are slightly modified from those in Ref. Bilandzic et al. (2011): $W_{\left\langle 2\right\rangle}=M(M-1)/2$, $W_{\left\langle 4\right\rangle}=M(M-1)(M-2)(M-3)/4!$ etc 333Since the correct weight should consider only unique $2k$-particle combinations instead of all permutations, the weight for 2k-particle correlation need to be divided by $(2k)!$. However since all events are affected by the same factor, this does not matter much in practice except for discussing the statistical power between different cumulant methods.. Each term has a simple interpretation. The two-particle correlator $\left\langle 2\right\rangle$ has two terms: the first term $Q_{n}^{2}$ contains $M^{2}$ pairs, and the second term corresponds to the contribution of $M$ duplicate pairs: $\phi_{1}=\phi_{2}$. Four-particle correlator $\left\langle 4\right\rangle$ has $P_{M,4}=M(M-1)(M-2)(M-3)$ quadruplets, expressed as $Q_{n}^{4}$ ($M^{4}$ quadruplets) minus contributions from quadruplets where the same particle appears more than once, e.g., combinations such as $\phi_{1}=\phi_{2}=\phi_{3}=\phi_{4}$, $\phi_{1}=\phi_{2}=\phi_{3}\neq\phi_{4}$, $\phi_{1}=\phi_{2}\neq\phi_{3}\neq\phi_{4}$ … etc. In total, there are $M^{4}-P_{M,4}$ quadruplets containing duplicated particles. The removal of duplicate combinations insures that the statistical fluctuation associated with finite particle multiplicity drops out from ${\bm{q}}_{n}$ after averaging over many events.

When only flow correlations are present, cumulants are fully determined from moments of the $p(v_{n})$:

[TABLE]

For any $p(v_{n})$ distribution, it is expected that $v_{n}\{2\}\geq v_{n}\{4\}$ since $v_{n}\{2\}^{4}-v_{n}\{4\}^{4}=\left\langle(v_{n}^{2}-\left\langle v_{n}^{2}\right\rangle)^{2}\right\rangle\geq 0$.

If the distribution of ${\bm{v}}_{n}$ is the Gaussian function defined in Eq. 2, flow harmonics defined by cumulants have a simple expression Voloshin et al. (2008):

[TABLE]

while the $v_{n}\{2\}$ also includes contribution from flow fluctuations, the higher-order cumulants $v_{n}\{4,6,..\}$ measure the component associated with average geometry.

In small systems, due to finite number of sources in the initial density distribution, flow fluctuations are expected to deviate from Gaussian, and were suggested to follow a power function Bzdak et al. (2014); Yan and Ollitrault (2014).

[TABLE]

where $\alpha$ is related to the number of sources for particle production in the initial state, and $\kappa$ is related to hydrodynamic response to the initial eccentricity $\epsilon_{n}$, $\kappa=v_{n}/\epsilon_{n}$. The higher-order cumulants $v_{n}\{4,6,..\}$ are expected to quickly converge to a non-zero value reflecting the influence of the non-Gaussian tail (controlled by $\alpha$) of the power distribution.

When both flow and non-flow are present and are un-correlated with each other, the generating function Eq. 7 becomes:

[TABLE]

and the measured cumulant is the sum of the separate contributions from cumulant of flow and cumulant of non-flow:

[TABLE]

Equation 19 provides a simple way to understand the influence of non-flow to $c_{n}\{2k\}$. In large collision systems, the number of non-flow sources is large and proportional to $M$. Since the orientations of these sources are weakly correlated, non-flow fluctuation $p(s_{n})$ is expected to approach Gaussian with a width that scales as $1/M$ (confirmed in HIJING simulation Jia and Mohapatra (2013)). In this case, the non-flow contributions to four-particle or higher-order cumulants are naturally suppressed, e.g., $c_{n}\{2k,s\}\approx 0$ for $k>1$. However, in small collision systems, the $p(s_{n})$ distribution is expected to be highly non-Gaussian since the number of sources is small but strongly fluctuates event to event, and the particle multiplicity in a single source can be large compared to $M$. In this case, the contribution from non-flow could be large, and the sign and magnitude of $c_{n}\{2k\}$ depend on the nature of $p(v_{n})$ and $p(s_{n})$.

IV Fluctuations of flow and non-flow and the sign of cumulants

The $c_{n}\{2k\}$ are more fundamental than the $v_{n}\{2k\}$; the latter is well motivated only when flow distribution is close to Gaussian. As we shall discuss below (also see Jia and Radhakrishnan (2015)), it is easy to find $p(v_{n})$ distributions for which $\left\langle v_{n}^{4}\right\rangle>2\left\langle v_{n}^{2}\right\rangle^{2}$ such that $c_{n}\{4\}>0$ and $v_{n}\{4\}$ become undefined 444This is also true for higher-order cumulants and Lee-Yang Zero method Bhalerao et al. (2003), whose existence is only proven for narrow fluctuations ($v_{n}^{\rm{max}}/v_{n}^{\rm{min}}\lesssim 2.3$) or nearly Gaussian fluctuations Bhalerao et al. (2003).. Furthermore due to the large power relating the two quantities, i.e., $c_{n}\{2k\}\propto(v_{n}\{2k\})^{2k}$, a large change in $c_{n}\{2k\}$ usually leads to a small change in $v_{n}\{2k\}$ for large $k$. A factor of two change in $c_{n}\{8\}$, for example, only leads to 9% change in $v_{n}\{8\}$. Given the large uncertainty (more than $10\%$) in current measurements of $v_{2}\{4\}$, $v_{2}\{6\}$ and $v_{2}\{8\}$ in small systems, they cannot yet place strong constraints on $p(v_{n})$ 555To distinguish between Gaussian and power distributions, for example, a precision of few percents is required for $v_{2}\{4\}$, $v_{2}\{6\}$, and $v_{2}\{8\}$ Yan and Ollitrault (2014)..

We consider a simple case where the data consist of two types of events: The first type has finite correlation signal $v_{0}$ with a probability of $a$, and the second type has zero signal with a probability of $1-a$.

[TABLE]

It can be shown that for certain range of $a$, the $c_{n}\{2k\}$ for $k>1$ always have the “wrong” sign such that $v_{n}\{2k\}$ are undefined. For example for $a=1/3$, $c_{n}\{4\}=\frac{1}{9}v_{0}^{4}$, $c_{n}\{6\}=-\frac{2}{9}v_{0}^{6}$, $c_{n}\{8\}=\frac{71}{9}v_{0}^{8}$. Note that we make no explicit distinction between flow and non-flow, but rather just require the events to have the same underlying p.d.f $\emph{P}(\phi)$. This can be achieved in a simulation where events are generated with fixed multiplicity $M$ and $a=1/3$. In the case of flow, one could generate $N$ events with the same $v_{0}$, but with a random phase event-to-event; these events are then combined with $2N$ events that are generated with zero signal. In the case of non-flow, one could generate $N$ events combined with the same dijet event but rotated randomly in azimuth. The multiplicity of the dijet is chosen such that its effective harmonic signal in the merged event is also $v_{0}$; these events are then combined with $2N$ events generated with zero signal. It is clear that both cases should give the same $c_{n}\{2k\}$ value.

Another important property of cumulants is that they are not additive (see Ref. Jia (2014)) due to the presence of non-linear term such as $\left\langle v_{n}^{2}\right\rangle^{2}$ in$c_{n}\{4\}$, i.e. the cumulants calculated for ensemble A+B $c_{n}\{2k,$A+B$\}$ is not equal to the average of the cumulants calculated separately for ensemble A $c_{n}\{2k,$A$\}$ and ensemble B $c_{n}\{2k,$B$\}$. To see this, consider a pure fluctuation driven scenario of Eq. 2 by setting $v_{n}^{0}=0$: $p(v_{n})\propto v_{n}e^{-v_{n}^{2}/(2\delta_{n}^{2})}$. The events for $p(v_{n})$ are divided into two equal halves, A and B, according to the magnitude of $v_{n}$ (this is possible in reality via event-shape engineering technique). Then we have $c_{n}\{2k\}=0$ for $k>1$, but both $c_{n}\{2k,$A$\}$ and $c_{n}\{2k,$B$\}$ are negative. As we will show later, this is one of the reason why multi-particle cumulants are sensitive to the event class definition.

These examples demonstrate that the values of $c_{n}\{2k\}$ depend only on the probability distribution of flow and non-flow via Eq. 19. Having the “correct sign” for $c_{n}\{2k\}$, i.e. negative for even-$k$ and positive for odd-$k$, in principle, is neither a necessary nor sufficient condition for the correlations being dominated by a global collectivity. On the other hand, it is entirely possible that Gaussian and power distributions could be a reasonable guess of the initial geometry in $pp$ collisions. In such cases, the “correct signs” of $c_{n}\{2k\}$ are still respected, but it is not yet proven they are the only models. ATLAS has shown that $c_{2}\{4\}$ calculated with the standard cumulant method in $pp$ collisions is sensitive to the event class definition Aaboud et al. (2017b). This is because, as shown below, $c_{2}\{4\}$ has a significant non-flow contribution, and the sensitivity of $c_{2}\{4\}$ to event class definition is simply due to the change of the underlying non-flow distribution when the event class is changed and the fact that multi-particle cumulants are not additive (discussed above). The main source of non-flow in $pp$ collisions is dijets, which are responsible for the large non-flow peak at \mbox{\Delta\phi}\sim\pi observed in two-particle correlations Aad et al. (2016). In the following we introduce a subevent cumulant method that further suppresses the non-flow correlations and as well as isolates genuine long-range collectivity that correlates particles in several distinct rapidity ranges.

V Subevent cumulants

In this paper we focus mainly on four-particle cumulants, but the generalization to higher-order is straightforward (formulae for some six- and eight-particle cumulants are provided in Appendixes B and C). The basic idea for subevents cumulants is shown in Figure 1a. The event is divided into three non-overlapping rapidity ranges, labeled as $a$,$b$ and $c$. Four-particle correlators are constructed by choosing two particles from subevent $a$ and one particle each from subevents $b$ and $c$. Dijets contributions, the main source for non-flow, are suppressed, since they can only produce particles in two subevents, except for the small chance of one jet falling on the boundary between two subevents. To increase the statistical power of the events, we include quadruplets obtained by permutations of $a$, $b$ and $c$. For comparison, we also consider cumulants based on two subevents as shown in Figure 1b. In this case, two particles each are chosen from $b$ and $c$, which effectively suppresses contributions from intra-jet correlations, but not inter-jet correlations when the two jets from a dijet land in different subevents.

V.1 Two-subevent cumulants

Two different types of four-particle correlators can be constructed from two subevents, expressed as $e^{in(\phi_{1}^{a}+\phi_{2}^{a}-\phi_{3}^{b}-\phi_{4}^{b})}$ or $e^{in(\phi_{1}^{a}+\phi_{2}^{b}-\phi_{3}^{a}-\phi_{4}^{b})}$. In the first type, particles in subevent $a$ are conjugated with particles in subevent $b$, and are denoted as “$(aa,b^{*}b^{*})$”, with the comma separating the subevents. We will use a short-handed notation “$a,a|b,b$” or simply “$2a|2b$”. In the second type, the two particles in each subevent are conjugated with each other, and are denoted as “$(aa^{*},bb^{*})$” or simply “$a,b|a,b$”.

Multi-particle correlators based on two subevents, each having $M_{a}$ and $M_{b}$ particles, can be derived by keeping track of the duplicates terms Bilandzic et al. (2011, 2014) similar to those in Eq. 12. The four-particle correlator of the first type and the related two-particle correlator, as well as their event-by-event weights have the following expressions in terms of flow vectors:

[TABLE]

For simplicity we have dropped the label “Re”, but it should be understood that only the real component of the $\left\langle 2k\right\rangle$ is kept. The superscripts $o=a,b$ denote the two subevents, and ${\bm{Q}}_{n,o}$ and ${\bm{q}}_{n,o}$ represent the flow vector and normalized flow vectors, respectively:

[TABLE]

The last three terms in the numerator of Eq. 22 account for three types of quadruplets containing duplicated particles, $\phi_{1}=\phi_{2}$, $\phi_{3}=\phi_{4}$, and $\phi_{1}=\phi_{2}$ & $\phi_{3}=\phi_{4}$, respectively. There are a total of $M_{a}^{2}M_{b}^{2}-M_{a}(M_{a}-1)\;M_{b}(M_{b}-1)$ such quadruplets.

Following the diagrammatic approach and the notation used in Ref. Borghini et al. (2001a), multi-particle correlators, after averaging over many events, can be expanded into contributions from flow and non-flow, assuming they are uncorrelated:

[TABLE]

where the “connected” terms denoted by $\left\langle\right\rangle_{c}$ represent genuine multi-particle correlations from non-flow, and we have kept the complex notation for flow vectors, e.g. $\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,b}^{*}\right\rangle=\left\langle v_{n,a}v_{n,b}\cos n(\Phi_{n}^{a}-\Phi_{n}^{b})\right\rangle$, to keep track of the twist of event-plane angles and fluctuations of flow magnitudes. The last two terms in the four-particle correlator with underline are non-flow terms involving only two particles, which can be removed if the four-particle cumulant is defined as:

[TABLE]

The remaining non-flow contribution, denoted by the first term, has to connect two particles in subevent $a$ with the other two particles in subevent $b$.

The second type of four-particle correlation has the following expressions:

[TABLE]

In this case, one can show that the four particle cumulants should be defined as:

[TABLE]

The second type of four-particle cumulant has a larger contribution from non-flow, since a single jet could correlate two particles in one subevent, as indicated by \left\llangle 2\right\rrangle_{a|a} and \left\llangle 2\right\rrangle_{b|b} that need to be subtracted from \left\llangle 4\right\rrangle_{a,b|a,b}.

V.2 Three-subevent cumulants

There are also two types of four-particle correlators constructed from three subevents, which can be expressed as either $e^{in(\phi_{1}^{a}+\phi_{2}^{a}-\phi_{3}^{b}-\phi_{4}^{c})}$ or $e^{in(\phi_{1}^{a}+\phi_{2}^{b}-\phi_{3}^{a}-\phi_{4}^{c})}$. In the first type, particles in $a$ are conjugated with particles in $b$ and $c$, and are denoted as “$(b^{*},aa,c^{*})$” or “$a,a|b,c$”. In the second correlator, the two particles from subevent $a$ are conjugated with each other, and is denoted as “$(b,aa^{*},c^{*})$” or “$a,b|a,c$”.

The first type of four-particle correlation can be written as:

[TABLE]

It can be expanded into contributions from flow and non-flow after averaging over many events:

[TABLE]

The second line contains non-flow involving only two particles, which can be removed by the following definition of the four-particle cumulant:

[TABLE]

The second type of four-particle correlation has the following expression:

[TABLE]

The corresponding four-particle cumulant should be defined as:

[TABLE]

This second type of four-particle cumulant has more contribution from non-flow, since a single jet correlates two particles in the same subevent $a$, as indicated by \left\llangle 2\right\rrangle_{a|a} that need to be subtracted from \left\llangle 4\right\rrangle_{a,b|a,c}.

V.3 Statistical power of different methods

The statistical power of the cumulant method is controlled by the number of unique combinations in each event. When the event multiplicity $M$ is large, the number of unique quadruplets for $\left\langle 4\right\rangle$ is expected to be $W\approx M^{4}/4!=M^{4}/24$ for standard cumulant, $W\approx(M/2)^{2}/2!(M/2)^{2}/2!=M^{4}/64$ for the two-subevent method and $W\approx 3\times(M/3)^{2}/2!(M/3)(M/3)=M^{4}/54$ for the three-subevent method, where the factor 3 takes into account the three permutations of $a$, $b$ and $c$. Therefore in comparison to the standard cumulant method, the two-subevent and three-subevent methods loose a factor of 2.7 and 2.3 statistics, respectively. It is also interesting to note that the three-subevent method has more statistical power than the two-subevent method. On the other hand, the number of unique pairs for $\left\langle 2\right\rangle$ is $M^{2}/2$, $(M/2)^{2}=M^{2}/4$ and $(M/3)^{2}=M^{2}/9$ for standard, two-subevent (see Eq. 27) and three-subevent (see Eq. 34) methods, respectively. In the end, the total statistical uncertainty on $c_{n}\{4\}$ is expected to be smallest for the standard method, but it is non-trivial to decide whether the two-subevent method or the three-subevent method has better statistical precision for $c_{n}\{4\}$.

VI Simulation setup

Particle production in $pp$ collisions is often described by QCD-inspired models implemented in Monte Carlo (MC) event generators such as PYTHIA Sjöstrand, Torbjörn and Mrenna, Stephen and Skands, Peter Z. (2008). The PYTHIA model contains significant non-flow correlations from jets, dijets and resonance decays but has no genuine long-range ridge correlations. In this paper, 200 million $pp$ collisions at $\sqrt{s}=13$ TeV are generated with PYTHIA8. Multi-particle cumulants based on standard method as well as subevent methods are calculated to quantify how they are biased by non-flow correlations as a function of charged particle multiplicity. Furthermore, flow signal is added to the generated event using a flow afterburner Masera et al. (2009), and the performance for recovering the input flow signal is studied.

In a typical azimuthal correlation analysis, multi-particle correlators $\left\langle 2k\right\rangle$ are calculated for particles passing selection criteria $X_{1}$, and are then averaged over many events with the same number of charged particles $N_{\mathrm{ch}}$ passing another selection criteria $X_{2}$ to obtain \left\llangle 2k\right\rrangle and then $c_{n}\{2k\}$ Aaboud et al. (2017b). In order to present the cumulants along a common $x$-axis, the $N_{\mathrm{ch}}$ based on criteria $X_{2}$ is often mapped to a common event activity measure, typically $N_{\mathrm{ch}}$ passing yet another criteria $X_{3}$. If non-flow is sufficiently suppressed, $c_{n}\{2k\}$ should measure the collectivity of the entire event and therefore should not depend on the intermediate criteria $X_{2}$ 666This was observed for $v_{n}$ obtained 2PC analysis in $p$+Pb collisions with peripheral subtraction Aad et al. (2014a): the $v_{n}$ obtained in bins of forward total energy $E_{\mathrm{T}}^{{\scriptscriptstyle\mathrm{Pb}}}$, but then mapped to $N_{\mathrm{ch}}$ based on correlation between $E_{\mathrm{T}}^{{\scriptscriptstyle\mathrm{Pb}}}$ and $N_{\mathrm{ch}}$, is found to agree well with $v_{n}$ obtained by binning events directly in $N_{\mathrm{ch}}$.. However, the ATLAS Collaboration observed that $c_{2}\{4\}$ calculated in both PYTHIA and with $pp$ data depends sensitively on criteria $X_{2}$ Aaboud et al. (2017b), and sensitivity was thought to arise from relative multiplicity fluctuation between $X_{2}$ and $X_{3}$. The real reason, as we shall argue below, is due to change in event-by-event non-flow distribution $p(s_{n})$ 777As we discussed before, the statistical component ${\bm{s}}_{\mathrm{n}}^{\rm{stat}}$ associated with multiplicity fluctuation does not influence the measured $v_{n}$.. We show that the three-subevent method nearly completely removes non-flow contributions, such that the resulting $c_{2}\{4\}$ is independent of $X_{2}$.

Table 1 summarizes the types of cumulants and charged particle selection criteria studied in PYTHIA8. For each method, cumulants are calculated in two $p_{\mathrm{T}}$ ranges ($X_{1}$) for event classes defined by the number of charged particles in four $p_{\mathrm{T}}$ ranges ($X_{2}$), denoted by $N_{\textrm{ch}}^{\textrm{sel}}$. For each combination, the $c_{n}\{2k\}$ are calculated for events with a fixed $N_{\textrm{ch}}^{\textrm{sel}}$ multiplicity, the results for one-particle-width bins in $N_{\textrm{ch}}^{\textrm{sel}}$ are then combined to broader $N_{\textrm{ch}}^{\textrm{sel}}$ intervals. When presenting the final results, the \left\langle\mbox{N_{\textrm{ch}}^{\textrm{sel}}}\right\rangle values for given $N_{\textrm{ch}}^{\textrm{sel}}$ interval are mapped to the average number of charged particles with $p_{\mathrm{T}}>0.4$ GeV ($X_{3}$), denoted by \left\langle\mbox{N_{\mathrm{ch}}}\right\rangle.

VII Results

Figure 2 shows the four-particle cumulants obtained with the standard method using charged particles in $0.3<p_{\mathrm{T}}<3$ GeV. The $N_{\textrm{ch}}^{\textrm{sel}}$ for calculating \left\llangle 2k\right\rrangle is defined in $0.3<p_{\mathrm{T}}<3$ GeV, $p_{\mathrm{T}}>0.2$ GeV, $p_{\mathrm{T}}>0.4$ GeV, or $p_{\mathrm{T}}>0.6$ GeV, which are then mapped to \left\langle\mbox{N_{\mathrm{ch}}}\right\rangle. The $c_{2}\{4\}$ values are found to differ significantly between the four choices of $N_{\textrm{ch}}^{\textrm{sel}}$, therefore confirming the observation made by the ATLAS Collaboration Aaboud et al. (2017b). The red and blue dashed lines indicate expected $c_{n}\{4\}$ values corresponding to 4% and 6% flow signal. Clearly non-flow contributions in the standard method are too large for a meaningful measurement of $v_{2}$. In fact, when $N_{\textrm{ch}}^{\textrm{sel}}$ is defined in the $p_{\mathrm{T}}$ range of $p_{\mathrm{T}}>0.6$ GeV, the sign of $c_{2}\{4\}$ is negative at small \left\langle\mbox{N_{\mathrm{ch}}}\right\rangle values. The non-flow contribution to $v_{3}$ is significantly smaller as dijets are not expected to contribute significantly to $v_{3}$; a $v_{3}\{4\}$ signal larger than 4% is measurable in the standard method.

Figure 3 shows the $c_{n}\{4\}$ obtained with the standard method using charged particles in a higher transverse momentum range $0.5<p_{\mathrm{T}}<5$ GeV. The $c_{2}\{4\}$ is much larger in all cases, which is consistent with a larger non-flow contribution expected from dijets. Figures 4 and 5 show $c_{n}\{4\}$ obtained from the two-subevent and three-subevent methods (the first type), respectively. The two-subevent method greatly reduces the values of $c_{n}\{4\}$ compared to the standard method, and the dependence on the choice of $N_{\textrm{ch}}^{\textrm{sel}}$ is also much smaller. This behavior is expected since the contributions from single jets or resonance decays are suppressed. Nevertheless, the magnitude of the residual non-flow signal is still comparable to a 4%–6% flow signal. In contrast, the three-subevent method almost completely suppresses the non-flow. For both $c_{2}\{4\}$ and $c_{3}\{4\}$, the magnitudes of the residual non-flow are less than 4% of flow signal. This means that the three-subevent method should be sensitive to genuine long-range flow signal as small as 4%.

Figure 6 shows a direct comparison between the standard method and the two- and three-subevent methods for $v_{2}\{4\}$ in two $p_{\mathrm{T}}$ ranges. The three-subevent method has the best performance in suppressing the non-flow effects.

To quantify the performance of the three methods for recovering the underlying flow signal, a flow afterburner Masera et al. (2009) is used to add a constant $v_{2}$ or $v_{3}$ signal to the generated PYTHIA events. Figure 7 shows the calculated $c_{2}\{4\}$ with 4% or 6% $v_{2}$ imposed on the generated events. In the case 4% input flow, only the three-subevent method can recover the input. In the case of 6% input flow, the two-subevent method can also recover the input flow for \left\langle\mbox{N_{\mathrm{ch}}}\right\rangle>80.

In the cumulant analysis, it has been argued that $c_{n}\{2k\}$ should be calculated for each fixed $N_{\textrm{ch}}^{\textrm{sel}}$ bin to minimize the event-by-event variation of particle multiplicity, then should be averaged over broader $N_{\textrm{ch}}^{\textrm{sel}}$ intervals Bilandzic et al. (2011); Chatrchyan et al. (2013b). Figure 8 shows $c_{2}\{4\}$ for charged particles in $0.3<p_{\mathrm{T}}<3$ GeV for events with 40\leq\mbox{N_{\mathrm{ch}}}<80, where $N_{\mathrm{ch}}$ is the number of charged particles also in $0.3<p_{\mathrm{T}}<3$ GeV. The $c_{2}\{4\}$ values are obtained with $N_{\mathrm{ch}}$ bin widths of 1, 5, 10, 20 and 40, respectively. The $c_{2}\{4\}$ values for each case are then averaged over the 40\leq\mbox{N_{\mathrm{ch}}}<80 interval to give a single $c_{2}\{4\}$ value. The difference of $c_{2}\{4\}$ for a given bin-width from $c_{2}\{4\}$ for unit-bin-width is then plotted. Clearly, the standard method is much more sensitive to the bin-width than the two-subevent and three-subevent methods. This sensitivity is due to the fact that the standard method still has significant non-flow, whose event-by-event fluctuation also has significant non-Gaussianity. Therefore the residual non-flow has significant dependence on the bin-width in the standard method 888Even if $c_{2}\{4\}$ is calculated in unit $N_{\mathrm{ch}}$ bin, there is still significant residual non-flow (as shown by Figs 3 and 4).. On the other hand, since non-flow contributions are significantly further suppressed in the subevent methods, the nature of the non-flow fluctuations no longer matter much for the subevent methods.

Finally, we also studied the $c_{n}\{4\}$ using a four subevent method, where $|\eta|<2.5$ range is divided into four equal $\eta$ ranges $a,b,c,d$ and the cumulant is calculated as c_{n}^{a,b|c,d}\{4\}\equiv\left\llangle 4\right\rrangle_{a,b|c,d}-2\left\llangle 2\right\rrangle_{a|c}\left\llangle 2\right\rrangle_{b|d}. Despite the significantly larger statistical uncertainties, the values of $c_{n}\{4\}$ converge to that from the three-subevent method 999Another test is also performed, where we increase the $\eta$ range to $|\eta|<5$ to include more particles, in that case, the four-subevent method agrees with the three-subevent method perfectly..

VIII Discussion and outlook

Ever since its discovery in small collision systems, the ridge phenomenon has been an area of insensitive study in the high-energy and heavy-ion physics community. The current debate is centered on the underlying multi-particle dynamics and collective nature of the long-range correlations, in particular whether the ridge is a global property of the event, and what is its connection to the collective flow in large collision systems Dusling et al. (2016); Schlichting and Tribedy (2016). The main challenge is to disentangle correlations involving all particles (flow or collectivity) from those involving a subset of the particles in restricted $\eta$ and $\phi$ space (non-flow). The two-particle correlation method based on peripheral subtraction or the template fit method Aad et al. (2014a, 2016); Aaboud et al. (2017a) can successfully remove the non-flow, but it does not address the multi-particle nature of the ridge. Results from multi-particle azimuthal cumulants, namely negative $c_{n}\{4\}$ and $c_{n}\{8\}$, positive $c_{n}\{6\}$, and the relation $v_{n}\{4\}\approx v_{n}\{6\}\approx v_{n}\{8\}$, have been used as a “definition” of collectivity Khachatryan et al. (2015a, 2017). However this is true only if the flow fluctuations are relatively narrow or close to Gaussian and non-flow contributions to $c_{2}\{2k\}$ are small. A recent study from ATLAS Aaboud et al. (2017b) shows that in $pp$ collisions, the standard cumulant method is contaminated by non-flow correlations over the entire measured $N_{\mathrm{ch}}$ range, which further complicates the relationship between different $v_{n}\{2k\}$s.

In this context, our paper addresses the issue of collectivity in small collision systems from two fronts. First, we clarify the statistical nature of multi-particle cumulants, and show that azimuthal cumulants are affected not only by the event-averaged flow and non-flow, but more importantly by the event-by-event fluctuations of flow and/or non-flow. We then propose alternative cumulant methods based on two or three subevents in non-overlapping $\eta$ ranges, and demonstrate that these methods can further suppress non-flow correlations.

When cumulant methods were first proposed in 2000 Borghini et al. (2001a), the importance of non-gaussian flow fluctuations was not recognized. If flow fluctuations or flow distribution $p(v_{n})$ are close to Gaussian, all higher-order cumulants should give the same answer $v_{n}\{4\}=v_{n}\{6\}...$, and the difference between $v_{n}\{2\}$ and $v_{n}\{4\}$ reflects the width of $p(v_{n})$ Voloshin et al. (2008). However when $p(v_{n})$ strongly deviates from Gaussian, $v_{n}\{2k\}$ (and also the Lee-Yang zero method Bhalerao et al. (2003)) may be mathematically undefined and one should instead rely on the original cumulant quantity $c_{n}\{2k\}$ Jia and Radhakrishnan (2015), which is always well defined. For example consider

[TABLE]

which is mathematically undefined only if $2\left\langle v_{n}^{2}\right\rangle^{2}<\left\langle v_{n}^{4}\right\rangle$. If this relation is violated due to the underlying $p(v_{n})$, the perceived hallmark relation for collectivity, $v_{n}\{4\}\approx v_{n}\{6\}\approx v_{n}\{8\}$, is no longer meaningful.

The contribution of non-flow to cumulants can be discussed following a similar probabilistic approach. The contributions of flow and non-flow to cumulants are additive:

[TABLE]

The sign of $c_{n}\{2k,\mbox{non-flow}\}$ can be either positive or negative depending on the shape of event-by-event fluctuations of non-flow $p(s_{n})$. The results from ATLAS and our follow-up study from PYTHIA8 suggest that the value and sign of $c_{n}\{4\}$ in $pp$ collisions are driven mainly by the non-flow component in the standard cumulant method. The values of $c_{n}\{2k,\mbox{non-flow}\}$ are found to be very sensitive to the particle selection criteria $p_{\mathrm{T}}>p_{\mathrm{T}}^{\rm{min}}$ used to define event for averaging: $c_{n}\{4\}$ is positive when $p_{\mathrm{T}}^{\rm{min}}=0.2$ and 0.4 GeV, but is negative for $p_{\mathrm{T}}^{\rm{min}}=0.6$ GeV in low multiplicity events. This suggests that the $p(s_{n})$ in PYTHIA8 is very sensitive to $p_{\mathrm{T}}$ of the particles, as expected from jets and dijets. It is also interesting to point out that by choosing the $p_{\mathrm{T}}$ range used to define $N_{\textrm{ch}}^{\textrm{sel}}$, one may be able to “dial” the $c_{n}\{2k,\mbox{non-flow}\}$ to be very close to zero, and thereby expose the genuine flow signal. Unfortunately one can not do this in a model independent way, therefore the observed negative $c_{2}\{4\}$ by CMS Collaboration Khachatryan et al. (2017) can be interpreted as a negative $c_{n}\{2k,\mbox{flow}\}$ from collectivity, if and only if the non-flow contribution $c_{n}\{2k,\mbox{non-flow}\}$ can be shown to be small 101010Indeed, the $c_{n}\{2k,\mbox{non-flow}\}$ from PYTHIA8 model is large and can be negative depending on the choice of $p_{\mathrm{T}}$ range. On the other hand, no model that include only non-flow can describe the observation in the $pp$ data..

Motivated by the discussions above, we have developed a cumulant method based on two- or three $\eta$-separated subevents. The two-subevent method suppresses the non-flow from single jet fragmentation, while the three-subevent method also suppresses correlations between the two jets in a dijet. The performance of the two methods is quantified using PYTHIA8 simulation with only non-flow, as well as with physical flow added with an afterburner. The three-subevent method is shown to give consistent results between different criteria used to define event classes for averaging, while the two-subevent method still has some sensitivity. The three-subevent method is able to recover flow signal as small as 4% for $v_{2}$ and $v_{3}$ for events with \left\langle\mbox{N_{\mathrm{ch}}}\right\rangle as low as 30. The subevent method also offers a more robust definition of collectivity in $pp$ or $p$+A collisions, i.e. the existence of long-range ridge signal that correlates between three or more distinct $\eta$ ranges. The method also provides a strong test of the gluon saturation models used to describe $c_{2}\{4\}$ obtained with the standard cumulant method, where it was argued that the sign change of $c_{2}\{4\}$ could be due to non-Gaussian correlations of the domains of strong QCD fields in the initial state Dumitru et al. (2015, 2015); Lappi et al. (2016). It would be interesting to see wether or not such color domains can simultaneously correlate three or more pseudorapidity ranges and contribute to subevent cumulants.

Looking into the future, there are several improvements or applications of the subevent cumulant method in both small as well as large collision systems:

• •

Rapidity gap between subevents: this can further suppress short-range sources falling on the boundary between two subevents. Furthermore, by requiring the same rapidity gap as in 2PC analysis, $v_{n}\{4\}$ can be compared directly with $v_{n}\{2\}$ obtained via a peripheral subtraction procedure.

• •

Four or higher number of subevents: this can further reduce non-flow and prove that the azimuthal collectivity indeed exists between particles in four or more distinct $\eta$ ranges. But this also reduce significantly the statistics.

• •

Subevents selected in more exclusive way: for example one could reject events containing jets or dijets. By doing so, the large tail in the non-flow distribution $p(s_{n})$ should be suppressed and $c_{n}\{2k,\mbox{non-flow}\}$ is reduced. One could also have one subevent containing all the reconstructed jets (several $\eta$ slices each one containing one or more jets), and the rest of the $\eta$ range is defined as the second subevent. This way the non-flow contributions from jets and dijets are contained in one subevent.

• •

Differential flow: the procedure for calculating $v_{n}\{4\}(p_{T},\eta)$ is straightforward in the three-subevent method. One just needs to restrict the particle in one subevent to a certain $p_{\mathrm{T}}$ or $\eta$ range.

• •

Flow longitudinal dynamics: for example, ignoring the non-flow terms, subevent cumulants such as Eqs. 27, 30, 33 or 34, are sensitive to decorrelation effects between different $\eta$ ranges:

[TABLE]

First, $\left\langle{\bm{v}}_{n,a}{\bm{v}}_{n,b}^{*}\right\rangle=\left\langle v_{n,a}v_{n,b}\cos n(\Phi_{n}^{a}-\Phi_{n}^{b})\right\rangle$ can be calculated from two-particle correlations (see Eq. 21), supplemented with peripheral subtraction and template fit Aad et al. (2014a, 2016); Aaboud et al. (2017a) to reject the non-flow contributions. $\left\langle v^{2}_{n,a}\right\rangle\left\langle v^{2}_{n,b}\right\rangle$ and $\left\langle{\bm{v}}^{2}_{n,a}{\bm{v}}^{*2}_{n,b}\right\rangle$ are then obtained from Eqs. 40 and 41. The study of flow fluctuations based on moments were discussed before in Ref. Bhalerao et al. (2015), but subevent cumulants allow us to directly measure these moments including longitudinal flow decorrelation effects, without relying on a reference detector as was usually done in experiments Khachatryan et al. (2015b).

• •

Generalization of subevent cumulants to symmetric cumulants or event-plane correlations Bilandzic et al. (2014): Some examples are discussed in Appendixes E and F; they provide access to the longitudinal dynamics of the initial state (e.g. correlation between $v_{2}$ and $v_{3}$ ) and final state mode-mixing (e.g. correlation between $v_{2}$ and $v_{4}$).

We thank A. Bilandzic for clarification on the framework of direct cumulants. We acknowledge Sooraj Radhakrishnan for help in the early preparation of this paper. We appreciate valuable comments and fruitful discussions with R. Lacey. This research is supported by the NSF Grant No. PHY-1613294.

Appendix

VIII.1 Performance of 2nd type of two- and three-subevent cumulants

We also calculate cumulants based on the second type of two-subevent and three-subevent methods mentioned in Section V. The results, as shown in Figs. 9 and 10, are found to be very sensitive to the particle $p_{\mathrm{T}}$ selection criteria used to define events for averaging, suggesting that these correlators have more contributions from non-flow. The larger non-flow in these correlators can be traced to the fact that the two particles from the same subevent are conjugated with each other as shown in Eqs. 30 and 36 (terms containing \left\llangle 2\right\rrangle_{a|a} or \left\llangle 2\right\rrangle_{b|b}).

VIII.2 Six- and eight-particle cumulants based on two subevents

The cumulants based on two subevents can be generated for more than four particles:

[TABLE]

The non-flow contribution to \left\llangle 6\right\rrangle_{3a|3b} is straightforward to derive, but the expression is quite lengthy, so we will not show it here. But one can show that only flow and intrinsic six-particle non-flow contribute to this definition of cumulants:

[TABLE]

Similarly the eight-particle cumulant can be derived as,

[TABLE]

The procedures for deriving $\left\langle 6\right\rangle_{3a|3b}$ and $\left\langle 8\right\rangle_{4a|4b}$ are provided in Appendix H. The correctness of these formulas has been cross-checked by comparing the results with exact calculation using nested-loops.

By combining information from subevent cumulants of different orders (Eqs.39–41), one can systematically calculate higher-order flow moments including the longitudinal dynamics, such as $\left\langle{\bm{v}}^{3}_{n,a}{\bm{v}}^{*3}_{n,b}\right\rangle$ and $\left\langle{\bm{v}}^{4}_{n,a}{\bm{v}}^{*4}_{n,b}\right\rangle$, with minimum contribution from non-flow effects.

VIII.3 Six-particle cumulants based on three subevents

We shall consider only six-particle correlations where two particles are chosen from each subevent, which leads to two types of correlators “$(bb^{*},aa^{*},cc^{*})$” and “$(bb,aa^{*},c^{*}c^{*})$”. In the first type, the two particles in each subevent are conjugated with each other. In the second type, particles in subevent $b$ are conjugated with particles in subevent $c$, while the two particles in subevent $a$ are conjugated with each other. These two types can also be denoted as “$a,b,c|a,b,c$” and “$a,b,b|a,c,c$”, respectively.

The expression for cumulants can be derived using the diagrammatic approach in Ref. Borghini et al. (2001a), except that we also need to keep track of the subevent indexes explicitly. The results for six-particle cumulants are:

[TABLE]

[TABLE]

The $c_{n}^{a,b,b|a,c,c}\{6\}$ has a smaller contribution from non-flow.

VIII.4 Cumulants based on four subevents

Four-particle correlations, where one particle is taken from each subevent can be written as Adler et al. (2002):

[TABLE]

One interesting six-particle correlation based on four subevents is $(cc,a,b*,d^{*}d^{*})$ or “$a,c,c|b,d,d$”:

[TABLE]

To minimize the residual non-flow, one can put subevents $a$ and $b$ in between subevents $c$ and $d$, and also require a small rapidity gap between $a$ and $b$.

VIII.5 Symmetric cumulants

The subevent cumulants can be generalized to symmetric cumulants first proposed in Ref. Bilandzic et al. (2014). There are three variants of two-subevent symmetric cumulants:

[TABLE]

where $n>m$ and the notation combines harmonic number and subevent index. These three types can also be schematically labeled as $(n_{a}m_{a},n_{b}^{*}m_{b}^{*})$, $(n_{a}m_{a}*,m_{b}n_{b}^{*})$, $(n_{a}n_{a}^{*},m_{b}m_{b}^{*})$, respectively. It is interesting to note that contributions of flow are expected to be nearly identical for Eq. 48 and Eq. 49, and they are sensitive to both flow magnitude fluctuations and event-plane twist effects. On the other hand, the correlator in Eq. 50 is only sensitive to the fluctuation of the flow magnitudes.

Since terms such as $\left\langle{\bm{v}}_{n,a}{\bm{v}}^{*}_{n,b}\right\rangle$, $\left\langle{\bm{v}}_{m,a}{\bm{v}}^{*}_{m,b}\right\rangle$, $\left\langle v_{n,a}^{2}\right\rangle$, $\left\langle v_{m,b}^{2}\right\rangle$ can be calculated from subevent cumulants involving harmonic flow of the same order via Eqs. 39–41, one can also calculate directly $\left\langle v_{n,a}^{2}v_{m,b}^{2}\right\rangle$ and $\left\langle{\bm{v}}_{n,a}{\bm{v}}^{*}_{n,b}{\bm{v}}_{m,a}{\bm{v}}^{*}_{m,b}\right\rangle$ by combining cumulants and symmetric cumulants, and therefore gain knowledge about the decorrelations between harmonics of different order.

Similarly, one can write down the formulas for symmetric cumulants from three subevents. There are four possibilities, $(n_{b}^{*},n_{a}m_{a},m_{c}^{*})$, $(n_{b}^{*},n_{a}m_{a}^{*},m_{c})$, $(m_{b},n_{a}n_{a}^{*},m_{c}^{*})$ and $(n_{b},m_{a}m_{a}^{*},n_{c}^{*})$, with the following formula:

[TABLE]

VIII.6 Asymmetric cumulants (event-plane correlators)

We consider two interesting cases, the correlation between $v_{2}$ and $v_{4}$

[TABLE]

and correlation between $v_{2}$, $v_{3}$ and $v_{5}$

[TABLE]

Depending on how the particles are arranged, removal of the duplicate terms (or autocorrelation) in subevent $a$ may play an important role. If the all ${\bm{q}}$ are placed in subevent $a$ and all ${\bm{q}}^{*}$ are placed in subevent $b$ (Eqs. 52 and 54), the autocorrelation effect is on the order of $\tau_{a}v_{n+m}/(v_{n}v_{m})\sim 1/M_{a}$ as first derived by Bhalerao, Ollitrault and Pal. Bhalerao et al. (2013), which is small when multiplicity in the subevent is large (e.g. in A+A collisions). On the other hand, when ${\bm{q}}$ and ${\bm{q}}^{*}$ appear in the same subevent, the autocorrelation contribution is on the order of $\tau_{a}v_{|n-m|}/(v_{n}v_{m})$, which could be large since $v_{|n-m|}\gg v_{n}v_{m}$. Therefore explicit removal of duplicate terms is important if the multiplicity in the subevent is small, for example if the $\eta$ range is narrow or in a small collision system, or if ${\bm{q}}$ and ${\bm{q}}^{*}$ appear in the same subevent.

Various combinations of subevents for the same correlators also contain different information on the non-linear mode-mixing effects and longitudinal decorrelation effects. It is well-established that the $v_{4}$ and $v_{5}$ in A+A collisions contain a linear contribution associated with initial conditions and a mode-mixing contribution due to non-linear hydrodynamic response Gardim et al. (2012); Teaney and Yan (2012); Aad et al. (2014b, 2015); Qiu and Heinz (2012):

[TABLE]

If the non-flow contribution can be ignored, we have

[TABLE]

These two correlators are sensitive to the second moments and first moment of the decorrelations between the two subevents, respectively. In central A+A collisions, where $v_{4}$ is dominated by the linear component $v_{4\rm{L}}$, the two correlators both measure $\left\langle{\bm{v}}_{2}^{2}{\bm{v}}^{*}_{\rm{4L}}\right\rangle$, which is sensitive to the correlations in the initial eccentricity, i.e.,

[TABLE]

In non-central collisions, where the mode-mixing component is important, one expects that the two correlators deviate from each other, with $\left\langle 3\right\rangle_{2_{a},2_{a}|4_{b}}<\left\langle 3\right\rangle_{2_{a},2_{b}|4_{a}}$. The deviation is expected to increase with the rapidity separation between subevents $a$ and $b$.

Note that the relations in Eqs. 58 can be independently probed by directly measuring two four-particle correlators \left\llangle 4\right\rrangle_{2a|2b}, \left\llangle 4\right\rrangle_{a,a|a,b} and their associated cumulants $c_{2}\{4\}^{2a|2b}$, $c_{2}\{4\}^{a,a|a,b}$ :

[TABLE]

This relation, together with Eqs. 39–41, can be used to calculate $\left\langle{\bm{v}}_{2,a}^{2}{\bm{v}}^{2*}_{2,b}\right\rangle$ and $\left\langle v^{2}_{2}{\bm{v}}_{2,a}{\bm{v}}^{*}_{2,b}\right\rangle$, and eventually the value of $\beta_{2,2}$.

Similarly, the three correlators involving $v_{5}$, $v_{2}$, and $v_{3}$ can be written as:

[TABLE]

The first correlator is sensitive to the decorrelation of both $v_{2}$ and $v_{3}$, while the second and the third correlators are sensitive to the decorrelation of $v_{3}$ only and $v_{2}$ only, respectively. It would be interesting to measure the centrality and $\eta$ dependence of these three correlators to extract the linear component which is directly sensitive to the correlation in the initial eccentricity:

[TABLE]

VIII.7 Particle weights

Often in data analysis, the weight $w_{i}$ needs to be applied for each particle. The general framework for deriving the expression for weighted cumulants is given in Refs. Bilandzic et al. (2011); Bilandzic (2012); Bilandzic et al. (2014); Di Francesco et al. (2017) or via a slightly simplified approach provided in Appendix H. All correlators can be conveniently expressed in terms of weighted per-particle flow vector:

[TABLE]

With this definition, the formulas for subevent correlators can be obtained by simply replacing $\tau^{k-1}{\bm{q}}_{kn}$ by $\tau_{k-1}{\bm{q}}_{kn;k}$ ($k\geq 2$). Here are some examples:

[TABLE]

where we have moved the subevent labels $a$ and $b$ outside the parentheses to compactify the formulas, i.e., $(1-\tau_{1})_{a}\equiv 1-\tau_{1,a}$. The advantage of using per-particle flow vector is that all quantities are smaller than 1, and they can be ordered in terms of $\tau_{k}\sim 1/M^{k}$ or $q_{kn;k}\sim q_{n}^{k}$.

Formulas for two-particle correlators and four-particle correlators based on three-subevents are

[TABLE]

For completeness, we also document the expression for standard cumulant, which again can be readily converted from the formulas given in Refs. Bilandzic et al. (2011); Bilandzic (2012); Bilandzic et al. (2014) or by using the procedure given in Appendix H (Eqs. 85 and 86). The situation is slightly more complicated since there are two or more particles that are complex-conjugated with each other.

[TABLE]

The event-by-event weight for any correlator is also straightforward to write out. It is simply the denominator times the $\sum_{i}w_{i}$ in each subevent raised to the power of number of particles chosen from that subevent. Here are some examples:

[TABLE]

If all particles have the same weights, $\left\langle 4\right\rangle$ and $\left\langle 6\right\rangle$ are simplified to:

[TABLE]

VIII.8 Iterative procedure for deriving the equations for multi-particle correlation

In the original paper on direct cumulants by A. Bilandzic, R. Sneilings and S. Voloshin  Bilandzic et al. (2011), the expression for multi-particle correlator $\left\langle 2k\right\rangle$ is obtained by expanding all lower-order correlators and solving many coupled linear equations to remove combinations with duplicated particle indexes. This was later improved into an iterative procedure in Ref. Bilandzic et al. (2014) for general multi-particle correlations including mixed harmonics and the algorithm for calculating these correlators has been provided Bilandzic et al. . Here we recast that general iterative procedure into a simplified form to handle two specific classes of correlators used in this paper: the correlation involving the whole event used in the standard method and two-subevent correlations.

VIII.8.1 unweighted case

We introduce the following notation for a multi-particle correlator in an event with $M$ particles:

[TABLE]

where, following convention in Ref. Bilandzic et al. (2011), $k*n$ and $-k*n$ are short-hand notations for $n,n,..,n$ and $-n,-n,..,-n$, each with $k$ indexes. With this notation, the correlators in the standard method and two-subevent method (Eq. LABEL:eq:af12) can be expressed as:

[TABLE]

Note that the number of permutations can be written as:

[TABLE]

To derive the expression for $\overline{k}_{k*n}$, we introduce an operator $\otimes$ and a set of expansion rules, which relate the number of permutations, $M(M-1)...(M-k+1)$, to the final expression in terms of flow vector ${\bm{Q}}_{n}$:

[TABLE]

The expansion rules explicitly keep track of the duplications as one goes from $k$-particle to $k+1$-particle correlation. The first few are:

[TABLE]

Therefore it is straightforward to write out the expression for $\overline{k}_{k*n}$,

[TABLE]

To derive the expression for $\overline{2k}_{k*n,-k*n}$ used in the standard method where all particle are chosen from the whole event, one just needs to extend the rules in Eq. 73 to include case where particle are conjugated with each other:

[TABLE]

where we have used superscripts + and - to indicate the sign of $\phi$ angle for the additional particle going from $k$-particle correlation to $k+1$-particle correlation. Some examples are:

[TABLE]

[TABLE]

The expressions for higher-order correlation are rather lengthy, therefore it is helpful to define new operator $\otimes_{2k\rightarrow 2k+2}$ that directly relates the $2k$- and $2k+2$-particle correlations:

[TABLE]

[TABLE]

where “$cc$” stands for complex conjugate.

VIII.8.2 with particle weights

All formula discussions above can be generalized to include particle weights, and we shall express all results using per-particle flow vectors:

[TABLE]

After replacing ${\bm{Q}}$ with ${\bm{q}}$, the two-subevent correlation and corresponding rules can be rewritten as:

[TABLE]

Note that $\overline{k}_{k*n}=\hat{k}_{k*n}\left(\sum_{i}w_{i}\right)^{k}$, and therefore the factor $\left(\sum_{i}w_{i}\right)^{k}$ cancels out between the numerator and denominator of $\left\langle 2k\right\rangle_{ka|kb}$. When particle weight is a constant, $\lim_{n\rightarrow 0}\hat{k}_{k*n}=(1-\tau)(1-2\tau)...(1-(k-1)\tau),\tau=1/M$.

Similarly the expansion rules for $\left\langle 2k\right\rangle_{k*n,-k*n}$ are given by:

[TABLE]

With these, we can rewrite Eqs. 77 and 78 as:

[TABLE]

[TABLE]

The expression for $\lim_{n\rightarrow 0}\widehat{2k}_{k*n,-k*n}$ can be obtained by setting ${\bm{q}}_{pn,-mn}=1$, for example,

[TABLE]

From Eqs .85 and 86, we obtain the formula for $\left\langle 4\right\rangle$ and $\left\langle 6\right\rangle$ shown in Eqs. 66 and 67, respectively.

Finally the event-by-event weights are:

[TABLE]