Examination of Flow and Non-Flow Factorization Methods in Small Collision Systems

TL;DR

This study evaluates methods for separating flow and non-flow effects in small collision systems using Monte Carlo simulations, revealing their effectiveness varies with collision energy and system size.

Contribution

The paper systematically tests flow and non-flow separation techniques in small collision systems using different Monte Carlo generators, highlighting their limitations and applicability.

Findings

Reasonable non-flow subtraction in p+p collisions at LHC energies.

Failure of techniques when applied to p+p collisions at RHIC energies.

Over-subtraction of non-flow effects in p+Au collisions at RHIC for pT > 1 GeV/c.

Abstract

Two particle correlations have been used extensively to study hydrodynamic flow patterns in heavy-ion collisions. In small collision systems, such as $p$$+$$p$ and $p$$+$$A$, where particle multiplicities are much smaller than in $A$$+$$A$ collisions, non-flow effects from jet correlations, momentum conservation, particle decays, etc. can be significant, even when imposing a large pseudorapidity gap between the particles. A number of techniques to subtract the non-flow contribution in two particle correlations have been developed by experiments at the Large Hadron Collider (LHC) and then used to measure particle flow in $p$$+$$p$ and $p$$+$Pb collisions. Recently, experiments at the Relativistic Heavy Ion Collider (RHIC) have explored the possibility of adopting these techniques for small collision systems at lower energies. In this paper, we test these techniques using Monte Carlo generators PYTHIA and HIJING, which do not include any collective flow, and AMPT, which does. We find that it is crucial to examine the results of such tests both for correlations integrated over particle transverse momentum $p_T$ and differentially as a function of $p_T$. Our results indicate reasonable non-flow subtraction for $p$$+$$p$ collisions at the highest LHC energies, while failing if applied to $p$$+$$p$ collisions at RHIC. In the case of $p$$+$Au collisions at RHIC, both HIJING and AMPT results indicate a substantial over-subtraction of non-flow for $p_{T}\gtrsim1~{\rm GeV}/c$ and hence an underestimate of elliptic flow.

Figures

Tables

Table 1. Table 1: Parameters used in ampt .

Parameter	Value
ISOFT	4
PARJ(41)	0.55
PARJ(42)	0.15
Parton screening mass	6.45d0 (0.75 mb)
$\alpha$ in parton cascade	0.47d0

Equations

$$f(\Delta\phi)=G\bigg{\{}1+2\sum_{n=1}^{\infty}a_{n}\cos(n\Delta\phi)\bigg{\}},$$

$$f(\Delta\phi)=G\bigg{\{}1+2\sum_{n=1}^{\infty}(c_{n}+d_{n})\cos(n\Delta\phi)\bigg{\}}.$$

$$J(\Delta\phi)=G\bigg{\{}2\sum_{n=1}^{\infty}d_{n}\cos(n\Delta\phi)\bigg{\}},$$

$$f(\Delta\phi)=J(\Delta\phi)+G\bigg{\{}1+2\sum_{n=1}^{\infty}c_{n}\cos(n\Delta\phi)\bigg{\}}.$$

$$\frac{J^{\mathrm{LM}(\mathrm{HM})}(\Delta \varphi )}{G^{\mathrm{LM}(\mathrm{HM})}}=2\sum _{n=1}^{\infty }{d}_n^{\mathrm{LM}(\mathrm{HM})}\cos (n\Delta \varphi ),$$

$$\frac{J^{\mathrm{HM}}(\Delta \varphi )}{G^{\mathrm{HM}}}=\frac{{\sum }_{n=1}^{\infty }{j}_n{d}_n^{\mathrm{LM}}\cos (n\Delta \varphi )}{{\sum }_{n=1}^{\infty }{d}_n^{\mathrm{LM}}\cos (n\Delta \varphi )}\frac{J^{\mathrm{LM}}(\Delta \varphi )}{G^{\mathrm{LM}}},$$

$$\frac{J^{\mathrm{HM}}(\Delta \varphi )}{G^{\mathrm{HM}}}=j^{\mathrm{ATLAS}}\frac{J^{\mathrm{LM}}(\Delta \varphi )}{G^{\mathrm{LM}}}.$$

$$f^{\mathrm{LM}}(\Delta \varphi )$$

$$f^{\mathrm{HM}}(\Delta \varphi )$$

$$\displaystyle=F^{\mathrm{temp}}f^{\mathrm{LM}}(\Delta\phi)+G^{\mathrm{temp}}\bigg{\{}1+2\sum_{n=2}^{\infty}{c_{n}^{\mathrm{temp}}}\cos(n\Delta\phi)\bigg{\}},$$

$$F^{\mathrm{temp}}$$

$$G^{\mathrm{temp}}$$

$$c_n^{\mathrm{temp}}$$

$$c_n^{\mathrm{HM}}=c_n^{\mathrm{temp}}-j^{\mathrm{ATLAS}}(c_n^{\mathrm{temp}}-c_n^{\mathrm{LM}}),$$

$$\frac{J^{\mathrm{HM}}(\Delta \varphi )}{G^{\mathrm{HM}}}=j^{\mathrm{CMS}}\frac{J^{\mathrm{LM}}(\Delta \varphi )}{G^{\mathrm{LM}}},$$

$$j^{\mathrm{CMS}}=\frac{G^{\mathrm{LM}}{\int }_{-1.2}^{1.2}(f_{\mathrm{SR}}^{\mathrm{HM}}(\Delta \varphi )-f_{\mathrm{LR}}^{\mathrm{HM}}(\Delta \varphi ))\mathrm{d}\Delta \varphi }{G^{\mathrm{HM}}{\int }_{-1.2}^{1.2}(f_{\mathrm{SR}}^{\mathrm{LM}}(\Delta \varphi )-f_{\mathrm{LR}}^{\mathrm{LM}}(\Delta \varphi ))\mathrm{d}\Delta \varphi },$$

$$c_n^{\mathrm{HM}}=a_n^{\mathrm{HM}}-j^{\mathrm{CMS}}{a}_n^{\mathrm{LM}}.$$

$$c_n^{\mathrm{HM}}=a_n^{\mathrm{HM}}-j^{\mathrm{CMS}}(a_n^{\mathrm{LM}}-c_n^{\mathrm{LM}}),$$

$$f(\Delta \varphi )$$

$$\displaystyle=\int B(0,0)\frac{S(\mbox{$\Delta\eta$},\mbox{$\Delta\phi$})}{B(\mbox{$\Delta\eta$},\mbox{$\Delta\phi$})}\mathrm{d}\mbox{$\Delta\eta$},$$

$$S(\mbox{$\Delta\eta$},\mbox{$\Delta\phi$})=\frac{1}{N^{\rm Trig}}\frac{\mathrm{d}^{2}N^{\rm Same}}{\mathrm{d}\mbox{$\Delta\eta$}\mathrm{d}\mbox{$\Delta\phi$}},$$

$$B(\mbox{$\Delta\eta$},\mbox{$\Delta\phi$})=\frac{1}{N^{\rm Trig}}\frac{\mathrm{d}^{2}N^{\rm Mix}}{\mathrm{d}\mbox{$\Delta\eta$}\mathrm{d}\mbox{$\Delta\phi$}}.$$

Full text

I Introduction

The standard evolution model of relativistic heavy ion collisions involves the translation of initial geometric anisotropies into final momentum correlations via nearly perfect hydrodynamic flow followed by hadronic rescattering Romatschke and Romatschke (2017); Heinz and Snellings (2013). Starting with first indications in $p+p$ collisions at the CERN Large Hadron Collider (LHC) Khachatryan et al. (2010), there is now a wealth of data indicating similar translation of geometry into flow in smaller collisions systems at the BNL Relativistic Heavy Ion Collider (RHIC) and the LHC Nagle and Zajc (2018). Here we focus on $p+p$ collisions at the LHC and $p+p$ and $p+\mathrm{Au}$ collisions at RHIC where one is pushing the limits of how small a system and how small a particle multiplicity still results in significant final flow signatures. It is imperative to consider both collision energies and systems since a simultaneous description within the perfect fluid paradigm currently gives the best data description Weller and Romatschke (2017); Aidala et al. (2019).

In these small systems, there are many contributions to the final particle correlations and some of these contributions increase in relative strength with decreasing multiplicity. There are many sources of correlations among only a subset of the final particles that result in an increased probability of particles to be nearby in momentum space (i.e., close in both pseudorapidity $\Delta\eta\approx 0$ and azimuthal angle $\Delta\phi\approx 0$). Examples include the decay of heavy resonances, correlations due to quantum statistics (HBT effects) Lisa et al. (2005), and the fragmentation of a high momentum quark or gluon—the last of these resulting in what are referred to as jet correlations. These correlations are effectively minimized by requiring the two particles to have a substantial gap in pseudorapidity, typically $|\Delta\eta|>2$. There are other correlations that survive such a cut including dijet correlations. In a leading-order hard scattering, a parton (quark or gluon) and partner parton will be nearly back-to-back in azimuth ($\Delta\phi\approx\pi$) and can be widely separated in pseudorapidity, since the incoming partons need not be longitudinally momentum balanced. As such, one may expect two-particle correlations with an enhancement near $\Delta\phi\approx\pi$ that extends long range in pseudorapidity.

As shown in Figure 1, one can imagine a toy scenario where a medium indicated by the blue region is created in a collision. In addition, a hard scattering results in two partons that largely fragment into hadrons outside of the medium. In a perfectly factorized picture, the medium particles may have angular correlations that reflect the flowing fluid while the jet particles are correlated with each other but are uncorrelated with the medium geometry. As one examines collisions with larger final particle multiplicities and larger volume, a larger fraction of the particles come from the medium and thus the influence of the jet correlations is reduced. In this simplest scenario the correlations from the jet particles are referred to as “nonflow” and those from the medium particles are referred to as “flow.” We highlight that it is impossible for this factorization to be perfectly true. The jet partons cannot hadronize without some color interaction with the medium or underlying event, and in some case many of the jet partons scatter with medium partons. For that reason there is no direct separability of the two, rather they are related through a convolution.

Another key source of nonflow correlations is simply global momentum conservation. The simplest case would be to imagine $N$ particles emitted isotropically but obeying momentum conservation. The root software package Brun and Rademakers (1997) has a class TGenPhaseSpace that allows for $N$-body decays simply following phase space filling and momentum conservation.

Of course, we do not expect that the $N$ particles are distributed simply according to phase space rules. Within a model such as pythia8 Sjostrand et al. (2008), one can think of every two-jet event as really being at least a four-jet event. As shown in panel (b) of Fig. 1, if two incoming partons have a large momentum transfer, the resulting beam remnants extend nearly longitudinal color strings (i.e., the other two jets). These strings can have a transverse momentum kick such that particles emitted from the upper (orange) string have a slight trend upward and from the lower (orange) string a slight trend downward. Thus, the partons and resulting hadrons will have a momentum conservation correlation that may be convolved with a flow correlation if these particles undergo additional final state scattering.

The $\Delta\phi$ correlation function with the requirement $|\Delta\eta|>2$ is shown in Fig. 2, where the $N$-body phase space calculation is shown in the left panel and the jet-type correlation from pythia8 is shown in the right panel. Most correlation analyses characterize the two-particle correlation in $\Delta\phi$ in terms of a Fourier decomposition Voloshin and Zhang (1996):

[TABLE]

where the coefficients $a_{n}$ are the Fourier coefficient at order $n$, $G$ is the normalization factor corresponding to the average number of associated particles per trigger particle in the sample, and all sine terms vanish due to symmetry. The coefficients $a_{n}$ are shown for the two nonflow contributions in Fig. 2. In the case of the $N$-body phase space, the dominant contribution is $c_{1}$ with a depletion near $\Delta\phi\approx 0$ and an enhancement of pairs near $\Delta\phi\approx\pi$. However, there are nonzero higher order contributions and their relative strengths will of course depend on the other correlations embedded along with momentum conservation. For the jet-type correlations, the dominant term is also $a_{1}$ due to the strong long-range away-side peak. However, the successive terms contribute significantly and with alternating signs in order to describe the nearly flat region around $\Delta\phi\approx 0$ and the peak around $\Delta\phi\approx\pi$.

In the following sections we define the nonflow subtraction methods and then detail the resulting tests of these methods using the Monte Carlo models pythia Sjostrand et al. (2008), hijing Gyulassy and Wang (1994), and ampt Lin et al. (2005) in different collisions systems and and at different energies.

II Definition of Methods

The most important assumption in all nonflow subtraction methods is the assumption that the correlation coefficients $a_{n}$ in Eq. (1) can be separated into two linearly additive contributions $a_{n}=c_{n}+d_{n}$, where $c_{n}$ is the flow coefficient quantifying correlations related to the initial geometry and $d_{n}$ is the nonflow coefficient of pair correlations. As discussed in Sec. I, this assumption may not be realistic, but this is the starting point for all nonflow subtraction methods, and so we begin our discussion here. With this assumption, Eq. (1) can be rewritten as

[TABLE]

In order to extract the true flow, $c_{n}$, it is necessary to remove the nonflow contributions, $d_{n}$. This is particularly important in small systems where the nonflow correlations can dominate the flow signal. We can define the nonflow contribution $J(\Delta\phi)$ to $f(\Delta\phi)$ as

[TABLE]

so that Eq. (2) becomes

[TABLE]

All techniques attempting to disentangle flow and nonflow operate by comparing correlations between two data selection samples: one at higher multiplicity (HM), where flow is expected to have a larger influence; and one at lower multiplicity (LM), where flow is expected to have a smaller influence. The flow in LM events is at times assumed to be negligible or non-existent; however, we need not consider such differences just yet. The above equations simply need to be trivially relabeled so that the LM and HM categories are distinct. Here we use LM (HM) as a superscript on all relevant quantities to indicate LM (HM) events.

The terms $c_{n}^{\mathrm{HM}}$ and $c_{n}^{\mathrm{LM}}$ are the flow correlation coefficients at low and high multiplicity, respectively; $d_{n}^{\mathrm{HM}}$ and $d_{n}^{\mathrm{LM}}$ are the nonflow correlation coefficients at low and high multiplicity, respectively.

Applying these labels and rewriting Eq. (3) we obtain

[TABLE]

meaning the $J^{\mathrm{LM}}(\Delta\phi)$ and $J^{\mathrm{HM}}(\Delta\phi)$ are related to each other as

[TABLE]

where the relational coefficients are $j_{n}=d_{n}^{\mathrm{HM}}/d_{n}^{\mathrm{LM}}$. By estimating the $j_{n}$ and the flow coefficients at low multiplicity, $c_{n}^{\mathrm{LM}}$, the nonflow contribution at high multiplicity can be determined and subtracted. There are a few different nonflow subtraction methods available on the market based on different assumptions regarding the $c_{n}^{\mathrm{LM}}$ and $j_{n}$ coefficients.

II.1 Method 1

The first method we detail has been developed by the ATLAS Collaboration, as originally applied in Ref. Aad et al. (2016). Here we detail the key assumptions made in using this method. The first assumption is that the nonflow correlation coefficient $j_{n}=d_{n}^{\mathrm{HM}}/d_{n}^{\mathrm{LM}}$ is independent of the harmonic number $n$, so that $j_{n}=j^{\mathrm{ATLAS}}$ for all $n$. This allows Eq. (6) to be rewritten as:

[TABLE]

This assumption follows from the idea that at the LHC the nonflow correlation is dominated by the jet contribution. The additional requirement is then that the shape of the jet correlation does not change with multiplicity event class.

A second assumption made is that the flow contribution at the lowest harmonic $c_{1}$ is negligible compared to the nonflow $d_{1}$, i.e., $c_{1}\ll d_{1}$. Hence $c_{1}$ is completely ignored; alternately one can think of it as being absorbed via redefinition of $d_{1}$. Following these assumptions, we can rewrite the angular correlation distributions in Eq. (4) for HM events in relation to LM events,

[TABLE]

and then combining the two using Eq. 7 to obtain

[TABLE]

where $F^{\mathrm{temp}}$, $G^{\mathrm{temp}}$, and $c_{n}^{\mathrm{temp}}$ are parameters defined by

[TABLE]

There parameters are obtained from fitting $f^{\mathrm{HM}}(\Delta\phi)$ with $f^{\mathrm{LM}}(\Delta\phi)$. Thus, this approach is often referred to as a “template-fitting” method because the LM correlation function serves as a template for the HM correlation function.

In the special case where the flow coefficients are identical at all orders between the low multiplicity and high multiplicity classes, i.e., $c_{n}^{\rm HM}=c_{n}^{\rm LM}$ for all $n$, then the extracted $c_{n}^{\mathrm{temp}}$ is exactly equal to $c_{n}^{\mathrm{HM}}$. This is a simplifying assumption that has no concise physics motivation. This special case was assumed in the ATLAS publications Aad et al. (2016); Aaboud et al. (2017).

Relaxing this special case, an additional correction should be applied to the fitted value to obtain $c_{n}^{\mathrm{HM}}$:

[TABLE]

where $j_{n}$ can be obtained using fit parameters $j^{\mathrm{ATLAS}}=G^{\mathrm{LM}}F^{\mathrm{temp}}/G^{\mathrm{HM}}$ and $c_{n}^{\mathrm{LM}}$ is usually estimated using the flow coefficient measured in the second lowest multiplicity sample. The ATLAS Collaboration applies this additional correction in Ref. Aaboud et al. (2019).

II.2 Method 2

The second method we detail has been developed by the CMS Collaboration, as originally applied in Ref. Khachatryan et al. (2017). Similar approaches have also been explored by the ATLAS and ALICE collaborations Aad et al. (2013); Abelev et al. (2013). Here we detail the key assumptions made in using this method.

Like the ATLAS method, this method also assumes that the nonflow correlation shape does not change with multiplicity class, hence $j_{n}=d_{n}^{\mathrm{HM}}/d_{n}^{\mathrm{LM}}$ is independent of the harmonic number $n$, so that $j_{n}=j^{\mathrm{CMS}}$ for all $n$. This allows Eq. (6) to be rewritten as

[TABLE]

in exact analogy to the ATLAS case. However, the method by which $j^{\mathrm{CMS}}$ is determined is different. They assume that the nonflow correlation coefficients, at least for $n>1$, are dominated by jet contributions. Since they assume the jet shape is independent of multiplicity class, they only need to determine the relative jet yields in the different multiplicity classes. They do this by measuring the near-side jet yield in the different multiplicity classes via the short range (SR) correlation, i.e., with $\Delta\phi\approx 0$ and $\Delta\eta\approx 0$. They isolate the near-side jet peak by taking the difference between the SR and long-range (LR) correlations as follows:

[TABLE]

where $f^{\mathrm{HM(LM)}}_{\mathrm{SR}}(\Delta\phi)$ is the short-range correlation with $|\Delta\eta|<1$ and $f^{\mathrm{HM(LM)}}_{\mathrm{LR}}(\Delta\phi)$ is the long-range correlation with $2<|\Delta\eta|<5$. The short-range and long-range correlations are usually normalized by number of trigger particles such that the integral of their difference over $-1.2<\Delta\phi<1.2$ is close to the number of particles produced per jet.

Next, direct Fourier fits are applied to the long-range correlations at low and high multiplicities, and the Fourier coefficient $a_{n}^{\mathrm{HM(LM)}}$ are extracted. This approach assumes that there is no flow contribution at all in the low multiplicity event selection, i.e., $c_{n}^{\mathrm{LM}}=0$. This is an extreme assumption that forces any extraction of flow coefficients as a function of multiplicity class to approach zero at low multiplicity, and is quite different from the ATLAS approach. Then at each order $n$ the flow coefficient at high multiplicity is given by

[TABLE]

As noted above, as measurements move towards low multiplicity, such that $c_{n}^{\mathrm{HM}}\approx c_{n}^{\mathrm{LM}}$, $c_{n}^{\mathrm{HM}}$ will converge to zero by construction. The CMS results using this method thus always trend to zero. However, if a nonzero flow coefficient in the low multiplicity sample is allowed, Eq. (18) could be rewritten as

[TABLE]

where $c_{n}^{\mathrm{LM}}$ can be estimated in the same way as in the ATLAS method. Once the generalized results as shown in Eq. (15) and Eq. (18) are used, one should be able to obtain similar results with the two methods.

One remaining difference between the methods would arise from the different ways of estimating the jet variable $j_{n}$. The nonflow correlation can have contributions from jets as well as overall momentum conservation. The momentum conservation contribution is predominantly in the $d_{1}$ coefficient. Thus, in the CMS method where they subtract the nonflow order by order, even though the $j^{\mathrm{CMS}}$ is determined from jets alone, it does not matter since they are not extracting a first-order $c_{1}$ at high multiplicity. In contrast, in the ATLAS template-fitting method, it assumes the nonflow shape (combining both jets and momentum conservation) scale in the same way. This seems unlikely given the findings by CMS via the short-range correlations. Since neither experiment extracts a $c_{1}$ at high multiplicity where they would get very different results, there is only the potential for a residual effect on the higher order $c_{n}$ via the ATLAS fitting procedure.

II.3 Jet Shape Assumption

Method 1 and 2 are derived from the same expression and heavily rely on the assumption that the jet shape is the same at low and high multiplicities. This assumption can be tested in pythia8, in which no flow contribution is expected, i.e., $c_{n}=0$. We perform a direct Fourier coefficient extraction in the low and high multiplicity event classes and then compute $j_{n}=d_{n}^{\mathrm{HM}}/d_{n}^{\mathrm{LM}}$. To examine if all $j_{n}$ are the same, we plot $(j_{n}-j_{2})/j_{2}$ as shown in Fig. 3. The points for $n=1$ and $n=3$ indicate that the assumption of a common $j_{n}$ is violated in pythia8 and that violation increases with higher multiplicity, reaching a level for $n=1$ ($n=3$) relative to $n=2$ of $-20\%$ (+40%). We note that the $j_{2}$ is decreasing with multiplicity, so the impact of this violation might be actually smaller at higher multiplicity. We also show the $j^{\mathrm{ATLAS}}$ and $j^{\mathrm{CMS}}$ extracted values. One can see that they closely agree with the $n=1$ case and deviate significantly from the $n=3$ case. Note that the $j^{\mathrm{ATLAS}}$ is actually identical to the $n=1$. This is because there is no first-order flow coefficient ($c_{1}=0$) by construction. The fit procedure results in $j^{\mathrm{ATLAS}}$ being determined completely by the first-order nonflow correlation coefficient ($d_{1}$), while the contribution from higher order nonflow correlation coefficients is absorbed into the higher order flow correlation coefficients.

We note that the $j_{n}$ assumption can also be sensitive to the particular pseudorapidity gap chosen and the low multiplicity reference selection. The sensitivity of the $j_{n}$ estimated by Method 1 to the choice of pseudorapidity gap and low multiplicity selection is shown in Fig. 4. The method has particular problems when the pseudorapidity gap is $1<|\Delta\eta|<1.5$, because there are contributions from both the near- and away-side jet. The method also is most sensitive when the low multiplicity selection is at its lowest.

All of these violations of the assumptions in Method 1 and 2 can only be gauged in terms of consequences on the extracted flow coefficients $c_{n}$ by testing the methods on various Monte Carlo physics models. In the following sections, we examine the results of such tests.

III Closure Tests with pythia8 and hijing

We now test these procedures on Monte Carlo generators. In Monte Carlo generators such as pythia8 Sjostrand et al. (2008) and hijing Gyulassy and Wang (1994), there is no collective flow in a sense that there are no final state interactions to translate a spatial geometry into momentum anisotropies. Thus, one expects that the application of a successful nonflow subtraction method should result in flow anisotropy coefficients of exactly zero in these cases. These are thus referred to as closure tests. By measuring the residual of these coefficients $c_{n}$, the level of closure can be quantified. We have applied the methods described in Sec. II to determine the level of closure in various Monte Carlo generators.

Before proceeding we want to define the nomenclature used in the following sections. In our studies we focus on the extraction of elliptic flow $n=2$, and show all results in terms of $v_{22}$ and $v_{22}^{\rm sub}$ extracted from the two-particle correlations. First we define \mbox{v_{22}}=a_{2}(p_{T,1},p_{T,2}), where the trigger and associated particles are in momentum selections around $p_{T,1}$ and $p_{T,2}$, respectively. Once the nonflow subtraction technique is applied and the $c_{n}$ coefficients are estimated, we define \mbox{v_{22}^{\rm sub}}=c_{2}(p_{T,1},p_{T,2}). These quantities are related but not equal to $v_{2}^{2}$. The standard differential $v_{2}$ as a function of $p_{T}$ would then be extracted as v_{2}(p_{T,1})=\mbox{v_{22}}/\sqrt{a_{2}(p_{T,2},p_{T,2})} and v_{2}^{\rm sub}(p_{T,1})=\mbox{v_{22}^{\rm sub}}/\sqrt{c_{2}(p_{T,2},p_{T,2})}. In the case where both $p_{T,1}$ and $p_{T,2}$ represent the same broad range in $p_{T}$, e.g., 0.5<\mbox{p_{T}}<5~{}{\rm GeV}/c, these expressions reduce to v_{2}=\sqrt{\mbox{v_{22}}}=\sqrt{a_{2}} and v_{2}^{\rm sub}=\sqrt{\mbox{v_{22}^{\rm sub}}}=\sqrt{c_{2}}, representing integral $v_{2}$ as a function of some event-level variable (multiplicity, centrality, etc.).

We also highlight that experiments have different techniques for selecting lower and higher multiplicity events. If the selection is based on charged particle multiplicity centered around midrapidity, i.e., in the same range as the two particles used for the correlations, we label the event categories in terms of $N^{\rm ch}$. In contrast, other measurements utilize multiplicity or energy in a forward or backward rapidity range, for example in the Pb-going direction in $p+\mathrm{Pb}$ collisions, and thus outside the range of the two particles used for the correlations. In this second case, we refer to the event selections by “centrality.”

III.1 LHC $p+p$ Case

These methods were developed particularly by the ATLAS and CMS Collaborations for use in the highest energy $p+p$ collisions, and we examine that collision system first. We note that the CMS Collaboration has presented the results of a pythia8 closure test for particles integrated over $p_{T}$ in their paper (see Fig. 4 of Ref. Khachatryan et al. (2017)), but such a study has not been published by ATLAS. We consider acceptance selection cuts similar to those in the experiments, though use one set so that all methods are compared apples-to-apples. Correlations are determined from all charged hadrons within $|\eta|<2.5$ from $10^{8}$ pythia8 $p+p$ events (SoftQCD:nonDiffractive = on) at $\sqrt{s}$ = 13 TeV. Charged hadrons with \mbox{p_{T}}>0.4~{}{\rm GeV}/c (\mbox{p_{T}}>0.5~{}{\rm GeV}/c) are used for event multiplicity categorization (particle correlations). We note that a slightly different $p_{T}$ selection (\mbox{p_{T}}>0.3~{}{\rm GeV}/c) has been used for data analysis by the CMS collaboration.

Figure 18 (shown in the Appendix for clarity) shows one dimensional $\Delta\phi$ two-particle correlation functions for short-range (|\mbox{\Delta\eta}|<1) and long-range (2<|\mbox{\Delta\eta}|<5) regions. Charged hadrons are required to satisfy \mbox{p_{T}}>0.5~{}{\rm GeV}/c and $|\eta|<2.5$. Each panel presents a different range of charged hadron multiplicity event selection $N^{\rm ch}$. In this case, there is no visually obvious shape variation of the $\Delta\phi$ correlation functions both in short and long ranges throughout the entire multiplicity range. This implies that the $j_{n}$ values are going to be approximately independent of $n$, though not exactly as demonstrated in the previous discussion. The dashed lines are fits to the distributions of long-range correlations to extract Fourier coefficients described in Eq. (1).

Figure 5 shows the second-order Fourier coefficients $v_{22}$ extracted directly from the correlation functions and after the nonflow subtraction technique is applied. Results are shown for charged hadrons with 0.5<\mbox{p_{T}}<5~{}{\rm GeV}/c as a function of event multiplicity $N^{\rm ch}$. Figure 6 show the results for the highest multiplicity selection and for trigger particles as a function of $p_{T}$. The left and right panels of both figures apply the Method 1 (ATLAS) and Method 2 (CMS) methods, respectively.

The nonzero value of the coefficients after subtraction indicate that the closure test is not perfectly satisfied. This is mainly due to remaining jet correlations from the near-side (\mbox{\Delta\phi}\approx 0), even with the large $\Delta\eta$ gap, and/or a small shape change from the away-side (\mbox{\Delta\phi}\approx\pi). The nonflow effect on $v_{22}$ is larger at low multiplicity events, and it becomes smaller as the event multiplicity increases. The two nonflow subtraction methods are applied with three different ranges of low multiplicity selection, and the results of the subtracted $v_{22}^{\rm sub}$ from each method are shown in each panel.

In Fig. 5 the gray bands correspond to $|v_{2}|=0.03$ (thus equivalent to \mbox{v_{22}}=0.03\times 0.03). There is some sensitivity to the low multiplicity selection. The $v_{22}^{\rm sub}$ with low multiplicity selection 10\leq\mbox{N^{\rm ch}}<20 and 20\leq\mbox{N^{\rm ch}}<30 are within the level of $|v_{2}|<0.03$ over the entire multiplicity range. Using events in 0\leq\mbox{N^{\rm ch}}<10 as a low multiplicity bin results in $v_{22}^{\rm sub}$ with a slightly larger deviation from zero in the lower multiplicity range, but it converges with the $v_{22}^{\rm sub}$ from the other cases at higher multiplicities.

The same test is done as a function of the $p_{T}$ of the trigger particle ($p_{T}^{\rm Trig}$) for the highest 5% multiplicity events (\mbox{N^{\rm ch}}\geq 85). Figure 19 in the Appendix shows the $\Delta\phi$ correlation functions for low (0\leq\mbox{N^{\rm ch}}<10) and high (\mbox{N^{\rm ch}}\geq 85) multiplicity events. Each panel shows a different $p_{T}^{\rm Trig}$ range, and the $p_{T}$ range of associated particles is 0.5<\mbox{p_{T}^{\rm Assoc}}<5~{}{\rm GeV}/c. The shape of the near-side peak in the short-range correlation function becomes narrower as $p_{T}^{\rm Trig}$ increases, and there is no significant shape difference between low and high multiplicity events at the same $p_{T}^{\rm Trig}$.

Figure 6 shows the $v_{22}$ of long-range $\Delta\phi$ correlations as a function of $p_{T}^{\rm Trig}$. Here, the $v_{22}$ increases with $p_{T}^{\rm Trig}$ indicating a stronger nonflow effect at higher $p_{T}^{\rm Trig}$. The lines represent $v_{22}^{\rm sub}$ using events from three different multiplicity ranges for the nonflow subtraction. As with the nonflow results as a function of multiplicity, the $v_{22}^{\rm sub}$ with reference events in 10\leq\mbox{N^{\rm ch}}<20 and 20\leq\mbox{N^{\rm ch}}<30 are quite consistent, and the $v_{22}^{\rm sub}$ values are smaller than 0.001 in 0.5<\mbox{p_{T}^{\rm Trig}}<5~{}{\rm GeV}/c. The $v_{22}^{\rm sub}$ with events of the lowest multiplicity range (0\leq\mbox{N^{\rm ch}}<10) is slightly different from the other two cases and shows a greater degree of nonclosure at higher $p_{T}$ using Method 2.

To summarize the closure test in $p+p$ collisions at $\sqrt{s}$ = 13 TeV from pythia8:

The $\Delta\phi$ two-particle correlation functions in pythia8 exhibit minor violations of the assumptions in the nonflow subtraction methods examined. 2. 2.

The pythia8 resulting $v_{22}^{\rm sub}$ using Method 1 and 2 pass the closure test much better than $|v_{2}|<0.03$ as long as one avoids the lowest multiplicity range for the reference selection. 3. 3.

The modest degree of nonclosure with pythia8 may be considered as a systematic uncertainty on the final extracted $v_{2}$ results.

III.2 RHIC $p+p$ Case

There are currently no measurements of extracted flow coefficients in $p+p$ collisions at RHIC. These are challenging measurements due to multiple collision pileup, much lower multiplicities compared with LHC collisions, and the more limited phase space acceptance of the RHIC experiments. However, we include this case for completeness and to inform future studies. For this test, we use all charged hadrons within $|\eta|<2.5$ from $10^{9}$ pythia8 $p+p$ events (SoftQCD:nonDiffractive = on) at $\sqrt{s}$ = 200 GeV, and charged hadrons in \mbox{p_{T}}>0.2~{}{\rm GeV}/c are used for event multiplicity categorization. For these studies, we have modeled a very large acceptance similar to that of the LHC experiments.

Two-particle $\Delta\phi$ and $\Delta\eta$ correlation functions are made with the same definition introduced in Sec. III.1. One-dimensional $\Delta\phi$ correlation functions in short (|\mbox{\Delta\eta}|<1) and long (2<|\mbox{\Delta\eta}|<5) ranges in various multiplicity selections are presented in Fig. 20 in Appendix A. Charged hadrons in 0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c and $|\eta|<2.5$ are used for the correlation function.

One obvious difference from the $\Delta\phi$ correlation functions of the LHC case (Figure 18) is that the shape of short-range $\Delta\phi$ correlation function significantly changes in the multiplicity range of 0\leq\mbox{N^{\rm ch}}<40 at $\sqrt{s}$ = 200 GeV. In the case of the lowest multiplicity bin, (0\leq\mbox{N^{\rm ch}}<10), the per-trigger-yield in the short-range correlation is smaller than that in the long-range correlation, so Method 2 is simply not applicable within this multiplicity bin because it results in a negative jet yield.

Figure 7 shows the $v_{22}$ from Fourier fits to the long-range $\Delta\phi$ correlation functions plotted versus multiplicity $N^{\rm ch}$. Additionally plotted as solid lines are the $v_{22}^{\rm sub}$ values using three different low multiplicity reference ranges. The gray bands correspond to a window of $|v_{2}|=0.03$. In results from Method 1 shown in the left panel. The $v_{22}^{\rm sub}$ clearly depends on the selection of low multiplicity reference, and the $v_{22}^{\rm sub}$ greatly deviate from zero in lower multiplicity ranges. This is because of the dramatic shape variation of the $\Delta\phi$ correlations in events of 0\leq\mbox{N^{\rm ch}}<40. When using events in the 20\leq\mbox{N^{\rm ch}}<30 multiplicity range as reference, the $v_{22}^{\rm sub}$ is within the window of $|v_{2}|=0.03$.

In the right panel of Fig. 7, the nonflow subtraction results using Method 2 are presented; however, the lowest multiplicity bin (0\leq\mbox{N^{\rm ch}}<10) is not included due to the negative jet yield as discussed earlier. The $v_{22}^{\rm sub}$ from Method 2 using the remaining reference multiplicity ranges are within the window of $|v_{2}|=0.03$. One possible reason for these results differing from Method 1 is an interplay of the shape variation of $\Delta\phi$ correlation functions at both short and long range. In Method 1, the particular multiplicity dependence of the long-range correlation function results in a large deviation of the $v_{22}^{\rm sub}$ from zero. However in Method 2, which also uses information about the scaling of the jet yield at short range, this effect is partially compensated and the deviation is significantly smaller. As shown in Fig. 20, the difference of minimum per-trigger yields in short-range and long-range correlations strongly changes with multiplicity, and the additional jet yield in low multiplicity events possibly results in a smaller scaling for $v_{22}$ of the low multiplicity events.

The outcome of the nonflow subtraction as a function of $p_{T}$ in $p+p$ collisions at $\sqrt{s}$ = 200 GeV has been studied as well. $\Delta\phi$ correlation functions with different $p_{T}^{\rm Trig}$ ranges at short and long ranges are presented in Fig. 21 in the Appendix. Charged hadrons in $|\eta|<2.5$ are used for the two-particle correlation, and the $p_{T}$ range of associated particles is 0.2<\mbox{p_{T}^{\rm Assoc}}<3~{}{\rm GeV}/c. Events of \mbox{N^{\rm ch}}\geq 50 corresponding to the highest 5% multiplicity are selected as high multiplicity events. The shape of short-range $\Delta\phi$ correlation functions in the lowest multiplicity range (0\leq\mbox{N^{\rm ch}}<10) show a large $p_{T}$ dependence in 0.2<\mbox{p_{T}^{\rm Trig}}<2~{}{\rm GeV}/c, and the short-range $\Delta\phi$ correlation function of 0.2<\mbox{p_{T}^{\rm Trig}}<0.5~{}{\rm GeV}/c in the high multiplicity bin also show a quite different shape compared to the $\Delta\phi$ correlation functions of other $p_{T}^{\rm Trig}$ bins.

The $v_{22}$ from Fourier fits as well as $v_{22}^{\rm sub}$ are presented in Fig. 8, and the left (right) panel shows the $v_{22}^{\rm sub}$ with Method 1 (Method 2). In the nonflow subtraction results with Method 1 using only long-range $\Delta\phi$ correlation functions, the $v_{22}^{\rm sub}$ strongly depends on $p_{T}^{\rm Trig}$, and the deviation from zero becomes larger as $p_{T}^{\rm Trig}$ increases. The difference is largest when using events in the lowest multiplicity range (0\leq\mbox{N^{\rm ch}}<10) as the reference. These $v_{22}^{\rm sub}$ results indicate that the shape of the long-range $\Delta\phi$ correlation function changes with multiplicity so that the scaled low multiplicity $\Delta\phi$ correlation does not perfectly describe those in high multiplicity events. The $v_{22}^{\rm sub}$ from Method 2 in the right panel shows a smaller deviation from zero than that of Method 1, and the nonflow subtraction result using the lowest multiplicity reference (0\leq\mbox{N^{\rm ch}}<10) is worse than the two other cases. It is notable that there is a much smaller shape variation when going from low to high multiplicity events when categorizing event activity at forward rapidity. This can be considered for future analyses with $p+p$ data at RHIC.

To summarize the closure test in $p+p$ collisions at $\sqrt{s}$ = 200 GeV from pythia8:

The pseudorapidity coverage $|\eta|<2.5$ was used to compare with the closure test results of the LHC (Sec. III.1). 2. 2.

The $\Delta\phi$ two-particle correlation functions exhibit a clear multiplicity and $p_{T}^{\rm Trig}$ dependent shape variation unlike the case in $p+p$ collisions at $\sqrt{s}$ = 13 TeV from pythia8. 3. 3.

This shape variation results in a significant dependence on the low multiplicity reference selection for both methods. 4. 4.

$v_{22}^{\rm sub}$ from Method 1 shows a much larger deviation from zero compared to the results in the LHC case. 5. 5.

$v_{22}^{\rm sub}$ from Method 2 gives a smaller deviation than that of Method 1, but Method 2 cannot be applied using the some low multiplicity references and $p_{T}^{\rm Trig}$ ranges due to a negative jet yield.

III.3 RHIC $p+\mathrm{Au}$ Case

There have been numerous extractions of flow coefficients at RHIC in $p+\mathrm{Au}$, $d+\mathrm{Au}$, and ${}^{3}\mathrm{He}+\mathrm{Au}$ collisions – highlighted by the PHENIX publications of elliptic $v_{2}$ and triangular $v_{3}$ flow in all three systems Aidala et al. (2019). The PHENIX results are shown with no subtraction of the nonflow and instead with asymmetric systematic uncertainties to estimate the possible contributions. These correlations have been checked with a pseudorapidity gap as large as $|\Delta\eta|>2.75$ using the PHENIX central arm tracks ($|\eta|<0.35$) and the Au-going Beam-Beam Counter ($-3.9<\eta<-3.1$).

Recently the STAR experiment has shown preliminary results using tracks in their Time Projection Chamber only, with pseudorapidity $|\eta|<0.9$ and $|\Delta\eta|>1.0$ Huang (2019). The much smaller pseudorapidity gap yields a much larger nonflow contribution with influences on both the near-side \mbox{\Delta\phi}\approx 0 and the away-side \mbox{\Delta\phi}\approx\pi. They have employed multiple of the above outlined subtraction techniques to extract preliminary flow coefficients and find smaller $v_{2}$ than the PHENIX results particularly for \mbox{p_{T}}>1.5~{}{\rm GeV}/c in high multiplicity $p+\mathrm{Au}$ and $d+\mathrm{Au}$ events Huang (2019).

Here we examine the nonflow subtraction in various kinematic ranges and multiplicity classifications in $p+\mathrm{Au}$ collisions. $1\times 10^{9}$ $p+\mathrm{Au}$ and $2\times 10^{9}$ $p+p$ events were generated with hijing Gyulassy and Wang (1994), and charged hadrons were selected for two-particle correlations. We explore a nonflow subtraction with $p+p$ events from hijing in addition to using low multiplicity or peripheral $p+\mathrm{Au}$ events. First, the case of a wide pseudorapidity coverage ($|\eta|<2.5$) similar to the LHC experiments has been studied to check any difference with the same kinematic range but in lower collision energy (as was done with the study in pythia8 $p+p$ collisions at $\sqrt{s}$ = 200 GeV). Then we detail a study modeling the more limited STAR acceptance.

Figure 22 in the Appendix shows two-particle $\Delta\phi$ correlation function in short (|\mbox{\Delta\eta}|<1) and long (2<|\mbox{\Delta\eta}|<5) range in various multiplicity bins in $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV from hijing. Charged hadrons in 0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c and $|\eta|<2.5$ are used for the correlation functions, and the multiplicity is defined as the number of charged hadrons in \mbox{p_{T}}>0.2~{}{\rm GeV}/c and $|\eta|<2.5$. Similar to the case of pythia8 $p+p$ in $\sqrt{s}$ = 200 GeV, the shape of the two-particle $\Delta\phi$ correlation function in the lowest multiplicity bin (0\leq\mbox{N^{\rm ch}}<10) is quite different from that in the higher multiplicity ranges.

The $v_{22}$ from direct Fourier fits to the two particle $\Delta\phi$ correlation in long range are presented in Fig. 9, and the lines represent the $v_{22}^{\rm sub}$ using low multiplicity events in three different ranges for nonflow subtraction. The $v_{22}^{\rm sub}$ with events in the lowest multiplicity range (0\leq\mbox{N^{\rm ch}}<10) show the largest deviation from zero in the results with both methods. The worse closure with Method 2 is probably related to the shape of near-side short-range correlation in low multiplicity events (0\leq\mbox{N^{\rm ch}}<10). The $v_{22}^{\rm sub}$ with events in the two other low multiplicity bins, 10\leq\mbox{N^{\rm ch}}<20 and 20\leq\mbox{N^{\rm ch}}<30, are within the window of $v_{2}=0.03$.

We further explore this closure test using kinematic ranges similar to the STAR experiment Huang (2019). Figure 23 in the Appendix show the two-particle $\Delta\phi$ correlation function in short (0<|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) range in various multiplicity bins in $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV from hijing. Charged hadrons in 0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c and $|\eta|<0.9$ are used for the two-particle correlations, and multiplicity is defined as the number of charged hadrons in \mbox{p_{T}}>0.2~{}{\rm GeV}/c and $|\eta|<0.9$. In this narrower pseudorapidity acceptance, the shape of two-particle $\Delta\phi$ correlation function in the lowest multiplicity bin (0\leq\mbox{N^{\rm ch}}<5) is different from that in higher multiplicity ranges, as was the case with a wider pseudorapidity acceptance ($|\eta|<2.5$). Another important thing to point out is that there is a clear peak shape at the near-side (\mbox{\Delta\phi}\approx 0) in long-range from jet correlations which is invisible in $\Delta\phi$ correlation functions with a larger $\Delta\eta$ gap (2<|\mbox{\Delta\eta}|<5) shown in Figure 22.

Figure 10 shows the $v_{22}$ from Fourier fits to the two-particle $\Delta\phi$ correlation in long-range (1<|\mbox{\Delta\eta}|<1.8), and the values are larger than those from the wider $\Delta\eta$ gap indicating a stronger nonflow effect with a smaller $\Delta\eta$ gap. The solid lines give $v_{22}^{\rm sub}$ using two different low multiplicity bins, 0\leq\mbox{N^{\rm ch}}<5 and 5\leq\mbox{N^{\rm ch}}<10. It is interesting that the $v_{22}^{\rm sub}$ using the lowest multiplicity bin is significantly different between the two methods. The positive $v_{22}^{\rm sub}$ with the Method 1 is due to the remaining jet correlation with a smaller $\Delta\eta$ gap in $\Delta\phi$ correlation functions of higher multiplicity bins which is barely seen in the $\Delta\phi$ correlations of the lowest multiplicity. Therefore, the scaled correlation function of the low multiplicity bin cannot describe the peak structure on the near-side, resulting in the positive $v_{22}^{\rm sub}$. The results from Method 2 show large negative values, and this is related to the different shape of short-range $\Delta\phi$ correlation function in the lowest multiplicity bin introducing a large scaling with jet yields. The $v_{22}^{\rm sub}$ with the next low multiplicity bin (5\leq\mbox{N^{\rm ch}}<10) show a better closure within the level of $|v_{2}|=0.03$.

Another way to categorize event activity in $p+\mathrm{Au}$ collisions is to use centrality defined with charged particle multiplicity in the Au-ion-going (backward) rapidity ($-5.0<\eta<-3.3$ in STAR and $-3.9<\eta<-3.1$ in PHENIX). Results from the PHENIX and STAR experiments in small systems are categorized in this manner. In $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV, the multiplicity correlation between mid and backward rapidity is weak, so the shape difference of the $\Delta\phi$ correlation functions seen in low multiplicity events defined at midrapidity may not appear in peripheral events defined at backward rapidity. Figure 24 in the Appendix shows two-particle $\Delta\phi$ correlation in short (|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) ranges in $p+p$ and various centrality bins of $p+\mathrm{Au}$ collisions from hijing. Charged hadrons in 0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c and $|\eta|<0.9$ are used for the correlation functions, and centrality is defined with charged hadrons in $-5.0<\eta<-3.3$. As discussed earlier, one thing notably different from the case of multiplicity categorized at midrapidity shown in Figure 23 is that the shape of two-particle $\Delta\phi$ correlation function is similar in $p+p$ and all centrality bins of $p+\mathrm{Au}$ collisions. At RHIC, selecting events based on multiplicity in the same kinematic range as the particles selected for correlations introduces a significant undesirable shape variation of correlation functions.

Figure 11 shows the $v_{22}$ from Fourier fits to the two particle $\Delta\phi$ correlation at long-range. Here, the nonflow effects become larger from peripheral events to central events. The lines are $v_{22}^{\rm sub}$ using events from $p+p$ and two different peripheral selections (50–85% and 85–100%) of $p+\mathrm{Au}$ events. Note that it is usually difficult to collect events of 85–100% centrality of $p+\mathrm{Au}$ collisions in experiments due to the difficulty of triggering on events with small multiplicity at forward and backward rapidity, so the next peripheral bin (50–85%) is also used for the nonflow subtraction. The $v_{22}^{\rm sub}$ from both methods with all three selections of low multiplicity events are within the level of $|v_{2}|=0.03$ shown as a gray band. This smaller $v_{22}^{\rm sub}$ with centrality compared to the results using multiplicity at midrapidity is mainly coming from the similar shape of the $\Delta\phi$ correlation functions from peripheral to central $p+\mathrm{Au}$ collisions.

As an additional test in $p+\mathrm{Au}$ collisions from hijing, we have checked the $p_{T}^{\rm Trig}$ dependence in 0–5% central collisions. Figure 25 in the Appendix shows $\Delta\phi$ correlations at short (|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) ranges for different $p_{T}^{\rm Trig}$ bins. Charged hadrons in $|\eta|<0.9$ are used for the two-particle correlation, and the $p_{T}$ range of associated particles is 0.2<\mbox{p_{T}^{\rm Assoc}}<3~{}{\rm GeV}/c. Centrality is defined with charged hadrons in $-5.0<\eta<-3.3$ (Au-going direction). The shape of the correlation function in 0.2<\mbox{p_{T}^{\rm Trig}}<0.5~{}{\rm GeV}/c is different from other $p_{T}^{\rm Trig}$ bins both in $p+p$ and 0–5% $p+\mathrm{Au}$ collisions, but the shapes in $p+p$ and 0–5% $p+\mathrm{Au}$ events at the same $p_{T}^{\rm Trig}$ bin look comparable.

Figure 12 shows the $v_{22}$ from Fourier fit to the long-range $\Delta\phi$ correlations as a function of $p_{T}^{\rm Trig}$ in 0–5% central $p+\mathrm{Au}$ collisions, and the lines represent $v_{22}^{\rm sub}$ with events in $p+p$ and two peripheral centrality bins of $p+\mathrm{Au}$ collisions. Note that the $y$-axis range is much larger than the previous plots where the $p_{T}$-integrated results are shown as a function of multiplicity and centrality, due to a much stronger nonflow contribution in the higher $p_{T}^{\rm Trig}$ range than that in the $p_{T}$-integrated case. $v_{22}^{\rm sub}$ shows a weak dependence on the selection of low multiplicity events for both methods. However, the results from Method 1 show a clear $p_{T}^{\rm Trig}$ dependence and large negative values indicating a significant over-subtraction similar to the case of $p+p$ collisions at $\sqrt{s}$ = 200 GeV shown in Fig. 8. Method 2 shows a smaller deviation from zero than Method 1, but the $v_{22}^{\rm sub}$ values are still larger than 0.001.

Figure 13 shows the $v_{33}$ from the Fourier fit and the subtraction methods. Unlike the $v_{22}$ case, $v_{33}$ from the Fourier fit yields negative values. With the subtraction procedures applied, $v_{33}^{\rm sub}$ from both methods yields positive values. Thus the methods result in an over subtraction, but with a negative sign, hence leading to a strong $p_{T}$-dependent enhancement.

The PHENIX results are measured via multiple detector systems covering a wide range in pseudorapidity. However, the midrapidity coverage is limited ($|\eta|<0.35$) and thus the short-range correlations used in Method 2 are not available. The PHENIX flow measurements generally use the event plane method, but there are comparisons using three sets of two-particle correlations, where one can then algebraically solve for the anisotropy at midrapidity. In principle one can apply Method 1 to each of the three sets of two-particle correlations. We reserve this more detailed study to a future publication.

To summarize the closure test in $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV from hijing:

There are very substantial jet shape modifications when selecting events based on charged particle multiplicity around midrapidity $N^{\rm ch}$. These lead to significant distortions in the nonflow subtraction in both Method 1 and 2 and for a larger and smaller acceptance. 2. 2.

Selection of event categories based on forward detectors away from midrapidity, as employed by PHENIX and STAR, significantly improve the results of the hijing closure test. 3. 3.

The $v_{22}^{\rm sub}$ in integrated $p_{T}$ as a function of centrality with the two methods have some dependence on the low multiplicity reference selection, but are generally within the $|v_{2}|<0.03$ level in both Method 1 and 2. 4. 4.

In contrast the $p_{T}$ dependent results indicate a significant over-subtraction of nonflow in both Method 1 and 2. The over-subtraction is very large in Method 1 for \mbox{p_{T}}>1~{}{\rm GeV}/c. 5. 5.

In case of the $v_{33}^{\rm sub}$ as a function of $p_{T}$, there is an indication of significant bias to increase $v_{33}^{\rm sub}$ in both Method 1 and 2.

IV Further Tests with ampt

The above tests have a significant limitation in that one is testing procedures to disentangle flow and nonflow on models that have only nonflow. ampt Lin et al. (2005) is a model that includes both contributions and can thus be further elucidating. We follow the previous ampt studies in Refs. Aidala et al. (2017) to obtain the “truth” flow with respect to the participant plane calculated by using initial-state coordinates of partons resulting from string melting. More details on the method to calculate the truth flow can be found in Ref. Aidala et al. (2017). We highlight that the finite number of partons used to define the geometry and possible factorization breaking make this a rough estimate of the truth. A key item to note is that previous ampt studies Nagle and Orjuela Koop (2019) indicate that flow and nonflow do not factorize, in part because partons from jets can rescatter with medium partons. For both of these reasons, we are not expecting a perfect closure but rather examining possible trends in the subtraction procedure.

For this test, we used ampt v2.26 with string melting, and a few important parameters are listed in Table 1. Note that the updated string parameters, PARJ(41) and PARJ(42), introduced in Ref. Lin (2014) have been found to be important in recovering the peak structure from jets in near-side short-range correlation functions.

We follow the same procedure applied above with hijing $p+\mathrm{Au}$ and $p+p$ events. Two-particle $\Delta\phi$ and $\Delta\eta$ correlations are made from charged hadrons from $6\times 10^{8}$ $p+\mathrm{Au}$ and $2\times 10^{8}$ $p+p$ ampt events. Figure 26 in the Appendix shows the two-dimensional correlation functions at short (0<|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) range and various multiplicity bins in $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV from ampt. Charged hadrons in 0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c and $|\eta|<0.9$ are used for the $\Delta\phi$ correlation functions, and multiplicity is defined as the number of charged hadrons in \mbox{p_{T}}>0.2~{}{\rm GeV}/c and $|\eta|<0.9$. One thing to point out from the comparison of correlations between ampt and hijing (shown in Fig. 23) is that there is a more pronounced peak structure at the near-side (\mbox{\Delta\phi}\approx 0) in the long-range correlations which has contributions from the truth flow in ampt. Another difference is that the shape of the short-range $\Delta\phi$ correlation function in the lowest multiplicity bin (0\leq\mbox{N^{\rm ch}}<5) is similar to the shape in higher multiplicity bins in ampt events so that the nonflow subtraction results may not depend as much on the selection of the low multiplicity bin.

Figure 14 shows the $v_{22}$ from Fourier fits to the two particle $\Delta\phi$ correlation functions in long range as a function of multiplicity at midrapidity. As with the previous studies, the solid lines are the $v_{22}^{\rm sub}$ with two different selections of low multiplicity events, and the dashed line is $v_{22}$ calculated as a product of truth flow ($v_{2}^{2}\{\rm PP\}$) with an assumption of factorization. The $v_{22}^{\rm sub}$ with both methods and $v_{2}^{2}\{\rm PP\}$ are close to zero for low multiplicity events and increase with multiplicity, but the $v_{22}^{\rm sub}$ is much larger than $v_{2}^{2}\{\rm PP\}$. One caveat is $v_{22}$, which is a root-mean-square of $v_{2}$ ($\langle v_{2}^{2}\rangle$), is higher than a product of average $v_{2}$ (${\langle v_{2}\rangle}^{2}$). We have estimated the difference via the relation $\langle v_{2}^{2}\rangle/{\langle v_{2}\rangle}^{2}\propto\langle\varepsilon_{2}^{2}\rangle/{\langle\varepsilon_{2}\rangle}^{2}$ is $\approx 20\%$, which is much smaller than the observed difference.

Another set of two-particle $\Delta\phi$ correlations in short (|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) ranges with charged hadrons in the same kinematic range but different event multiplicity categorization defined at backward rapidity are shown in Fig. 27 in the Appendix. From $p+p$ collisions to various centrality ranges of $p+\mathrm{Au}$ collisions, the shape of the correlations are quite similar except for the near-side peak structure at long range, possibly due to the truth flow in ampt. Figure 15 shows the $v_{22}$ from Fourier fits in long-range as a function of centrality. The solid lines are the $v_{22}^{\rm sub}$ with two different selections of low multiplicity events, and the dashed lines represent $v_{22}$ from the truth flow, $v_{2}^{2}\{\rm PP\}$. The $v_{22}^{\rm sub}$ with the Method 1 in central $p+\mathrm{Au}$ collisions is larger than the $v_{2}^{2}\{\rm PP\}$ like the $v_{22}^{\rm sub}$ as a function of multiplicity at midrapidity. The $v_{22}^{\rm sub}$ using Method 2 is smaller than the $v_{22}^{\rm sub}$ from Method 1 and is consistent with the $v_{2}^{2}\{\rm PP\}$.

An additional test with ampt events is done for different $p_{T}^{\rm Trig}$ bins in 0–5% central $p+\mathrm{Au}$ collisions. Figure 28 in the Appendix shows the two-particle $\Delta\phi$ correlations in short (|\mbox{\Delta\eta}|<0.5) and long (1<|\mbox{\Delta\eta}|<1.8) range for various $p_{T}^{\rm Trig}$ bins in $p+p$ and 0–5% central $p+\mathrm{Au}$ collisions from ampt. The correlations are made from charged hadrons in $|\eta|<0.9$, and the $p_{T}$ range of associated particles is 0.2<\mbox{p_{T}^{\rm Assoc}}<3~{}{\rm GeV}/c. Centrality is defined with charged hadron multiplicity in $-5.0<\eta<-3.3$ in the direction of the Au ion. As the case of hijing $p+p$ and $p+\mathrm{Au}$ events, the shape of the two-particle $\Delta\phi$ correlations in the lowest $p_{T}^{\rm Trig}$ bin is different from that in other $p_{T}^{\rm Trig}$ bins both in $p+p$ and 0–5% $p+\mathrm{Au}$ collisions, but the shape at the same $p_{T}^{\rm Trig}$ bin of $p+p$ events and 0–5% $p+\mathrm{Au}$ events looks comparable.

Figure 16 shows the $v_{22}$ from Fourier fits to the long-range two-particle $\Delta\phi$ correlation as a function of $p_{T}^{\rm Trig}$ in 0–5% central $p+\mathrm{Au}$ collisions. The solid lines are $v_{22}^{\rm sub}$ with low multiplicity selections from $p+p$ and two peripheral $p+\mathrm{Au}$ centrality bins, and the dashed line represents the $v_{22}$ from the truth flow, $v_{2}^{\rm Trig}\{{\rm PP}\}\times v_{2}^{\rm Assoc}\{{\rm PP}\}$. Similar to the case of hijing $p+\mathrm{Au}$ events shown in Fig. 12, the nonflow effects increase with $p_{T}^{\rm Trig}$. The $v_{22}^{\rm sub}$ with both methods are below zero in \mbox{p_{T}^{\rm Trig}}>2.5~{}{\rm GeV}/c using Method 1 and \mbox{p_{T}^{\rm Trig}}>1~{}{\rm GeV}/c using Method 2. This clearly indicates an over-subtraction in these $p_{T}^{\rm Trig}$ ranges, because the $v_{22}$ from truth flow in ampt should give positive values even though the truth $v_{22}$ from two-particle correlation may not be exactly same as $v_{2}^{\rm Trig}\{{\rm PP}\}\times v_{2}^{\rm Assoc}\{{\rm PP}\}$, shown as the dashed line. The deviation of the $v_{22}^{\rm sub}$ from zero is larger with Method 2 which is opposite to the closure test results with hijing events shown in Fig. 12. This is possibly due to the different shape of two particle $\Delta\phi$ correlation function between hijing and ampt both at short and long ranges presented in Figs. 25 and 28.

Figure 17 shows the $v_{33}$ from the Fourier fit and the subtraction methods. Note that the $v_{33}$ from the truth flow, $v_{3}^{\rm Trig}\{{\rm PP}\}\times v_{3}^{\rm Assoc}\{{\rm PP}\}$, are non zero but have very small values. Similar to the hijing study, $v_{33}$ from the Fourier fit yields negative values, but the magnitude is much smaller in ampt, possibly due to the truth flow contribution in ampt. In the subtraction results, $v_{33}^{\rm sub}$ from both methods yields positive values that are significantly larger than the truth flow. As was the case in the hijing study, there is an over-subtraction, again with a negative sign, and thus a strong $p_{T}$-dependent enhancement.

The summary of the ampt studies is as follows:

Selection of $p+\mathrm{Au}$ event categories based on midrapidity multiplicity $N^{\rm ch}$ result in significant undersubtraction of nonflow contributions in ampt. 2. 2.

Selection of $p+\mathrm{Au}$ event categories based on forward detectors away from midrapidity, as employed by PHENIX and STAR, improve the results of the ampt test. The results as a function of centrality for the $v_{22}^{\rm sub}$ values in the two methods match the truth extracted ampt flow qualitatively. 3. 3.

In contrast the $p_{T}$ dependent results indicate a significant over-subtraction of nonflow in both Method 1 and 2. The over-subtraction is very large in Method 2 for \mbox{p_{T}}>1~{}{\rm GeV}/c. 4. 4.

In case of the $v_{33}^{\rm sub}$ as a function of $p_{T}$, there is an indication of significant bias to increase $v_{33}^{\rm sub}$ in both Method 1 and 2.

V Conclusions

We examined closure tests of the nonflow subtraction methods used with two-particle correlations developed by ATLAS and CMS Collaborations to study the collective behavior of particle production in small collision systems at the LHC. Monte Carlo event generators (pythia and hijing) including no collective flow are use to quantify a level of closure. In the test for $p+p$ collisions at $\sqrt{s}$ = 13 TeV at the LHC where a large pseudorapidity gap (|\mbox{\Delta\eta}|>2) can be applied, the shape of the $\Delta\phi$ correlation functions in short and long range is relatively stable throughout the event multiplicity and $p_{T}$ ranges for which the resulting second-order Fourier coefficients after nonflow subtraction ($v_{22}^{\rm sub}$) are less than 0.001 in these ranges.

We also tested the nonflow subtraction methods for $p+p$ collisions at $\sqrt{s}$ = 200 GeV, and it is observed that the shape of the $\Delta\phi$ correlation functions changes significantly with event multiplicity and $p_{T}$ ranges. Therefore, the nonflow subtraction results are very sensitive to the selection of the low multiplicity reference used for the subtraction procedure, indicating that one should be extremely careful applying the same nonflow subtraction techniques to the $p+p$ data at RHIC in future analyses.

In the test with hijing for $p+\mathrm{Au}$ collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV, the large sensitivity to the low multiplicity event selection that was observed in $p+p$ collisions at $\sqrt{s}$ = 200 GeV is still apparent when event activity is categorized at the same rapidity range as that of the two-particle correlations. The test with a smaller $\Delta\eta$ gap (1<|\mbox{\Delta\eta}|<1.8) shows a significant difference between the methods, because it starts to be affected by short-range jet correlation contributions. Although the event activity categorization at backward rapidity helps to reduce the dependence of the low multiplicity event selection on the results of integrated $p_{T}$ (0.2<\mbox{p_{T}}<3~{}{\rm GeV}/c), the test results as a function of $p_{T}$ show a clear over-subtraction in \mbox{p_{T}}>1~{}{\rm GeV}/c with both methods for $v_{2}$ and $v_{3}$.

We extended this study with ampt, which includes both flow and nonflow effects, and the nonflow subtraction results were compared to the truth flow with respect to the participant plane. The $v_{22}^{\rm sub}$ and $v_{33}^{\rm sub}$ are inconsistent with the truth flow in ampt, and the $p_{T}$ dependent results show negative values at higher $p_{T}$ with both methods. Thus in $p+\mathrm{Au}$ collisions, both hijing and ampt studies indicate that both Method 1 and 2 result in an over-subtraction of nonflow and thus an underestimate of real elliptic flow and overestimate of real triangular flow.

Acknowledgments

We acknowledge Jiangyong Jia and Wei Li for useful discussions and a careful reading of the manuscript. We acknowledge Julia Velkovska for useful discussions and the suggestion to extend this study to third order flow coefficients. SHL, KKH, JLN acknowledge support from the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Contract No. DE-FG02-00ER41152. QH and DVP acknowledge support from the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Contract No. DE-SC0018117.

Appendix A Two-particle $\Delta\phi$ correlation functions

Here we present two-particle $\Delta\phi$ correlation functions described as Eq. (1) in Sec. I with charged hadrons from event generators, and these correlation functions are used for $v_{22}$ and $v_{22}^{\rm sub}$ calculations (see Figs. 18–28). We specifically follow the formula for per-trigger yields as a function of $\Delta\phi$ within a certain $\Delta\eta$ range described in Ref. Khachatryan et al. (2017),

[TABLE]

where $\Delta\eta$ and $\Delta\phi$ are the differences of $\eta$ and $\phi$ between trigger and associated particles. S(\mbox{\Delta\eta},\mbox{\Delta\phi}) and B(\mbox{\Delta\eta},\mbox{\Delta\phi}) are the yields of pairs normalized by the number of trigger particles, $N^{\mathrm{Trig}}$, in the same and mixed events respectively,

[TABLE]

[TABLE]

The purpose of normalizing by the mixed event pair distributions is to account for the geometric $\Delta\eta$-dependent pair acceptance effect.