Azimuthal correlation anisotropies in p + p collisions simulated using Pythia

One of the most interesting aspects of heavy ion collisions is the azimuthal angle distribution of particles in the final state. This is because the overlap geometry of these collisions is anisotropic, similar to an almond shape. This anisotropic geometry is imprinted on the final state particle momentum distribution because of interactions among the particles before they freeze out [1]. These interactions drive the collective expansion of the overlap collision system, a phenomenon called collective flow [2]. A comparison between the momentum anisotropy and initial geometry anisotropy can provide information about the interactions in a hot and dense medium, the quark-gluon plasma (QGP), created in these collisions.

To study azimuthal anisotropies, we often express particle distributions in Fourier harmonic series in terms of the particle azimuthal angle (ϕ) relative to the reaction plane (Ψ) [3]:

The coefficients $ v_n $ characterize the harmonic anisotropies. The most straightforward way to calculate the anisotropies is through two-particle correlations [4]:

where $ \Delta\phi = \phi_1-\phi_2 $ is the pair azimuthal opening angle. The Fourier coefficients characterizing the two-particle cumulants are given by

If the azimuthal anisotropy in $ \Delta\phi $ results all from global correlations of all particles with respect to the reaction plane (i.e. flow), then $ v_n = \sqrt{V_n} $.

However, two-particle correlations can result from processes such as resonance decays and jets. These few-body correlations contribute to $ V_n $ and are often referred to as nonflow [5−7], in contrast to the global flow correlations. Thus, $ V_n $ contains both flow and nonflow contributions. Because nonflow arises from two- and few-body particle correlations, their effect is diluted by particle multiplicity, resulting in a typical inverse multiplicity dependence in the two-particle cumulants,

Because of the large multiplicity in mid-central and central heavy-ion collisions, nonflow contamination is minor compared with collective flow [8]. Nonflow can dominate in peripheral collisions where multiplicity is small.

In studies on heavy ion collisions, proton-proton ($ p+p $) and proton-nucleus ($ p+A $) collisions are often taken as references to infer the properties of the QGP or nuclear effects in hot quantum chromodynamics (QCD). Researchers have long considered that the particle azimuthal distributions in those small-system ($ p+p $ and $ p+A $) collisions should be uniform (no preferred azimuthal direction exists), and the two-particle anisotropy $ V_n $ in these small-system collisions must all be nonflow. However, recent experiments suggest that collective flow may exist in $ p+A $ and $ p+p $ collisions, particularly in high-multiplicity events [9−21]; a strong indication is that the four- and multi-particle cumulants, essentially devoid of nonflow, are finite [22−24]. In two-particle cumulant measurements of Eq. (4), the nonflow effect can be severe because the event multiplicity, and thus the multiplicity dilution, is small. This is the case even in high-multiplicity events of these small-system collisions, where the multiplicity is comparable to that in peripheral heavy ion collisions. Many nonflow mitigation methods have been used [8, 25, 26]. In small-system collisions, particularly at LHC energies, finite positive $ v_n^2 $ have been observed to remain after nonflow mitigation [26]. This led to the conclusion that a small QGP droplet may be formed in small-system collisions [26].

On the other hand, the overlap geometry in small-system collisions is lumpy owing to fluctuations, and final-state momentum anisotropy can result if there is some degree of interaction [27−30], e.g., multi-parton interactions [31−34], in these collisions. The magnitude of these interactions, gauged by the magnitude of the final-state momentum anisotropy relative to the initial geometry anisotropy, would then discriminate whether the interactions occur in a low-density environment or high-density QGP.

In addition to these physics problems, the technical question of how well nonflow has been mitigated is not settled and is under extensive investigation. Therefore, anisotropies in model simulations of $ p+p $ collisions are interesting to examine, particularly using the Pythia event generator [35−38]. Pythia is naively not expected to generate single-particle azimuthal anisotropy because it is not a cascade model with final-state interactions. However, a recent implementation of multi-parton interactions (color reconnection) was perceived to mimic some degree of final-state interactions; thus, we cannot completely exclude collective flow [36, 39, 40]. Therefore, studies on anisotropies using Pythia would be valuable in at least two fronts. One is to aid in studying nonflow subtraction techniques and associated uncertainties [8, 41, 42], and the other is to discern any collective flow in $ p+p $ interactions to guide the interpretations of experimental anisotropy measurements in small-system collisions [25, 26].

We hereby report a study on azimuthal correlation anisotropies in $ p+p $ interactions simulated using Pythia version 8.303 [38], which includes multi-parton interactions and color reconnection [39, 40].

We have used Pythia 8.303 [38] to simulate $ p+p $ collisions at the center-of-mass energy of $\sqrt{s_{NN}} = 200$ GeV. We generate $ 3\times10^8 $ events each for diffractive and non-diffractive interactions with soft processes (SoftQCD) [36]. The soft processes in Pythia represent the entire cross-section ($ \sigma_{\rm{tot}} $) of the hadron-hadron interaction, i.e., minimum bias (MB) interaction. Rare processes, such as jet physics at all scales, occur as part of the total cross-section. For $ p+p $ collisions in Pythia, the cross section is calculated using the 1992 Donnachie-Landshoff parameterization [36]. This total cross-section is the sum of several components: $ \sigma_{\rm{tot}} = \sigma_{\rm{sd}} + \sigma_{\rm{dd}} + \sigma_{\rm{nd}} + \sigma_{\rm{el}} $ [43]. The diffractive part is divided into the single (sd) and double (dd) components, where the former has the$ pp \rightarrow pX $ processes, whereas the latter has the $ pp \rightarrow ppX $ processes. The diffractive events in our simulation include both single-diffractive and double-diffractive processes. The non-diffractive (nd) part includes all other inelastic processes $ pp \rightarrow X $. Because the elastic (el) process, $ pp \rightarrow pp $, does not produce particles, it is not included in our simulation. The comparison between diffractive and non-diffractive interactions is expected to provide insight into initial- vs. final-state effects, because diffractive interactions should be more closely related to initial state of the proton, whereas non-diffractive interactions could be influenced by final-state effects, if any.

Table 1 lists the majority of particle species and their average compositions in non-diffractive events. In our simulations, we keep the electromagnetic- and weak-decay hadrons stable.

Particle	Multiplicity	Particle	Multiplicity
$\pi^+$	12.140	$\pi^-$	11.422
$\pi^0$	13.289
$K^+$	1.631	$K^-$	1.459
$K^0$	2.263	$\bar{K}^0$	1.059
p	4.162	$\bar{p}$	0.700
n	1.768	$\bar{n}$	0.699
Λ	0.370	$\bar{\Lambda}$	0.217

Table 1.  Average full-phase space multiplicities of various particle species in non-diffractive events generated using Pythia 8.303. A total of approximately $3 \times 10^8$ events are generated.

We include only $ \pi^\pm, K^\pm, p, \bar{p} $ in our study, referred to as particles of interest (POI). Because experiments (such as STAR [44] at RHIC and ALICE [45] at the LHC) typically measure charged particles in the limited midrapidity region, we focus on charged hadrons within the pseudorapidity range $ |\eta|<1 $ and transverse momentum range $ 0.2< p_\perp<2 $ GeV/c. Figure 1(a) shows the diffractive and non-diffractive event probability distributions in POI multiplicity within $ |\eta|<1 $ and $ 0.2< p_\perp<2 $ GeV/c, referred to as $ N_{{\rm{ch}}} $ hereinafter. As reference, the distribution in the total number of final-state particles ($ N_{{\rm{tot}}} $), including both charged and neutral particles over all phase space, is shown in Fig. 1(b). Diffractive events have fewer produced particles than non-diffractive events, as expected.

Figure 1.  (color online) Probability distributions of diffractive and non-diffractive $p+p$ collisions simulated using Pythia 8.303 as a function of (a) $N_{\rm{ch}}$, the midrapidity (within $|\eta|<1$ and $0.2<{p_ \bot }<2$ GeV/c) charged hadron ($\pi^\pm,K^\pm,p,\bar{p}$) multiplicity, and (b) $N_{\rm{tot}}$, the total final-state particle multiplicity (including all charged and neutral particles, over $4\pi$ solid angle).

Events of differing $ N_{{\rm{ch}}} $ are expected to have varying physics, particularly in non-diffractive events, influences of final-state effects such as multi-parton interactions and color reconnection on possibly collective behavior. Thus, we divide the events into the following $ N_{{\rm{ch}}} $ bins: 1–3, 4–6, 7–10, 11–14, 15–20. We study anisotropies of POI in $ N_{{\rm{ch}}} $ bins and present our results as functions of the mean multiplicity $ \langle { N_{{\rm{ch}}}} \rangle $ of each $ N_{{\rm{ch}}} $ bin.

Figure 2 shows $ \sqrt{V_n} $ as a function of the total multiplicity $ N_{{\rm{tot}}} $ for both diffractive and non-diffractive events. We observe that $ \sqrt{V_2} $ decreases with $ N_{{\rm{tot}}} $ and appears to saturate at large $ N_{{\rm{tot}}} $. Motivated by Eq. (4), we fit the data using $ \sqrt{p_0^2+p_1/ N_{{\rm{tot}}}} $. The fitted $ p_0 $ parameter is non-zero for both diffractive and non-diffractive events, suggesting that both types of events may have flow-like global anisotropies. The magnitudes of the $ p_0 $ parameter are comparable between diffractive and non-diffractive events. This is rather surprising given the vastly different physics processes in diffractive and non-diffractive interactions.

Figure 2.  (color online) Two-particle cumulant parameter, $\sqrt{V_2}$, as a function of the total multiplicity $N_{\rm tot}$. Curves are fit using $\sqrt{V_2}=\sqrt{p_0^2+p_1/N_{\rm tot}}$, where $p_0$ and $p_1$ are fit parameters.

Note that the fit quality is poor, indicating that the simple relation in Eq. (4) may not hold. For example, the nonflow correlation strength such as jet correlations may increase faster than linearly in multiplicity, which would result in a slower decrease of $ V_n $ than $ 1/ N_{{\rm{tot}}} $. The global flow-like correlation of the underlying event may also change with the event multiplicity such that a constant $ p_0 $ parameter may not be adequate. It is clear that a more rigorous investigation is required to identify the structure of $ p+p $ events from Pythia simulations.

Two-particle $ ( \Delta\eta, \Delta\phi) $ correlations are widely used to identify sources of correlations, where $ \Delta\eta\equiv\eta_1-\eta_2 $ and $ \Delta\phi\equiv\phi_1-\phi_2 $ are the two-particle distances in the pseudorapidity and azimuthal angle [46]. For example, resonance decays manifest as a near-side peak at $ ( \Delta\eta, \Delta\phi)\approx (0,0) $; so do intra-jet correlations. The effects of global momentum conservation and dijet correlations manifest as away-side ridge peaked at $ \Delta\phi = \pi $ and roughly uniform in $ \Delta\eta $. Global momentum conservation may be modeled using a negative $ V_1 = \langle {\cos \Delta\phi} \rangle $ term [47], and others are often modeled using Gaussians. In our simulation of soft QCD interactions (i.e. MB events), jet contributions are minimal. However, high-multiplicity events are known to be biased towards stronger jet production [48].

Figure 3 shows $ ( \Delta\eta, \Delta\phi) $ correlations of POI from diffractive and non-diffractive events. The two-particle acceptance, approximately a triangle in $ \Delta\eta $ as discussed below, has not been corrected in these plots. The shapes are similar between diffractive and non-diffractive events, whereas the magnitude is significantly smaller in the former because of the smaller number of particles produced in diffractive events.

Figure 3.  (color online) Two-particle ($\Delta \eta $, $\Delta \phi$) distributions from diffractive (a) and non-diffractive (b) interactions. The two-particle ''triangle'' acceptance in $\Delta \eta $ is not corrected. The z axis (color coordinate) is the number of pairs per unit of pseudorapidity per radian.

Before continuing, we want to correct for an important aspect, namely, the two-particle acceptance. Considering a particle in a detector of limited acceptance in η, a second particle a distance $ | \Delta\eta| $ away may fall outside the acceptance on one side. Correcting for such a loss in acceptance would be helpful. This acceptance is approximately a triangle in $ \Delta\eta $ (for approximately uniform single-particle distribution in η), and is often referred to as the ''triangle'' acceptance [49]. We have assumed solenoidal detectors of full acceptance in ϕ (which is important for azimuthal anisotropy measurements), and because of the azimuthal periodicity, the two-particle $ \Delta\phi $ acceptance is uniform. We obtain the two-particle acceptance factor using the mixed-event technique, taking one particle from one event and the other from another event of the same multiplicity value as the first event. The mixed-event correlations are obtained separately for diffractive and non-diffractive events. The mixed-event correlation function is normalized such that the magnitude is unity at $ \Delta\eta = 0 $ because, with one particle within the η acceptance (in our case $ |\eta|<1 $), the probability for the second particle to be inside the acceptance is 100% if $ \Delta\eta = 0 $. The normalized mixed-event correlation is divided to obtain the acceptance-corrected two-particle correlations.

Figure 4 shows the ''triangle'' acceptance corrected $ ( \Delta\eta, \Delta\phi) $ correlations in non-diffractive interactions for three multiplicity ranges. The correlations feature a near-side peak $ ( \Delta\eta\approx0, \Delta\phi\approx0) $ (most prominent in high multiplicity events) and a broad away-side ridge at $ \Delta\phi\approx\pi $. These features are characteristic of intra-jet and inter-jet (dijet) correlations, as well as resonance decays and global momentum conservation. Correlations in diffractive interactions are similar. Note that the near-side peak is not visible in low-multiplicity events. However, resonance decays (and similarly jet correlations) must also be present in these events. As shown later, the negative $ V_1 $ component is large in these low-multiplicity events, such that the near-side correlation peak is submerged by the negative dip [14].

Figure 4.  (color online) Two-particle correlations in non-diffractive interactions, as in Fig. 3(b), in various multiplicity ranges of $1\leq N_{\rm{ch}}\leq3$ (a), $4\leq N_{\rm{ch}}\leq6$ (b), and $15\leq N_{\rm{ch}}\leq20$ (c). The two-particle “triangle” acceptance in $\Delta \eta $ is corrected by mixed events of the corresponding multiplicity range. The z axis (color coordinate) is the number of pairs per unit of pseudorapidity per radian. For example, the $\Delta \eta $-acceptance corrected total number of pairs is approximately 20 in panel b), giving an average amplitude of approximately 0.8.

Far more resonances decay into daughters with an opposite-sign (OS) charge pair than a same-sign (SS) charge pair. Intra-jet correlated hadrons have more opposite-charge pairs than same-charge pairs because jets are mostly neutral gluon jets (i.e., charge-ordering phenomenon in jets). Thus, one has stronger correlation in OS pairs than in SS pairs, and we can gain more insights by analyzing OS and SS correlations separately.

Figure 5 shows the acceptance-corrected correlations, separately for OS and SS pairs, in non-diffractive $ p+p $ collisions. One multiplicity bin is shown as an example; the results are similar in shape for other multiplicity bins. The near-side correlations are more prominent in OS correlations than in SS correlations, reflecting, presumably, resonance decays and charge ordering in (mini-)jet fragmentation. The away-side correlations ($ \Delta\phi\approx\pi $) are similar between OS and SS correlations, indicating no charge-dependent correlations on the away side.

Figure 5.  (color online) Acceptance-corrected two-particle $(\Delta \eta ,\Delta \phi)$ correlations in non-diffractive collisions with $4\leq N_{\rm{ch}}\leq6$ as in Fig. 4(b) but separately for OS pairs (a) and SS pairs (b). The correlations are similar in shape for other multiplicity bins and diffractive events.

Figure 6 shows the $ \Delta\eta $ and $ \Delta\phi $ projections of the OS and SS 2D correlation functions for non-diffractive events with $ 4\leq N_{{\rm{ch}}}\leq6 $ from Fig. 5. Other multiplicity bins are similar in the correlation shapes. The $ \Delta\eta $ correlation is characteristic of a Gaussian centered at $ \Delta\eta = 0 $. The $ \Delta\phi $ correlations are dominated by the away-side peak at $ \Delta\phi = \pm\pi $, with essentially no near-side peak at $ \Delta\phi = 0 $ visible at the scale of the plots because of the aforementioned large negative dipole.

Figure 6.  (color online) Projections of the 2D correlation functions from Fig. 5 onto $\Delta \eta $ (a,c) and $\Delta \phi$ (b,d) for OS (a,b) and SS (c,d) pairs in non-diffractive events with $4\leq N_{\rm{ch}}\leq6$. The correlations are normalized per event and by bin size.

We know that the OS correlation contains more contributions from resonance decays and jets than the SS correlation because of contributions from neutral resonances and charge ordering in jet fragmentation. The number of OS and SS pairs from the underlying event also differ significantly, as shown in Fig. 6. Careful analysis is required to study these charge-dependent nonflow correlations. We leave it to future research. In this study, we are interested in the possible flow-like anisotropies in $ p+p $ events simulated using Pythia. Because the SS pairs are less contaminated by nonflow, we focus on SS correlations in the following sections.

Although to a less degree than the OS pairs, SS pairs are also affected by jet and dijet correlations. Intra-jet correlations are concentrated on the near side $ \Delta\phi\sim0 $. The dijet correlations are back-to-back and concentrated on the away side $ \Delta\phi\sim\pm\pi $. Because the underlying parton-parton interaction kinematics are not necessarily symmetric about midrapidity, the away-side correlations are approximately $ \Delta\eta $ independent. Other back-to-back correlations, such as global momentum conservation and possibly flow-like harmonics, are also mostly independent of $ \Delta\eta $. Thus, we may take the difference between small- and large-$ \Delta\eta $ correlations to assess the near-side nonflow correlations.

We define small $ \Delta \eta $ to be $ |\Delta \eta| < 0.5 $ and large $ \Delta \eta $ to be $ 1.5 < |\Delta \eta| < 2 $ (the upper limit is automatically satisfied because of our single-particle acceptance requirement of $ |\eta|<1 $). Note that we have defined an equal $ | \Delta\eta| $ size for the two regions, but the specific $ \Delta\eta $ values for the definitions are less important as long as the large-$ \Delta\eta $ region contains minimal near-side contribution. Figure 7(a) shows the $ \Delta\phi $ correlations at small and large $ \Delta\eta $ for the medium multiplicity range, as an example. The $ \Delta\eta $ acceptance has already been corrected in those correlations. The small- and large-$ \Delta\eta $ correlations are similar in amplitude, suggesting a more or less uniform two-particle distributions in $ \Delta\eta $. The correlations are fit with the functional form of two Gaussians atop a constant pedestal, with one Gaussian centered at $ \Delta\phi = 0 $ and the other at $ \Delta\phi = \pm\pi $. The fitted pedestal values are approximately equal between them. Figure 7(b) shows the ratio of the correlations between small and large $ \Delta\eta $ values. The ratio is fit with the same functional form as above, and the fitted pedestal value is $ 0.96\pm 0.06 $.

Figure 7.  (color online) (a) SS correlations in $\Delta \phi$ within small $|\Delta \eta |<0.5$ and large $1.5<|\Delta\eta|<2$ for non-diffractive events with $4\leq N_{\rm{ch}}\leq6$. The correlations are obtained from projections of the acceptance-corrected $(\Delta \eta ,\Delta \phi)$ correlations in Fig. 5 in corresponding $\Delta \eta $ ranges and are then divided by $\left\langle {{N_{{\rm{ch}}}}} \right\rangle $; therefore, they are SS dihadron correlations per hadron. Each superimposed curve is a fit to the corresponding data with the sum of a pedestal and two Gaussians, one centered at $\Delta \phi=0$ and the other at $\Delta \phi=\pm\pi$. The fitted pedestals are $0.0708\pm 0.0007$ and $0.0704\pm 0.0001$ for small and large $\Delta \eta $, respectively. The central Gaussian amplitude for the large-$\Delta \eta $ correlation is negative. (b) Ratio of small-$\Delta \eta $ correlation over a large-$\Delta \eta $ one. The superimposed curve is a fit to the same functional form as in panel a). The fitted pedestal value is $0.96\pm 0.06$.

Figure 8 shows the difference between small- and large-$ \Delta\eta $ dihadron correlations, with the large-$ \Delta\eta $ correlation first multiplied by the fitted pedestal to the ratio in Fig. 7(b). Because of this choice of normalization, the correlation difference is not necessarily zero outside the near-side region. If we normalized the small- and large-$ \Delta\eta $correlations by matching their away side, the correlation differences would all be averaged at zero outside the near-side region. For the low-multiplicity bin, there is no near-side peak because of the large negative dipole and likely minimal intra-jet correlations. The near-side peak begins to show up in the mid-multiplicity bin and is predominant for the high-multiplicity bin, suggesting strong contributions from intra-jet correlations. Resonance decay contributions are likely minimal in SS pairs.

Figure 8.  (color online) Differences between small- and large-$\Delta\eta$ dihadron correlations for three multiplicity ranges of $1\leq N_{\rm{ch}}\leq3$, $4\leq N_{\rm{ch}}\leq6$, and $15\leq N_{\rm{ch}}\leq20$, where the large-$\Delta\eta$ correlation is first scaled by a normalization factor. Those correlations for the medium multiplicity range are shown in Fig. 7(a). The normalization factor is taken to be the pedestal value of the corresponding ''pedestal + two-Gaussian'' fit, such as the one in the lower panel of Fig. 7(b); the pedestal values for $1\leq N_{\rm{ch}}\leq3$ and $15\leq N_{\rm{ch}}\leq20$ (that are not shown in Fig. 7) are $0.95 \pm 0.01$ and $0.96 \pm 0.07$, respectively.

A large-$ \Delta\eta $ is not significantly contaminated by intra-jet correlations; therefore, it is cleaner to extract flow anisotropies from large-$\Delta\eta $, although dijet correlations are still present on the away side. Figure 9 depicts the large-$ \Delta\eta $ per-trigger dihadron correlations for the three selected multiplicity bins. For easy display, the correlations are chopped off a constant pedestal such that the near-side region is zero for all three multiplicity ranges. The pedestal is obtained from a ''pedestal + Gaussian'' fit where the Gaussian describes the away-side correlation about $ \Delta\phi = \pm\pi $.

Figure 9.  (color online) SS per-trigger dihadron $\Delta \phi$ correlations at large $\Delta \eta $ ($1.5<|\Delta \eta |<2$) in non-diffractive p+p collisions simulated using Pythia 8.303 for three multiplicity ranges: $1 \leq N_{\rm{ch}} \leq 3$ (low mult.), $4\leq N_{\rm{ch}}\leq6$ (medium mult.), and $ 15 \leq N_{\rm{ch}} \leq $20 (high mult.). Each distribution is subtracted by the constant pedestal from the “pedestal + Gaussian” fit, similar to that in Fig. 7 for the medium-multiplicity range. The subtracted pedestals are $0.035 \pm 0.010$, $0.07 \pm 0.0001$, and $0.342 \pm 0.001$ for the three multiplicity ranges, respectively.

The correlation shapes are similar, with a prominent away-side peak and minimal near-side correlations. The away-side peak is primarily of a negative dipole shape ($ -\cos \Delta\phi $); this coresponds with the expectation from global momentum conservation [47] and with the expectation that back-to-back dijet correlations at low $ p_\perp $values are also primarily of a negative dipole shape. The away-side peak amplitude increases with the event multiplicity. This appears to suggest that the back-to-back dijet contribution increases with the event multiplicity, as global momentum conservation should not depend on multiplicity. This is likely owing to multiplicity bias–when high-multiplicity is demanded, we likely select events with relatively stronger jet contributions [50−52].

With the understanding that the SS correlations are composed of intra-jet and dijet correlations atop an underlying event with possible harmonic oscillations, we proceed to fit the mixed-event acceptance corrected $ ( \Delta\eta, \Delta\phi) $ correlations. We use the following fit function:

The first term on the r.h.s. of Eq. (5) describes the near-side 2D Gaussian in $ ( \Delta\eta, \Delta\phi) $; the repeated $ \Delta\phi $-Gaussian at $ \pm2\pi $ ensures azimuthal periodicity with those beyond $ 2\pi $ safely neglected. The second term describes the $ \Delta\phi $ ridge that is present as a Gaussian at $ \Delta\eta = 0 $; such a ridge could be the result of the aforementioned beam jet fragmentation [53]. The third term attempts to describe a feature present in the correlation functions that bends up towards large $ | \Delta\eta| $, which can be attributed to back-to-back emission of particles in η due to, e.g., local momentum conservation, enabling $ \eta_{1}+\eta_{2} $ to fluctuate about zero in the Gaussian, and integrating it within the acceptance boundaries of $ |\eta|<1 $results in the marginal distribution of an error function in$ \Delta\eta $. The fourth term is the underlying event contribution with possible harmonic oscillations; we keep only the first four harmonics in the fit, where the $ V_n $ $ (n = 1,2,3,4) $ quantify their amplitudes. Note that these harmonic anisotropy amplitudes are with respect to the underlying event multiplicity quantified by the baseline parameter C, excluding contributions from the nonflow (intra-jet and dijet) particles. We find that no additional structure is required on the away side to describe the away-side jet. This is because the away-side jet contribution, in our $ p_\perp $-inclusive study (dominated by low $ p_\perp $ particles), can presumably be described by a negative $ V_1 $ dipole and blended in the $ V_1 $ harmonic term.

We fit the $ ( \Delta\eta, \Delta\phi) $ SS correlation in each multiplicity bin using Eq. (5). Figure 10 shows, as an example, the fit result for the $ 7\leq N_{{\rm{ch}}}\leq10 $ multiplicity bin of the non-diffractive events. The $ \chi^2/{\rm{NDF}} $ values are 1.24, 1.54, 1.12, 1.06, and 0.99 for the multiplicity bins of 1–3, 4–6, 7–10, 11–14, and 15–20 of the non-diffractive events, respectively. For diffractive events, the $ \chi^2/{\rm{NDF}} $ values are 1.78, 1.80, 1.08, and 1.01 for the multiplicity bins of 1–3, 4–6, 7–10, and 11–14, respectively. The number of degrees of freedom (NDF) is 3988 for all correlation functions. We note that the $ \chi^2/{\rm{NDF}} $ values are not optimal, suggesting that the fit function of Eq. (5) is not ideal. Nevertheless, the fitted $ V_n $ parameters may still provide a reasonable description of the underlying modulation. The fit quality can be improved by introducing more components (parameters) into the fit function; however, such improvements may be at the expense of clarity.

Figure 10.  (color online) Example of the 2D model fit by Eq. (5) to the SS pair correlations for the multiplicity range $11\leq N_{\rm{ch}}\leq14$.

We are most interested in the $ V_n $ parameters. Multiple methods can be used to obtain $ V_n $, whose physical interpretation depends on the method used [8]. It can be obtained from the two-particle cumulant using the Fourier coefficient of Eq. (3) and the fitting method of Eq. (5). Figure 11 shows $ V_2 $ obtained using various methods from the SS pair correlations for a selected multiplicity bin of diffractive events and two multiplicity bins of non-diffractive events: the leftmost data point is the Fourier coefficient of SS correlations within $ | \Delta\eta|<0.5 $, which provides a large value because of nonflow contribution from the near-side peak, and the second data point from the left is that without any $ \Delta\eta $ cut; the third point is that with $ | \Delta\eta|>1.0 $, and the fourth one is that with $ | \Delta\eta|>1.5 $, both of which yield a significantly reduced value because of the exclusion of the main part of nonflow at small $ | \Delta\eta| $. The rightmost data point is the $ V_2 $ obtained from the 2D fit, which is the smallest for the results shown, suggesting a large removal of nonflow effects.

Figure 11.  (color online) $V_2$values obtained from SS correlation functions using various methods. The results shown are for diffractive events of one multiplicity bin and non-diffractive events of two multiplicity bins. The rightmost set of data points are from 2D fits. The rest are Fourier coefficients of correlation functions with various $\Delta \eta $ cuts.

Figure 12 shows $ V_1 $, $ V_2 $, and $ V_3 $ values obtained from the Fourier coefficients of the SS correlations at $ | \Delta\eta|>1.5 $ (such as the one in Fig. 7(a)) and from the 2D fit to the SS $ ( \Delta\eta, \Delta\phi) $ correlations. These two methods should have the smallest nonflow contamination. The values are plotted as functions of $ N_{{\rm{ch}}} $ for both non-diffractive (open black points) and diffractive events (filled red points). We can infer from Fig. 12(a) that the negative dipole $ V_1 $ amplitude decreases with increasing multiplicity. However, the multiplicity dependence is not $ 1/ N_{{\rm{ch}}} $, which would be characteristic of global momentum conservation [47]. Thus, the results suggest that the $ V_1 $ component in these $ p+p $ interactions is not of the nature of global momentum conservation. Note that the $ V_1 $ amplitude is larger in diffractive events than in non-diffractive events at similar multiplicities. The $ V_1 $ amplitude from the 2D fit method is generally larger than the Fourier coefficient magnitude from the large-$ \Delta\eta $ correlations.

Figure 12.  (color online) $V_n$ parameters obtained from the 2D fit (circles, connected by lines), together with the Fourier coefficients calculated from the SS correlations at $1.5<|\Delta \eta |<2$ (squares). Results are shown for diffractive (open black points) and non-diffractive (filled red points) events as functions of multiplicity.

The $ V_2 $ coefficient shown in Fig. 12(b) is positive and exhibits a decreasing trend with multiplicity. The trends may be suggestive of nonflow, which would be inversely proportional to multiplicity (e.g., nonflow contribution from the near-side peak). It is possible that the 2D fits are not adequate to separate the correlations into flow and nonflow in these $ p+p $ interactions, such that the fitted harmonics still reflect some nonflow characteristics. An interesting observation is that the 2D fit $ V_2 $ values in the two lowest multiplicity bins of non-diffractive events are higher than the Fourier coefficients calculated at large $ | \Delta\eta|>1.5 $, whereas they are close to each other for the three higher multiplicity bins. This may reflect the difficulty in fitting the nonflow components in these low multiplicity events. The reason for this difficulty might be the near-side peak not having a large presence, i.e., it is barely noticeable for the two low multiplicity bins. Additionally, the fitted $ V_1 $ and $ V_2 $ parameters are anti-correlated, and this may contribute to the lowest multiplicity bin of the non-diffractive events. It may be at play for the diffractive events, where the three low multiplicity bins reveal larger fit $ V_2 $ and more negative fit $ V_1 $ than the Fourier counterparts.

Finally, the $ V_3 $ coefficient shown in Fig. 12(c) is close to zero with large uncertainties. The overall magnitude may be negative, suggestive of the influence of the away-side peak. Note that no qualitative differences seem to appear in the extracted $ V_n $ harmonics between diffractive and non-diffractive events despite distinctive physics processes.

One method to mitigate nonflow contributions is to subtract correlations in low-multiplicity events from those in high-multiplicity events [8]. This is possible under the assumption that the nonflow correlation shape does not change with multiplicity and the abundance of the nonflow correlation sources is proportional to multiplicity [8]. The latter also assumes, implicitly, that the selection of events according to event multiplicity has no biases to nonflow correlations. Under these assumptions, the ''flow'' harmonics in high-multiplicity events would be

where the superscripts "high" and "low" denote those quantities in high- and low-multiplicity events, respectively. Figure 13 shows the $ V_2^{\rm{sub}} $ as functions of multiplicity for non-diffractive events. The results are shown for two cases: one uses SS correlations without any $ \Delta\eta $ cut, and the other uses those with $ | \Delta\eta|>1.5 $. The results differ, and the interpretation can be that they may either (i) reflect collective flow–which can be possible in $ p+p $ interactions because of multi-parton interactions or color reconnection–and a strong $ \Delta\eta $decorrelation effect occurs or (ii) not reflect collective flow–which is naively expected to be $ \Delta\eta $ independent owing to its origin in early dynamics–but nonflow, as nonflow is not independent of multiplicity or strong multiplicity biases exist in event selection. For comparison, the 2D fit $ V_2 $ values in the three high multiplicity bins are also plotted in Fig. 13; those for the lower multiplicity bins are significantly higher as shown in Fig. 12. No systematics are apparent from the comparison.

Figure 13.  (color online) Low-multiplicity subtracted $V_2$ values as functions of multiplicity. Inverse multiplicity scaled subtraction (blue points) is shown for no $\Delta \eta $ cut and $|\Delta \eta |>1.5$. Dipole-scaled subtraction (pink points) is shown for $|\Delta \eta |>1.5$. Also shown are $V_2$ values from 2D fit (red squares). The results are obtained from SS correlation functions of non-diffractive events.

The $ V_2^{\rm{sub}} $ results shown in Fig. 13 from the various methods suggest the single-particle azimuthal anisotropy $ v_2 = \sqrt{V_2} $ of the order of 5%–7%. This is comparable to experimental findings at the LHC [22−24], despite the lower RHIC energy used in our simulation. This is particularly true for the $ V_2^{\rm{sub}} $ results extracted from correlations in the large-$ \Delta\eta $ region of $ | \Delta\eta|>1.5 $, the approximate region used in LHC experiments. It would be interesting to further investigate the possible anisotropic flows in Pythia simulations of $ p+p $ interactions at LHC energies [40].

Another recipe for nonflow subtraction is to assume that the negative $ V_1 $ dipole is all nonflow and the nonflow correlation shape does not change with multiplicity, i.e., the nonflow $ V_n $ ($ n>1 $) is proportional to $ V_1 $. This would lead to

The $ V_2^{\rm{sub}} $ values obtained from this dipole-scaled subtraction are indicated by pink points in Fig. 13 for SS correlations at $ | \Delta\eta|>1.5 $. They are all consistent with zero, implying that the $ V_1 $ and $ V_2 $ harmonics of the large-$ \Delta\eta $ SS correlations scale with each other; we confirm this by comparing Fig. 12(a) and Fig. 12(b).

Azimuthal anisotropies have been observed in high-energy small-system (proton-proton and proton/deuteron/ $ ^3{\rm He}$-nucleus) collisions, similar to those measured in relativistic heavy ion collisions. This raises an intriguing possibility of collective flow behavior and the formation of a quark-gluon plasma droplet in these small-system collisions and has aroused intense research interest. Because of severe contamination of nonflow (genuine few-body) correlations in these small systems, analysis and interpretation of the two-particle azimuthal cumulant anisotropies are challenging. Several nonflow mitigation methods are available with various assumptions and degrees of validity. At least three interpretations are possible: (i) the measured anisotropies are part of nonflow correlations that differ from low- to high-multiplicity events; (ii) the measured anisotropies are a result of multi-parton interactions and color connections, which are different from hydrodynamic-like interactions in the final state of heavy ion collisions; and (iii) an appreciable contribution comes from hydrodynamic-like final-state interactions and thus collective flow in these small-system collisions.

To elucidate further, we studied two-particle angular correlations in pseudorapidity and azimuthal differences via simulated $ p+p $ interactions using the Pythia 8.303 event generator. The model includes multi-parton interactions and color reconnection but no final-state interactions or hydrodynamic expansion. We present our study pedagogically and report azimuthal anisotropies extracted using several methods. These methods have different assumptions and sensitivities to nonflow correlations; thus, our results represent a fair range in analysis method dependency and degree of validity. We observe finite azimuthal anisotropies from several but not all of the employed methods. We observe no qualitative difference between diffractive and non-diffractive events, despite their distinctive physics processes. The results are qualitatively, and in some cases semi-quantitatively, similar to those observed in experimental data. Therefore, our study provides an important benchmark that can aid in improving data analysis and interpreting experimental findings.

Our study can be built up on and further enhanced in several ways. One is to turn off multi-parton interactions and color reconnection to examine the effects on the final azimuthal anisotropies. Another is to simulate $ p+p $ interactions at the higher LHC energies and expand to proton-nucleus and nucleus-nucleus collisions using the same Pythia model and investigate the systematics of the azimuthal anisotropies on system size and collision energy. A thorough comparison of experimental data to Pythia simulations with variations in physics processes, at the level of both the raw correlation functions and the extracted azimuthal anisotropies, should be performed. Some of such studies have already been reported, but a systematic study would be invaluable.

MST acknowledges support from the Purdue-Colombia Exchange and Partnership Program.