Belle Preprint 2022-11 KEK Preprint 2022-10

Y.-C. Chen ,

# Two-particle angular correlations in $e^{+}e^{-}$ collisions to hadronic
final states in two reference coordinates at Belle

0000-0002-9038-5324 Y.-J. Lee  0000-0003-2593-7767 P. Chang
0000-0003-4064-388X I. Adachi  0000-0003-2287-0173 H. Aihara ,
0000-0002-1907-5964 S. Al Said  0000-0002-4895-3869 D. M. Asner ,
0000-0002-1586-5790 H. Atmacan  0000-0003-2435-501X V. Aulchenko
0000-0002-5394-4406 T. Aushev  0000-0002-6347-7055 R. Ayad
0000-0003-3466-9290 V. Babu ,  0000-0003-0419-6912 P. Behera ,
0000-0002-1527-2266 K. Belous  0000-0003-0014-2589 J. Bennett
0000-0002-5440-2668 M. Bessner ,  0000-0003-1776-0439 V. Bhardwaj ,
0000-0001-8857-8621 B. Bhuyan  0000-0001-6254-3594 T. Bilka
0000-0003-1449-6986 D. Bodrov ,  0000-0001-5279-4787 G. Bonvicini ,
0000-0003-4861-7918 J. Borah  0000-0003-2990-1913 A. Bozek
0000-0002-5915-1319 M. Bračko  0000-0002-2495-0524 P. Branchini
0000-0002-2270-9673 A. Budano ,  0000-0002-0856-1131 M. Campajola
0000-0003-2518-7134 D. Červenkov ,  0000-0002-1865-741X M.-C. Chang
0000-0002-8650-6058 V. Chekelian  0000-0001-8860-8288 B. G. Cheon ,
0000-0002-8803-4429 K. Chilikin  0000-0001-7620-2053 H. E. Cho
0000-0002-7008-3759 K. Cho  0000-0003-1705-7399 S.-J. Cho  0000-0002-1673-5664
S.-K. Choi  0000-0003-2747-8277 Y. Choi  0000-0003-3499-7948 S. Choudhury
0000-0001-9841-0216 D. Cinabro ,  0000-0001-7347-6585 S. Das ,
0000-0001-6857-966X G. De Pietro ,  0000-0001-8442-107X R. Dhamija
0000-0001-7052-3163 F. Di Capua  0000-0001-9076-5936 J. Dingfelder
0000-0001-5767-2121 Z. Doležal  0000-0002-5662-3675 T. V. Dong
0000-0003-3043-1939 D. Epifanov ,  0000-0001-8656-2693 T. Ferber ,
0000-0002-6849-0427 D. Ferlewicz  0000-0002-4374-1234 B. G. Fulsom ,
0000-0002-5862-9739 R. Garg  0000-0002-7406-4707 V. Gaur  0000-0002-8880-6134
N. Gabyshev ,  0000-0002-8593-6857 A. Garmash  0000-0003-2599-1405 A. Giri
0000-0002-8895-0128 P. Goldenzweig  0000-0001-8785-847X E. Graziani ,
0000-0001-8602-5652 T. Gu ,  0000-0002-1470-6536 C. Hadjivasiliou ,
0000-0002-2234-0001 K. Hayasaka ,  0000-0002-6347-433X H. Hayashii
0000-0002-5138-5903 M. T. Hedges ,  0000-0001-6504-1872 D. Herrmann
0000-0001-9772-9989 W.-S. Hou ,  0000-0002-4260-5118 C.-L. Hsu
0000-0002-1641-430X K. Inami ,  0000-0003-2765-7072 N. Ipsita ,
0000-0002-2927-3366 A. Ishikawa  0000-0002-3561-5633 R. Itoh
0000-0003-1590-0266 M. Iwasaki  0000-0002-9402-7559 W. W. Jacobs ,
0000-0002-9996-6336 E.-J. Jang  0000-0002-1935-9887 S. Jia ,
0000-0001-8176-8545 Y. Jin  0000-0002-7323-0830 K. K. Joo  0000-0002-5515-0087
K. H. Kang ,  0000-0002-6816-0751 G. Karyan ,  0000-0001-5365-3716 T. Kawasaki
,  0000-0002-4089-5238 C. Kiesling ,  0000-0002-2209-535X C. H. Kim
0000-0002-5743-7698 D. Y. Kim  0000-0001-8125-9070 K.-H. Kim ,
0000-0002-4659-1112 Y.-K. Kim ,  0000-0002-9695-8103 P. Kodyš ,
0000-0002-8644-2349 A. Korobov ,  0000-0001-5959-8172 S. Korpar
0000-0003-0971-0968 E. Kovalenko  0000-0001-8084-1931 P. Križan
0000-0002-4967-7675 P. Krokovny ,  0000-0002-1236-4667 M. Kumar ,
0000-0002-6627-9708 R. Kumar  0000-0002-6277-2626 K. Kumara
0000-0003-1572-5365 Y.-J. Kwon ,  0000-0001-9448-5691 T. Lam ,
0000-0001-9128-6806 J. S. Lange  0000-0003-0234-0474 S. C. Lee ,
0000-0002-9835-1006 J. Li ,  0000-0001-5520-5394 L. K. Li  0000-0002-7366-1307
Y. Li ,  0000-0002-4413-6247 Y. B. Li  0000-0002-9909-2851 L. Li Gioi ,
0000-0003-2024-5649 J. Libby  0000-0002-1219-3247 C.-W. Lin
0000-0003-2187-8237 K. Lieret  0000-0003-2792-7511 D. Liventsev ,
0000-0003-3416-0056 M. Masuda ,  0000-0002-7109-5583 T. Matsuda
0000-0003-4673-570X S. K. Maurya ,  0000-0002-7764-5777 F. Meier
0000-0002-6088-0412 M. Merola  0000-0002-7082-8108 K. Miyabayashi ,
0000-0003-4352-734X R. Mizuk ,  0000-0002-2209-6969 G. B. Mohanty
0000-0001-6850-7666 M. Mrvar ,  0000-0001-6388-3005 R. Mussa ,
0000-0002-0294-9071 M. Nakao ,  0000-0001-8424-7075 Z. Natkaniec ,
0000-0003-0486-9291 A. Natochii  0000-0002-1076-814X M. Nayak ,
0000-0002-2572-4692 N. K. Nisar ,  0000-0001-9562-1253 S. Nishida
0000-0001-6373-2346 S. Ogawa ,  0000-0002-7310-5079 H. Ono ,
0000-0003-4486-0064 Y. Onuki ,  0000-0002-1646-6847 P. Oskin
0000-0002-7524-0936 G. Pakhlova ,  0000-0001-7518-3022 S. Pardi ,
0000-0001-7994-0537 H. Park  0000-0001-6087-2052 S.-H. Park
0000-0001-6019-6218 S. Patra ,  0000-0002-4114-1091 S. Paul
0000-0002-8813-0437 T. K. Pedlar  0000-0001-9839-7373 R. Pestotnik
0000-0003-1804-9470 L. E. Piilonen  0000-0001-6836-0748 T. Podobnik
0000-0002-6131-819X E. Prencipe ,  0000-0002-9465-2493 M. T. Prim
0000-0002-1407-7450 N. Rout ,  0000-0002-4310-3638 G. Russo ,
0000-0001-5823-4393 S. Sandilya ,  0000-0002-4199-4369 A. Sangal
0000-0001-5853-349X L. Santelj  0000-0003-3904-2956 V. Savinov ,
0000-0002-9184-2830 G. Schnell ,  0000-0002-7336-3246 C. Schwanda ,
0000-0003-4844-5028 R. Seidl ,  0000-0002-6552-6973 Y. Seino
0000-0002-8378-4255 K. Senyo  0000-0002-1615-9118 M. E. Sevior ,
0000-0002-4824-101X M. Shapkin  0000-0002-4098-9592 C. Sharma ,
0000-0002-1312-0429 C. P. Shen ,  0000-0002-9012-4618 J.-G. Shiu ,
0000-0002-8478-5639 E. Solovieva ,  0000-0002-5735-4059 M. Starič ,
0000-0001-8751-5944 Z. S. Stottler  0000-0002-1898-5333 J. F. Strube
0000-0001-7470-9301 M. Sumihama ,  0000-0002-8954-0585 T. Sumiyoshi ,
0000-0002-0486-3896 M. Takizawa ,  0000-0001-8225-3973 U. Tamponi ,
0000-0001-6651-0706 K. Tanida ,  0000-0002-8255-3746 F. Tenchini ,
0000-0003-3469-9377 K. Trabelsi ,  0000-0001-6567-3036 M. Uchida ,
0000-0003-4904-6168 T. Uglov ,  0000-0002-4944-1830 Y. Unno
0000-0003-3355-765X K. Uno  0000-0002-2209-8198 S. Uno  0000-0002-3401-0480 P.
Urquijo  0000-0002-0887-7953 R. van Tonder ,  0000-0002-7448-4816 G. Varner ,
0000-0002-0302-8151 A. Vinokurova ,  0000-0003-4220-8056 A. Vossen ,
0000-0003-0983-4936 E. Waheed  0000-0001-7774-0363 E. Wang ,
0000-0001-6391-5118 M.-Z. Wang  0000-0002-0979-8341 X. L. Wang
0000-0001-5805-1255 M. Watanabe  0000-0001-6917-6694 S. Watanuki ,
0000-0002-5241-6628 E. Won ,  0000-0002-4245-7442 W. Yan ,
0000-0003-0713-0871 S. B. Yang ,  0000-0002-9543-7971 J. H. Yin ,
0000-0002-1479-9349 C. Z. Yuan ,  0000-0002-1652-6686 Y. Yusa
0000-0002-4001-9748 Y. Zhai  0000-0001-7207-5122 Z. P. Zhang ,
0000-0001-6140-2044 V. Zhilich  0000-0002-0907-5565 V. Zhukova
0000-0002-8253-641X

###### Abstract

We present the analysis of two-particle angular correlations using coordinate
systems defined with the conventional beam axis and the event thrust axis. We
propose the latter as a good representation for the correlation structure
interpretation in the $e^{+}e^{-}$ collision system. The $e^{+}e^{-}$
collisions to hadronic final states at center-of-mass energies of
$\sqrt{s}=10.52$ GeV and 10.58 GeV are recorded by the Belle detector at KEKB.
In this paper, results on the first dataset are supplementary to the previous
Belle publication Belle:2022fvl . At the same time, the latter is the first
two-particle correlation measurement at collision energy on the $\mathchar
28935\relax(4S)$ resonance and is sensitive to its decay products.
Measurements are reported as a function of the charged-particle multiplicity.
Finally, a qualitative understanding of the correlation structure is discussed
using a combination of Monte Carlo simulations and experimental data.

##  1 Introduction

In nucleon-nucleon and heavy-ion collision experiments, two-particle angular
correlations STAR:2005ryu ; STAR:2009ngv ; PHOBOS:2009sau ; Chatrchyan:2012wg
; Aamodt:2011by ; Adam:2019woz  are extracted for the study of the Quark-Gluon
Plasma (QGP) Busza:2018rrf , color reconnection OrtizVelasquez:2013ofg , color
rope interaction Bierlich:2020naj  and the search for initial-state
correlation effects such as the Color Glass Condensate Dumitru:2010iy . In
these measurements, a long-range angular correlation, the ridge-like structure
STAR:2009ngv ; PHOBOS:2009sau , has been observed in various collision systems
and at various collision energies. Since the beginning of the LHC operation, a
ridge structure in the near-side ($\Delta\phi\approx 0$ between particles) has
been observed also in high-multiplicity proton-proton collisions by the CMS
collaboration Khachatryan:2010gv  and confirmed by experiments at LHC and RHIC
using smaller collision systems such as proton-ion CMS:2012qk ; Aaij:2015qcq
collisions. The studies have later been extended to the away-side
($\Delta\phi\approx\pi$) in proton-proton Aad:2015gqa , proton-ion
Abelev:2012ola ; ATLAS:2012cix  and deuteron-ion PHENIX:2013ktj  collisions.
In heavy-ion collisions, the ridge structure is associated with the
hydrodynamic expansion of a strongly interacting and expanding system
Ollitrault:1992bk  produced from collisions at different impact parameters and
the fluctuating initial state Alver:2010gr . However, the physical origin of
the ridge structure in small systems is still under debate Dusling:2013qoz ;
Bozek:2011if ; He:2015hfa ; Nagle:2018nvi . Recently, there has been growing
interest in measuring two-particle correlations in even smaller collision
systems of $ep$ ZEUS:2019jya  and $e^{+}e^{-}$ Nagle:2017sjv ; Badea:2019vey ;
Belle:2022fvl , as they serve as counterparts complementary to the results in
large collision systems and can be used to find the minimal conditions for
collective behavior Nagle:2017sjv . The use of an electron beam removes
complications such as multiple parton interactions and initial-state
correlations. Those measurements have provided additional insights on the
ridge signal Altinoluk:2020wpf ; Bierlich:2020naj ; Castorina:2020iia ;
Agostini:2021xca . In this work, two-particle correlations calculated using a
coordinate system with the beam pipe direction assigned as the $z$ axis are
reported in full detail. However, as pointed out in the study of ALEPH
$e^{+}e^{-}$ collisions via hadronic $Z$ decay data Badea:2019vey , this
conventional coordinate system might not be the optimal configuration in the
$e^{+}e^{-}$ collision system. In the $e^{+}e^{-}$ system, back-to-back dijets
fragmenting from the $q\bar{q}$ pair production dominate the event topology,
and the directions of such fragmentation are not fixed but emitted with a
$1+\cos^{2}\theta$ angular distribution relative to the beam pipe. Thus, pair-
wise correlations of the global dijet topology contribute differently to the
two-particle angular correlation functions from one event to another,
depending on the orientation of the dijet in each event. Such smearing effects
can degrade the detection sensitivity of the anomalous particle collectivity
originating from the possibly thermalized medium, which is, if any,
intuitively assumed to take place between the two out-going quarks. In the
$e^{+}e^{-}$ collision system, without the complexity of parton distribution
functions and large underlying events, one can mitigate the smearing effect to
allow for an understanding of finer correlation structures by defining the
angular correlations relative to the event thrust axis, which is a well-
defined surrogate of the hard process quark-pair direction using final-state
particles. On the other hand, it is hard to conclude the overall direction of
the parton-level interaction with a single axis in hadron-hadron collision. In
this paper, we present a detailed description of the analysis methodology and
results for the thrust-axis two-particle correlation function measurement.
Measurement with the thrust axis analysis provides a different view from the
conventional beam axis analysis. For the $e^{+}e^{-}$ annihilation process,
the thrust axis analysis is less descriptive of the leading dijet correlation
while more sensitive to higher-order soft emissions. We report on the
measurement of two-particle angular correlation functions in high-multiplicity
$e^{+}e^{-}$ annihilation events at $\sqrt{s}=10.52$ GeV and 10.58 GeV data
taken by the Belle experiment. The center-of-mass energy of the second dataset
is set on the $\mathchar 28935\relax(4S)$ resonance, with about one-fourth the
events decaying from $\mathchar 28935\relax(4S)$ states into $B\bar{B}$
mesons. The two-particle correlation of the latter dataset is studied for the
first time. This measurement is sensitive to anisotropic correlations in the
azimuthal angle due to the resonance decays in the dataset, which we further
discuss aided by Monte Carlo simulations. The other dataset is 60 MeV below
the resonance, and mostly composed of quark pairs fragmentating into hadrons.
Those are tests of the existence of ridge-like signals without the
complications from initial states. This study could also provide significant
new input to the phenomenological fragmentation models at the low-energy
regime of the Belle data. Taking advantage of the high statistics of the data,
an extreme charged-particle multiplicity reach of around 14 particles may be
achieved, corresponding to 0.04% of all hadronic $e^{+}e^{-}$ annihilation
events and being about twice the average multiplicity. In addition, we use
simulations to better understand the correlation structure and a comparison of
the expectations from different event generators, offering insights to the
fragmentation models and a deeper look into the explanation of two-particle
correlation. The results are partially reported in ref. Belle:2022fvl ; this
paper puts more emphasis on introducing the thrust-axis two-particle
correlations analysis methodology and understanding the measured results by
means of simulations. The presentation is organized as follows: The data
sample and corresponding event thrust distributions are reported in section 2.
Descriptions of the thrust-axis two-particle correlation function and the
corrections applied are given in sections 3 and 4. Two-particle correlations
and the projected azimuthal differential associated yields are shown in
section 5, as well as an estimate of the ridge-signal upper limits. With Monte
Carlo (MC) event generators, we compare and study correlation functions in
order to further understand the $e^{+}e^{-}$ collision at the low-energy
regime in section 6.

##  2 Experimental setup, data sample and the event thrust

We perform two-particle correlation measurements with Belle $e^{+}e^{-}$
collision data. Two datasets, one at $\sqrt{s}=10.52$ GeV and the other at
10.58 GeV are analyzed, corresponding to about 89.5 $\mbox{\,fb}^{-1}$ and
332.2 $\mbox{\,fb}^{-1}$ of integrated luminosity, respectively. The
$\sqrt{s}=10.58$ GeV collision data is targeted at the $\mathchar
28935\relax{(4S)}$ resonance and referred to as the “on-resonance data”, while
the other is called “off-resonance data”. Partial results for the off-
resonance dataset are released Belle:2022fvl  with the intention of providing
a complementary understanding of two-particle correlations measured in hadron
collisions with results from high-energy $e^{+}e^{-}$ collision. The full
measurement and the analysis details are documented in this paper. The Belle
experiment is configured with asymmetric 8 GeV electron and 3.5 GeV positron
beams, customized to study the decay-time-dependent $CP$ violation in $B$
decays. Data are collected with the Belle detector Abashian:2000cg ;
Belle:2012iwr , which is a large-solid-angle magnetic spectrometer consisting
of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an
array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement
of time-of-flight scintillation counters (TOF), and an electromagnetic
calorimeter (ECL) comprising CsI(Tl) crystals, all located inside a
superconducting solenoid coil that provides a 1.5 T magnetic field. An iron
flux-return surrounding the coil is instrumented to detect $K_{L}^{0}$ mesons
and muons (KLM). The evtgen Ryd:2005zz \- and pythia6 2001CoPhC.135..238S \-
based Belle MC sample is used as the simulation of $e^{+}e^{-}$ annihilation
events, including hadronic $q\bar{q}$ ($q=u,d,s$ and $c$) fragmentation,
$\mathchar 28935\relax{(4S)}$ decays, radiative Bhabha events, low
multiplicity $e^{+}e^{-}\rightarrow l^{+}l^{-}$ ($l=e$ or $\mu$) and two-
photon processes. This MC dataset is used to study reconstruction
inefficiencies and thereby deriving efficiency correction factors for the
data. In the following, results with the Belle MC sample are also provided
accompanying data results. Primary charged particles (table 1) in hadronic
events Belle:2001jqo  are used for the two-particle correlation study. We
require those particles to be within the detector’s angular coverage, and
apply a loose requirement on the track’s distance with respect to the
interaction point to ensure the selected tracks are from collisions. We focus
on high-momentum tracks with the transverse momentum $p_{\rm T}>0.2$ GeV$/c$.
Reconstruction effects such as duplicate low-$p_{\rm T}$ tracks and electron-
positron pairs from photon conversions are also studied, with specific removal
treatments listed in table 1. After all selections, the subsequent
reconstructed-track multiplicity (${\rm N}_{\rm trk}^{\rm rec}$) is used for
studying the correlation function’s multiplicity dependence. Considering the
computational cost, events with ${\rm N}_{\rm trk}^{\rm rec}<12$ are analyzed
with 11.5 $\mbox{\,fb}^{-1}$ and 11.7 $\mbox{\,fb}^{-1}$ partial samples for
the off-resonance and on-resonance dataset, respectively, whereas the full
datasets are used for analyzing the rarer high-multiplicity events. Figure 1
shows the reconstructed-track multiplicity for the off-resonance and on-
resonance datasets. The multiplicity classes used in this study, their
corresponding fraction of data, and the mapping of average reconstructed
multiplicities $\langle{\rm N}_{\rm trk}^{\rm rec}\rangle$ to average
multiplicities after efficiency correction $\langle{\rm N}_{\rm trk}^{\rm
corr}\rangle$ are listed in table 2. The efficiency correction procedure is
later introduced in section 4.1.

Table 1: Summary table for particle selections. Track selections have been
elucidated in ref. Belle:2022fvl . Neutral particles are selected mainly for
the event thrust calculation. Charged particles
---
Primary tracks |  • [reconstructed] decay promptly or from the long-lifetime
| particle whose decay vertex from the interaction point is
| less than 1 cm in the transverse plane
|  Long-lifetime particle candidates:
|  • $K_{S}^{0}$ ($0.480$-$0.516$ GeV$/c^{2}$)
|  • $\mathchar 28931\relax^{0}/\bar{\mathchar 28931\relax}^{0}$
($1.111$-$1.121$ GeV$/c^{2}$)
|  • [MC-truth] decay promptly or from particles with proper
| lifetime $\tau\leq 1$ cm/$c$
Acceptance | particles’ polar angle (defined with respect to the opposite
| direction of the $e^{+}$ beam): $17^{\circ}\leq\theta\leq 150^{\circ}$
High quality tracks |  $p_{\rm T}\geq 0.2$ GeV$/c$
Impact parameter | (tracks’ radial and $z$\- distance from the interaction point)
|  $|\Delta r|<2$ cm, $|\Delta z|<5$ cm
Duplicate track removal | veto the softer track of a low-momentum pair ($p_{\rm T}<0.4$ GeV$/c$)
| travelling with a small opening angle $\delta$:
|  1\. same-sign charges with $\cos\delta>0.95$
|  2\. opposite-sign charges with $\cos\delta<-0.95$
Photon conversion veto | veto track pairs which can form common vertices ($V^{0}$ objects)
| with
|  1\. $z$ distance between two tracks $<10$ cm
|  2\. reconstructed $V^{0}$’s mass $<0.25$ GeV$/c^{2}$
|  3\. decay-vertex radius $>1.5$ cm
Neutral particles (for thrust calculation)
Cluster selection | No association from tracks’ extrapolations
Acceptance | particles’ polar angle (defined with respect to the opposite
| direction of the $e^{+}$ beam): $17^{\circ}\leq\theta\leq 150^{\circ}$
Energy requirement | Forward endcap: $E\geq 0.10$ GeV
| Backward endcap: $E\geq 0.15$ GeV
| Barrel: $E\geq 0.05$ GeV

(a) off-resonance

(b) on-resonance

Figure 1: The reconstructed-track multiplicity (${\rm N}_{\rm trk}^{\rm rec}$) distributions for (a) off-resonance and (b) on-resonance samples. The corresponding integrated luminosity used for plotting is about 1.3 $\rm fb^{-1}$ for each. Results of the data (black dots) and the pythia6-based Belle MC (histograms) are overlaid. The colored bins in gold, blue, red, brown and gray are the events with ${\rm N}_{\rm trk}^{\rm rec}$ in $[0,6)$, $[6,10)$, $[10,12)$, $[12,14)$ and $[14,\infty)$ intervals. Table 2: Average multiplicities and corrected multiplicities with high-quality tracks for different ${\rm N}_{\rm trk}^{\rm rec}$ intervals for the off-resonance and on-resonance datasets. These tracks pass requirements in table 1: $p_{T}\geq 0.2$ GeV/$c$ and $17^{\circ}\leq\theta\leq 150^{\circ}$, hence the high efficiency. ${\rm N}_{\rm trk}^{\rm rec}$ interval | Off-resonance | On-resonance
---|---|---
|  Fraction (%) |  $\left<{\rm N}_{\rm trk}^{\rm rec}\right>$ |  $\left<{\rm N}_{\rm trk}^{\rm corr}\right>$ |  Fraction (%) |  $\left<{\rm N}_{\rm trk}^{\rm rec}\right>$ |  $\left<{\rm N}_{\rm trk}^{\rm corr}\right>$
$[6,10)$ |  44.33 |  6.98 |  7.05 |  48.80 |  7.06 |  7.18
$[10,12)$ |  2.65 |  10.26 |  10.12 |  3.93 |  10.29 |  10.20
$[12,14)$ |  0.29 |  12.20 |  11.90 |  0.53 |  12.23 |  11.99
$[14,\infty)$ |  0.02 |  14.22 |  14.24 |  0.04 |  14.26 |  13.93

Since the electron and positron beam energies delivered by KEK-B are
different, events are analyzed after being boosted to the $e^{+}e^{-}$
collision center-of-mass frame. The event thrust axis $\hat{n}$
PhysRevLett.39.1587  is obtained through maximizing the sum of the projected
particle momenta on itself in the center-of-mass frame, formulated as

$T~{}~{}{\accentset{\rm
max}{=}}~{}~{}\frac{\sum\nolimits_{i}\left|\vec{p_{i}}\cdot\hat{n}\right|}{\sum\nolimits_{i}\left|\vec{p_{i}}\right|},$
(1)

where $T$ is the resulting thrust value and $\vec{p_{i}}$ is the momentum of
the $i$-th particle, with the sum extending over charged, neutral particles
and the missing momentum. Neutral particles are ECL clusters without an
association from tracks’ extrapolations. Requirements on the ECL acceptance
and a cluster energy greater than 50–150 MeV depending on detector regions are
also applied. Detailed selections are summarized in table 1. The invisible
objects, such as neutrinos, and non-reconstructed particles at the
reconstruction level are taken into accounts as the missing momentum in the
calculation of the event thrust, given by

$\vec{p}_{\rm MET}=-\sum\nolimits_{\rm neu,chg}\vec{p}.$ (2)

Figure 2 shows the event thrust distributions for the off-resonance and on-
resonance dataset in bins of ${\rm N}_{\rm trk}^{\rm rec}$, with the
comparisons to the MC. The differences between data and MC are addressed by
reweighting studies which are discussed in Section 4.2.

(a) off-resonance

(b) on-resonance

Figure 2:  The thrust distributions for (a) off-resonance and (b) on-resonance
samples using the same partial samples as figure 1. Results of the data
(points) and the pythia6-based Belle MC (hatched areas) are overlaid. The
multi-panel histograms show the inclusive distributions and events with ${\rm
N}_{\rm trk}^{\rm rec}$ in $[0,6)$, $[6,10)$, $[10,12)$, $[12,14)$ and
$[14,\infty)$ intervals in black, gold, blue, red, brown and gray,
respectively.

##  3 Two-particle correlation function in the thrust axis coordinate

To measure two-particle correlation functions, we first perform a boost back
into the $e^{+}e^{-}$ center-of-mass frame. Two-particle correlations are then
calculated by using methods previously employed in heavy ion collision and
hadron collision measurements, which have been introduced in ref.
Belle:2022fvl  and formulated as

$\displaystyle\frac{1}{{\rm N}_{\rm trk}^{\rm corr}}\frac{d^{2}{\rm N}^{\rm
corr}_{\rm pair}}{d\Delta\eta d\Delta\phi}$
$\displaystyle=B(0,0)\times\frac{S(\Delta\eta,\Delta\phi)}{B(\Delta\eta,\Delta\phi)},$
(3)

where ${\rm N}_{\rm trk}^{\rm corr}$ denotes the number of tracks and ${\rm
N}^{\rm corr}_{\rm pair}$ the pairing yield, after efficiency correction in
the event. The correlation is analyzed in terms of the track-pair’s
pseudorapidity difference: $\Delta\eta=\pm(\eta_{i}-\eta_{j})$ and azimuthal
angle difference: $\Delta\phi=\pm(\phi_{i}-\phi_{j})$. We count four entries
of angular differences per track pair by assuming that the correlation
function is symmetric about the origin $(\Delta\eta,\Delta\phi)=(0,0)$. The
right-hand side of the equation indicates that the correlation is obtained by
factoring out the baseline background correlation function
$B(\Delta\eta,\Delta\phi)$ from the signal correlation function
$S(\Delta\eta,\Delta\phi)$, and scaled by the $B(0,0)$ factor. The signal and
background correlation, formulated as

$\displaystyle S(\Delta\eta,\Delta\phi)$ $\displaystyle=\frac{1}{{\rm N}_{\rm
trk}^{\rm corr}}\frac{d^{2}{\rm N}^{\rm same}}{d\Delta\eta d\Delta\phi},$ (4)
$\displaystyle B(\Delta\eta,\Delta\phi)$ $\displaystyle=\frac{1}{{\rm N}_{\rm
trk}^{\rm corr}}\frac{d^{2}{\rm N}^{\rm mix}}{d\Delta\eta d\Delta\phi},$

count the per-trigger-particle 111To describe the pairing between two
particles $i$-th and $j$-th, this paper uses the term “trigger particle” for
the $i$-th particle and the “associated particle” for the $j$-th one for
convenience. In other two-particle correlation analyses, there can be
different selection criteria for trigger and associated particles, but in this
paper, both are from the same collection if not further specified. pairing
yields within the same event and with the “mixed event” CMS:2011cqy ;
Chatrchyan:2012wg ; CMS:2012qk ; CMS:2013jlh ; CMS:2014und . A mixed event in
this analysis is formed by pairing events with similar ${\rm N}_{\rm trk}^{\rm
rec}$. In the analysis, trigger particles from a signal event are matched with
associated tracks from three other events separately, in order to increase the
statistical accuracy of the mixed event sample. The pairing yield in the
background correlation function (${\rm N}^{\rm mix}$) is obtained by counting
the number of associated tracks in the mixed event for each trigger particle
in the signal event. Since the associated track and the trigger track come
from different events, the pair is totally uncorrelated. The scale factor
$B(0,0)$ is the bin value at the origin of the background correlation
function, accounting for the normalization of the artificially constructed
background function. In this work, we label the reference coordinate system of
this classic two-particle correlation formulation, referring to the beam pipe
as its $z$-axis for the calculation of tracks’ $\eta$ and $\phi$ coordinates,
as the “beam axis coordinates.” However, the dominant energy flow of the
$e^{+}e^{-}$ collision events follow the directions of the out-going quark and
anti-quark produced in the electron-position annihilation, which differ event-
by-event. Consequently, in the fixed beam axis coordinate system, such
dominant dijet-like topology can be smeared out, affecting the observation of
any smaller correlation structure. More specifically, the leading and the
subleading jets, which make up the main back-to-back correlation feature in
$e^{+}e^{-}$ collision data, will contribute with arbitrary pseudorapidity
differences (i.e., not necessarily at large $\Delta\eta$) depending on the
outgoing dijet direction. In addition to the conventional beam axis analysis
construction, this paper documents the methodology of novel two-particle
correlations measured in the “thrust axis coordinates.” This analysis method
is devised initially by ref. Badea:2019vey  for $e^{+}e^{-}$ collision data,
taking the event thrust axis, which approximates the outgoing direction of
$e^{+}e^{-}$ collision events at the parton level, as the reference $z$-axis.
The difference in the thrust axis analysis is that tracks’ coordinates
($\eta,\phi$) are recalculated with respect to the new reference $z$-axis, the
thrust axis $\hat{n}$. The reference axis for $\phi=0$ (or new reference
$x$-axis) is assigned with $\hat{n}\times(\hat{n}\times\hat{z})$, where
$\hat{n}$ and $\hat{z}$ are the unit vectors of the event thrust and the beam
direction. This defines a new set of coordinate values in the thrust axis
reference frame ($\eta_{T},\phi_{T}$), which can be used to calculate the
angle separations and to form the signal and background correlation
distributions in eq. (3). For the background correlation function calculation,
we use the same event-mixing method as the typical beam axis analysis by
combining different events to form a mixed event, which is then used to
account for the random pairing baseline correlation in the thrust axis
analysis. However, since the thrust axis of each paired event is randomly
orientated, a thrust-mixing reweighting correction on particles’ thrust-
referenced angular distribution ($\eta_{T},\phi_{T}$) is adopted. Such a
correction is needed since the variation of the thrust axis direction causes
detector acceptance in the thrust coordinates to vary on an event-by-event
basis. The reweighting factor is obtained from a histogram division of the
physical spectrum by the mixed-event spectrum, written out as

$m(|\eta_{T}|,\phi_{T},{\rm N}_{\rm trk}^{\rm rec})=\left[\frac{d^{2}{\rm
N}^{\rm phys}}{d|\eta_{T}|d\phi_{T}}/\frac{d^{2}{\rm N}^{\rm
mix}}{d|\eta_{T}|d\phi_{T}}\right]_{{\rm N}_{\rm trk}^{\rm rec}},$ (5)

which can force the mixed-event single-particle angular distribution to mimic
the physical one. The reweighting factor is further applied as track weights
on the background two-particle correlation function
$B(\Delta\eta,\Delta\phi)$. This ensures that the event topology of mixed
events conform to that of physical events, and leaves the pairing of the
trigger particle with tracks from the mixed event uncorrelated at the same
time. Since the particle’s pseudorapidity evolves non-linearly towards
divergence when getting close to the reference axis, the two-particle
correlation function observable defined on the $\Delta\eta$-$\Delta\phi$ plane
monitors mainly final states traveling in the transverse direction to the
reference axis. In other words, particle pairs with very large $\eta$ values
are beyond the scope of interest for the $\Delta\eta$ window. Thus, when
measuring with the correlation defined with respect to the beam axis for heavy
ion collisions, one is looking at the mid-rapidity region where final states
expand from a possible QGP, while backgrounds from beam remnants in the
extreme forward and backward regions are out of the pseudorapidity coverage of
the observable defined. For the thrust axis coordinates, the mid-rapidity
region probes the physics of soft emission while the leading-order quark-
initiated jet correlations with large pseudorapidity differences are omitted.
A projection of correlations to the particle pairs’ azimuthal angle difference
$\Delta\phi$ is also studied, given by integrating and averaging two-
dimensional correlation functions over $\Delta\eta$ from $\Delta\eta_{\rm
min}$ to $\Delta\eta_{\rm max}$,

$Y(\Delta\phi)=\frac{1}{{\rm N}_{\rm trk}^{\rm corr}}\frac{d{\rm N}^{\rm
corr}_{\rm pair}}{d\Delta\phi}=\frac{1}{\Delta\eta_{\rm max}-\Delta\eta_{\rm
min}}\int\limits_{\Delta\eta_{\rm min}}^{\Delta\eta_{\rm max}}\frac{1}{{\rm
N}_{\rm trk}^{\rm corr}}\frac{d^{2}{\rm N}^{\rm corr}_{\rm pair}}{d\Delta\eta
d\Delta\phi}d\Delta\eta.$ (6)

We explore this azimuthal differential associated yield for the full range of
$\Delta\phi$ and in three $\Delta\eta$ regions: short range
($0\leq|\Delta\eta|<1$), middle range ($1\leq|\Delta\eta|<1.5$) and long range
($1.5\leq|\Delta\eta|<3.0$). The long-range azimuthal differential associated
yield is denoted as $Y_{l}(\Delta\phi)$ in the following content.

##  4 Corrections

###  4.1 Tracking efficiency correction

With the Belle MC sample as introduced in section 2, we study the
reconstruction effects of nonuniform detection efficiency and
misreconstruction bias. A reweighting factor is applied to the reconstructed
tracks by the inverse of the tracking efficiency, bringing the reconstructed-
track spectra closer to reflecting the geometrical acceptance of the detector.
The efficiency is given by

$\varepsilon(p_{\rm T},\theta,\phi,{\rm N}_{\rm trk}^{\rm
rec})=\left[\frac{d^{3}{\rm N}^{\rm reco}}{dp_{\rm T}d\theta
d\phi}/\frac{d^{3}{\rm N}^{\rm gen}}{{dp_{\rm T}d\theta d\phi}}\right]_{{\rm
N}_{\rm trk}^{\rm rec}},$ (7)

where ${\rm N}^{\rm gen}$ denotes the number of primary charged particles (as
defined in table 1) counted at the generator level, and ${\rm N}^{\rm reco}$
denotes those at the reconstruction level. The tracking efficiency correction
addresses reconstruction effects such as the inclusion of secondary particles
and detector effects.

###  4.2 MC reweighting

There are potential differences between data and simulation: for instance,
some processes are not taken into account, or the proportion of each process
is not estimated accurately, etc. The MC sample used to model the
reconstruction effects is thus reweighted to correct for the imperfection in
MC simulation. A “MC reweighting factor $r$” is applied as the event weight to
adjust the discrepancy. This factor is obtained by a histogram ratio,

$r({\rm N}_{\rm trk},T)=\langle\frac{d^{2}{\rm N}^{\rm data}}{d{\rm N}_{\rm
trk}dT}\rangle/\langle\frac{d^{2}{\rm N}^{\rm MC}}{d{\rm N}_{\rm
trk}dT}\rangle,$ (8)

where the numerator and the denominator are the data and MC event multiplicity
(${\rm N}_{\rm trk}$) and thrust ($T$) distributions normalized by number of
events, respectively. For events with the number of tracks less than $12$, the
ratios are close to 1 ($0.98$–$1.08$), while the ratios are larger in high-
multiplicity events, shooting up to $1.14$–$1.34$. The MC reweighting factor
$r$ is incorporated into the calculation of the efficiency correction factor
(see eq. (7)) to derive the new factor $\varepsilon^{\prime}$, given by

$\varepsilon^{\prime}(p_{\rm T},\theta,\phi,{\rm N}_{\rm trk}^{\rm
rec})=\left[\frac{d^{3}[r({\rm N}_{\rm trk},T){\rm N}^{\rm reco}]}{dp_{\rm
T}d\theta d\phi}/\frac{d^{3}[r({\rm N}_{\rm trk},T){\rm N}^{\rm
gen}]}{{dp_{\rm T}d\theta d\phi}}\right]_{{\rm N}_{\rm trk}^{\rm rec}}.$ (9)

###  4.3 Bin-size effect correction — $B(0,0)$ normalization extrapolation and
long-range correlations scaling

In the calculation of two-particle correlations (see eq. (3)), we account for
the normalization factor $B(0,0)$ by extrapolation. This procedure is used to
resolve any artificial bin-width effect222 The correlation function taken as
the number of pairing yields per $\Delta\eta(\Delta\phi)$ is an expression of
the density function ($d^{2}N/\text{(bin width)}^{2}$), hence there is no
trivial bin-width dependency. However, different $\Delta\eta(\Delta\phi)$
window sizes also affect the pairing yields average; this second-order bin
size effect is estimated in this work.  in the conventional method, where
$B(0,0)$ is given by the yield in the zeroth bin. The effects introduced by
using finite-bin histogramming to approximate the correlation function are
studied independently for $\Delta\eta$ and $\Delta\phi$. The impact of the
variation of the $\Delta\phi$ bin width is checked to have a small effect on
the $B(0,0)$ value compared to that caused by the variation of the
$\Delta\eta$ bin width. Therefore, we focus only on the correction of the
$\Delta\eta$ bin choice dependency on the raw $B(0,0)$ value, leaving the
uncertainty regarding the $\Delta\phi$ bin choice as one of the sources of
systematic uncertainty. The trend of the $B(0,0)$ value as a function of the
$\Delta\eta$ bin-width configuration is modeled with a second-order polynomial
asymptotically approaching a constant as the $\Delta\eta$ bin width goes to
zero. We use as the extrapolation of $B(0,0)$, the value at the zero
$\Delta\eta$ bin width. For the beam axis analysis, a residual bin-size effect
on the normalization of the long-range azimuthal differential yield
$Y_{l}(\Delta\phi)$ is found. Hence, we calibrate the long-range correlation
magnitude by the same approach as that to correct the $B(0,0)$ value. The
final beam-axis $Y_{l}(\Delta\phi)$ is multiplied by a scaling factor
$\frac{\int Y_{l}(\Delta\phi)d\Delta\phi]_{ext}}{\int
Y_{l}(\Delta\phi)d\Delta\phi]_{hist}}$, where the numerator is the total
associated yield in the long-range region using a second-order-polynomial
extrapolation, and the denominator is the raw histogram value.

###  4.4 Residual MC correction

After the above corrections, remaining possible reconstruction effects are
handled with the bin-by-bin correction method Choudalakis:2011rr . The
correction factor is derived from the ratio of MC correlation functions at the
reconstruction and generator level as ${\rm
C}(\Delta\phi)=\frac{Y(\Delta\phi)_{{\rm gen},i_{g}}}{Y(\Delta\phi)_{{\rm
reco},i_{r}}}$, where indices $i_{g}$ and $i_{r}$ are ${\rm N}_{\rm trk}^{\rm
rec}$ bins counted at the generator and reconstruction level, respectively.
The factors correct $2$-$5\%$ discrepancies between the reconstruction and
generator level. Final data correlation results are obtained from the
multiplication of the original correlation function with the bin-by-bin
correction factor.

##  5 Results

Partial results of the two-particle correlations measured using the off-
resonance dataset are published in ref. Belle:2022fvl . In figure 3, beam-axis
two-particle correlations for the off-resonance and on-resonance data are
shown as a function of ${\rm N}_{\rm trk}^{\rm rec}$. The correlation
structure in view of the beam axis analysis can be understood in the following
decomposition. The near-side peak correlation, residing near
$(\Delta\eta,\Delta\phi)=(0,0)$, has contributions mainly from track pairs
that originate in the same jet. The away-side correlation, elongated along
$\Delta\phi\approx\pi$, results from back-to-back momentum balancing. This
correlation in $e^{+}e^{-}$ annihilation has a feature that tilts upwards to
large $\Delta\eta$. Since the annihilation products of the $e^{+}e^{-}$
collision do not happen to align with the direction of the beamline, this
increases the probability of finding back-to-back emitted final states at
larger pseudorapidity difference. The typical ridge correlation is the
structure with $\Delta\phi\approx 0$ and persisting over a wide range of
$\Delta\eta$ Abelev:2012ola . This is not visually observed in the Belle data.

Figure 3: Beam-axis two-particle correlation functions for the off-resonance
(top) and on-resonance (bottom) data versus the reconstructed-track
multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ in $[0,6)$, $[6,10)$, $[10,12)$ and
$[14,\infty)$ intervals from left to right.

Zoom-in displays of beam-axis correlation functions in different $\Delta\eta$
ranges are shown by azimuthal differential associated yield distributions in
figures 4 and 5 with the off-resonance and on-resonance data, respectively.
Results for data (black points with error bars) and MC (blue curves with error
bars) are overlaid in the short range (left), middle range (middle) and long
range (right) region for comparison. Plot labels indicate the correlated
systematic uncertainties over $\Delta\phi$ bins; uncorrelated systematic
uncertainties depending on the $\Delta\phi$ value are plotted with gray boxes.
There are some discrepancies in the near-side ($\Delta\phi\approx 0$) and the
away-side ($\Delta\phi\approx\pi$) peak values between the on-resonance data
and MC results, while for the $e^{+}e^{-}\rightarrow q\bar{{q}}$ off-resonance
sample, the MC provides a better match with the data.

Figure 4: Beam-axis azimuthal differential associated yields as a function of
the reconstructed-track multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ for the
off-resonance sample in the short-range ($0\leq|\Delta\eta|<1$), middle-range
($1\leq|\Delta\eta|<1.5$) and long-range ($1.5\leq|\Delta\eta|<3.0$) regions
in columns from left to right, respectively. Figure 5: Beam-axis azimuthal
differential associated yields as a function of the reconstructed-track
multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ for the on-resonance sample in the
short-range ($0\leq|\Delta\eta|<1$), middle-range ($1\leq|\Delta\eta|<1.5$)
and long-range ($1.5\leq|\Delta\eta|<3.0$) regions in columns from left to
right, respectively.

In figure 6, thrust-axis two-particle correlations for the off-resonance and
on-resonance data are shown as a function of the reconstructed-track
multiplicity. The magnitude of correlations decreases significantly compared
with that of the beam axis analysis, due to the exclusion of the dominant
back-to-back correlation at the extreme pseudorapidity difference. The main
correlation structure is the bump at the away-side region; compared with the
$e^{+}e^{-}$ collisions at the $Z$-pole energy Badea:2019vey , a sizable near-
side peak correlation is lacking. We investigate this latter effect in section
6 by simulations, finding that the near-side peak correlation of $e^{+}e^{-}$
annihilation depends strongly on the collision energy. In general,
qualitatively similar results are observed for the off-resonance and on-
resonance data.

Figure 6: Thrust-axis two-particle correlation functions for the off-resonance
(top) and on-resonance (bottom) data versus the reconstructed-track
multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ in $[0,6)$, $[6,10)$, $[10,12)$ and
$[14,\infty)$ intervals from left to right.

Figures 7 and 8 show details of the data and MC comparison for thrust-axis
coordinate azimuthal differential associated yields with the off-resonance and
on-resonance data, respectively. In the small $\Delta\eta$ region, there are
mild-level magnitude differences between data and MC correlation. In the long
range region, a qualitatively better agreement in correlation function shapes
between the data and MC is observed for off-resonance results; however, there
is up to a 5% larger discrepancy seen in the correlation function magnitude
for the on-resonance sample. For the search of the ridge signal, the
$e^{+}e^{-}\rightarrow q\bar{{q}}$ off-resonance dataset reveals no evidence
of special enhancement, as reported by the measurement with ALEPH data
Badea:2019vey  and the recent Belle Letter Belle:2022fvl . For the on-
resonance data, it is worth noting that there is an enhanced long-range near-
side correlation. The enhancement is visible for events with multiplicity
greater than 10 in the rightmost columns of figure 8. However, it doesn’t
resemble the typical ridge structure observed widely in the large collision
systems, which is an enhanced correlation region over a broad range of
$\Delta\eta$ with small $\Delta\phi$. We also observe a similar enhancement
with MC, and provide further investigation of the enhancement using MC in
section 6.

Figure 7: Thrust-axis azimuthal differential associated yields as a function
of the reconstructed-track multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ for the
off-resonance sample in the short-range ($0\leq|\Delta\eta|<1$), middle-range
($1\leq|\Delta\eta|<1.5$) and long-range ($1.5\leq|\Delta\eta|<3.0$) regions
in columns from left to right, respectively. Figure 8: Thrust-axis azimuthal
differential associated yields as a function of the reconstructed-track
multiplicity ${\rm N}_{\rm trk}^{\rm rec}$ for the on-resonance sample in the
short-range ($0\leq|\Delta\eta|<1$), middle-range ($1\leq|\Delta\eta|<1.5$)
and long-range ($1.5\leq|\Delta\eta|<3.0$) regions in columns from left to
right, respectively.

The systematic uncertainties due to selection and correction operations are
calculated with respect to the long-range associated yield
$Y_{l}(\Delta\phi)$, and listed as percentages in tables 3 and 4 for the beam-
axis and thrust-axis analyses, respectively. In the tables, correlated
systematic uncertainties for all $\Delta\phi$ bins are reported. The
associated yield variation from each systematic source discussed below is
observed in general to be uniform along $\Delta\phi$, while a mild level of
non-uniform systematic uncertainty on the associated yield is observed in the
variation of the thrust mixing reweighting factor parametrization (see eq.
(5)). The effect is estimated to be of order $0.01\%$. The uncorrelated
uncertainties depending on $\Delta\phi$ values are shown in the azimuthal
differential associated yield plots (figures 4, 5, 7 and 8).

Table 3: Beam axis analysis systematic uncertainties for the off-resonance and on-resonance data as a function of the reconstructed-track multiplicity ${\rm N}_{\rm trk}^{\rm rec}$. All values are reported as percentages of the long-range differential associated yield. ${\rm N}_{\rm trk}^{\rm rec}$ |  Primary particle selection |  Tracking efficiency |  Event selection |  $\Delta\phi$ bin width |  $B(0,0)$ extrapolation |  Long-range scaling |  MC reweighting
---|---|---|---|---|---|---|---
off-resonance
$[6,10)$ |  0.42 |  0.35 |  0.38 |  0.01 |  0.01 |  0.01 |  0.01
$[10,12)$ |  0.30 |  0.35 |  $<$0.01 |  $<$0.01 |  0.03 |  0.04 |  0.03
$[12,14)$ |  0.39 |  0.35 |  $<$0.01 |  0.02 |  0.04 |  0.08 |  0.04
$[14,\infty)$ |  6.27 |  0.35 |  $<$0.01 |  0.02 |  0.14 |  0.30 |  0.20
on-resonance
$[6,10)$ |  0.25 |  0.35 |  0.25 |  0.02 |  0.01 |  0.01 |  0.02
$[10,12)$ |  0.02 |  0.35 |  $<$0.01 |  $<$0.01 |  0.02 |  0.03 |  $<$0.01
$[12,14)$ |  0.04 |  0.35 |  $<$0.01 |  0.03 |  0.01 |  0.01 |  0.04
$[14,\infty)$ |  1.02 |  0.35 |  $<$0.01 |  0.02 |  0.02 |  0.05 |  0.08
Table 4: Thrust axis analysis systematic uncertainties for the off-resonance and on-resonance data as a function of the reconstructed-track multiplicity ${\rm N}_{\rm trk}^{\rm rec}$. All values are reported as percentages of the long-range differential associated yield. ${\rm N}_{\rm trk}^{\rm rec}$ |  Primary particle selection |  Tracking efficiency |  Event selection |  $\Delta\phi$ bin width |  $B(0,0)$ extrapolation |  MC reweighting |  Thrust mixing reweighting
---|---|---|---|---|---|---|---
off-resonance
$[6,10)$ |  0.15 |  0.35 |  0.04 |  0.05 |  0.02 |  $<$0.01 |  0.16
$[10,12)$ |  0.01 |  0.35 |  0.26 |  0.06 |  0.04 |  0.02 |  0.03
$[12,14)$ |  0.71 |  0.35 |  1.32 |  0.03 |  0.06 |  0.02 |  0.09
$[14,\infty)$ |  3.81 |  0.35 |  0.15 |  0.23 |  0.21 |  0.12 |  0.25
on-resonance
$[6,10)$ |  0.03 |  0.35 |  0.07 |  0.04 |  0.01 |  0.03 |  $<$0.01
$[10,12)$ |  0.27 |  0.35 |  0.06 |  0.05 |  0.03 |  0.01 |  0.04
$[12,14)$ |  0.47 |  0.35 |  0.02 |  0.02 |  0.01 |  0.03 |  0.02
$[14,\infty)$ |  0.12 |  0.35 |  1.13 |  0.01 |  0.03 |  0.06 |  0.45

The primary particle selection systematic uncertainty is estimated by changing
the generator definition of the proper lifetime requirement $\tau<1$ cm/$c$ to
the vertex requirement $V_{r}<1$ cm. This variation of different truth-level
definitions enters in the correction factor calculation for the tracking
efficiency. A 0.35% uncertainty on the reconstruction efficiency for high-
momentum tracks with $p_{\rm T}>200$ ${\mathrm{\,Me\kern-1.00006ptV\\!/}c}$
BaBar:2014omp  propagates directly to the calculation of the per-trigger-
particle associated yield magnitude. The sample in this analysis passes
hadronic event selections Belle:2001jqo , which requires the energy sum
($E_{\rm sum}$) in the ECL to be greater than $0.18\sqrt{s}$. We examine the
systematic effect of hadronic event selection by tightening it to
$0.23\sqrt{s}$. The variation of the event selection has an impact in the low-
multiplicity region primarily, while the effect is not significant for high
multiplicity events. We correct for the histogramming bin-size effect of the
two-particle correlation function by extrapolating to the zero $\Delta\eta$
bin width discussed in section 4.3 and consider the $\Delta\phi$ bin choice as
a source of systematic uncertainties. The nominal $\Delta\phi$ bin width is
0.157; we recalculate with an alternative bin choice, corresponding to a bin
width of 0.174. Either are enough to describe the correlation function with
perfect granularity. The resulting difference is quoted as one of the
systematics. We only use this alternative configuration to evaluate the
$\Delta\phi$ bin-width systematics, and assume similar changes in the case for
decreasing bin-width, given that the size of this uncertainty is roughly of
order $0.01\%$ and therefore small in comparison to the dominating
uncertainties. For the $B(0,0)$ normalization factor, we calibrate the bin-
size effect by extrapolation with a second-order fitting curve. The fitting
error is propagated to the $Y_{l}(\Delta\phi)$ and is considered as one of the
systematic uncertainties. The uncertainty on the scaling factor of the long-
range correlation magnitude is also given by the propagation of the fitting
error. The systematic uncertainties of reweighting factors are studied by
changing their parametrization. For the MC reweighting factor (see eq. (8)),
we change the dependence from the event multiplicity and thrust to
multiplicity only; for the thrust mixing reweighting factor (see eq. (5)), we
change the dependence from ($|\eta|$, $\phi$) to $|\eta|$ only. The
corresponding difference on the correlation observable is quoted as the
systematic uncertainty. The systematic uncertainties on the residual MC
correction factors are assessed by changing the correction factors’
parametrization to two dimensions $(\Delta\eta,\Delta\phi)$; hence, the
alternative correction factors are determined using the ratio of two-
dimensional correlation functions. Uncertainties in the correction factors are
found to be no greater than $\mathcal{O}(0.01\%)$, which is a negligible
effect compared with other sources of systematics. Systematic uncertainties
from selections (including primary particles, tracks and the hadronic events),
totalling from approximately $0.4\%$–$6.3\%$, contribute as dominant
systematic uncertainties in this measurement. Larger uncertainties appear in
high-multiplicity bins, which are due to the limitation in statistics to
acquire the precise track efficiency correction factor, the MC reweighting
factor and the thrust-mixing reweighting factor. In the beam axis analysis,
the variation of the event selection criteria has little effect on the change
of associated yields in events with multiplicity greater than 10, while in the
thrust axis analysis, since an additional operation of the thrust-mixing
reweighting is applied, the uncertainties are greater. Other uncertainties are
comparably small, contributing from $0.01\%$ up to $0.5\%$ in the difference
of long-range associated yields. A common measure for the structure of the
associated yield function uses the “zero yield at minimum” (ZYAM) method
Ajitanand:2005jj . The minimum of the $Y(\Delta\phi)$ associated yield is
shifted to the zero yield after a flat correlation contribution, denoted as a
constant $C_{\rm ZYAM}$, is subtracted off. Since correlation functions are
symmetric about $\Delta\phi=0$, we consider three even-function fit templates

$f(\Delta\phi)=\left\\{\begin{array}[]{l}v_{0}+2\sum\limits_{n=1}^{3}v_{n}\cos(n\Delta\phi),\\\
a_{0}+a_{1}(\Delta\phi)^{2}+a_{2}\cos(2\Delta\phi),\\\
a_{0}+a_{1}(\Delta\phi)^{2}+a_{2}(\Delta\phi)^{4},\end{array}\right.$ (10)

for the determination of the constant $C_{\rm ZYAM}$ correlation and the $x$
coordinate of the minimum ($\Delta\phi_{\rm min}$). Thereafter, the ridge
yields are quantified by integrating ZYAM-subtracted azimuthal differential
associated yields over the long-range near-side
($0\leq\Delta\phi\leq\Delta\phi_{\rm min}$) region,

$Y_{\rm ridge}=\int\limits_{0}^{\Delta\phi_{\rm min}}[Y_{l}(\Delta\phi)-C_{\rm
ZYAM}]d\Delta\phi.$ (11)

For the measurement shows no obvious ridge-like structure, the ridge yield
upper limit is reported using the bootstrap method Efron:1979bst . The
measured azimuthal differential associated yields are varied according to
their statistical and systematic uncertainties and re-evaluated with the ZYAM
method to obtain an alternative ridge yield. The procedure is repeated two
million times to form the bootstrapped ridge yield datasets using the three
fit templates (see eq. (10)), where the most conservative 95% upper limits are
reported. We report 95% upper limits of the ridge yield or the confidence
level of ridge-signal exclusion at $10^{-7}$ as a function of the averaged
corrected multiplicity $\langle{\rm N}_{\rm trk}^{\rm corr}\rangle$. With the
large Belle data size, strong $5\sigma$ ridge-signal exclusions at $10^{-7}$
are set in the beam axis analysis in most multiplicity ranges for both the
off-resonance and on-resonance data, except the highest multiplicity bin
(${\rm N}_{\rm trk}^{\rm corr}\geq 14$) with the least statistics in the off-
resonance sample, a $97\%$ confidence level of $10^{-7}$ exclusion is set. The
upper limit results of the thrust axis analysis are summarized in figure 9.
The off-resonance upper limits are less than $O(10^{-3})$ level for the
majority of events. The smallness of on-resonance long-range near-side
enhancement, as pointed out in figure 8, is quantified at about $O(10^{-3})$
to $O(10^{-2})$. Upper limits are also provided considering the quantitative
significance level of the measured central values.

Figure 9:  Central values (square markers for on-resonance) or upper limits of
the ridge yield as a function of $\langle{\rm N}_{\rm trk}^{\rm corr}\rangle$
for the off-resonance (red) and on-resonance (black) datasets for the thrust
axis analyses. The associated uncertainties with the measured ridge-yield
central values are smaller than the symbols. The label “$>5\sigma$” indicates
a $5\sigma$ confidence level upper limit.

##  6 Studies with event generators

Complementary to ref. Belle:2022fvl , we provide the off-resonance and on-
resonance results in full ${\rm N}_{\rm trk}^{\rm rec}$ analysis ranges
comparing data with pythia6-based Belle MC, pythia 8.240 Sjostrand:2014zea ,
herwig 7.2.2 Bellm:2015jjp  and sherpa 2.2.5 Gleisberg:2008ta  event
generators for long-range azimuthal differential associated yields after
performing the ZYAM method. In this sub-section, we will call the
pythia6-based Belle MC by its version number pythia6, in comparison with other
event generators. Figure 10 displays the results. Generators (pythia8, herwig
and sherpa) compared here are configured based on their default settings. More
specifically, in pythia8 event generation, the Monash 2013 tune Skands:2014pea
($e^{+}e^{-}$ collision’s default tune) is used. The $e^{+}e^{-}\rightarrow
q\bar{q}$ annihilation process is simulated with the full interference between
photon and $Z$-boson propagators (option WeakSingleBoson:ffbar2gmZ), and the
on-resonance production is generated by decaying the $\mathchar
28935\relax(4S)$ particle evenly into charged and neutral $B$-meson pairs,
where the mass and decay width of $\mathchar 28935\relax(4S)$ are set to the
corresponding world-average values PDG . In the herwig event generation, we
require the order of electroweak coupling to be 2 (OrderInAlphaEW=2) and the
order of strong coupling to be 0 (OrderInAlphaS=0) in the hard process
simulation. For the $\mathchar 28935\relax(4S)$ resonance production, the
built-in matrix element for $e^{+}e^{-}$ colliding into a spin-$1$ vector
meson (MEee2VectorMeson) is used. Decays of $B$ mesons are further done by the
evtgen/pythia8 interface. The sherpa MC samples are generated with
$\alpha_{\rm S}(M_{Z})=0.1188$, and a two-loop correction of the running of
$\alpha_{\rm S}$ (ORDER_ALPHAS=1) is specified. Full mass effects of charm and
beauty quarks are considered. For the $\mathchar 28935\relax(4S)$ resonance
production, charged and neutral $B$-meson pairs are generated with the same
probabilities. In general, there is no significant indication of multiplicity-
dependent trends in the comparisons of data and MC long-range correlations,
while minor inconsistencies across different multiplicity ranges could come
from imperfections of the simulated multiplicity distributions. As pointed out
in ref. Belle:2022fvl  using high-multiplicity results, pythia6 shows good
agreements with data in the off-resonance results (figure 10(a)). herwig and
sherpa generators are also capable of describing the data correlations in the
near-side region, but have discrepancies in the away-side correlation
magnitudes. In the thrust-axis analysis, there are differences in data and MC
correlation distribution shapes. MCs have flatter correlations in the near-
side region, while their away-side peaks are steeper than that of data.
Besides this general behavior, a notable feature is that herwig overshoots
data as $\Delta\phi<1$ and $\Delta\phi\approx\pi$; sherpa shows similar but
smaller deviations. Similarly, disagreements in the away-side correlation
magnitude also appear in the on-resonance beam-axis results for pythia8,
herwig, and sherpa generators in figure 10(b). But unlike the off-resonance
case, there are over-predictions in the beam-axis away-side magnitude from
pythia6. There is mild near-side enhancement in the on-resonance thrust-axis
correlations. This feature is also found in generator results, with pythia6,
herwig, and sherpa overshooting the data points, while pythia8 undershoots
them.

(a) off-resonance

(b) on-resonance

Figure 10: Comparison of the ZYAM-subtracted $Y(\Delta\phi)$ in the range
$1.5\leq|\Delta\eta|<3.0$ for beam (top) and thrust (bottom) axis analyses as
a function of the charged-particle multiplicity. The colored bands show
simulation predictions from pythia 6 (green), pythia8 (yellow), herwig (blue)
and sherpa (violet). Reconstruction effects on data (data points) are
corrected, with error bars representing associated statistical uncertainties,
and the gray boxes the systematic uncertainties. For visual purposes, the
minimal statistical uncertainties of the MC correlation colored bands are set
to be $0.4\%$ of the plotting ranges, and residual panels have minimum
thresholds of $3\%$ of the plotting ranges.

###  6.1 The thrust-axis origin-peak jet correlations with simulation

As pointed out in section 5, the result with Belle data in the thrust axis
analysis lacks a significant origin-peak correlation. We study this behavior
with the Sherpa event generator, simulating $e^{+}e^{-}$ annihilation events
with collision energies ranging from 10.52 $\mathrm{\,Ge\kern-1.00006ptV}$ up
to the $Z$-pole energy. The corresponding results are shown in figure 11.
Initial state radiation is turned off in this simulation. For different
collision energies, radiation photons’ momentum kicks have different levels of
recoil effects on $q\bar{q}$ fragmentation events, which is an additional
unwanted complex factor for the correlation trend study as a function of
fragmenting quarks’ energy. However, a consistent evolving trend of the
magnitude of the origin-peak correlation is seen with increasing collision
energy for every multiplicity interval. Under the picture in the thrust axis
analysis, the main contribution to the origin-peak correlation is contributed
from jets fragmenting in the transverse direction of the event thrust axis.
Increasing the collision energy can also make these jets more energetic.

Figure 11: Thrust-axis two-particle correlation functions and the short-range
($0\leq|\Delta\eta|<1$) associated yield projections (left two columns) for
$e^{+}e^{-}$ annihilation simulated by sherpa for different collision
energies. The increasing trend of origin-peak correlations relative to the
collision energy is observed regardless of the event multiplicity, and ${\rm
N}_{\rm trk}^{\rm rec}\in[6,10)$ is shown as representative here. For
comparison, beam-axis two-particle correlations under each energy benchmark
are shown in the right two columns.

###  6.2 The thrust-axis long-range near-side enhancement of on-resonance data

We observe an enhancement of long-range near-side correlations with the on-
resonance data, as reported in section 5. MC simulation also reproduces a
similar order of long-range near-side excess; hence, we study this effect with
the on-resonance Belle MC sample to further understand the impact on two-
particle correlation functions arising from $B$-decay events. We show for
consistency the two-particle correlations and ZYAM-subtracted associated
yields with separated $B$-decay and $q\bar{{q}}$ fragmenting processes in
figure 12, based on the truth level MC. As per expectation, the excess is
mainly originated from $B$-decay. However, as mentioned in section 5, the
enhanced structures seen for these $B$-decay events (or for that matter the
on-resonance sample) have an intrinsic difference in shape from those reported
with heavy ions collisions. Recall for the latter, the typical ridge structure
is an elongated feature spanning a wide pseudorapidity range, but the
enhancement seen in the $B$-decay event sample are lumps only located at the
long range regions. On the other hand, an anisotropy in azimuthal angle for
charged particle distributions is also not anticipated in the $B$ decay
system, which is believed to be well described by the decay relations and
energy-momentum conservation. To understand the long-range near-side
enhancement in the $B$-decay events under the thrust axis analysis, a
decomposition study of the correlation contribution among two $B$ mesons is
performed. In the following, we examine the correlation function with track
pairs selected specifically from:

  1. 1.

Single $B$ meson decay products:
The associated particle and the trigger particle should be descendent from the
same $B$ meson.

  2. 2.

Different $B$ meson decay products:
The associated particle and the trigger particle should be descendant from
opposite $B$ mesons.

Figure 12: Two-particle correlation functions and long-range
($1.5\leq|\Delta\eta|<3.0$) ZYAM-subtracted azimuthal differential associated
yields for separated MC sources with the Belle MC. Panels (1a)-(1d) are
results with $B$-meson decays and panels (2a)-(2d) are those from $q\bar{q}$
fragmentation, in which (a) and (b) show the lower multiplicity results with
$6\leq N_{\rm trk}<10$; (c) and (d) show the higher multiplicity ones with
$N_{\rm trk}\geq 12$. The long-range near-side enhancements are seen
specifically in the $B$ decays MC sample in both the low- and high-
multiplicity cases.

In figure 13, distinct patterns in the two-particle correlation functions and
the long-range ZYAM-subtracted yields are observed under these two
circumstances. The sizable enhancement in the long-range near-side region is
seen exclusively in the correlation function calculated by track pairs from
different $B$ mesons. A complementary exercise with a toy sample is checked in
figure 13(3a)-13(3b). The toy sample is blended with two $B$ mesons picked up
from uncorrelated events, and the event thrust axis is recalculated with the
new mixture of constituents. The calculation of the correlation function with
track pairs from different $B$ mesons’ decay products using this toy sample
dataset guarantees that there is no special physical correlation across these
two $B$ meson’s decay products. The only correlation in the blended toy sample
is the momentum balance for the individual $B$ mesons’ decays. Therefore, the
result concerns only the decay topology and the alignment of the thrust axis
with respect to two $B$ decay systems. The blended toy sample replicates the
correlation function structure seen in crossing contribution in physical
$B$-decay events (figure 13(2a)-13(2b)). We provide a schematic explanation of
enhanced near-side correlations under such $B$-meson-pair-decay events in
figure 14. The derived event thrust axis for the two-$B$-decay event has a
large probability to lie in the plane formed by two $B$ mesons’ main decay
directions. The special circumstances further makes most of the constituents
sit on either $\phi$ or $\pi+\phi$, yielding aggregated correlation at
$\Delta\phi\approx 0$ or $\Delta\phi\approx\pi$. It is concluded that the
sculpted ridge signal in the $B$-decay events (or the on-resonance sample) is
a consequence of the special event topology and thrust axis alignment in two
$B$ meson decay systems.

Figure 13: Two-particle correlation functions and long-range
($1.5\leq|\Delta\eta|<3.0$) ZYAM-subtracted azimuthal differential associated
yields for tracks (1) from the same or (2) across different $B$ mesons for
every $\mathchar 28935\relax(4S)\rightarrow B\bar{B}$ event with the Belle MC.
Panel (3) shows the pairing yields counted from different $B$ mesons from the
toy sample mixed from uncorrelated events. Figure 14:  A schematic diagram
indicating how a special arrangement of two $B$ decay system can have a long-
range near-side correlation. The orange cones represent the decay products
from one of the $B$ mesons, and the blue cones are from the other $B$ decay.
Two $B$ decay directions are expressed with notations $\hat{e}_{B,1}$,
$\hat{e}_{B,2}$. The event thrust direction $\hat{n}$ is prone to lie on the
plane spanned by $\hat{e}_{B,1}$, $\hat{e}_{B,2}$, hence dividing the decay
products into two parts, aggregating at azimuthal coordinates around $\phi$ or
$\pi+\phi$, respectively. In this configuration, the decay products from the
opposite $B$ mesons (the orange and blue cones) easily contribute to the long-
range near-side correlation.

##  7 Conclusions

This paper details the methodology and results of the two-particle
correlations analysis in thrust axis coordinates. The method has been first
utilized in the measurement with ALEPH data Badea:2019vey  and subsequently
adopted in the Belle measurement Belle:2022fvl . To express particle
correlations for a randomly oriented event topology, such as arising in the
$e^{+}e^{-}$ collision system, we suggest that this thrust axis analysis is
the more appropriate coordinate system compared to conventional fixed beam
axis coordinates. The thrust axis provides an estimate of the
$e^{+}e^{-}\rightarrow q\bar{{q}}$ direction in the hard process. Constructed
upon the thrust-axis, two-particle correlations for $e^{+}e^{-}\rightarrow
q\bar{{q}}$ events are an observable sensitive to higher-order soft emission
transverse to the leading dijet axis. In a broader view, we report full
results of two-particle correlations in both the beam and thrust axis
coordinate systems performed on datasets of $e^{+}e^{-}$ to hadronic final
states at $\sqrt{s}=10.52$ GeV and $\sqrt{s}$ =10.58 GeV. This amplifies the
preceding Belle publication Belle:2022fvl , which focuses on the pure
$e^{+}e^{-}\rightarrow q\bar{{q}}$ fragmentation dataset at $\sqrt{s}=10.52$
GeV. As already pointed out in that Letter, the thrust-axis correlation
function distribution is different from measurements in hadron collisions and
high-energy $e^{+}e^{-}$ collisions in its shape, with the significant near-
side peak structure, interpreted as intra-jet correlation, not being seen. In
this paper, we research qualitatively the correlation structure taking
advantage of simulations. By studying with the sherpa event generator, we
characterize the evolution of the near-side-peak correlation magnitude from
low-energy collisions towards the high collision energies. The two-particle
correlation measurement with the $\mathchar 28935\relax(4S)$ on-resonance
dataset is presented for the first time. A low-scaled long-range near-side
enhancement is observed in data for the thrust-axis two-particle correlations.
However, this enhancement is different from the more familiar elongated ridge
structure over a broad $\Delta\eta$ range as reported in $pp$ collisions and
heavy ion experiments. Monte Carlo simulations are able to qualitatively
replicate these new features of the on-resonance correlation data. It is
concluded that such a special topological arrangement within a decay system
can give rise to such correlations.

## Acknowledgements

This work, based on data collected using the Belle detector, which was
operated until June 2010, was supported by the Ministry of Education, Culture,
Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the
Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of
Nagoya University; the Australian Research Council including grants
DP180102629, DP170102389, DP170102204, DE220100462, DP150103061, FT130100303;
Austrian Federal Ministry of Education, Science and Research (FWF) and FWF
Austrian Science Fund No. P 31361-N36; the National Natural Science Foundation
of China under Contracts No. 11675166, No. 11705209; No. 11975076; No.
12135005; No. 12175041; No. 12161141008; Key Research Program of Frontier
Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the
Ministry of Education, Youth and Sports of the Czech Republic under Contract
No. LTT17020; the Czech Science Foundation Grant No. 22-18469S; Horizon 2020
ERC Advanced Grant No. 884719 and ERC Starting Grant No. 947006 “InterLeptons”
(European Union); the Carl Zeiss Foundation, the Deutsche
Forschungsgemeinschaft, the Excellence Cluster Universe, and the
VolkswagenStiftung; the Department of Atomic Energy (Project Identification
No. RTI 4002) and the Department of Science and Technology of India; the
Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation
(NRF) of Korea Grant Nos. 2016R1D1A1B02012900, 2018R1A2B3003643,
2018R1A6A1A06024970, RS202200197659, 2019R1I1A3A01058933, 2021R1A6A1A03043957,
2021R1F1A1060423, 2021R1F1A1064008, 2022R1A2C1003993; Radiation Science
Research Institute, Foreign Large-size Research Facility Application
Supporting project, the Global Science Experimental Data Hub Center of the
Korea Institute of Science and Technology Information and KREONET/GLORIAD; the
Polish Ministry of Science and Higher Education and the National Science
Center; the Ministry of Science and Higher Education of the Russian
Federation, Agreement 14.W03.31.0026, and the HSE University Basic Research
Program, Moscow; University of Tabuk research grants S-1440-0321, S-0256-1438,
and S-0280-1439 (Saudi Arabia); the Slovenian Research Agency Grant Nos.
J1-9124 and P1-0135; Ikerbasque, Basque Foundation for Science, Spain; the
Swiss National Science Foundation; the Ministry of Education and the Ministry
of Science and Technology of Taiwan; and the United States Department of
Energy and the National Science Foundation. These acknowledgements are not to
be interpreted as an endorsement of any statement made by any of our
institutes, funding agencies, governments, or their representatives. We thank
the KEKB group for the excellent operation of the accelerator; the KEK
cryogenics group for the efficient operation of the solenoid; and the KEK
computer group and the Pacific Northwest National Laboratory (PNNL)
Environmental Molecular Sciences Laboratory (EMSL) computing group for strong
computing support; and the National Institute of Informatics, and Science
Information NETwork 6 (SINET6) for valuable network support.

## References

  * (1) Belle collaboration, Y. C. Chen et al., _Measurement of Two-Particle Correlations of Hadrons in $e^{+}e^{-}$ Collisions at Belle_, _Phys. Rev. Lett._ 128 (2022) 142005, [2201.01694].
  * (2) STAR collaboration, J. Adams et al., _Distributions of charged hadrons associated with high transverse momentum particles in $pp$ and $\mathrm{Au}+\mathrm{Au}$ collisions at $\sqrt{{s}_{NN}}=200\text{ }\mathrm{GeV}$_, _Phys. Rev. Lett._ 95 (2005) 152301, [nucl-ex/0501016].
  * (3) STAR collaboration, B. I. Abelev et al., _Long range rapidity correlations and jet production in high energy nuclear collisions_ , _Phys. Rev. C_ 80 (2009) 064912, [0909.0191].
  * (4) PHOBOS collaboration, B. Alver et al., _High transverse momentum triggered correlations over a large pseudorapidity acceptance in $\mathrm{Au}+\mathrm{Au}$ collisions at $\sqrt{{s}_{NN}}=200\text{ }\mathrm{GeV}$_, _Phys. Rev. Lett._ 104 (2010) 062301, [0903.2811].
  * (5) CMS collaboration, S. Chatrchyan et al., _Centrality dependence of dihadron correlations and azimuthal anisotropy harmonics in $\mathrm{PbPb}$ collisions at $\sqrt{s_{NN}}=2.76\text{ }\mathrm{TeV}$_, _Eur. Phys. J._ C72 (2012) 2012, [1201.3158].
  * (6) ALICE collaboration, K. Aamodt et al., _Harmonic decomposition of two-particle angular correlations in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76\text{ }\mathrm{TeV}$_, _Phys. Lett. B_ 708 (2012) 249–264, [1109.2501].
  * (7) STAR collaboration, J. Adam et al., _Azimuthal Harmonics in Small and Large Collision Systems at RHIC Top Energies_ , _Phys. Rev. Lett._ 122 (2019) 172301, [1901.08155].
  * (8) W. Busza, K. Rajagopal and W. van der Schee, _Heavy Ion Collisions: The Big Picture, and the Big Questions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 339–376, [1802.04801].
  * (9) A. Ortiz Velasquez, P. Christiansen, E. Cuautle Flores, I. Maldonado Cervantes and G. Paić, _Color Reconnection and Flowlike Patterns in $pp$ Collisions_, _Phys. Rev. Lett._ 111 (2013) 042001, [1303.6326].
  * (10) C. Bierlich, S. Chakraborty, G. Gustafson and L. Lönnblad, _Setting the string shoving picture in a new frame_ , _JHEP_ 03 (2021) 270, [2010.07595].
  * (11) A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian, T. Lappi and R. Venugopalan, _The Ridge in proton-proton collisions at the LHC_ , _Phys. Lett. B_ 697 (2011) 21–25, [1009.5295].
  * (12) CMS collaboration, V. Khachatryan et al., _Observation of Long-Range Near-Side Angular Correlations in Proton-Proton Collisions at the LHC_ , _JHEP_ 09 (2010) 091, [1009.4122].
  * (13) CMS collaboration, S. Chatrchyan et al., _Observation of long-range near-side angular correlations in proton-lead collisions at the LHC_ , _Phys. Lett._ B718 (2013) 795–814, [1210.5482].
  * (14) LHCb collaboration, R. Aaij et al., _Measurements of long-range near-side angular correlations in $\sqrt{s_{\text{NN}}}=5\text{ }\mathrm{TeV}$ proton-lead collisions in the forward region_, _Phys. Lett._ B762 (2016) 473–483, [1512.00439].
  * (15) ATLAS collaboration, G. Aad et al., _Observation of Long-Range Elliptic Azimuthal Anisotropies in $\sqrt{s}=$13 and 2.76 TeV $pp$ Collisions with the ATLAS Detector_, _Phys. Rev. Lett._ 116 (2016) 172301, [1509.04776].
  * (16) ALICE collaboration, B. Abelev et al., _Long-range angular correlations on the near and away side in $p$-Pb collisions at $\sqrt{s_{NN}}=5.02\text{ }\mathrm{TeV}$_, _Phys. Lett._ B719 (2013) 29–41, [1212.2001].
  * (17) ATLAS collaboration, G. Aad et al., _Observation of Associated Near-Side and Away-Side Long-Range Correlations in $\sqrt{s_{NN}}$=5.02 TeV Proton-Lead Collisions with the ATLAS Detector_, _Phys. Rev. Lett._ 110 (2013) 182302, [1212.5198].
  * (18) PHENIX collaboration, A. Adare et al., _Quadrupole Anisotropy in Dihadron Azimuthal Correlations in Central $d+$Au Collisions at $\sqrt{s_{{}_{NN}}}$=200 GeV_, _Phys. Rev. Lett._ 111 (2013) 212301, [1303.1794].
  * (19) J.-Y. Ollitrault, _Anisotropy as a signature of transverse collective flow_ , _Phys. Rev._ D46 (1992) 229–245.
  * (20) B. Alver and G. Roland, _Collision geometry fluctuations and triangular flow in heavy-ion collisions_ , _Phys. Rev._ C81 (2010) 054905, [1003.0194].
  * (21) K. Dusling and R. Venugopalan, _Comparison of the color glass condensate to dihadron correlations in proton-proton and proton-nucleus collisions_ , _Phys. Rev._ D87 (2013) 094034, [1302.7018].
  * (22) P. Bozek, _Collective flow in p-Pb and d-Pd collisions at TeV energies_ , _Phys. Rev._ C85 (2012) 014911, [1112.0915].
  * (23) L. He, T. Edmonds, Z.-W. Lin, F. Liu, D. Molnar and F. Wang, _Anisotropic parton escape is the dominant source of azimuthal anisotropy in transport models_ , _Phys. Lett._ B753 (2016) 506–510, [1502.05572].
  * (24) J. L. Nagle and W. A. Zajc, _Small System Collectivity in Relativistic Hadronic and Nuclear Collisions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 211–235, [1801.03477].
  * (25) ZEUS collaboration, I. Abt et al., _Two-particle azimuthal correlations as a probe of collective behaviour in deep inelastic $ep$ scattering at HERA_, _JHEP_ 04 (2020) 070, [1912.07431].
  * (26) J. L. Nagle, R. Belmont, K. Hill, J. Orjuela Koop, D. V. Perepelitsa, P. Yin et al., _Minimal conditions for collectivity in $e^{+}e^{-}$ and $p+p$ collisions_, _Phys. Rev._ C97 (2018) 024909, [1707.02307].
  * (27) A. Badea, A. Baty, P. Chang, G. M. Innocenti, M. Maggi, C. Mcginn et al., _Measurements of two-particle correlations in $e^{+}e^{-}$ collisions at 91 GeV with ALEPH archived data_, _Phys. Rev. Lett._ 123 (2019) 212002, [1906.00489].
  * (28) T. Altinoluk and N. Armesto, _Particle correlations from the initial state_ , _Eur. Phys. J. A_ 56 (2020) 215, [2004.08185].
  * (29) P. Castorina, D. Lanteri and H. Satz, _Strangeness enhancement and flow-like effects in $e^{+}e^{-}$ annihilation at high parton density_, _Eur. Phys. J. A_ 57 (2021) 111, [2011.06966].
  * (30) P. Agostini, T. Altinoluk and N. Armesto, _Multi-particle production in proton–nucleus collisions in the color glass condensate_ , _Eur. Phys. J. C_ 81 (2021) 760, [2103.08485].
  * (31) A. Abashian et al., _The Belle Detector_ , _Nucl. Instrum. Meth. A_ 479 (2002) 117–232.
  * (32) Belle collaboration, J. Brodzicka et al., _Physics Achievements from the Belle Experiment_ , _PTEP_ 2012 (2012) 04D001, [1212.5342].
  * (33) A. Ryd, D. Lange, N. Kuznetsova, S. Versille, M. Rotondo, D. P. Kirkby et al., _EvtGen: A Monte Carlo Generator for B-Physics_ , EVTGEN-V00-11-07.
  * (34) T. Sjöstrand, P. Edén, C. Friberg, L. Lönnblad, G. Miu, S. Mrenna et al., _High-energy-physics event generation with PYTHIA 6.1_ , _Comput. Phys. Commun._ 135 (Apr., 2001) 238–259, [hep-ph/0010017].
  * (35) Belle collaboration, K. Abe et al., _Measurement of inclusive production of neutral pions from $\mathchar 28935\relax(4S)$ decays_, _Phys. Rev. D_ 64 (2001) 072001, [hep-ex/0103041].
  * (36) E. Farhi, _Quantum chromodynamics test for jets_ , _Phys. Rev. Lett._ 39 (Dec, 1977) 1587–1588.
  * (37) CMS collaboration, S. Chatrchyan et al., _Long-range and short-range dihadron angular correlations in central PbPb collisions at a nucleon-nucleon center of mass energy of 2.76 TeV_ , _JHEP_ 07 (2011) 076, [1105.2438].
  * (38) CMS collaboration, S. Chatrchyan et al., _Multiplicity and Transverse Momentum Dependence of Two- and Four-Particle Correlations in pPb and PbPb Collisions_ , _Phys. Lett. B_ 724 (2013) 213–240, [1305.0609].
  * (39) CMS collaboration, V. Khachatryan et al., _Long-range two-particle correlations of strange hadrons with charged particles in pPb and PbPb collisions at LHC energies_ , _Phys. Lett. B_ 742 (2015) 200–224, [1409.3392].
  * (40) ATLAS collaboration, G. Choudalakis, _Unfolding in ATLAS_ , in _PHYSTAT 2011_ , (Geneva), CERN, 4, 2011, 1104.2962, DOI.
  * (41) Ed. A.J. Bevan, B. Golob, Th. Mannel, S. Prell, and B.D. Yabsley, _The Physics of the B Factories_ , _Eur. Phys. J. C_ 74 (2014) 3026, [1406.6311].
  * (42) N. N. Ajitanand, J. M. Alexander, P. Chung, W. G. Holzmann, M. Issah, R. A. Lacey et al., _Decomposition of harmonic and jet contributions to particle-pair correlations at ultra-relativistic energies_ , _Phys. Rev._ C72 (2005) 011902, [nucl-ex/0501025].
  * (43) B. Efron, _Bootstrap methods: Another look at the jackknife_ , _Ann. Statist._ 7 (01, 1979) 1–26.
  * (44) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159–177, [1410.3012].
  * (45) J. Bellm et al., _Herwig 7.0/Herwig++ 3.0 release note_ , _Eur. Phys. J. C_ 76 (2016) 196, [1512.01178].
  * (46) T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert et al., _Event generation with SHERPA 1.1_ , _JHEP_ 02 (2009) 007, [0811.4622].
  * (47) P. Skands, S. Carrazza and J. Rojo, _Tuning PYTHIA 8.1: the Monash 2013 Tune_ , _Eur. Phys. J. C_ 74 (2014) 3024, [1404.5630].
  * (48) Particle Data Group collaboration, P. A. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.