Non-flow effects in correlation between harmonic flow and transverse momentum in nuclear collisions

A large anti-correlation signal between elliptic flow $v_2$ and average transverse momentum $[p_{\mathrm{T}}]$ was recently measured in small collision systems, consistent with a final-state hydrodynamic response to the initial geometry. This negative $v_2$-$[p_{\mathrm{T}}]$ correlation was predicted to change to positive correlation for events with very small charged particle multiplicity $N_{\mathrm{ch}}$ due to initial-state momentum anisotropies of the gluon saturation effects. However, the role of non-flow correlations is expected to be important in these systems, which is not yet studied. We estimate the non-flow effects in $pp$, $p$Pb and peripheral PbPb collisions using {\tt Pythia} and {\tt Hijing} models, and compare them with the experimental data. We show that the non-flow effects are largely suppressed using the rapidity-separated subevent cumulant method (details of the cumulant framework are also provided). The magnitude of the residual non-flow is much less than the experimental observation in the higher $N_{\mathrm{ch}}$ region, supporting the final-state response interpretation. In the very low $N_{\mathrm{ch}}$ region, however, the sign and magnitude of the residual non-flow depend on the model details. Therefore, it is unclear at this moment whether the sign change of $v_2$-$[p_{\mathrm{T}}]$ can serve as evidence for initial state momentum anisotropies predicted by the gluon saturation.


I. INTRODUCTION
In high-energy hadronic collisions, particle correlations are an important tool to study the multi-parton dynamics of QCD in the strongly coupled non-perturbative regime [1]. Measurements of azimuthal correlations in small collision systems, such as pp and p+A collisions [2-6], have revealed a strong harmonic modulation of particle densities dN dφ ∝ 1 + 2 ∑ ∞ n=1 v n cos n(φ − Φ n ). Measurement of v n and their event-by-event fluctuations have been performed as a function of charged particle multiplicity N ch in pp and p+A collisions. It is found that the azimuthal correlations involve all particles over a wide pseudorapidity range. A key question is whether this multi-particle collectivity reflects initial momentum correlation from gluon saturation effects (ISM) [7], or a final-state hydrodynamic response to the initial transverse collision geometry (FSM) [8].
Recently, the correlation between v n and [p T ], the average transverse momentum of particles in each event, was proposed to be a sensitive observable to distinguish between the initial-state and final-state effects [9]. The lowest order of such correlation is characterized by the covariance cov(v 2 n , [p T ]) ≡ ⟨v 2 n [p T ]⟩ − ⟨v 2 n ⟩ ⟨[p T ]⟩ [10] with the average carried over events, which have been measured at the LHC [11,12]. In the final-state dominated scenario, the flow harmonics are diven by the initial spatial eccentricity ε n , v n ∝ ε n , while the [p T ] is related to the transverse size of the overlap region: events with similar total energy but smaller transverse size in the initial state are expected to have a stronger radial expansion and therefore larger [p T ] [13]. Hydrodynamic model calculations with negligible initial transverse momentum predict a positive cov(v 2 n , [p T ]) at large N ch which changes to negative cov(v 2 n , [p T ]) towards small N ch region [14][15][16], whereas at small enough Nch, initial momentum anisotropy can, in fact, dominate. In a gluon saturation picture, these correlations are expected to give a positive contribution to cov(v 2 n , [p T ]) [9]. Therefore, the N ch dependence of cov(v 2 2 , [p T ]), after considering both initial and final-state effects, is predicted to exhibit a double sign change as a function N ch . The experimental observation of such sign change was further argued to provide a strong evidence for the gluon saturation physics [9].
On the other hand, momentum correlations could also arise from "non-flow" effects from resonance-decays, jets and dijets [17]. Such non-flow correlations usually involve a few particles from one or two localized pseudorapidity regions, in contrast to the initial momentum correlation from gluon saturation, which spans continuously over a large rapidity range similar to hydrodynamic flow. The non-flow effects are often suppressed by correlating particles from two or more subevents separated in pseudorapidity. This so-called subevent cumulant method [18] has been validated for several multi-particle correlators involving flow harmonics of same or different orders [18][19][20], such as four-particle cumulants c n {4} = ⟨v 4 n ⟩ − 2 ⟨v 2 n ⟩ 2 , four-particle symmetric cumulants ⟨v 2 n v 2 m ⟩ − ⟨v 2 n ⟩ ⟨v 2 m ⟩ and three particle asymmetric cumulants ⟨v n v m v n+m cos (nΦ n + mΦ m − (n + m)Φ n+m )⟩. It is found that results from the standard cumulant method are contaminated by non-flow correlations in pp, pA and peripheral AA collisions, while they are largely suppressed in the subevent method that requires three or more subevents [21,21,22]. Since covariance cov(v 2 n , [p T ]) is a threeparticle correlator, it can be measured with two subevent or three-subevent methods, which suppress the non-flow while keeping the genuine long range multi-particle correlations associated with ISM and FSM.
In this paper, we study the influence of the non-flow correlations to covariance cov(v 2 n , [p T ]) in pp, pPb and PbPb collisions using Pythia8 [23] A2 tune and Hijing v1.37 [24] models in the standard and subevent methods. We find that the non-flow correlations give a positive contribution to cov(v 2 2 , [p T ]), which are strongly suppressed in the three-subevent method, but not completely eliminated. The sign and magnitude of the residual non-flow are model dependent. Therefore, the mere observation of change of cov(v 2 2 , [p T ]) from negative to positive towards low N ch in the experimental results may not serve as evidence for the presence of gluon saturation.

II. METHODOLOGY AND MODEL SETUP
The covariance cov(v 2 n , [p T ]) is a three-particle correlator, which is obtained by averaging over unique triplets in each event, and then over all events in an event class [10,14]: where the indices i, j and k loop over distinct charged particles to account for all unique triplets, the particle weight w i is constructed to correct for detector effects, and the ⟨⟩ denotes average over events. In order to reduce short-range "non-flow" correlations, pseudorapidity gaps are often explicitly required between the particles in each triplet. This analysis uses the so-called standard, two-subevent and three-subevent methods [18] to explore the influence of non-flow correlations as detailed below. The choices of η ranges for the subevents are identical to those used by the ATLAS experiment [11,12]. In the standard method, all charged particles within η < 2.5 are used. In the two-subevent method, triplets are constructed by combining particles from two subevents labeled as a and c with a ∆η gap in between to reduce non-flow effects: −2.5 < η a < −0.75 , 0.75 < η c < 2.5. The two particles contributing to the flow vector are chosen as one particle each from a and c, while the third particle providing the p T weight is taken from either a or c. In the three-subevent method, three non-overlapping subevents a, b and c are chosen: −2.5 < η a < −0.75 , η b < 0.5 , 0.75 < η c < 2.5. The particles contributing to flow are chosen from subevents a and c while the third particle is taken from subevent b.
A direct calculation of the nested-loop in Eq. (1) is computationally expensive. Instead, it can be expanded algebraically within the multi-particle cumulant framework [18,25] into a polynomial function of vectors and scalars: where the sum runs over particles in a given event or subevent and "k" and "m" are natural integer powers. It is straightforward to show that expansion of Eq. (1) in the three methods gives: R[(q n;1 p 1;1 − τ 1 o n;2 ) a (q * n;1 ) c + (q n;1 p 1;1 − τ 1 o n;2 ) c (q * n;1 ) a ] 1 − (τ 1 ) a + 1 − (τ 1 ) c ⟩ cov(v 2 n , [p T ]) 3sub = ⟨R[(q n;1 ) a (q * n;1 ) c ](p 1;1 ) b ⟩ where the R denotes the real component of the complex number. Experimentally, the v n -[p T ] correlation is aften presented in normalized form known as Pearson's correlation coefficient [10], where the var([p T ]) and var(v 2 n ) are variances of p T fluctuations and v 2 n fluctuations, respectively. The var([p T ]) is obtained using all the pairs in the full event η < 2.5, The dynamical variance var(v 2 n ) are calculated in terms of two-particle cumulant c n {2} and four particle cumulants The c n {4}, being a four-particle correlator, is known to be relatively insensitive to non-flow correlations but usually has poor statistical precision. Therefore it is obtained from the standard cumulant method using the full event. On the other hand, the two particle cumulants c n {2} is more susceptible to non-flow correlations and therefore is calculated from the two-subevent method with the η choices discussed above. This definition is mostly free of non-flow in large collision systems. But in small systems, this definition could still be biased by non-flow effects as we discussed in Appendix A.
To evaluate the influence of non-flow correlations to cov(v 2 n , [p T ]) and ρ(v 2 n , [p T ]), the Pythia8 A2 tune [23] and Hijing v1.37 [24] models are used to generate pp events at √ s = 13 GeV, pPb and peripheral PbPb events at √ s NN = 5.02 TeV, respectively. These models contain significant non-flow correlations from jets, dijets, and resonance decays and can be used to quantify the efficacy of non-flow suppression in these methods. In these simulations, the particle weight are set to be unity, w i = 1 and events are classified by N ch , the number of charged particles in η < 2.5 with p T > 0.1 GeV. The cov(v 2 n , [p T ]) are calculated in three p T ranges using the standard and subevent methods: 0.2 < p T < 2 GeV, 0.5 < p T < 2 GeV, and 0.5 < p T < 5 GeV. They are presented as a function of charged particle density at mid-rapidity dN ch dη, which is assumed to be 1/5 of N ch , N ch ≈ 5dN ch dη. Figure 1 compares the results of cov(v 2 n , [p T ]) from the standard and subevent methods in pp collisions from Pythia8 model. The values from the standard method are positive for all harmonics. This is because the correlations are dominated by the jet fragmentations, which produce clusters of particles with larger p T and enhanced azimuthal correlations at ∆φ ∼ 0, and therefore tend to simultaneously increase the v 2 n and [p T ]. The values from the twosubevent method are positive for even harmonics and negative for odd harmonics, consistent with the dominance of correlations from away-side jet fragments: the away-side correlations are expected to give a more negative v 2 3 and larger [p T ], and therefore a negative value of cov(v 2 3 , [p T ]). For the three-subevent method, the values of cov(v 2 2 , [p T ]) are positive at dN ch dη ≲ 10 and are slightly negative for dN ch dη > 10. The magnitudes of cov(v 2 n , [p T ]) are largest for the standard method, and smallest for the three-subevent method. Similar ordering among the three methods are observed in all three collision systems and all p T selections, and the magnitudes of signal from three-subevent method are always the smallest, suggesting that this method is least affected by non-flow. For the remaining discussion, we focus on discussing results from the three-subevent method.  To further investigate the origin of the sign change of cov(v 2 2 , [p T ]) in the low dN ch dη region, Figure 3 compares the pp results from Pythia8 with those obtained from the Hijing model. The results are in good quantitative agreement for dN ch dη > 20. In the dN ch dη < 30 range and towards lower dN ch dη, the Hijing results show a stronger decrease compared to the Pythia8 results. The Hijing results start to increase at dN ch dη < 10 similar to Pythia8, but except for the lowest p T range of 0.2 < p T < 2 GeV, the increase is not enough for the cov(v 2 2 , [p T ]) to change sign. The results from pp collisions at √ s = 5 TeV are also shown in Figure 3. The values are more negative than those for the √ s = 13 TeV results for dN ch dη < 20, suggesting that residual non-flow is larger at lower √ s at the same dN ch dη.  ) between pp, pPb and PbPb collisions, separately in three p T ranges. The pPb and PbPb values are negative at low dN ch dη region, whose magnitudes increase with p T . This is different from the pp results, which are positive at dN ch dη ≲ 8 region. In the dN ch dη > 10 region, the pp values are negative and lower than those for the pPb and PbPb collisions. The values for pPb collisions are close to but consistently lower than those in PbPb collisions, suggesting a slightly larger residual non-flow in pPb collisions.

III. RESULTS
In order to estimate the non-flow effects on the ρ(v 2 n , [p T ]), we need to choose an appropriate normalization in Eq. (6). The var(v 2 n ) directly obtained from these models should not be used, because they only contain non-flow. Instead, we estimate var(v 2 n ) from the previous published measurements of v n {2} and v n {4} in these three collision systems [6,21,26] as: The v n,tmp {2} were measured using the two-particle correlation and improved template method from Ref. [26] that explicitly subtracts the non-flow correlations. The p T dependence of the v n,tmp {2} are taken from Ref. [26]. The values of v 2 {4} v 2 {2} are taken from Ref. [21] for pp and pPb and from Ref. [6] for PbPb, which are found to be in the range of 0.71-0.74 as a function of dN ch dη, and they are assumed to be independent of p T . The v 2 {4} term ) as a function of dN ch dη from the three-subevent method compared between three collision systems for 0.2 < p T < 2 GeV (left), 0.5 < p T < 2 GeV (middle), and 0.5 < p T < 5 GeV (right).
leads to a 28% reduction to var(v 2 2 ). For third-order harmonics, the values of v 3 {4} v 3 {2} have been found to be very small Ref. [22] and therefore is neglected in this study, i.e. we assume var(v 2 3 ) = v 3,tmp {2} 4 . Examples of the dN ch dη dependence of var(v 2 2 ) and var(v 2 3 ) are given in Figure 7 of Appendix A.
, current statistical uncertainties do not provide a precise lower limit for the non-flow contributions in pPb and PbPb collisions. But in pp collisions and higher p T , the non-flow effects could lead to ρ(v 2 3 , [p T ]) values significantly larger than one. Equipped with these detailed knowledge of non-flow, we are ready to discuss its impact on the interpretation of the v n -[p T ] correlation in terms of ISM and FSM. The top panels of Figure 6 compare the non-flow expectation of cov(v 2 2 , [p T ]) with the ATLAS data [11] 1 . The strength of the non-flow correlations is much smaller than the experimental data in the PbPb collisions (which covers dN ch dη > 20 region), but could be significant in pPb collisions in 0.5 < p T < 2 GeV, reaching a level of around 30-40% of the experimental values at dN ch dη ∼ 20. The results are also compared to the CGC-hydro model for pPb [9] that includes both ISM and FSM but without non-flow. In the dN ch dη > 20 region where the FSM dominates, the model over-predicts the experimental data. In the dN ch dη < 10 region, the CGC-hydro model is dominated by a positive ISM signal, which seems to be smaller in magnitude than the expected non-flow contribution. It might be that the combined non-flow and ISM would still remain negative for the 0.5 < p T < 2 GeV range. For the 0.2 < p T < 2 GeV range where the non-flow contribution is smaller, the combined signal could be slightly positive around dN ch dη ∼ 5 − 10, but would still remain negative at dN ch dη ∼ 5. 2 < p T < 2 GeV (left) and 0.5 < p T < 2 GeV (right). The results are obtained from the three-subevent method and compared with experimental data from ATLAS and CGC-hydro model calculations [9] that include the initial-state momentum correlations. The data points with dN ch dη < 10 in the left panels have been rescaled by the factors in order to fit into the y-ranges.
The bottom panels of Figure 6 show the same comparison in terms of ρ(v 2 2 , [p T ]). The qualitative behaviors are largely the same, with a few important quantitative differences from cov(v 2 2 , [p T ]). The non-flow contributions relative to the experimental data are larger, especially in the pPb collisions, reaching more than 50% of the experimental values at dN ch dη ∼ 20 in 0.5 < p T < 2 GeV. This is due to the fact that the values of var([p T ]) in HIJING are about a factor of 2 smaller than the experimental values, leading to a more negative ρ(v 2 2 , [p T ]) closer to the data. We also caution that the ATLAS var(v 2 n ) data, calculated via Eq. (8), are still biased by non-flow contributions (see Appendix A), which reduce ρ(v 2 2 , [p T ]) slightly further. The main message of Figure 6 is that the interpretation of the cov(v 2 2 , [p T ]) at low dN ch dη region is rather complicated. Firstly, the non-flow contributions from our model studies are negative and could account for some of the observed negative signal in the low dN ch dη range that are also associated with the FSM. Secondly, the negative non-flow contributions compete with the ISM and may eliminate the sign-change in the actual measurement. Thirdly, the fact that Pythia8 model shows a positive ρ(v 2 2 , [p T ]) at dN ch dη < 10 ( Figure 5) suggests that the sign of non-flow contributions is model-dependent and could also be positive. In the latter case, even if experiments observe a positive ρ(v 2 2 , [p T ]), one could not easily interpret this signal as generated by the ISM.

IV. SUMMARY
The influences of non-flow effects to the three-particle correlation between harmonic flow v n and event-by-event average transverse momentum [p T ], cov(v 2 n , [p T ]), are studied in pp, pPb and peripheral PbPb collisions for n = 2 − 4. This study is performed using Pythia8 and Hijing event generators, which contain only non-flow correlations such as fragmentation of jet and dijets and resonance decays, but have no genuine long-range multi-particle correlations from the initial-state or the final-state evolution. The efficacy of non-flow suppression via the rapidity separated three-subevent method has been tested, and is observed to give smallest cov(v 2 n , [p T ]) values in comparison to the standard and two-subevent methods for all harmonics, collision systems and p T ranges investigated in this paper. The values of cov(v 2 2 , [p T ]) from the three-subevent method are negative in the region dN ch dη > 20 and approach zero towards higher dN ch dη. The magnitudes of the cov(v 2 2 , [p T ]) are much smaller than the experimentally measured values in the pPb and PbPb collisions, suggesting that the measured cov(v 2 2 , [p T ]) values in dN ch dη > 20 reflect genuine correlations arising from the final-state interactions. In the region dN ch dη < 20, the values of cov(v 2 2 , [p T ]) decrease toward more negative values in Hijing simulations of pPb and PbPb collisions, but increases in Hijing and Pythia8 simulations of pp collisions. They reach a maximum (positive for Pythia8 but is negative in Hijing) at around dN ch dη ∼ 20 before decreasing again for dN ch dη < 5. The differences between Hijing and Pythia8 suggest that the non-flow contributions in dN ch dη < 20 region are highly model-dependent. The predicted sign change of dN ch dη from initial-state momentum correlation due to gluon saturation physics may not be observed if the non-flow contributions are negative, or unambiguous if the non-flow contributions are positive. Further detailed quantitative model investigation of these different sources are required.
We appreciate comments from G. Giacalone, S. Huang, J. Nagle, and B. Schenke. We thank B. Schenke for sharing the CGC-hydro model calculation, and J. Nagle'suggestion to compare pp collisions between Pythia8 and Hijing. This work is supported by DEFG0287ER40331 and PHY-1913138. In the ATLAS measurement [12], the var(v 2 n ) was calculated using Eq. (8). In the low dN ch dη region, the c n {2} 2sub and the resulting ρ(v 2 n , [p T ]) could be strongly biased by the non-flow correlations. Figure 7 compares var(v 2 n ) from Eq. (8) with those estimated via Eq. (9) based on published v n data in three collision systems. They are presented in terms of 4 var(v 2 n ) in order to be shown in the familiar scale as the single-particle v n values. In pPb and PbPb collisions, the non-flow is sub-dominant for dN ch dη > 20 but can be larger than the genuine flow signal at lower dN ch dη values. In the pp collisions, the non-flow contribution is comparable or larger than the genuine flow signal over the full dN ch dη range.
To estimate the possible bias of the non-flow, we add the var(v 2 n ) from flow and non-flow of Figure 7 in quadrature sum: var(v 2 n ) mod = var(v 2 n ) 2 flow + var(v 2 n ) 2 non−flow . The var(v 2 n ) mod is then used to obtain a modified form of Pearson coefficient ρ(v 2 n , [p T ]) mod . The results are shown in Figure 8. Comparing to the original unbiased results in Figure 5, the magnitudes of the ρ(v 2 n , [p T ]) mod are much reduced in the low dN ch dη region due to the large non-flow bias to var(v 2 n ). The differences between the three systems are also artificially reduced. Therefore, it is important to use a var(v 2 n ) that is free of non-flow effects by following the procedure given in Eq. (9).
[  3 ) (bottom) as a function of dN ch dη compared between three collision systems for 0.2 < p T < 2 GeV (left), 0.5 < p T < 2 GeV (middle), and 0.5 < p T < 5 GeV (right). The data points are calculated from Pythia8 and Hijing models via Eq. (8) and the lines are estimated via Eq. (9) from published vn data.