This is the accepted manuscript made available via CHORUS. The article has been
published as:

Multiparticle azimuthal correlations for extracting event-byevent elliptic and triangular flow in Au + Au collisions at
sqrt[s_{NN}]=200 GeV
A. Adare et al. (PHENIX Collaboration)
Phys. Rev. C 99, 024903 — Published 6 February 2019
DOI: 10.1103/PhysRevC.99.024903

1
2

Multi-particle azimuthal correlations for extracting event-by-event elliptic and√triangular flow in Au+Au collisions at sN N = 200 GeV

A. Adare,13 C. Aidala,41 N.N. Ajitanand,58, ∗ Y. Akiba,53, 54, † M. Alfred,23 N. Apadula,28, 59 H. Asano,34, 53
4B. Azmoun,7 V. Babintsev,24 A. Bagoly,17 M. Bai,6 N.S. Bandara,40 B. Bannier,59 K.N. Barish,8 S. Bathe,5, 54
5A. Bazilevsky,7 M. Beaumier,8 S. Beckman,13 R. Belmont,13, 41, 47 A. Berdnikov,56 Y. Berdnikov,56 D.S. Blau,33, 44
6
M. Boer,36 J.S. Bok,46 K. Boyle,54 M.L. Brooks,36 J. Bryslawskyj,5, 8 V. Bumazhnov,24 S. Campbell,14, 28
7
V. Canoa Roman,59 C.-H. Chen,54 C.Y. Chi,14 M. Chiu,7 I.J. Choi,25 J.B. Choi,10, ∗ T. Chujo,62 Z. Citron,64
8
M. Connors,21, 54 M. Csanád,17 T. Csörgő,18, 65 T.W. Danley,48 A. Datta,45 M.S. Daugherity,1 G. David,7, 59
9
K. DeBlasio,45 K. Dehmelt,59 A. Denisov,24 A. Deshpande,54, 59 E.J. Desmond,7 A. Dion,59 P.B. Diss,39
10
J.H. Do,66 A. Drees,59 K.A. Drees,6 J.M. Durham,36 A. Durum,24 A. Enokizono,53, 55 S. Esumi,62 B. Fadem,42
11
W. Fan,59 N. Feege,59 D.E. Fields,45 M. Finger,9 M. Finger, Jr.,9 S.L. Fokin,33 J.E. Frantz,48 A. Franz,7
12
A.D. Frawley,20 C. Gal,59 P. Gallus,15 P. Garg,3, 59 H. Ge,59 F. Giordano,25 A. Glenn,35 Y. Goto,53, 54 N. Grau,2
13
S.V. Greene,63 M. Grosse Perdekamp,25 T. Gunji,12 T. Hachiya,43, 53, 54 J.S. Haggerty,7 K.I. Hahn,19
14
H. Hamagaki,12 H.F. Hamilton,1 S.Y. Han,19 J. Hanks,59 S. Hasegawa,29 T.O.S. Haseler,21 K. Hashimoto,53, 55
15
X. He,21 T.K. Hemmick,59 J.C. Hill,28 K. Hill,13 A. Hodges,21 R.S. Hollis,8 K. Homma,22 B. Hong,32 T. Hoshino,22
16
N. Hotvedt,28 J. Huang,7 S. Huang,63 K. Imai,29 M. Inaba,62 A. Iordanova,8 D. Isenhower,1 D. Ivanishchev,52
17
B.V. Jacak,59 M. Jezghani,21 Z. Ji,59 J. Jia,7, 58 X. Jiang,36 B.M. Johnson,7, 21 D. Jouan,50 D.S. Jumper,25
18
S. Kanda,12 J.H. Kang,66 D. Kawall,40 A.V. Kazantsev,33 J.A. Key,45 V. Khachatryan,59 A. Khanzadeev,52
19
C. Kim,32 D.J. Kim,30 E.-J. Kim,10 G.W. Kim,19 M. Kim,57 B. Kimelman,42 D. Kincses,17 E. Kistenev,720
R. Kitamura,12 J. Klatsky,20 D. Kleinjan,8 P. Kline,59 T. Koblesky,13 B. Komkov,52 D. Kotov,52, 56 B. Kurgyis,17
21
K. Kurita,55 M. Kurosawa,53, 54 Y. Kwon,66 R. Lacey,58 J.G. Lajoie,28 A. Lebedev,28 S. Lee,66 S.H. Lee,28, 59
22
M.J. Leitch,36 Y.H. Leung,59 N.A. Lewis,41 X. Li,11 X. Li,36 S.H. Lim,36, 66 M.X. Liu,36 S. Lökös,17, 18 D. Lynch,723
T. Majoros,16 Y.I. Makdisi,6 M. Makek,67 A. Manion,59 V.I. Manko,33 E. Mannel,7 M. McCumber,36
24
P.L. McGaughey,36 D. McGlinchey,13, 36 C. McKinney,25 A. Meles,46 M. Mendoza,8 A.C. Mignerey,39
25
D.E. Mihalik,59 A. Milov,64 D.K. Mishra,4 J.T. Mitchell,7 G. Mitsuka,31, 54 S. Miyasaka,53, 61 S. Mizuno,53, 62
26
A.K. Mohanty,4 P. Montuenga,25 T. Moon,66 D.P. Morrison,7 S.I. Morrow,63 T.V. Moukhanova,33
27
T. Murakami,34, 53 J. Murata,53, 55 A. Mwai,58 K. Nagashima,22 J.L. Nagle,13 M.I. Nagy,17 I. Nakagawa,53, 54
28
H. Nakagomi,53, 62 K. Nakano,53, 61 C. Nattrass,60 P.K. Netrakanti,4 T. Niida,62 S. Nishimura,12 R. Nouicer,7, 54
29
T. Novák,18, 65 N. Novitzky,30, 59 A.S. Nyanin,33 E. O’Brien,7 C.A. Ogilvie,28 J.D. Orjuela Koop,13 J.D. Osborn,41
37
30 A. Oskarsson,
K. Ozawa,31, 62 R. Pak,7 V. Pantuev,26 V. Papavassiliou,46 J.S. Park,57 S. Park,53, 57, 59 S.F. Pate,46
28
31 M. Patel,
J.-C. Peng,25 W. Peng,63 D.V. Perepelitsa,7, 13 G.D.N. Perera,46 D.Yu. Peressounko,33 C.E. PerezLara,5928
32
J. Perry, R. Petti,7, 59 C. Pinkenburg,7 R. Pinson,1 R.P. Pisani,7 M.L. Purschke,7 P.V. Radzevich,56 J. Rak,30
33
B.J. Ramson,41 I. Ravinovich,64 K.F. Read,49, 60 D. Reynolds,58 V. Riabov,44, 52 Y. Riabov,52, 56 D. Richford,534
T. Rinn,28 S.D. Rolnick,8 M. Rosati,28 Z. Rowan,5 J.G. Rubin,41 J. Runchey,28 B. Sahlmueller,59 N. Saito,31
35
T. Sakaguchi,7 H. Sako,29 V. Samsonov,44, 52 M. Sarsour,21 S. Sato,29 B. Schaefer,63 B.K. Schmoll,60 K. Sedgwick,853, 5436 R. Seidl,
A. Sen,28, 60 R. Seto,8 P. Sett,4 A. Sexton,39 D. Sharma,59 I. Shein,24 T.-A. Shibata,53, 61 K. Shigaki,22
28, 4337 M. Shimomura,
P. Shukla,4 A. Sickles,7, 25 C.L. Silva,36 D. Silvermyr,37, 49 B.K. Singh,3 C.P. Singh,3 V. Singh,341
38
M.J. Skoby, M. Slunečka,9 M. Snowball,36 R.A. Soltz,35 W.E. Sondheim,36 S.P. Sorensen,60 I.V. Sourikova,739
P.W. Stankus,49 M. Stepanov,40, ∗ S.P. Stoll,7 T. Sugitate,22 A. Sukhanov,7 T. Sumita,53 J. Sun,59 Z. Sun,1640
J. Sziklai,65 A. Taketani,53, 54 K. Tanida,29, 54, 57 M.J. Tannenbaum,7 S. Tarafdar,63, 64 A. Taranenko,44, 58
41
R. Tieulent,21, 38 A. Timilsina,28 T. Todoroki,53, 54, 62 M. Tomášek,15 C.L. Towell,1 R. Towell,1 R.S. Towell,164
42 I. Tserruya,
Y. Ueda,22 B. Ujvari,16 H.W. van Hecke,36 J. Velkovska,63 M. Virius,15 V. Vrba,15, 27 X.R. Wang,46, 54
43
Y. Watanabe,53, 54 Y.S. Watanabe,12, 31 F. Wei,46 A.S. White,41 C.P. Wong,21 C.L. Woody,7 M. Wysocki,49
44
B. Xia,48 C. Xu,46 Q. Xu,63 L. Xue,21 S. Yalcin,59 Y.L. Yamaguchi,12, 54, 59 A. Yanovich,24 J.H. Yoo,32
45
I. Yoon,57 H. Yu,46, 51 I.E. Yushmanov,33 W.A. Zajc,14 A. Zelenski,6 S. Zharko,56 S. Zhou,11 and L. Zou8
46
(PHENIX Collaboration)
3

47
48
49
50
51
52
53

1

Abilene Christian University, Abilene, Texas 79699, USA
Department of Physics, Augustana University, Sioux Falls, South Dakota 57197, USA
3
Department of Physics, Banaras Hindu University, Varanasi 221005, India
4
Bhabha Atomic Research Centre, Bombay 400 085, India
5
Baruch College, City University of New York, New York, New York, 10010 USA
6
Collider-Accelerator Department, Brookhaven National Laboratory, Upton, New York 11973-5000, USA
7
Physics Department, Brookhaven National Laboratory, Upton, New York 11973-5000, USA228

University of California-Riverside, Riverside, California 92521, USA
Charles University, Ovocný trh 5, Praha 1, 116 36, Prague, Czech Republic
1056
Chonbuk National University, Jeonju, 561-756, Korea
11
57
Science and Technology on Nuclear Data Laboratory, China Institute58
of Atomic Energy, Beijing 102413, People’s Republic of China
1259
Center for Nuclear Study, Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan
13
60
University of Colorado, Boulder, Colorado 80309, USA
1461
Columbia University, New York, New York 10027 and Nevis Laboratories, Irvington, New York 10533, USA
15
62
Czech Technical University, Zikova 4, 166 36 Prague 6, Czech Republic
16
63
Debrecen University, H-4010 Debrecen, Egyetem tér 1, Hungary
17
64
ELTE, Eötvös Loránd University, H-1117 Budapest, Pázmány P. s. 1/A, Hungary
18
65
Eszterházy Károly University, Károly Róbert Campus, H-3200 Gyöngyös, Mátrai út 36, Hungary
1966
Ewha Womans University, Seoul 120-750, Korea
20
67
Florida State University, Tallahassee, Florida 32306, USA
21
68
Georgia State University, Atlanta, Georgia 30303, USA
22
69
Hiroshima University, Kagamiyama, Higashi-Hiroshima 739-8526, Japan
23
70
Department of Physics and Astronomy, Howard University, Washington, DC 20059, USA
24
71
IHEP Protvino, State Research Center of Russian Federation, Institute for High Energy Physics, Protvino, 142281, Russia
25
72
University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA
26
73
Institute for Nuclear Research of the Russian Academy of Sciences, prospekt 60-letiya Oktyabrya 7a, Moscow 117312, Russia
27
74Institute of Physics, Academy of Sciences of the Czech Republic, Na Slovance 2, 182 21 Prague 8, Czech Republic
28
75
Iowa State University, Ames, Iowa 50011, USA
29
76
Advanced Science Research Center, Japan Atomic Energy Agency, 2-4
77Shirakata Shirane, Tokai-mura, Naka-gun, Ibaraki-ken 319-1195, Japan
30
78
Helsinki Institute of Physics and University of Jyväskylä, P.O.Box 35, FI-40014 Jyväskylä, Finland
31
79
KEK, High Energy Accelerator Research Organization, Tsukuba, Ibaraki 305-0801, Japan
32
80
Korea University, Seoul 02841, Korea
33
81
National Research Center “Kurchatov Institute”, Moscow, 123098 Russia
34
82
Kyoto University, Kyoto 606-8502, Japan
35
83
Lawrence Livermore National Laboratory, Livermore, California 94550, USA
36
84
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
37
85
Department of Physics, Lund University, Box 118, SE-221 00 Lund, Sweden
38
86
IPNL, CNRS/IN2P3, Univ Lyon, Universit Lyon 1, F-69622, Villeurbanne, France
39
87
University of Maryland, College Park, Maryland 20742, USA
40
88
Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003-9337, USA
41
89
Department of Physics, University of Michigan, Ann Arbor, Michigan 48109-1040, USA4290
Muhlenberg College, Allentown, Pennsylvania 18104-5586, USA
43
91
Nara Women’s University, Kita-uoya Nishi-machi Nara 630-8506, Japan
44
92
National Research Nuclear University, MEPhI, Moscow Engineering Physics Institute, Moscow, 115409, Russia
45
93
University of New Mexico, Albuquerque, New Mexico 87131, USA
46
94
New Mexico State University, Las Cruces, New Mexico 88003, USA
47
95
Physics and Astronomy Department, University of North Carolina at Greensboro, Greensboro, North Carolina 27412, USA
48
96
Department of Physics and Astronomy, Ohio University, Athens, Ohio 45701, USA
49
97
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
50
98
IPN-Orsay, Univ. Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, BP1, F-91406, Orsay, France
51
99
Peking University, Beijing 100871, People’s Republic of China
52
100
PNPI, Petersburg Nuclear Physics Institute, Gatchina, Leningrad region, 188300, Russia
53
101
RIKEN Nishina Center for Accelerator-Based Science, Wako, Saitama 351-0198, Japan
54
102
RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973-5000, USA
55
103
Physics Department, Rikkyo University, 3-34-1 Nishi-Ikebukuro, Toshima, Tokyo 171-8501, Japan
56
104
Saint Petersburg State Polytechnic University, St. Petersburg, 195251 Russia
57
105
Department of Physics and Astronomy, Seoul National University, Seoul 151-742, Korea
58
106
Chemistry Department, Stony Brook University, SUNY, Stony Brook, New York 11794-3400, USA
59
107
Department of Physics and Astronomy, Stony Brook University, SUNY, Stony Brook, New York 11794-3800, USA
60
108
University of Tennessee, Knoxville, Tennessee 37996, USA
61
109
Department of Physics, Tokyo Institute of Technology, Oh-okayama, Meguro, Tokyo 152-8551, Japan
62
110
Tomonaga Center for the History of the Universe, University of Tsukuba, Tsukuba, Ibaraki 305, Japan
63
111
Vanderbilt University, Nashville, Tennessee 37235, USA
64
112
Weizmann Institute, Rehovot 76100, Israel
65
113
Institute for Particle and Nuclear Physics, Wigner Research Centre for Physics, Hungarian
114
Academy of Sciences (Wigner RCP, RMKI) H-1525 Budapest 114, POBox 49, Budapest, Hungary
66
115
Yonsei University, IPAP, Seoul 120-749, Korea
67
116
Department of Physics, Faculty of Science, University of Zagreb, Bijenička c. 32 HR-10002 Zagreb, Croatia
117(Dated: December 7, 2018)
54

55

9

3
We present measurements of elliptic and triangular azimuthal anisotropy of charged particles√
detected at forward rapidity 1 < |η| < 3 in Au+Au collisions at sN N = 200 GeV, as a functionof centrality. The multiparticle cumulant technique is used to obtain the elliptic flow coefficientsv2 {2}, v2 {4}, v2 {6}, and v2 {8}, and triangular flow coefficients v3 {2} and v3 {4}. Using the smallvariance limit, we estimate the mean and variance of the event-by-event v2 distribution from v2 {2}and v2 {4}. In a complementary analysis, we also use a folding procedure to study the distributionsof v2 and v3 directly, extracting both the mean and variance. Implications for initial geometricalfluctuations and their translation into the final state momentum distributions are discussed.

∗ Deceased† PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov

4
I.

118

INTRODUCTION

Collisions of heavy nuclei at ultra-relativistic energies are believed to create a state of matter called the strongly
coupled quark-gluon plasma, as first observed at the Relativistic Heavy Ion Collider (RHIC) [1–4]. The quark-gluon
121 plasma evolves hydrodynamically as a nearly perfect liquid as evinced by the wealth of experimental measurements
122 and theoretical predictions (or descriptions) of the azimuthal anisotropy of the produced particles. [5]. Multi-particle
123 correlations are generally taken as strong evidence of hydrodynamical flow, which necessarily affects most or all
124 particles in the event [6]. This is different from mimic correlations (generically called nonflow) that are not related to
125 the hydrodynamical evolution and typically involve only a few particles.
126
Multi-particle correlations are also interesting because they have different sensitivities to the underlying event-by127 event fluctuations, which can provide additional insights into the initial geometry and its translation into final state
128 particle distributions [7, 8].
129
Recently, experimental and theoretical efforts have been directed towards measuring the fluctuations directly, using
130 event-by-event unfolding techniques. In principle, the multi-particle correlations and unfolding techniques provide
131 the same information about the underlying fluctuations, though in practice with different sensitivities [9].
The
132 techniques used at the Large Hadron Collider (LHC) are experimentally very different and provide complementary
133 information [10, 11].
134
In this manuscript we present measurements of 2-, 4-, 6-, and 8-particle correlations as well as event-by-event
135 measurements of the azimuthal anisotropy parameters corresponding to elliptic v2 and triangular v3 flow. We estimate
136 the relationship between the mean and variance with both techniques and discuss the implications for understanding
137 the detailed shape of the v2 and v3 distributions. These measurements, while the first of their kind at forward rapidity,
138 are consistent with previous measurements at midrapidity by STAR [12] and PHOBOS [13].
119
120

139

II.

EXPERIMENTAL SETUP

In 2014, the PHENIX experiment [14] at RHIC collected nearly 2 × 1010 minimum bias (MB) events of Au+Au√
141 collisions at a nucleon-nucleon center-of-mass energysN N = 200 GeV. The present analysis makes use of a subset
9
142 (≈ 10
events) of the total 2014 data sample. The PHENIX beam-beam counters (BBC) are used for triggering and
143 centrality determination. The BBCs [15] are located ± 144 cm from the nominal interaction point and cover the
144 full azimuth and 3.1 < |η| < 3.9 in pseudorapidity. By convention, the north side is forward rapidity (η > 0) and
145 the south side is backward rapidity (η < 0). Each BBC comprises an array of 64 phototubes with a fused quartz
146 Čerenkov radiator on the front. Charged particles impinging on the radiator produce Čerenkov light which is then147 amplified and detected by the phototube. The PHENIX MB trigger for the 2014 data sample of Au+Au collisions at√
148
sN N = 200 GeV was defined by at least two phototubes in each side of the BBC having signal above threshold and
149 an online z-vertex within ± 10 cm of the nominal interaction point. Additionally, PHENIX has a set of zero-degree
150 calorimeters (ZDC) that measure spectator neutrons from each incoming nucleus [15]. We require a minimum energy
151 in both ZDCs to remove beam related background present at the highest luminosities.
152
The centrality definition is based on the combined signal in the north and south BBCs. The charge distribution is
153 fitted using a Monte Carlo (MC) Glauber [16] simulation to estimate the number of participating nucleons (Npart ) and
154 a negative binomial distribution to describe the BBC signal for fixed Npart . All quantities in the present manuscript
155 are reported as a function of centrality and the corresponding Npart values are shown in Table I.
156
The main detector used in the analysis is the forward silicon vertex detector (FVTX). The FVTX [17] is a silicon
157 strip detector comprising two arms, north and south, covering 1 < |η| < 3. In Au+Au collisions there is a strong
158 correlation between the total signal in the BBCs and the total number of tracks in the FVTX. To remove beam related
159 background, we apply an additional event selection on the correlation between the total BBC signal and the number
160 of tracks in the FVTX.
161
Each FVTX arm has four layers. In the track reconstruction software, a minimum of three hits is required to
162 reconstruct a track. However, it is possible for there to be hit sharing with the central rapidity detector (VTX), so
163 that one or two of the three required hits can be in the VTX. We select tracks using a stricter requirement of at least
164 three hits in FVTX, irrespective of the number of hits in the VTX. We further require that the track reconstruction
2
165 algorithm have a goodness of fit of χ /d.o.f. < 5 for each track. Lastly, we require that each track has a distance
166 of closest approach (DCA) of less than 2 cm. The DCA is defined as the distance between the event vertex and the
167 straight-line extrapolation point of the FVTX track onto a plane which is perpendicular to the z-axis and contains
168 the event vertex. A 2 cm cut selects the FVTX tracks that likely originate from the event vertex, and is conservative
169 in accepting the nonzero DCA tail that stems from the uncertainty in the determination of the vertex position and
170 the bending of the actual track in the experimental magnetic field.
Due to the orientation of the FVTX strips
171 relative to the magnetic field, momentum determination is not possible using the tracks in the FVTX alone. However,140

5
TABLE I. Npart values for various centrality categories.
hNpart i
350.8 ± 3.1
301.7 ± 4.7
236.1 ± 3.8
167.6 ± 5.5
115.5 ± 5.8
76.1 ± 5.5
47.0 ± 4.7
26.7 ± 3.6
13.6 ± 2.4
6.1 ± 1.3

Centrality
0%–5%
5%–10%
10%–20%20%–30%30%–40%40%–50%50%–60%60%–70%70%–80%80%–93%

geant-4 [18] simulations have determined that the tracking efficiency is relatively independent of momentum for
pT & 0.3 GeV/c. Figure 1 shows the pT dependence of the FVTX tracking efficiency averaged over 1 < |η| < 3.
174 Figure 2 shows the tracking efficiency as a function of η in the FVTX for two different z-vertex selections. The single
175 particle tracking efficiency has a maximum value of 98.6% as a function of η. When averaging over 1 < |η| < 3, the176 maximum value of the pT -dependent efficiency is 67.9%, and pT = 0.3 GeV/c the efficiency is at 75% of its maximum
177 value.

172
173

acceptance × efficiency

0.7
0.6
0.5
0.4
0.3
|z| < 10 cm, 1 < |η| < 3

0.2
0.1
0
0

1

2

3

4

5

6

7

8 9 10p (GeV/c)T

FIG. 1. Tracking efficiency and acceptance in the FVTX as a function of pT . At pT = 0.3 GeV/c the efficiency is 75% of its
asymptotic value.

178

III.

ANALYSIS METHODS

The azimuthal distribution of particles in an event can be represented by a Fourier series [19]:XdN∝1+2vn cos(n(φ − ψn )),
dφ
n

(1)

efficiency

61-10 cm < z < -9 cm
9 cm < z < 10 cm

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
−3

−2

−1

0

1

2

3
η

FIG. 2. Tracking efficiency in the FVTX as a function of η for two different z-vertex selections.

where n is the harmonic number, φ is the azimuthal angle of some particle, ψn is the symmetry plane, and vn =hcos(n(φ − ψn ))i. There are many experimental techniques for estimating the vn coefficients, some of which we discuss
181 in this section.
The main ingredient in the present analysis, for both the cumulant results and the folding results, is the Q-vector.
The Q-vector is a complex number Qn = Qn,x + iQn,y with the components defined as

179

180

RQn = Qn,x =NX

cos nφi ,(2)

sin nφi ,(3)i

IQn = Qn,y =NXiwhere φi is the azimuthal angle of some particle and N is the number of particles in some event or subevent—a subevent
is a subset of a whole event, usually selected based on some kinematic selection, e.g. pseudorapidity. Because the
184 PHENIX FVTX detector subsystem is split into two separate arms, north (1 < η < 3) and south (−3 < η < −1), it
185 is natural to use tracks in the two arms as separate subevents for some calculations. In other calculations, all tracks
186 from north and south will be combined into a single event.
Additional corrections to the data are needed to account for any nonuniformity in the azimuthal acceptance. In thecase of uniform azimuthal acceptance, the event average of the Q-vector components is zero: hQn,x i = hQn,y i = 0. In
the case of nonuniform acceptance, there can be a systematic bias such that this relation does not hold. In the case of
few-particle correlations, i.e. 2- and 4-particle correlations, the bias can be corrected analytically in a straightforwardmanner [20]. In the case of correlations with a larger number of particles, however, this becomes impractical. The
total number of terms in a k-particle cumulant calculation without the assumption of azimuthal uniformity is givenby the Bell sequence: 1, 2, 5, 15, 52, 203, 877, 4140, ...—that is, the number of terms for 2- and 4-particle correlations
is a rather manageable 2 and 15, respectively; contrariwise, the number of terms for 6- and 8-particle correlations isa rather unmanageable 203 and 4140 terms, respectively. For that reason, the only practicable choice is to perform
calculations on corrected Q-vectors. The present analysis makes use of Q-vector re-centering [21]. In this procedure
one has the relation

182
183

average
Qcorrected= Qraw,nn − Qn

(4)

7
where
Qaverage= N hcos nφi + iN hsin nφi,n(5)hcos nφi = hQn,x /N i,hsin nφi = hQn,y /N i.

(6)
(7)

and

In the present analysis, we perform the Q-vector re-centering procedure for each FVTX arm separately and as
FVTX
a function of Ntracks. To assess the associated systematic uncertainty, we perform the Q-vector re-centering as a
189 function of centrality instead, as a function of additional secondary variables (event vertex and operational time
190 period), and for combined arms instead of separate.
187

188

A.

191

Cumulants

The cumulant method for flow analysis was first proposed in Ref. [22]. In the present analysis, we use the recursion
algorithm developed in Ref. [23], which is a generalization of the direct calculations using Q-vector algebra first
derived in Ref. [20]. We consider 2-, 4-, 6-, and 8- particle correlations. The multi-particle correlations are denoted
hki for k-particle correlations and are estimators for the k-th moment of vn , i.e. hki = hvnk i. In terms of the angular
relationships between different particles, they are
h2i = hcos(n(φ1 − φ2 ))i,h4i = hcos(n(φ1 + φ2 − φ3 − φ4 ))i,h6i = hcos(n(φ1 + φ2 + φ3 − φ4 − φ5 − φ6 ))i,h8i = hcos(n(φ1 + φ2 + φ3 + φ4 − φ5 − φ6 − φ7 − φ8 ))i,
192

(8)
(9)
(10)
(11)

where φ1,...,8 represent the azimuthal angles of different particles in the event.
The k-particle cumulants, denoted cn {k} are constructed in such a way that potential contributions from lower
order correlations are removed. Because the cumulants mix various terms that are of equal powers of vn , the cumulant
method vn , denoted vn {k}, is proportional to the k-th root of the cumulant. The cn {k} are constructed as follows:cn {2} = h2i,

(12)2cn {4} = h4i − 2h2i ,(13)3

cn {6} = h6i − 9h4ih2i + 12h2i ,

(14)

2

2

4

cn {8} = h8i − 16h6ih2i − 18h4i + 144h4ih2i − 144h2i ,(15)and the vn {k} arevn {2} = (cn {2})1/2 ,
1/4

(16)

vn {4} = (−cn {4})

,

(17)

vn {6} = (cn {6}/4)

1/6

(18)

,
1/8

vn {8} = (−cn {8}/33)

.

(19)

It is also possible to construct cumulants in two or more subevents, though in the present analysis we will only
concern ourselves with two subevents. For 2-particle correlations, rather than h2i = hcos(n(φ1 − φ2 ))i where φ1 andab
195 φ2 are from the same subevent, one has instead h2ia|b = hcos(n(φ1 − φ2 ))i where a, b denote two different subevents.
2
196 The cumulant and vn have the same relationship as in the single event case, i.e. vn {2}a|b = cn {2}a|b = h2ia|b . The
197 two-subevent 2-particle cumulant is also known as the scalar product method [24].
Subevent cumulants for correlations with four or more particles were first proposed in Ref. [25]. For two-subevent
4-particle correlations, there are two possibilities:

193

194

h4iab|ab = hcos(n(φa1 + φb2 − φa3 − φb4 ))i,(20)

h4iaa|bb = hcos(n(φa1 + φa2 − φb3 − φb4 ))i,

(21)

8
where the former allows 2-particle correlations within single subevents and the latter excludes them. The latter istherefore less susceptible to nonflow than the former, although both are less susceptible to nonflow than single event
4-particle correlations. The cumulants take the form
cn {4}ab|ab = h4iab|ab − h2ia|a h2ib|b − h2i2a|b ,(22)

cn {4}aa|bb = h4iaa|bb − 2h2i2a|b ,(23)

and the vn {4} values have the same relationship to the cumulants as in the single particle case, i.e. vn {4}ab|ab =(−cn {4}ab|ab )1/4 and vn {4}aa|bb = (−cn {4}aa|bb )1/4 .
200
To determine systematic uncertainties associated with event and track selection for the cumulant analysis, we vary
201 the event and track selection criteria and assess the variation on the final analysis results. The z-vertex selection
202 is modified from ± 10 cm to ± 5 cm.
The track selections are independently modified to have a goodness of fit
2
203 requirement χ /d.o.f. < 3. These changes move the cumulant results by an almost common multiplicative value and204 thus we quote the systematic uncertainties as a global scale factor uncertainty for each result.
198
199

205

B.

Folding

Here we describe an alternative approach where one utilizes the event-by-event Qn distribution to extract theevent-by-event vn distribution via an unfolding procedure. For our analysis we attempt a procedure similar to that
208 used by ATLAS as described in Ref. [10].
In brief, ATLAS successfully carries out the unfold in Pb+Pb collisions at 2.76 TeV and finds that the event-by-event
probability distribution for elliptic flow p(v2 ) is reasonably described by a Bessel-Gaussian function
206
207

vn −
p(vn ) = 2 eδ vn

2 +(v RP )2 )
(vn
n
2
2δv
n


I0

vn vnRP
δv2n


,

(24)

where vnRP and δv2n are function parameters that are related but not equal to the mean and variance of the distribution,
210 respectively. Because flow is a vector quantity, it has both a magnitude and a phase. When measuring vn one measures
211 the modulus of the complex number, meaning there is a reduction in the number of dimensions from two to one. If the212 fluctuations in each dimension are Gaussian, one then expects the final distribution of values to be Bessel-Gaussian.
213 Recently the CMS experiment has carried out a similar flow unfolding and observes small deviations from the Bessel214 Gaussian form, favoring the elliptic power distribution [26].
215
For the unfolding, ATLAS determines the response matrix in a data driven way. The smearing in the response matrix
216 is modest as Pb+Pb collisions have a high multiplicity and the ATLAS detector has large phase space coverage for
217 tracks −2.5 < η < +2.5. In our case, the multiplicity of Au+Au collisions is lower in comparison with the multiplicity
218 in Pb+Pb collisions and the phase space coverage of the FVTX detector is significantly smaller. Hence, the smearing
219 as encoded in the response matrix is significantly larger and the unfolding is more challenging.
220
To estimate the response function of the detector, we follow the procedure from ATLAS [10], which is to examine
221 the difference between two subevents for both Qx and Qy . We compare the Q-vector determined in the south arm ofsouth
222 the FVTX, Q, to the Q-vector determined in the north arm of the FVTX, Qnorth . Figure 3 shows an example
223 of this procedure for the n = 2 case. Figure 3 (a) shows the 2-dimensional distribution for the 20%–30% centrality
224 selection; Fig 3 (b) shows the one-dimensional projection of this onto the x-axis (i.e. the one-dimensional distributionnorth
225 of Qx); Fig 3 (c) shows the one-dimensional projection of this onto the y-axis (i.e. the one-dimensional− Qsouthx
north
226 distribution of Qy). These distributions in all centrality selections are Gaussian over four orders of− Qsouthy
227 magnitude, and we characterize them via their Gaussian widths δ2SE which are given in Table II. It is notable that
228 these widths are more than a factor of two larger than those quoted by ATLAS in Pb+Pb collisions, for example
229 δ2SE = 0.050 for Pb+Pb 20%–25% central events [10].
230
If there is a modest longitudinal decorrelation between the two subevents, it will manifest as a slight increase in the
231 δ2SE parameter. The final vn is averaged over that decorrelation. This effect, as in previous unfolding analyses [10],232 is neglected.
233
We highlight that even in the case of a perfect detector with perfect acceptance, there remains a smearing due
234 to the finite particle number in each event. This raises a question regarding the meaning of a true vn that is being
235 unfolded back to. In a hydrodynamic description, there is a continuous fluid from which one can calculate a single
236 true anisotropy vn for each event. If one then has the fluid breakup into a finite number of particles N , e.g. via
237 Cooper-Frye freeze-out, the anisotropy of those N particles will fluctuate around the true fluid value. However, in
238 a parton scattering description, for example ampt [27], the time evolution is described in terms of a finite number
239 of particles N . In this sense there is no separating of a true vn from that encoded in the N particles themselves.
209Qy (North) - Qy (South)

91

PHENIX Au+Au sNN=200 GeV
n=2, 20-30% Central

0.5

0

700
600

−1
−1

3

10

400

102

102

10

10

100
−0.5

0
0.5
1
Qx (North) - Qx (South)

104

10

200

(a)

δ2SE,x = 0.1209

500

300
−0.5

104

0

1
−1

3

1

(b)
−0.5

0
0.5
1
Qx (North) - Qx (South)δ2SE,y = 0.1200

(c)

−1

−0.5

0
0.5
1
Qy (North) - Qy (South)√
FIG. 3. Example distribution of Qnorth − Qsouth for the n = 2 case corresponding to Au+Au collisions at sN N = 200 GeV
and centrality 20–30%. (a) The two-dimensional distribution. (b) The projection onto Qx . (c) The projection onto Qy . Shownfor (b) and (c) are the extracted Gaussian widths δ2SE —the χ2 /d.o.f. values of the fits are 1.02 and 1.18 for (b) and (c),respectively.

TABLE II. Resolution parameter δ2SE values for n = 2 and n = 3 for various centrality categories.
Centrality
0%–5%
5%–10%
10%–20%20%–30%30%–40%40%–50%240241δ2SE (n = 2)
0.117
0.115
0.115
0.121
0.133
0.154

δ2SE (n = 3)
0.115
0.113
0.113
0.119
0.130
0.151

Regardless, one can still mathematically apply the unfolding and compare experiment and theory as manipulated
through the same algorithm.
As noted before, the one-dimensional radial projection of a two-dimensional Gaussian is the so-called BesselGaussian function. In this case it means that the conditional probability to measure a value vnobs given a true value
vn has the following Bessel-Gaussian form: obs obs )2 +v 2
(vn
n
vn vn
p(vnobs |vn ) ∝ vnobs e− 2δ2
I0
,
(25)
δ2

where δ is the smearing parameter characterizing the response due to the finite particle number (including from thedetector efficiency and acceptance), and I0 is a modified Bessel function of the first kind. The smearing parameter
244 δ uses the combination of the two FVTX arms and is related to the result from the difference by δ = δ2SE /2. We
245 highlight that the Bessel-Gaussian in Eqn. 25 is different from the Bessel-Gaussian in Eqn. 24, though both arise from246 a similar dimensional reduction.
247
We have employed the iterative Bayes unfold method as encoded in RooUnfold [28] and our own implementation of
248 a Singular Value Decomposition (SVD) unfold method [29, 30]. In both cases, the response matrix is populated using
249 Eqn. 25 using the data-determined smearing parameters. The FVTX-determined event-by-event Qn distributions are
250 used as input to the unfold. Figure 4 shows this dimensional reduction for the n = 2 and n = 3 case, respectively.
251 Figure 4 (a) and (c) show the two-dimensional distribution of Qn and Figure 4 (b) and (d) show the one-dimensional
252 distribution of |Qn | for the 20%–30% centrality range.
253
The Au+Au 20%–30% centrality class is expected to provide the best conditions, in terms of the predicted hv2 i
254 and resolution δ, to determine p(v2 ) via unfolding. However, because it is quite challenging to unfold the measured
255 distribution directly, we constructed a test version of the problem to illustrate the procedure using the SVD method.
256 The details of this test are given in Appendix V, but the end result is that the unfolding procedure is inherently
257 unstable and therefore fails to converge for the resolution parameters in the present analysis.
258
As a result, instead of inverting the response matrix, we can make an ansatz that the probability distributions,
242

243

10

Qy

1

105

PHENIX Au+Au sNN=200 GeV
Data n=2, 20-30% Central
Bessel-Gaussian Truth
Bessel-Gaussian Truth +Forward Fold

0.9

0.4

0.8
0.7

0.2

0.6
0.5

0

0.4
0.3

−0.2

0.2
−0.4

(a)
−0.4

−0.2

0

0.2

0.4
Qx

0.1
0

Qy

1
0.9

0.4

0.8
0.7

0.2

0.6
0.5

0

0.4
0.3

−0.2

0.2

(c)

−0.4
−0.4

−0.2

0

0.2

0.4
Qx

0.1
0

104
103
102
10

(b)

1
0

106
5

10

104

0.2

0.4
Q = Q2x + Q2y

PHENIX Au+Au sNN=200 GeV
Data n=3, 20-30% Central
Bessel-Gaussian Truth
Bessel-Gaussian Truth +Forward Fold

103
102
10
1
0

(d)
0.1

0.2

0.3
Q = Q2x + Q2y

FIG. 4. Example distribution of Q for the (a), (b) elliptic n = 2 case and (c), (d) triangular n = 3 case (lower). (a) and (c) show
the 2-dimensional distribution. (b) and (d) show the 1-dimensional distribution |Q|. The distribution corresponds to Au+Au√
collisions at sN N = 200 GeV and centrality 20–30%. Also shown in (b), (d) is the best fit for this run Bessel-Gaussian truth
distribution and the corresponding forward folded result (i.e. pushing the truth distribution through the response matrix).

p(v2 ) and p(v3 ) are exactly Bessel-Gaussian in form. Under this restrictive assumption, because the Bessel-Gaussian
form has only two parameters as shown in Eqn. 24, we can simply evaluate a large grid of parameter combinations as
261 guesses for the truth distribution and forward fold them, i.e. passing them through the response matrix to compare
262 to the observed Qn distribution. We have carried out such a “forward fold” procedure with over 10,000 parameter
263 combinations. We then determine the statistical best fit parameters and their statistical uncertainties based on a2
264 χ
mapping. We consider only the Gaussian statistical uncertainties here, and detail our treatment of systematic2
265 uncertainties in Section IV B. Some advantages of this procedure are that we explore the full χspace and have no266 sensitivity to an unfolding prior, regularization scheme, and number of iterations. The disadvantage of course is the
267 ansatz that the distribution is precisely Bessel-Gaussian.
268
Examples of this χ2 forward fold mapping are shown in Figure 5. It is striking that for the v2 case in the Au+Au
269 20%–30% central bin, the forward folding reveals a tight constraint on the Bessel-Gaussian parameters. In contrast,
270 for the same centrality bin and v3 , there is a band of parameter combinations providing a roughly equally good match
271 to the experimental Q3 distribution. Shown in Figure 4 for v2 (upper) and v3 (lower) are the best Bessel-Gaussian
272 fit distributions (red) and their forward folded results (blue). Both cases show good agreement between the forward
2
273 folded results and the measured experimental distribution. The corresponding best χ
min values indicate a good match.
259
260

11
It is notable that in the v3 case, the χ2min values are slightly worse in all cases and significantly worse in the most
275 peripheral selection. The more peripheral data have a slightly larger tail at high Q3 values which could indicate an
276 incompatibility with the Bessel-Gaussian ansatz.
274

25

PHENIX Au+Au sNN=200 GeV 20-30%

0.025

20

T

EbE v2: χ2-χ2

min

10

0.01

5

(a)

0.03 0.035 0.04 0.045 0.05 0.055 0.06
vRP
2

0

20

h+/- 1<|η|<3 , p >0.3 GeV/c
T

0.015

EbE v3: χ2-χ2

min

15

0.015

0.005

25
PHENIX Au+Au sNN=200 GeV 20-30 %h+/- 1<|η|<3 , p >0.3 GeV/c

0.02

0.02

δ3

δ2

0.03

15

0.01

10

0.005

5

(b)
0.005

0.01

0.015

0.02
vRP
3

0

FIG. 5. Two dimensional color plot of χ2 − χ2min values as a function of the two Bessel-Gaussian parameters, vnRP and δvn forthe 20%–30% centrality bin. (a) shows the second harmonic and (b) shows the third harmonic. Only χ2 values up to χ2min + 25
are shown.

277

278
279

280

IV.

RESULTS AND DISCUSSION

Here we detail the full set of results for the elliptic and triangular flow moments and distributions in Au+Au√
collisions at sN N = 200 GeV. We start by detailing the cumulant results.

A.

Cumulants Results

First we show in Figure 6 (a) and (b) the centrality dependence of v2 {2} and v2 {4}, respectively. The statisticaluncertainties are shown as vertical lines and the systematic uncertainty is quoted as a global factor uncertainty. (a)283 shows a dramatic difference for centrality larger than 40% between the red points, obtained without requiring a
284 pseudorapidity gap in the particle pair, and the magenta points, which have a pseudorapidity gap of |∆η| > 2.0. This
285 is due to the fact that the pseudorapidity gap removes a large amount of nonflow, especially in the peripheral collisions
286 where nonflow is combinatorially less suppressed relative to central collisions. Contrariwise, (b) shows no difference
287 between the black points (no pseudorapidity gap) and two different 2-subevent cumulants, one where short-range pairs
288 are allowed (blue points) and one where they are not (red points). The absence of any effect here indicates that the289 4-particle correlation sufficiently suppresses nonflow combinatorially such that the kinematic separation of particles
290 provides no additional benefit.
Note that this is not necessarily the case in smaller collision systems—subevent
291 cumulants have been shown to significantly reduce nonflow in p+p/Pb collisions at the LHC [31], and are of potential
3
292 interest in p/d/ He+Au collisions at RHIC.
293
Figure 7 shows the centrality dependence of multi-particle v2 , with 2, 4, 6, and 8 particles. The 4-, 6-, and 8294 particles v2 values are consistent with each other, as expected from the small-variance limit [7]. When accounting
295 for the η-dependence of v2 as measured by PHOBOS [32], which indicates that v2 at 1 < |η| < 3 is about 1.25 times
296 lower than it is at |η| < 1, the 2-, 4-, and 6-particle cumulant v2 are in good agreement with the STAR results [12].
ppConsidering that vn {2} = vn2 + σv2n and that in the small variance limit vn {4} ≈ vn2 − σv2n [8], one can estimate
the relative fluctuations ass
σ vn(vn {2})2 − (vn {4})2
≈
.
(26)
hvn i
(vn {2})2 + (vn {4})2
281
282

0.1

0.09
0.08

v2{2}
v2{2,|∆η|>2}

Au+Au sNN = 200 GeV
h± 1<|η|<3
PHENIX

0.08

v2

v2

12

0.07
0.06

0.07
0.06

0.05

0.05

0.04

0.04

0.03

0.03

v2{4}
v2{4}
ab|ab
v2{4}
aa|bb

0.02

0.02

Sys. Uncert. 6%

0.01
0
0

Au+Au sNN = 200 GeV
h± 1<|η|<3
PHENIX

10

20

30

40

50

0.01

(a)
60

70

0
0

80 90 100Centrality (%)

Sys. Uncert. 6%
10

20

30

40

50

(b)
60

70

80 90 100Centrality (%)

v2

FIG. 6. Centrality dependence of (a) v2 {2} and (b) v2 {4}. (a) The red points indicate no pseudorapidity gap whereas themagenta points indicate a pseudorapidity gap of |∆η| > 2.0. (b) The black points indicate v2 {4} with no pseudorapidity gap,
the blue points indicate a two-subevent method with |∆η| > 2.0 but where some short-range pairs are allowed, and the redpoints indicate a two-subevent method with |∆η| > 2.0 where no short-range pairs are allowed.

0.08
0.07

v2{2,|∆η|>2}
v2{4}
v2{6}
v2{8}

Au+Au sNN = 200 GeV
h± 1<|η|<3
PHENIX

0.06
0.05
0.04
0.03
0.02
0.01
0
0

Sys. Uncert. 6%
10

20

30

40

50

60

70

80 90 100Centrality (%)√
FIG. 7. Multi-particle v2 as a function of centrality in Au+Au collisions at sN N = 200 GeV. The magenta open diamonds
indicate the v2 {SP }, the blue open squares indicate v2 {4}, the black open circles indicate v2 {6}, and the green filled diamonds
indicate v2 {8}.

Figure 8 shows the centrality dependence of this cumulant estimate of σv2 /hv2 i. Despite the difference in the rapidity
region where the data are measured, they are in good agreement with STAR [12] and PHOBOS [13]. Also shown is a
299 comparison with ampt analyzed via cumulants in the same way as the experimental data. There is good agreement
300 between the two, indicating that the Monte Carlo Glauber initial conditions in ampt and their fluctuations capture
301 the key event-by-event varying ingredients. We can also calculate the event-by-event variations in the initial conditions
302 directly via Monte Carlo Glauber. In this case we utilize the event-by-event spatial eccentricity εn distributions. If
303 there is a linear mapping between initial spatial eccentricity and final momentum anisotropy (εn ∝ vn ), we should
304 expect a good match between σε / hεn i and σv / hvn i. Also show in Fig. 8 is the Monte Carlo Glauber result viann
305 the calculation of cumulants (solid blue line), as well as the direct calculation of the variance and mean from the full
306 εn distribution (dashed blue line). One sees that in midcentral 10%–50% collisions, the data and both theory curves
307 agree reasonably. For more central collisions, the Monte Carlo Glauber data-style calculation shows the same trend
297
298

13
as the data whereas the Monte Carlo Glauber direct calculation is significantly lower. This is due to the fact thatthe small-variance limit is not a valid approximation in central collisions. In peripheral collisions, both Monte Carlo
310 Glauber curves under-predict the data. This has been attributed to the nonlinear response in hydrodynamics [33].
308309

σv2/ 〈v2〉1.2
PHENIX

Au+Au sNN = 200 GeV
h± 1<|η|<3
MC Glauber, cumulant based estimate
MC Glauber, direct calculation
AMPT
Data

1

0.8
0.6
0.4
0.2
0
0

Sys. Uncert. 2%
10

20

30

40

50

60

70

80 90 100Centrality (%)√
FIG. 8. Cumulant method estimate of σv2 /hv2 i as a function of centrality in Au+Au collisions at sN N = 200 GeV. The
data are shown as black open squares. The same calculation as done in data is done in ampt, shown as a solid green line.
Calculations of σε2 /hε2 i performed in the Monte Carlo Glauber model are shown as blue lines. The solid blue line is theMonte Carlo Glauber calculation done using the same estimate as the data, the dashed blue line is the direct calculation of the
moments of the MC Glauber ε2 distribution.

Now we consider the v3 case. Figure 9 shows v3 {2, |∆η| > 2} as a function of centrality. The centrality dependenceof v3 is much smaller than that of v2 , which is expected because triangular flow is generated dominantly through313 fluctuations.
311

v3

312

0.02

0.018
0.016

v3{2,|∆η|>2}

Au+Au sNN = 200 GeV
h± 1<|η|<3
PHENIX〈v23〉 from folding

0.014
0.012
0.01
0.008
0.006
0.004
0.002
0
0

10

20

30

40

50
60
70
Centrality (%)√
FIG. 9. Centrality dependence of v3 {2, |∆η| > 2} in Au+Au collisions at sN N = 200 GeV, shown as magenta diamonds. Thepsystematic uncertainty is indicated as a shaded magenta band. Also shown as black squares are hv32 i as determined from thefolding analysis, which is shown in the next section.
314

Figure 10 (a) shows the results for c3 {4} as a function of centrality. The results are always positive within the

14systematic uncertainties and shows a trend towards even larger positive values as one moves away from the mostcentral collisions. Since v3 {4} = (−c3 {4})1/4 the v3 {4} are complex-valued.
317
Recently positive valued c2 {4} has been observed in p+Au collisions at RHIC [34] and p+p collisions at the
318 LHC [31, 35] and has been interpreted as arising from short range nonflow contributions.
The use of subevents,
319 especially when requiring the particles in the cumulant to be separated in rapidity, significantly reduces nonflow
320 contributions and yields negative values of c2 {4} where standard cumulant analysis does not [25]. In the Au+Au
321 analysis presented here, due to the FVTX acceptance, one includes both short range particle combinations (some or
322 all particles in a single FVTX arm) and long-range combinations.
323
We explore the potential influence of such short range nonflow contributions as well as the opposite effect from
324 long-range decorrelations by changing the FVTX arm requirements of the particle combinations. The most extreme
325 is requiring all particles in a single arm, shown in Figure 10 (b), and the result is an even larger positive c3 {4}—i.e. in
326 the direction expected from increased short range nonflow and opposite to the expectation of long-range decorrelations
327 causing the positive c3 {4}. We can also consider combinations of two subevents, with two particles in each FVTX
328 arm. One case, labeled ab|ab, has some short range correlations though fewer than the standard, whereas the other329 case, labeled aa|bb, doesn’t allow any. One sees a consistent behavior emerge: c3 {4}aa|bb < c3 {4}ab|ab < c3 {4} <330 c3 {4}singlearm All of these results go in the direction of a large nonflow influence which may be exacerbated by the
331 very small v3 flow signal particularly, at forward rapidity.
√
332
The STAR experiment has also measured c3 {4} in Au+Au collisions at sN N = 200 GeV, though at midrapidity
333 |η| < 1.0 [36]. Their results, also shown in Figure 10 (a), are consistent with zero and fluctuate between positive and334 negative c3 {4} values. The difference between the STAR and PHENIX data points likely stems from the different
335 acceptance in pseudorapidity (the STAR points are measured over |η| < 1 while the PHENIX points are measured
336 over 1 < |η| < 3 as discussed above). Differences in nonflow, event plane decorrelations, and the relative contribution
337 from fluctuations as a function of pseudorapidity may all contribute to these observations.
√
These results seem to indicate that the small-variance limit is not applicable to v3 in Au+Au collisions at sN N
338
339 = 200 GeV for any centrality. Regardless, the measurement of these 2- and 4-particle cumulants is insufficient to340 constrain the mean and variance of the triangular flow event-by-event distribution.
315

316

−6

−6

c3{4}

c3{4}

0.15 ×10
Au+Au sNN = 200 GeV
h± 1<|η|<3
PHENIX
0.1

0.5 ×10
Au+Au sNN = 200 GeV
h± 1<|η|<3
0.4
PHENIX

Combined arms
Single arm

0.3
0.05
0.2
0
−0.05
(a)
−0.1
0

0.1

c3{4}
c3{4}
ab|ab
c3{4}
aa|bb
STAR |η|<1
10

20

30

40

0
(b)

50
60
Centrality (%)

−0.1
0

10

20

30

40

50
60
Centrality (%)√
FIG. 10. Centrality dependence of c3 {4} for Au+Au collisions at sN N = 200 GeV. (a) Calculations using both arms: c3 {4}(black circles), c3 {4}ab|ab (blue diamonds), c3 {4}aa|bb (red squares), and comparison to STAR [36] (black stars). (b) Comparisonof c3 {4} determined using both arms (open symbols) and a single arm (closed symbols). Note that the open black circles are
the same in (a) and (b).

341

B.

Folding Results

Now we turn to the results from the event-by-event forward fold. As detailed in Section III B, in the v2 case theBessel-Gaussian parameters are well-constrained apart from the most central events. In the v3 case, however, the344 Bessel-Gaussian parameters are not well-constrained for any centrality class. However, despite the broad range of
345 possible δv
and v3RP values, these correspond to a rather small range for the real mean hv3 i and root-mean-square or3
346 variance σvof the distributions. This means that despite the lack of constraint on the parameters, the first (v3 ) and3
347 second (σv ) moments of the distribution are nevertheless well-constrained.
3
342
343

15
We can quantify hvn i and σvn by varying the Bessel-Gaussian parameters within the one- and two- standard
deviation statistical constraints. In addition, we determine the systematic uncertainties on these quantities by varying
350 the z-vertex and analyzing loose and tight cuts (as described for the cumulants analysis). An additional systematic
351 uncertainty on the response matrix is estimated by splitting the data sample into two subsets, one with higher
352 extracted δ and one with lower, forward folding the two data sets separately, and then assessing the difference.
353
Figure 11 (a) shows the extracted first moment hv2 i, Fig. 11 (b) shows the extracted second moment σv2 , and
354 Fig. 11 (c) shows the relative fluctuations σv /hv2 i, each as determined from the folding method and as a function of2
355 centrality. Likewise, Fig. 12 (a) shows the extracted hv3 i, Fig. 12 (b) shows the extracted σv , and Fig. 12 (c) shows3
356 the relative fluctuations σv /hv3 i. The colored bands indicate the statistical uncertainties at the 68.27% confidence
32357 level (red) and the 95.45% confidence level (green) from the χ
analysis. The thin black lines indicate the systematic
358 uncertainties. Also shown in 11 as blue squares are results from the cumulant based calculation as discussed in the359 previous section.
The hv2 i values are in excellent agreement for all centralities, and the σv2 and σv2 /hv2 i are in
360 reasonable agreement for 10%–50% centrality, where the small-variance limit holds. Figure 9 shows a comparisonp
361 between the cumulant result v3 {2, |∆η| > 2|} and the folding analysis result
hv32 i (calculated from the results in
362 Fig. 12). These results are consistent within the systematic uncertainties.
348

0.03

(a)

0
0

(b)

0.02

0.04

0.02

σv2 /<v 2> = RMS / MEAN

0.06

σv2 = RMS

<v2> = MEAN

349

PHENIX Au+Au 200 GeV
Bessel-Gaussian p(v2)± 1 σ (stat)± 2 σ (stat)
Cumulant Method
Sys. Uncertainty

10

20

30

40
50
Centrality

0.01

0
0

(c)
0.6

0.4

0.2

10

20

30

40
50
Centrality

0
0

10

20

30

40
50
Centrality

FIG. 11. Folding results for (a) v2 , (b) σv2 , and (c) σv2 /hv2 i. The black lines above and below the points indicate the systematic
uncertainties. The red (green) boxes indicate the statistical uncertainties at the 68.27% (95.45%) confidence level. In the caseof hv2 i, the statistical uncertainties at the 68.27% confidence level are too small to be seen, and the uncertainties at the 95.45%
confidence level are visible but noticeably smaller than the marker size. Shown as blue squares are the same quantities as
determined using the cumulant based calculation—these points are slightly offset in the x-coordinate to improve visibility.

We highlight that the σv2 /hv2 i values agree well with those determined from the cumulant method as shown in
Figure 8, except in the most central and peripheral Au+Au events. The most central 0%–5% events are exactly
365 where the Monte Carlo Glauber results in Figure 8 indicate a breakdown in the small-variance approximation. This366 is a good validation of the forward folding procedure and another confirmation that the event-by-event elliptic flow√
367 fluctuations in Au+Au collisions at
sN N = 200 GeV are dominated by initial geometry fluctuations.
368
Intriguingly, whereas the values of σv2 /hv2 i vary significantly as a function of centrality, the values of σv3 /hv3 i
369 are almost precisely 0.52 independent of centrality. To understand this better, we need to consider a rather peculiar
370 feature of the Bessel-Gaussian Function. Figure 13 shows the σv /hvn i of the Bessel-Gaussian as a function of the rationRPRP371 δ/vn . For values of δ > vn , the observed σv /hvn i saturates at a value of about 0.52. Thus, any Bessel-Gaussiann
372 in the large variance limit will have a σv /hvn i of the same value.
n373This observation can, in fact, help shed light on the observed discrepancy between the CMS [11] and ATLAS [10]√374 data on σv /hv3 i. Figure 14 shows σv /hv2 i and σv /hv3 i as a function of centrality in Pb+Pb collisions at
sN N
3
2
3
375 = 2.76 TeV from CMS and ATLAS. The CMS results are obtained using the cumulant method assuming the small376 variance limit. In contrast the ATLAS results are obtained via an event-by-event unfolding and calculating the exact
377 mean and variance of the distribution.
378
The σv2 /hv2 i values are in very good agreement, which appears to validate the small variance approximation (as√379 was also validated in the Au+Au at
sN N = 200 GeV case in this analysis). In contrast, there is a large difference in
380 the σv /hv3 i between the different methods. The ATLAS σv /hv3 i values are all very close to 0.52, exactly as observed
3
3
363
364

(a)

0.01

0.015

0.01

σv3 /<v 3> = RMS / MEAN

0.02

σv3 = RMS

<v3> = MEAN

16

(b)

0.005
PHENIX Au+Au 200 GeV

(c)

0.4

0.2

Bessel-Gaussian p(v3)

0.005

0.6

± 1 σ (stat)± 2 σ (stat)00

Sys. Uncertainty

10

20

30

0
0

40
50
Centrality

10

20

30

0
0

40
50
Centrality

10

20

30

40
50
Centrality

σvn/v n = RMS / MEAN

FIG. 12. Folding results for (a) hv3 i, (b) σv3 , and (c) σv3 /hv3 i. The black lines above and below the points indicate thesystematic uncertainties. The red (green) boxes indicate the statistical uncertainties at the 68.27% (95.45%) confidence level.
The σv3 /hv3 i values are all ≈ 0.52, the apparent limiting value of this quantity for the Bessel-Gaussian distribution.

0.6

0.4

0.2

0
0

0.5

1

1.5

2

2.5

3
δvn/v RPnFIG. 13. The observed ratio of the mean to the standard deviation (i.e. σvn /hvn i) of the Bessel-Gaussian as a function of theratio of the two free parameters δvn /vnRP . For values of δvn > vnRP , the observed σvn /hvn i saturates at 0.52.

above in the present Au+Au data and as found to be a limiting case for the Bessel-Gaussian function. To better
understand the σv3 /hv3 i, we also show σε3 /hε3 i as determined from MC Glauber calculations. The dashed red-line uses
383 the small-variance limit estimate with cumulants, as is done for the CMS data, and the agreement is quite reasonable.
384 The solid red line is calculated from the moments of the ε3 distribution directly, and shows good agreement with the
385 ATLAS data. This represents a quantitative confirmation of the event-by-event fluctuations and the breakdown in
386 the small variance approximation. The v3 {4} at forward rapidity at RHIC is found to be complex-valued, which may
387 be the result of a very small flow v3 and significant nonflow contributions.

381

382

388

389
390

V.

SUMMARY AND CONCLUSIONSIn summary, we have presented measurements of elliptic and triangular flow in Au+Au collisions at 200 GeV for
charged hadrons at forward rapidity 1 < |η| < 3. In particular, we compare flow cumulants (v2 {2}, v2 {4}, v2 {6},17σvn/ 〈vn〉

1.2
1

0.8

n=2 Cumulants (CMS)
n=2 E-by-E (ATLAS)
n=3 Cumulants (CMS)
n=3 E-by-E (ATLAS)
Glauber σ/< ε3> Cumulants
Glauber σ/< ε3> E-by-E

0.6
0.4
0.2
0
0

10 20 30 40 50 60 70 80 90 100
Centrality (%)

FIG. 14. The observed ratio of the standard deviation to the mean (i.e. σvn /hvn i) for n = 2 and n = 3 as a function of√centrality in Pb+Pb collisions at sN N = 2.76 TeV [10, 11]. Also shown are σε3 /ε3 values from MC Glauber estimated using
the small-variance limit cumulant technique (shown as the dashed line) and the direct calculation from the moments (shownas the solid line).

v2 {8} and v3 {2}, v3 {4}) and the mean and variance of the v2 and v3 event-by-event distributions using a forward-foldprocedure with a Bessel-Gaussian ansatz. These measurements are complementary in terms of sensitivity to initial
393 state geometry fluctuations and additional fluctuations from the evolution of the medium, for example via dissipative
394 hydrodynamics.
395
In the small-variance limit, where the event-by-event flow fluctuations are small compared to the average flow value
396 i.e. σv /hvn i < 1, we expect the cumulants extraction and the forward-fold results to agree. This is the case forn
397 elliptic flow in Au+Au collisions from 10%–50% central and both results agree with event-by-event fluctuations in the
398 initial geometry as calculated via Monte Carlo Glauber.
399
In contrast, we find that the small-variance limit fails for triangular flow for all centralities at RHIC and the LHC.
400 For LHC Pb+Pb results, the large-variance result for the cumulants can be described purely via Monte Carlo Glauber
401 initial geometry fluctuations. However, for RHIC Au+Au collisions the complex values of v3 {4} indicate that there
402 may be additional nonflow influences as well as sources of fluctuations in the translation of initial geometry into
403 final state momentum triangular anisotropies. Detailed comparisons with event-by-event hydrodynamic calculations
404 should be elucidating to understand the nature of these fluctuations.
391
392

405

ACKNOWLEDGMENTS

We thank the staff of the Collider-Accelerator and Physics Departments at Brookhaven National Laboratory andthe staff of the other PHENIX participating institutions for their vital contributions. We acknowledge support from
408 the Office of Nuclear Physics in the Office of Science of the Department of Energy, the National Science Foundation,
409 Abilene Christian University Research Council, Research Foundation of SUNY, and Dean of the College of Arts and
410 Sciences, Vanderbilt University (U.S.A), Ministry of Education, Culture, Sports, Science, and Technology and the
411 Japan Society for the Promotion of Science (Japan), Conselho Nacional de Desenvolvimento Cientı́fico e Tecnológico
412 and Fundação de Amparo à Pesquisa do Estado de São Paulo (Brazil), Natural Science Foundation of China (People’s
413 Republic of China), Croatian Science Foundation and Ministry of Science and Education (Croatia), Ministry of Edu414 cation, Youth and Sports (Czech Republic), Centre National de la Recherche Scientifique, Commissariat à l’Énergie
415 Atomique, and Institut National de Physique Nucléaire et de Physique des Particules (France), Bundesministerium
416 für Bildung und Forschung, Deutscher Akademischer Austausch Dienst, and Alexander von Humboldt Stiftung (Ger406
40718many), J. Bolyai Research Scholarship, EFOP, the New National Excellence Program (ÚNKP), NKFIH, and OTKA
(Hungary), Department of Atomic Energy and Department of Science and Technology (India), Israel Science Founda419 tion (Israel), Basic Science Research and SRC(CENuM) Programs through NRF funded by the Ministry of Education
420 and the Ministry of Science and ICT (Korea), Physics Department, Lahore University of Management Sciences (Pak421 istan), Ministry of Education and Science, Russian Academy of Sciences, Federal Agency of Atomic Energy (Russia),
422 VR and Wallenberg Foundation (Sweden), the U.S. Civilian Research and Development Foundation for the Indepen423 dent States of the Former Soviet Union, the Hungarian American Enterprise Scholarship Fund, the US-Hungarian
424 Fulbright Foundation, and the US-Israel Binational Science Foundation.
417
418

425

APPENDIX: TEST CASE FOR FULL UNFOLD

For this test case, the response matrix Â is shown in Fig. 15 (a) and is identical to that for the real data 20%–30%~ true = Q~ obs , where Qobs has been obtained in thecentrality class. We then attempt to solve the inverse problem ÂQ2
2
2
428 limit of infinite statistical precision, assuming a truth-level distribution with parameters such that the smearing, asT
429 encoded in Â, yields a distribution similar to that measured in data. The singular value factorization Â = Û Σ̂Ŵ
430 of the matrix is obtained, where Û and Ŵ are unitary matrices whose column vectors, ui and wi , are the left- and431 right- singular vectors of Â, respectively, and Σ̂ is a diagonal matrix, whose nonzero entries σi are its singular values.
432 Figure 15 (b) shows a few selected right-singular vectors wi . Notice that some vectors, namely those corresponding to
433 the largest singular values, are harmonic, whereas those corresponding to the smallest singular values are essentially
434 noise.
Because the response matrix is singular, we use the SVD decomposition to construct the solution of the inverse
problem as a linear combination of all right-singular vectors, as follows:!
Dim(A)
X~ obs~uTi · Q
2
~
w
~ i.
(27)
ϕi
Q2 =
σi
i=1
426

427

The damping factors ϕi = σi2 /(σi2 + λ2 ), for some λ ∈ R, are introduced to attenuate the contribution of the noisy
436 singular vectors to the sum. It is important to point out that in most implementations of SVD used in high-energy
437 physics, including RooUnfold, the above sum is simply truncated to include only a subset of the harmonic singular
438 vectors, potentially leading to loss of information.
439
To determine which singular vectors contribute to the solution in a meaningful manner, it is useful to examine the440 Picard plot [37] for the problem at hand, shown in Fig. 15 (c), which displays the singular values σi of Â, as well as~obs onto the singular vectors ~uT · Q~ obs , and the solution coefficients ~uT · Q~ obs /σi . Notice that the
441 the projection of Q
2
2
i
i
2
442 singular values and the Fourier coefficients drop sharply many orders of magnitude before leveling off, yet in such a
443 way that their ratio is roughly constant. The implication is then that all singular vectors appear to contribute equally
444 to the solution, which is clearly problematic given the noisy nature of most of them. In general, it is desirable for
445 Fourier coefficients to drop off faster than the singular values (to fulfill the so-called discrete Picard condition), such
446 that the Picard plot will reveal the appropriate set of terms to include in the solution, as identified by a sharp drop
447 in the solution coefficients.
448
Given that our problem does not satisfy the Picard condition, we introduce the attenuation factors ϕi in Eqn. 27.
true
449 The resulting unfolded Q2 is shown in Fig. 15 (d), along with the true Q2, and smeared Qobs
2 . We observe that the
450 unfolding works well, yielding a good description of the true distribution shape, with uncertainties associated with
451 varying the regularization parameter λ.
452However, in this case the unfolding procedure constitutes an ill-posed inverse problem, such that small perturbations
obs
453 in the input vector—that is, Q2—translate to very large errors in the solution, compounded by the fact that the
454 Picard condition is violated. In particular, we have verified with our test problem that the statistical fluctuations
obs
455 in Q2
when sampling a finite number of events, comparable to those recorded in data, indeed limit the number of
456 available harmonic singular vectors, thus causing the solution to be dominated by noise.
457
We now examine the application of the above unfolding method to data. Fig. 16 (a) shows an ansatz for Qtrue
2458 assuming a Bessel-Gaussian form, and the corresponding refolded smeared distribution. It compares very well to the459 data, as shown in the ratio plot in Fig. 16 (b). In principle, given the good quality of the fit, one would expect the
460 unfolding procedure to work with the data as input. However, the statistical fluctuations apparent in the ratio plot
461 perturb the solution in such a way that the noisy nonharmonic singular vectors are enhanced even more than in the
462 test problem, as shown in Fig. 16 (c). As a result, the number of available harmonic singular vectors is reduced,
463 and the problem has no satisfactory solution, even when regularization is applied. Thus, to be explicit, the unfolding
464 procedure fails. We note that if we apply our test example with a significantly better resolution, i.e. as in the ATLAS
465 Pb+Pb case, the method does converge as expected.
435

19
0

0.035
0.03

0.05

vobs
2

0.1
0.15
0.2

(a)
0

10

10
10

0

10
10

0.05

0.1

0.15

0.2

v2

0.25

w 13

0.1

0.015

0

w 20

-0.1
-0.2

(b)

-0.3
-0.4

0

4

10

5

0.02

0

w4

0.3
0.2

0.005

w1

0.4

0.025

0.01

0.25

0.5

20

40

60

80

Vector Entry

105

Unfolded
Truth
Smeared

3.5
3
2.5

i

uTi Q obs

2

uTi Q obs / i

-5

1.5

-10

1
0.5

10-15
10-20

100

(c)
0

0
20

40

60

Singular Value Index

80

100

-0.5

(d)
0.05

0.1

Q=

0.15

0.2

0.25

Q2x+Q 2y

FIG. 15. (a) Response matrix for test unfolding problem; (b) Selected right-singular values of the response matrix; (c) Picard~ true~ obsplot for inverse problem ÂQ=Q
2
2 , see text for details; (d) True, smeared, and unfolded Q2 as determined using SVD withattenuation factors.
Data/Refold

1.03

(a)

10 1

/ndf = 64.35 / 73 = 0.88

2

(b)

1.02

10 2

1.01
10

3

10

Au+Au 20-30% n=2

6

10 0

0.98

p(v ) truth
p(v
0.05

2
obs
2

0.1

) smear
0.15

0.2

0.25

0.3

0.35

0.4
Q2

0.970

10

20

10

0.99

PHENIX Data

5

30

(c)

i

uT Q obs
iuTi

Q obs/

i

1010

1

10 4

10

0.05

0.1 0.15 0.2 0.25

0.3 0.35

0.4
Q2

0

10

-10

10

-20

0

20

40

60

Singular Value Index

80

100

√
FIG. 16. (a) Data distribution for Qmeas
in Au+Au collisions at sN N = 200 GeV in the 20%–30% centrality class. Also
2
shown is an assumed Bessel-Gaussian truth distribution and its resolution smeared result. (b) Data divided by the resolution
smeared solution showing a good agreement within statistical uncertainties. (c) Picard plot for inverse problem with data
~ true
~ obs
ÂQ
=Q
2
2 , see text for details;

20[1] K. Adcox et al. (PHENIX Collaboration), “Formation of dense partonic matter in relativistic nucleus-nucleus collisions at
RHIC: Experimental evaluation by the PHENIX collaboration,” Nucl. Phys. A 757, 184 (2005).
468[2] J. Adams et al. (STAR Collaboration), “Experimental and theoretical challenges in the search for the quark gluon plasma:
469
The STAR Collaboration’s critical assessment of the evidence from RHIC collisions,” Nucl. Phys. A 757, 102 (2005).
470[3] B. B. Back et al. (PHOBOS Collaboration), “The PHOBOS perspective on discoveries at RHIC,” Nucl. Phys. A 757, 28471(2005).
472[4] I. Arsene et al. (BRAHMS Collaboration), “Quark gluon plasma and color glass condensate at RHIC? The Perspective473
from the BRAHMS experiment,” Nucl. Phys. A 757, 1 (2005).
474[5] P. Romatschke and U. Romatschke, “Relativistic Fluid Dynamics Out of Equilibrium,” ArXiv:1712.05815.
475[6] U. Heinz and R. Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci.
476
63, 123 (2013).
477[7] S. A. Voloshin, A. M. Poskanzer, A. Tang, and G. Wang, “Elliptic flow in the Gaussian model of eccentricity fluctuations,”
478
Phys. Lett. B 659, 537 (2008).
479[8] J.-Y. Ollitrault, A. M. Poskanzer, and S. A. Voloshin, “Effect of flow fluctuations and nonflow on elliptic flow methods,”
480
Phys. Rev. C 80, 014904 (2009).
481[9] J. Jia, “Event-shape fluctuations and flow correlations in ultra-relativistic heavy-ion collisions,” J. Phys. G 41, 124003
482(2014).
483 [10] G. Aad et al. (ATLAS Collaboration), “Measurement of the distributions of event-by-event flow harmonics in lead-lead484
collisions at = 2.76 TeV with the ATLAS detector at the LHC,” J. High Energy Phys. 11 (2013) 183.
485 [11] S. Chatrchyan et al. (CMS Collaboration), “Measurement of higher-order harmonic azimuthal anisotropy in PbPb collisions√486
at sN N = 2.76 TeV,” Phys. Rev. C 89, 044906 (2014).
√
487 [12] J. Adams et al. (STAR Collaboration), “Azimuthal anisotropy in Au+Au collisions at
sN N = 200 GeV,” Phys. Rev. C
488
72, 014904 (2005).
489 [13] B. Alver et al. (PHOBOS Collaboration), “Non-flow correlations and elliptic flow fluctuations in gold-gold collisions at√
490
sN N = 200 GeV,” Phys. Rev. C 81, 034915 (2010).
491 [14] K. Adcox et al. (PHENIX Collaboration), “PHENIX detector overview,” Nucl. Instrum. Methods Phys. Res., Sec. A 499,492
469 (2003).
493 [15] M. Allen et al. (PHENIX Collaboration), “PHENIX inner detectors,” Nucl. Instrum. Methods Phys. Res., Sec. A 499, 549
494(2003).
495 [16] C. Loizides, J. L. Nagle, and P. Steinberg, “Improved version of the PHOBOS Glauber Monte Carlo,” SoftwareX 1-2, 13496(2015).
497 [17] C. Aidala et al., “The PHENIX Forward Silicon Vertex Detector,” Nucl. Instrum. Methods Phys. Res., Sec. A 755, 44498(2014).
499 [18] S. Agostinelli et al. (GEANT4 Collaboration), “GEANT4: A Simulation toolkit,” Nucl. Instrum. Methods Phys. Res., Sec.
500
A 506, 250 (2003).
501 [19] S. A. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by Fourier expansion of Azimuthal particle
502
distributions,” Z. Phys. C 70, 665 (1996).
503 [20] A. Bilandzic, R. Snellings, and S. Voloshin, “Flow analysis with cumulants: Direct calculations,” Phys. Rev. C 83, 044913
504(2011).
505 [21] A. M. Poskanzer and S. A. Voloshin, “Methods for analyzing anisotropic flow in relativistic nuclear collisions,” Phys. Rev.
506
C 58, 1671 (1998).
507 [22] N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, “A New method for measuring azimuthal distributions in nucleus-nucleus
508
collisions,” Phys. Rev. C 63, 054906 (2001).
509 [23] A. Bilandzic, C. H. Christensen, K. Gulbrandsen, A. Hansen,
and Y. Zhou, “Generic framework for anisotropic flow
510
analyses with multiparticle azimuthal correlations,” Phys. Rev. C 89, 064904 (2014).
511 [24] C. Adler et al. (STAR Collaboration), “Elliptic flow from two and four particle correlations in Au+Au collisions at√
512
sN N = 200 GeV,” Phys. Rev. C 66, 034904 (2002).
513 [25] J. Jia, M. Zhou, and A. Trzupek, “Revealing long-range multiparticle collectivity in small collision systems via subevent514
cumulants,” Phys. Rev. C 96, 034906 (2017).
√
515 [26] A. M Sirunyan et al. (CMS Collaboration), “Non-Gaussian elliptic-flow fluctuations in PbPb collisions at
sNN = 5.02
516
TeV,” ArXiv:1711.05594.
517 [27] Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, “A Multi-phase transport model for relativistic heavy ion collisions,”
518
Phys. Rev. C 72, 064901 (2005).
519 [28] T. Adye, “Unfolding algorithms and tests using RooUnfold,” in Proceedings, PHYSTAT 2011 Workshop on Statistical
520
Issues Related to Discovery Claims in Search Experiments and Unfolding, CERN,Geneva, Switzerland 17-20 January521
2011 , CERN (CERN, Geneva, 2011) p. 313, arXiv:1105.1160 [physics.data-an].
522 [29] P. C. Hansen, “Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numer523
ical rank,” SIAM J. Sci. Stat. Comput. 11, 503 (1990).
524 [30] A. Hocker and V. Kartvelishvili, “SVD approach to data unfolding,” Nucl. Instrum. Methods Phys. Res., Sec. A 372, 469
525(1996).
526 [31] M. Aaboud et al. (ATLAS Collaboration), “Measurement of long-range multiparticle azimuthal correlations with the
466

467

21
subevent cumulant method in pp and p+Pb collisions with the ATLAS detector at the CERN Large Hadron Collider,”Phys. Rev. C 97, 024904 (2018).
529 [32] B. B. Back et al. (PHOBOS Collaboration), “Centrality and pseudorapidity dependence of elliptic flow for charged hadrons√
530
in Au+Au collisions at sN N = 200 GeV,” Phys. Rev. C 72, 051901 (2005).
531 [33] J. Noronha-Hostler, L. Yan, F. G. Gardim, and J.-Y. Ollitrault, “Linear and cubic response to the initial eccentricity in
532
heavy-ion collisions,” Phys. Rev. C 93, 014909 (2016).
533 [34] C. Aidala et al. (PHENIX Collaboration), “Measurements of Multiparticle Correlations in d + Au Collisions at 200, 62.4,
534
39, and 19.6 GeV and p + Au Collisions at 200 GeV and Implications for Collective Behavior,” Phys. Rev. Lett. 120,535
062302 (2018).
536 [35] V. Khachatryan et al. (CMS Collaboration), “Evidence for collectivity in pp collisions at the LHC,” Phys. Lett. B 765,537
193 (2017).
538 [36] L. Adamczyk et al. (STAR Collaboration), “Third Harmonic Flow of Charged Particles in Au+Au Collisions at sqrtsNN
539= 200 GeV,” Phys. Rev. C 88, 014904 (2013).
540 [37] P. C. Hansen, “The discrete picard condition for discrete ill-posed problems,” BIT Numerical Mathematics 30, 658 (1990).

527
528