The azimuth structure of nuclear collisions -- I
Version 1.6
The azimuth structure of nuclear collisions – I
Thomas A. Trainor and David T. Kettler
CENPA 354290, University of Washington, Seattle, WA 98195
(Dated: November 5, 2018)
We describe azimuth structure commonly associated with elliptic and directed flow in the context
of 2D angular autocorrelations for the purpose of precise separation of so-called nonflow (mainly
minijets) from flow. We extend the Fourier-transform description of azimuth structure to include
power spectra and autocorrelations related by the Wiener-Khintchine theorem. We analyze several
examples of conventional flow analysis in that context and question the relevance of reaction plane
estimation to flow analysis. We introduce the 2D angular autocorrelation with examples from
data analysis and describe a simulation exercise which demonstrates precise separation of flow and
nonflow using the 2D autocorrelation method. We show that an alternative correlation measure
based on Pearson’s normalized covariance provides a more intuitive measure of azimuth structure.
PACS numbers: 13.66.Bc, 13.87.-a, 13.87.Fh, 12.38.Qk, 25.40.Ep, 25.75.-q, 25.75.Gz
I. INTRODUCTION
A major goal of the RHIC is production of color-
deconfined or QCD matter in heavy ion (HI) collisions, a
bulk QCD medium extending over a nontrivial space-
time volume which is in some sense thermalized and
whose dynamics are dominated in some sense by quarks
and gluons as the dominant degrees of freedom [1]. “Mat-
ter” in this context means an aggregate of constituents
in an equilibrium state, at least locally in space-time,
such that thermodynamic state variables provide a nearly
complete description of the system. Demonstration of
thermalization is seen by many as a necessary part of the
observation of QCD matter.
A. Global variables
One method proposed to demonstrate the existence of
QCD matter is to measure trends of global event vari-
ables, statistical measures formulated by analogy with
macroscopic thermodynamic quantities and based on in-
tegrals of particle yields over kinematically accessible mo-
mentum space. E.g., temperature analogs include spec-
trum inverse slope parameter T and ensemble-mean pt
p̂t. Chemical analogs include particle yields and their
ratios, such as the ensemble-mean K/π ratio [2]. Corre-
sponding fluctuation measures have been formulated for
the event-wise mean pt 〈pt〉 (“temperature” fluctuations)
and K/π ratio (chemical or flavor fluctuations) [3, 4, 5].
Arguments by analogy are less appropriate when deal-
ing with small systems (‘small’ in particle number, space
and/or time) where large deviations from macroscopic
thermodynamics may be encountered.
B. Flow analysis
One such global feature is the large-scale angular struc-
ture of the event-wise particle distribution. The compo-
nents of angular structure described by low-order spheri-
cal or cylindrical harmonics are conventionally described
as “flows.” The basic assumption is that such structure
represents collective motion of a thermalized medium,
and hydrodynamics is therefore an appropriate descrip-
tion. Observation of larger flow amplitudes is therefore
interpreted by many to provide direct evidence for event-
wise thermalization in heavy ion collisions [6]. Given
those assumptions each collision event is treated sepa-
rately. Event-wise angular distributions are fitted with
model functions associated with collective dynamics. The
model parameters are interpreted physically in a thermo-
dynamic (i.e., collective, thermalized) context.
However, collective flow in many-body physics is a
complex topic with longstanding open issues. There is
conflict in the description of nuclear collisions between
continuum hydrodynamics and discrete multiparticle sys-
tems which echoes the state of physics prior to the study
of Brownian motion by Einstein and Perrin [7, 8]. Be-
yond the classical dichotomy between discrete and con-
tinuous dynamics there is the still-uncertain contribu-
tion of quantum mechanics to the early stages of nuclear
collisions. Quantum transitions may play a major role
in phenomena perceived to be “collective.” Premature
imposition of hydrodynamic (hydro) models on collision
data may hinder full understanding.
C. Nonflow and multiplicity distortions
A major concern for conventional flow analysis is the
presence of “nonflow,” non-sinusoidal contributions to
azimuth structure often comparable in amplitude to sinu-
soid amplitudes (flows). Nonflow is treated as a system-
atic error in flow analysis, reduced to varying degrees by
analysis strategies. Another significant systematic issue
is “multiplicity distortions” associated with small event
multiplicities, also minimized to some degree by analysis
strategies. Despite corrections nonflow and small multi-
plicities remain a major limitation to conventional flow
measurements. For those reasons flow measurements in
peripheral heavy ion collisions are typically omitted. The
http://arxiv.org/abs/0704.1674v1
opportunity is then lost to connect the pQCD physics of
elementary collisions to nonperturbative, possibly collec-
tive dynamics in heavy ion collisions.
D. Minijets
A series of recent experiments has demonstrated that
the nonsinusoidal components of angular correlations at
full RHIC energy are dominated by fragments from low-
Q2 partons or minijets [9, 10, 11, 12]. Minijets in RHIC
p-p and A-A collisions have been studied extensively via
fluctuations [9, 10] and two-particle correlations [11, 12].
Minijets may dominate the production mechanism for the
QCD medium [13], and may also provide the best probe
of medium properties, including the extent of thermal-
ization. Comparison of minijets in elementary and heavy
ion collisions may relate medium properties and collec-
tive motion to a theoretical QCD context.
Ironically, demonstrating the existence and properties
of collective flows and of jets (collective hadron motion
from parton collisions and fragmentation) is formally
equivalent. Identifying a partonic “reaction plane” and
determining a nucleus-nucleus reaction plane require sim-
ilar techniques. For example, sphericity has been used to
obtain evidence for deformation of particle/pt/Et angu-
lar distributions due to parton collisions [14] and collec-
tive nucleon flow [15]. At RHIC we should ask whether
final-state angular correlations depend on the geometry
of parton collisions (minijets), N-N collisions or nucleus-
nucleus collisions (flows), or all three. The analogy is
important because to sustain a flow interpretation one
has to prove that there is a difference: e.g., to what ex-
tent do parton collisions contribute to flow correlations
or mimic them? The phenomena coexist on a continuum.
To resolve such ambiguities we require analysis meth-
ods which treat flow and minijets on an equal footing
and facilitate their comparison, methods which do not
impose the hypothesis to be tested on the measurement
scheme. To that end we should: 1) develop a consistent
set of neutral symbols; 2) manipulate random variables
with minimal approximations; 3) introduce proper statis-
tical references so that nonstatistical correlations of any
origin can be isolated unambiguously; 4) treat azimuth
structure ab initio in a model-independent manner using
standard mathematical methods (e.g., standard Fourier
analysis); and 5) include what is known about minijets
(“nonflow”) and “flow” in a more general analysis based
on two-particle correlations.
E. Structure of this paper
An underlying theme of this paper is the formal re-
lation between event-wise azimuth structure in nuclear
collisions and Brownian motion, and how that relation
can inform our study of heavy ion collisions. We begin
with a review of Fourier transform theory and the rela-
tion between power spectra and autocorrelations. That
material forms a basis for analysis of sinusoidal compo-
nents of angular correlations in nuclear collisions which
is well-established in standard mathematics.
We then review the conventional methods of flow anal-
ysis from Bevalac to RHIC. Five papers are discussed in
the context of Fourier transforms, power spectra and au-
tocorrelations. To facilitate a more general description of
angular asymmetries we set aside flow terminology (ex-
cept as required to make connections with the existing lit-
erature) and move to a model-independent description in
terms of spherical and cylindrical multipole moments. We
emphasize the relation of “flows” to multipole moments
as model-independent correlation measures. Physical in-
terpretation of multipole moments is an open question.
We then consider whether event-wise estimation of the
reaction plane is necessary for “flow” studies. The con-
ventional method of flow analysis is based on such esti-
mation, assuming that event-wise statistics are required
to demonstrate collectivity, and hence thermalization, in
heavy ion collisions.
We define the 2D joint angular autocorrelation and de-
scribe its properties. The autocorrelation is fundamental
to time-series analysis, the Brownian motion problem and
its generalizations and astrophysics, among many other
fields. It is shown to be a powerful tool for separating
“flow” from “nonflow.” The autocorrelation eliminates
biases in conventional flow analysis stemming from finite
multiplicities, and makes possible bias-free study of cen-
trality variations in A-A collisions down to N-N collisions.
“Nonflow” is dominated by minijets (minimum-bias
parton fragments, mainly from low-Q2 partons) which
can be regarded as Brownian probe particles for the QCD
medium, offering the possibility to explore small-scale
medium properties. Minijet systematics provide strong
constraints on “nonflow” in the conventional flow con-
text. Interaction of minijets with the medium, and par-
ticularly its collective motion, is the subject of paper II
of this series.
Finally, we consider examples from RHIC data of au-
tocorrelation structure. We show the relation between
“flow” and minijets, how conventional flow analysis is
biased by the presence of minijets, and how the autocor-
relation method eliminates that bias and insures accurate
separation of different collision dynamics.
We include several appendices. In App. A we review
Brownian motion and its formal connection to azimuth
correlations in nuclear collisions. In App. B we review
the algebra of random variables in relation to conven-
tional flow analysis techniques. We make no approxima-
tions in statistical analysis and invoke proper correlation
references to obtain a minimally-biased, self-consistent
analysis system in which flow and nonflow are precisely
distinguished. In App. C we review the mathematics
of spherical and cylindrical multipoles and sphericity. In
App. D we review subevents, scalar products and event-
plane resolution. In App. E we summarize some A-A
centrality issues related to azimuth multipoles and mini-
jets.
II. FOURIER ANALYSIS
The azimuth structure of nuclear collisions is part of
a larger problem: angular correlations of number, pt and
Et on angular subspace (η1, η2, φ1, φ2). There is a for-
mal similarity between event-wise particle distributions
on angle and the time series of displacements of a parti-
cle in Brownian motion. In either case the distribution
is discrete, combining a large random component with
the possibility of a smaller deterministic component. The
mathematical description of Brownian motion includes as
a key element the autocorrelation density, related to the
Fourier power spectrum through the Wiener-Khintchine
theorem (cf. App. A).
The Fourier series describes arbitrary distributions on
bounded angular interval 2π or distributions periodic on
an unbounded interval. The azimuth particle distribu-
tion from a nuclear collision is drawn from (samples)
a combination of sinusoids nearly invariant on rapidity
near midrapidity, conventionally described as “flows,”
and other azimuth structure localized on rapidity and
conventionally described as “nonflow.” The two contri-
butions are typically comparable in amplitude.
We first consider the mathematics of the Fourier trans-
form and power spectrum and their role in conventional
flow analysis [16]. We assume for simplicity that the
only angular structure in the data is represented by a few
lowest-order Fourier terms. In conventional flow analysis
the azimuth distribution is described solely by a Fourier
series, and corrections are applied in an attempt to com-
pensate for “nonflow” as a systematic error. We later re-
turn to the more general angular correlation problem and
consider non-sinusoidal (nonflow) structure described by
non-Fourier model functions in the larger context of 2D
(joint) angular autocorrelations. Precise description of
the composite structure requires a hybrid mathematical
model.
Event-wise random variables are denoted by a tilde.
Variables without tildes are ensemble averages, indicated
in some cases explicitly by overlines. Event-wise averages
are indicated by angle brackets. The algebra of random
variables is discussed in App. B. Where possible we em-
ploy notation consistent with conventional flow analysis.
A. Azimuth densities
The event-wise azimuth density (particle, pt or Et) is
a set of n samples from a parent density integrated over
some (pt, η) acceptance, a sum over Dirac delta functions
(particle positions)
ρ̃(φ) =
ri δ(φ− φi) (1)
The ri are weights (1, pt or Et) appropriate to a given
physical context. We assume integration over one unit of
pseudorapidity, so multiplicity n estimates dn/dη (simi-
larly for pt and Et). The continuum parent density, not
directly observable, is the object of analysis. Fixed parts
of the parent density are estimated by a histogram aver-
aged over an event ensemble. The discrete nature of the
sample distribution and its statistical character present
analysis challenges which are one theme of this paper.
The correlation structure of the single-particle azimuth
density is manifested in ensemble-averaged multiparticle
(two-particle, etc.) densities. Accessing that structure
by projection of multiparticle spaces to 2D or 1D with
minimal distortion is the object of correlation analysis.
In each event the two-particle density is the Cartesian
product ρ̃(φ1, φ2) = ρ̃(φ1) ρ̃(φ2)
ρ̃(φ1, φ2) =
r2i δ(φ1 − φi)δ(φ2 − φi) (2)
n,n−1
rirj δ(φ1 − φi)δ(φ2 − φj),
where the first term represents self pairs. Ensemble-
averaged two-particle distribution ρ(φ1, φ2) with correla-
tions is not generally factorizable. By comparing the av-
eraged two-particle distribution to a factorized (or mixed-
pair) statistical reference two-particle correlations are re-
vealed. In a later section we compare multipole moments
from two-particle correlation analysis on azimuth to re-
sults from conventional flow analysis methods.
B. Fourier transforms on azimuth
A Fourier series is an efficient representation if azimuth
structure approximates a constant plus a few sinusoids
whose wavelengths are integral fractions of 2π. A Fourier
representation of a peaked distribution (e.g., jet cone)
is not an efficient representation. We assume a simple
combination of the lowest few Fourier terms.
The Fourier forward transform (FT) is
ρ̃(φ) =
exp(imφ) (3)
cos(m[φ−Ψm]),
where boldface Q̃m is an event-wise complex amplitude,
Q̃m is its magnitude and Ψm its phase angle. The second
line arises because ρ̃(φ) is a real function, and the prac-
tical upper limit on m is particle number n (wavelength
∼ mean interparticle spacing). Q̃m/2π = ρ̃m is the am-
plitude of the density variation associated with the mth
sinusoid. With ri → 1 Q̃m is the corresponding num-
ber of particles in 2π if that density were uniform. Ψm
is event-wise by definition and does not require a tilde.
The reverse transform (RT) is
Q̃m =
dφ ρ̃(φ) exp(−imφ) (4)
ri exp(−imφi)
= Q̃m exp(imΨm).
For the discrete transform φ ∈ [−π, π] is partitioned
into M equal bins with bin contents r̃l and bin centers
φl. We multiply Eq. (3) by bin width δφ = 2π/M , and
the Fourier transform pair becomes
r̃l =
cos(m[φl −Ψm]) (5)
Q̃m =
l=−M/2
r̃l exp(−imφl),
where the r̃l are also random variables. The upper limit
M/2 in the first line is a manifestation of the Nyquist
sampling theorem [17]. With ri → 1 r̃l → ñl and Qm/M
is the maximum number of particles in a bin associated
with the mth sinusoid.
C. Autocorrelations and power spectra
The azimuth autocorrelation density ρA(φ∆) is a pro-
jection by averaging of the pair density on two-particle
azimuth space (φ1, φ2) onto difference axis φ∆ = φ1−φ2.
The autocorrelation concept is not restricted to peri-
odic or bounded distributions or discrete Fourier trans-
forms [16]. In what follows ρ̃A includes self pairs. The
autocorrelation density is defined as [18]
ρ̃A(φ∆) ≡
dφ ρ̃(φ) ρ̃(φ+ φ∆) (6)
i,j=1
dφ δ(φ− φi) δ(φ− φj + φ∆)
i,j=1
rirj δ(φi − φj + φ∆).
Using Eq. (3) we obtain the FT as
ρ̃A(φ∆) =
dφ × (7)
exp(imφ)×
m′=−∞
Q̃∗m′
exp(−im′[φ+ φ∆])
[2π]2
exp(−i mφ∆)
[2π]2
[2π]2
cos(mφ∆),
and the RT as
Q̃2m = n
2〈r cos(m[φ−Ψm])〉2 (8)
dφ∆ ρ̃A(φ∆) cos(mφ∆)
i,j=1
rirj cos(m[φi − φj ])
= n〈r2〉+ n(n− 1)〈r2 cos(mφ∆)〉.
= n〈r2〉+ n(n− 1)〈r2 cos2(m[φ−Ψr])〉.
Phase angle Ψm has been eliminated, and Ψr will be
identified with the reaction plane angle.
We can write the same relations for ensemble-averaged
quantities because the terms are positive-definite
ρA(φ∆) ≡
[2π]2
[2π]2
cos(mφ∆), (9)
with RT
Q2m = 2π
dφ∆ ρA(φ∆) cos(mφ∆) (10)
i,j=1
rirj cos(m[φi − φj ])
= n〈r2〉+ n(n− 1)〈r2 cos(mφ∆)〉.
We have adopted the convention Q2m = Q̃
m to lighten
the notation. Coefficients Q2m are power-spectrum ele-
ments on wave-number indexm. That FT transform pair
expresses the Wiener-Khintchine theorem which relates
power-spectrum elements Q2m to autocorrelation density
ρA(φ∆). The autocorrelation provides precise access to
two-particle correlations given enough collision events, no
matter how small the event-wise multiplicities.
D. Autocorrelation structure
The autocorrelation concept was developed in response
to the Brownian motion problem and the Langevin equa-
tion, a differential equation describing Brownian motion
which contains a stochastic term. The concept is al-
ready apparent in Einstein’s first paper on the subject
(cf. App. A). The large-scale, possibly-deterministic mo-
tion of the Brownian probe particle must be separated
from its small-scale random motion due to thermal col-
lisions with molecules. Similarly, we want to extract az-
imuth correlation structure persisting in some sense over
an event ensemble from event-wise random variations.
The autocorrelation technique is designed for that pur-
pose.
The statistical reference for a power spectrum is the
white-noise background representing an uncorrelated sys-
tem. The reference is typically uniform up to large
wave number or frequency (hence white noise), an in-
evitable part of any power spectrum from a discrete pro-
cess (point distribution). The “signal” is typically limited
to a bounded region (signal bandwidth) of the spectrum
at smaller frequencies or wave numbers.
From Eq. (8) the power-spectrum elements are
Q̃2m =
r2i +
n,n−1
rirj cos(m[φi − φj ]) (11)
= n〈r2〉+ n(n− 1)〈r2 cos(mφ∆)〉
The first term in Eq. (11) is Q̃2ref , the white-noise back-
ground component of the power spectrum common to
all spectrum elements. The second term, which we de-
note Ṽ 2m, represents true two-particle azimuth correla-
tions. Note that Q̃20 = n
2〈r2〉, whereas Ṽ 20 = n(n−1)〈r2〉.
In terms of complex (or vector) amplitudes we can write
Q̃m = Q̃ref + Ṽm, (12)
where Q̃ref represents a random walker. There is no
cross term in Eq. (11) because Qref and Vm are uncor-
related.
Inserting the power-spectrum elements into Eq. (9) we
obtain the ensemble-averaged autocorrelation density
ρA(φ∆) =
n〈r2〉
δ(φ∆) +
n(n− 1)〈r2〉
[2π]2
[2π]2
cos(mφ∆).
The first term is the self-pair or statistical noise term,
which can be excluded from ρA by definition simply
by excluding self pairs. The second term, with V 20 =
n(n− 1)〈r2〉, is a uniform component, and the third term
is the sinusoidal correlation structure. The self-pair term
is referred to in conventional flow analysis as the “auto-
correlation,” in the sense of a bias or systematic error,
but that is a notional misuse of standard mathematical
terminology. The true autocorrelation density is the en-
tirety of Eq. (13), including (in this simplified case) the
self-pair term, the uniform component and the sinusoidal
two-particle correlations.
In general, the single-particle ensemble-averaged dis-
tribution ρ0 may be structured on (η, φ). We want to
subtract the corresponding reference structure from the
two-particle distribution to isolate the true correlations.
In what follows we assume ri = 1 for simplicity, therefore
describing number correlations. We subtract factorized
reference autocorrelation ρA,ref (φ1, φ2) = ρ0(φ1) ρ0(φ2)
representing a system with no correlations, with ρ0 =
n̄/2π ≃ d2n/dηdφ in this simple example, to obtain the
difference autocorrelation
∆ρA(φ∆) = ρA − ρA,ref (14)
σ2n − n̄
[2π]2
[2π]2
cos(mφ∆).
The first term measures excess (non-Poisson) multiplicity
fluctuations in the full (pt, η, φ) acceptance. The second
term is a sum over cylindrical multipoles. We now divide
the autocorrelation difference by
ρA,ref = ρ0 = n̄/2π
to form the density ratio
ρA,ref
σ2n − n̄
2π n̄
cos(mφ∆) (15)
≡ ∆ρA[0]√
ρA,ref
∆ρA[m]√
ρA,ref
cos(mφ∆),
The first term ∆ρA[0]/
ρA,ref is the density ratio av-
eraged over acceptance (∆η,∆φ). Its integral at full
acceptance is normalized variance difference ∆σ2
(σ2n− n̄)/n̄, which we divide by acceptance integral 1×2π
to obtain the mean of the 2D autocorrelation density
∆ρA[0]/
ρA,ref = ∆σ
(∆η,∆φ)/∆η∆φ. The sinusoid
amplitudes are ∆ρA[m]/
ρA,ref = n(n− 1) ṽ2m/(2πn̄) ≡
m, defining unbiased vm. The event-wise ṽ
〈cos(mφ∆)〉 are related to conventional flow measures
vm, but may be numerically quite different for small
multiplicities due to bias in the latter. Different mean-
value definitions result in different measured quantities
(cf. App. B). Whereas the Q2m are power-spectrum ele-
ments which include the white-noise reference, the V 2m (∝
squares of cylindrical multipole moments) represent the
true azimuth correlation signal. This important result
combines several measures of fluctuations and correla-
tions within a comprehensive system.
E. Azimuth vectors
Power-spectrum elements Q2m are derived from com-
plex Fourier amplitudes Qm. In an alternate representa-
tion the complex Qm can be replaced by azimuth vectors
~Qm. The ~Qm are conventionally referred to as “flow vec-
tors” [22], but they include a statistical reference as well
as a “flow” sinusoid. The ~Vm defined in this paper are
more properly termed “flow vectors,” to the extent that
such terminology is appropriate. We refer to the ~Qm by
the model-neutral term azimuth vector and define them
by the following argument.
The cosine of an angle difference—cos(m[φ1 − φ2]) =
cos(mφ1) cos(mφ2) + sin(mφ1) sin(mφ2)—can be repre-
sented in two ways, with complex unit vectors u(mφ) ≡
exp(imφ) or with real unit vectors ~u(mφ) ≡ cos(mφ)̂ı+
sin(mφ)̂ [complex plane (ℜz,ℑz) vs real plane (x, y)].
Thus,
cos(m[φ1 − φ2]) = ℜ{u(mφ1)u∗(mφ2)} (16)
= ~u(mφ1) · ~u(mφ2)).
If an analysis is reducible to terms in cos(m[φ1 − φ2])
the same results are obtained with either representation.
Thus, we can rewrite the first line of Eq. (3) as
ρ̃(φ) =
· ~u(mφ) (17)
cos(m[φ−Ψm]),
in which case
~̃Qm =
dφ ρ̃(φ) ~u(mφ) =
ri~u(mφi) (18)
= n〈r[cos(mφ), sin(mφ)]〉
= Q̃m~u(mΨm),
an event-wise random real vector.
III. CONVENTIONAL METHODS
We now use the formalism in Sec. II to review con-
ventional flow analysis methods in a common framework.
We consider five significant papers in chronological order.
The measurement of angular correlations to detect col-
lective dynamics (e.g., parton fragmentation and/or hy-
drodynamic flow) proceeds from directivity (1983) at the
Bevalac (1 - 2 GeV/u fixed target) to transverse spheric-
ity predictions (1992) for the SPS/RHIC (
sNN = 17 -
200 GeV), then to Fourier analysis of azimuth distribu-
tions (1994, 1998) and v2 centrality trends (1999). We
pay special attention to the manipulation of random vari-
ables (RVs). RVs do not follow the algebra of ordinary
variables, and the differences are especially important for
small sample numbers (multiplicities, cf. App. B).
A. Directivity at the Bevalac
An important goal of the Bevalac HI program was ob-
servation of the collective response of projectile nucleons
to compression during nucleus-nucleus collisions, called
directed flow, which might indicate system memory of
the initial impact parameter as opposed to isotropic ther-
mal emission. Because of finite (small) multiplicities and
large fluctuations relative to the measured quantity the
geometry of a given collision may be poorly defined, but
the final state may still contain nontrivial collective infor-
mation. The analysis goal becomes separating possible
collective signals from statistical noise.
An initial search for collective effects was based on the
3D sphericity tensor S̃ =
i ~pi~pi [14, 15] described in
App. C 3. Alternatively, the directivity vector was defined
in the transverse plane [19]. In the notation of Sec. II,
including weights ri → wipti, directivity is
~̃Q1 ≡
wi~pti =
wipti~u(φi) (19)
= Q̃1~u(Ψ1),
azimuth vector ~̃Qm with m = 1 (corresponding to di-
rected flow). Event-plane (EP) angle Ψ1 estimates true
reaction-plane (RP) angle Ψr. To maintain correspon-
dence with SPS and RHIC analysis we simplify the de-
scription in [19] to the case that n single nucleons are
detected and no multi-nucleon clusters; thus a → 1 and
A→ n.
The terms in ~̃Q1 are weighted by wi = w(yi) = ±1 cor-
responding to the forward or backward hemisphere rela-
tive to the CM rapidity, with a region [−δ, δ] about mid-
rapidity excluded from the sum (wi = 0). Q̃1 then ap-
proximates quadrupole moment q̃21 derived from spher-
ical harmonic ℜY 12 ∝ sin(2θ) cos(φ), as illustrated in
Fig. 1 (left panel, dashed lines). In effect, a rotated
quadrupole is modeled by two opposed dipoles, point-
symmetric about the CM in the reaction plane. It is ini-
tially assumed that EP angle Ψ1 of vector ~Q1 estimates
RP angle Ψr, and magnitude Q1 measures directed flow.
The first part of the analysis was based on subevents—
nominally equivalent but independent parts of each
event. The dot product ~Q1A · ~Q1B for subevents A
and B of each event was used to establish the exis-
tence of a significant flow phenomenon and the angu-
lar resolution of the RP estimation via the distribu-
tion on Ψ1A − Ψ1B. The EP resolution was defined by
cos(Ψ1 − Ψr) = 2
cos(Ψ1A −Ψ1B). The magnitude of
~̃Q1 was then related to an estimate of the mean trans-
verse momentum in the RP. Integrating over rapidity
with weights wi we obtain event-wise quantities
Q̃21 =
w2i p
n,n−1
wi wj ~pti · ~ptj (20)
≃ n〈p2t 〉+ n(n− 1)〈p2t cos(φ∆)〉
≡ Q̃2ref + Ṽ 21 .
The last line makes the correspondence with the notation
of this paper.
The initial analysis in [19] used Q̃21 − Q̃2ref = Ṽ 21 =
n(n−1)〈p2x〉, assuming that x̂ is contained in the RP and
there are no non-flow correlations. Note that n(n− 1) ≡
∑n,n−1
i6=j |wiwj | contains weights wiwj implicitly. Since
wi ∼ sin[2 θ(yi)] (cf. Fig. 1 – left panel), what is actu-
ally calculated in [19] for the single-particle (no multi-
nucleon clusters) case is the pt-weighted r.m.s. mean
of the spherical harmonic ℜY 12 (θ, φ) ∝ quadrupole mo-
ment q̃21, thereby connecting rank-1 tensor ~Q1 (with
weights on rapidity) to rank-two sphericity tensor S (cf.
App. C 3). The mean-square quantity calculated is
p2x ≡
n(n− 1)
≃ sin
2(2θ)p2t cos
2(φ −Ψr)
sin2(2θ)
, (21)
a minimally-biased statistical measure as discussed in
App. B, from which we obtain px =
p2x estimating the
transverse momentum per particle in the reaction plane.
The second part of the analysis sought to obtain px(y),
the weighted-mean transverse momentum in the RP as a
function of rapidity. It was decided to determine pxi =
pti cos(φi−Ψr) for the ith particle relative to the RP, but
with Ψr estimated by EP angle Ψ1. The initial attempt
was based on
xi ≡ wi~pti · ~u(Ψ1) = wi~pti ·
j wj ~ptj
k wk ~ptk|
. (22)
Summing wi ~pti over all particles in a y bin gives
〈p′x〉 =
∑n,n−1
i6=j wiwj ~pti · ~ptj
l |wl| |
k wk ~ptk|
n〈p2t 〉+ Ṽ 21
n Q̃1
= Q̃1/n =
〈p2t 〉+
〈p2x〉,
from which we obtain ensemble mean p′x =
˜〈p′x〉. That
result can be compared directly with the linear speed in-
ferred from a random walker trajectory, which is ‘infinite’
in the limit of zero time interval (cf. App. A). The first
term in the radicand is said to produce “multiplicity dis-
tortions” in conventional flow terminology. The second
term contains the unbiased quantity.
In contrast to the first part of the analysis the sec-
ond method retains the statistical reference within Q̃1
as part of the result, so that Q̃1/n ∼
〈p2t 〉/n for small
multiplicities and/or flow magnitudes, a false signal com-
parable in magnitude to the true flow signal. The unbi-
ased directed flow px was said to be “distorted” by the
presence of the statistical reference (called unwanted self-
correlations or “autocorrelations”) to the strongly-biased
value p′x dominated by the statistical reference.
An attempt was made to remove statistical distortions
arising from self pairs by redefining ~̃Q1 → ~̃Q1i, a vector
complementary to each particle i with that particle omit-
ted from the sum. The estimator of Ψr for particle i is
then Ψ1i in ~̃Q1i ≡
j 6=i wj ~ptj = Q̃1i ~u(Ψ1i) and
xi = wi~pti · ~u(Ψ1i) = wi~pti ·
j 6=i wj ~ptj
k 6=i wk ~ptk|
.(24)
0 1 2
vm/σ = √(2nvm
2 ) ~ √(2Vm
2 /n)
n = 5
n = 10, ... , 50
√(n-1)Vm
0 1 2 3
FIG. 1: Left panel: Comparison of directed flow data from
Fig. 3(a) of the event-plane analysis and the ℜY 12 ∝
sin(2θ[ylab]) spherical harmonic, with amplitude 95 MeV/c
obtained from V 21 . Weights in the form w(ylab) are denoted
by dashed lines. The correspondence with sin(2θ[ylab]) (solid
curve) is apparent. Right panel: The EP resolution obtained
from [22] (dashed curve) and from ratio
(n− 1)V 2
/nQ′2
defined in this paper (solid and dotted curves for several n
values).
Summing over i within a rapidity bin one has
〈p′′x〉 =
l |wl|
j 6=i wiwj ~pti · ~ptj
k 6=i wk ~ptk|
〈p2x〉
(n− 1)〈p2x〉
〈p2t 〉+ (n− 2)〈p2x〉
〈p2x〉 〈cos(Ψ′m −Ψr)〉
with Q̃′1 ≡
(n− 1)〈p2t 〉+ (n− 1)(n− 2)〈p2t cos(φ∆)〉.
Since V 21 = n(n− 1)〈p2x〉 = n(n− 1)〈p2t cos(φ∆)〉 one
sees that the division by Q̃′1 is incorrect, even though
it seems to follow the chain of argument based on RP
estimation and EP resolution with correction. The
correct (minimally-biased) quantity is px ≡
p2x =
/n(n− 1). The new EP definition removes the ref-
erence term from the numerator, but Q̃′1 in the denom-
inator retains the statistical reference in p′′x. There is
the additional issue that
x2 6= x̄. Two different mean
values are represented by 〈px〉 and px =
p2x. The dif-
ference can be large for small event multiplicities.
The remaining bias was attributed to the EP resolu-
tion. The resolution correction factor derived from the
initial subevent analysis (cf. App. D) was applied to
〈p′′x〉 to further reduce bias. In Fig. 1 (right panel) we
compare the EP resolution correction from [22] (dashed
curve) with factor
(n− 1)Ṽ 2
/nQ̃′2
required to convert
〈p′′x〉 from Eq. (25) to px from Eq. (21). The agreement
is very good. Eq. (21) is the least biased and most direct
way to obtain px ≡
p2x, both globally over the detec-
tor acceptance and locally in rapidity bins, without EP
determination or resolution corrections.
In Fig. 1 (left panel) we show the data (points) for
px(ylab) from the EP-corrected analysis and the solid
curve px sin[2θ(ylab)] ∝ q21Y21(θ(y), 0), where px =
/n(n− 1) = 95 MeV/c. The agreement is good,
and the similarity of sin[2θ(ylab)] (solid curve) to weights
w(ylab) (dashed lines) noted above is apparent. Loca-
tion of the sin[2θ(ylab)] extrema near the kinematic limits
(vertical lines) is an accident of the collision energy and
nucleon mass. These results are for Ecm = 1.32 GeV ∼√
2. In the notation of this paper V 21 = 4.7
(GeV/c)2,
/n(n− 1) = 0.095 GeV/c = wpx/a, and
Qx ≡ n̄
/n(n− 1) = 2.17 GeV/c ≃
(not Q1).
By direct and indirect means (directivity and RP esti-
mation) quadrupole moment q21 ∝ V1 was measured.
B. Transverse sphericity at higher energies
The arguments and techniques in [20] suggest a smooth
transition from relativistic Bevalac and AGS energies
(collective nucleon and resonance flow) to intermediate
SPS and ultra-relativistic RHIC energies (possible trans-
verse flow, possibly QCD matter, collectivity manifested
by correlations of produced hadrons, mainly pions). For
all heavy ion collisions thermalization is a key issue.
Clear evidence of thermalization is sought, and collective
flow is expected to provide that evidence.
Two limiting cases are presented in [20] for SPS flow
measurements: 1) linear N-N superposition with no col-
lective behavior (no flow); 2) thermal equilibrium – col-
lective pressure in the reaction plane – fluid dynamics
leading to “elliptic” flow. In a hydro scenario the initial
space eccentricity transitions to momentum eccentricity
through thermalization and early pressure. The paper
considers flow measurement techniques appropriate for
the SPS and RHIC, and in particular newly defines trans-
verse sphericity St.
According to [20] the 3D sphericity tensor introduced
at lower energies [15] can be simplified in ultra-relativistic
heavy ion collisions to a 2D transverse sphericity tensor.
Sphericity is transformed to 2D by ~p → ~pt, omitting the
momentum ẑ component near mid-rapidity. Transverse
sphericity (in dyadic notation) is
2S̃t ≡ 2
~pti~pti (26)
p2ti {I + C(φi)}
≡ n〈p2t 〉 {I + α̃1 C(Ψ2)}
defining α̃1 and Ψ2 in the tensor context, with
C(φ) ≡
cos(2φ) sin(2φ)
sin(2φ) − cos(2φ)
. (27)
This α̃ definition corresponds to Eq. (3.1) of [20].
We next form the contraction of S̃t with itself
2S̃t : S̃t = 2
(~pti · ~ptj)2 (28)
p4ti + 2
p2tip
tj cos
2(φi − φj)
= 2n〈p4t 〉+ 2n(n− 1)〈p4t cos2(φ∆)〉
= n(n+ 1)〈p4t 〉+ n(n− 1)〈p4t cos(2φ∆)〉
using the dyadic contraction notation A : B ≡ AabBab,
with the usual summation convention. That self-
contraction of a rank-2 tensor can be compared to the
more familiar self-contraction of a rank-1 tensor ~̃Q2 · ~̃Q2 =
Q̃22 = n〈p2t 〉 + n(n − 1)〈p2t cos(2φ∆)〉. The quantity
2[S̃t : S̃t]ref = n(n + 1)〈p4t 〉 is the (uncorrelated) refer-
ence for the rank-2 contraction, whereas Q̃2
2,ref = n〈p2t 〉
is the reference for the rank-1 contraction. Subtracting
the rank-2 reference contraction gives
S̃t : S̃t − [S̃t : S̃t]ref
= n(n− 1)〈p4t cos(2φ∆)〉
≃ 〈p2t 〉Ṽ 22 , (29)
which relates transverse sphericity to two-particle corre-
lations in the form Ṽ 22 = Q̃
2 − Q̃22,ref , thus establishing
the exact correspondence between S̃t and ~̃Q2.
From the definition of α1 in Eq. 26 above and Eq. (3.1)
of [20] we also have
2 S̃t : S̃t = n2〈p2t 〉2(1 + α̃21) (30)
which implies
n2〈p2t 〉2α̃21 = n2 σ̃2p2t + n〈p
t 〉 (31)
+ n(n− 1)〈p4t cos(2φ∆)〉
≃ n2 σ̃2p2t + 〈p
t 〉Q̃22. (32)
The first relation is exact, given the definition of α̃1,
but produces a complex statistical object containing the
event-wise variance of p2t in its numerator and random
variable n2 in its denominator. For 1/n → 0 it is true
that α̃1 → 〈cos(2[φ−Ψr])〉 (since Ψ2 → Ψr also), but α̃1
is a strongly biased statistic for finite n.
The definition
α̃2 =
〈p2x − p2y〉
〈p2x + p2y〉
from Eq. (2.5) of [20] seems to imply α̃2 = 〈p2t cos(2[φ−
Ψr])〉/〈p2t 〉 → 〈cos(2[φ−Ψr])〉, assuming that x̂ lies in the
RP. However, the latter relation fails for finite multiplic-
ity [the effect of the statistical reference or self-pair term
n in Eq. (31)] because each of event-wise 〈p2x〉 and 〈p2y〉 is
a random variable, and their independent random varia-
tions do not cancel in the numerator. The exact relation
is the first line of
n2〈p2t 〉2α̃22 = n2〈p2t cos(2[φ−Ψ2])〉2 (34)
≃ n2〈p2t 〉〈pt cos(2[φ−Ψ2])〉2
= 〈p2t 〉 Q̃22
α̃2 ≃ Q̃2/Q̃0 6= Ṽ2/Q̃0
The second line is an approximation which indicates that
α̃2 is more directly related to Q̃2 than is α̃1. But Q
2 is a
poor substitute for V 22 which represents true two-particle
azimuth correlations in a minimally-biased way by incor-
porating a proper statistical reference. The effect of the
reference contribution is termed a ‘distortion’ in [20].
C. Fourier series I
Application of Fourier series to azimuth particle dis-
tributions was introduced in [21]. Fourier analysis is de-
scribed as model independent, providing variables which
are “easy to work with and have clear physical interpreta-
tions.” Sinusoids or harmonics are associated with trans-
verse collective flow, the model-dependent language fol-
lowing [20] closely. To facilitate comparisons we convert
notation in [21] to that used in this paper: r(φ) → ρ(φ),
(xm, ym) → ~Qm, ψm → Ψm, vm → Qm and ṽm → Vm.
According to the proposed method density ρ̃(φ) repre-
sents, within some (pt, η) acceptance, an event-wise par-
ticle distribution on azimuth φ including weights ri =
1, pti or Eti. The FT in terms of azimuth vectors is
ρ̃(φ) =
riδ(φ− φi) (35)
· ~u(mφ)
cos(m[φ−Ψm]),
and the RT is
~̃Qm =
ri~u(mφi) ≡ Q̃m~u(mΨm), (36)
forming a conventional Fourier transform pair [cf.
Eqs. (3) and (4)]. Scalar amplitude Q̃m =
i ri~u(mφi) ·
~u(mΨm) = n〈r cos(m[φ − Ψm])〉 is the proposed flow-
analysis quantity. Q̃m is said to measure the flow mag-
nitude, and Ψm estimates reaction-plane angle Ψr. It
is proposed that Q̃m(η) evaluated within bins on η may
characterize complex “event shapes” [densities on (η, φ)].
As with directivity and sphericity, multiplicity fluctu-
ations are seen as a major obstacle to flow analysis with
Fourier series. Finite multiplicity is described as a source
of ‘bias’ which must be suppressed. It is stated that
(Ṽm,Ψr) are the “parameter[s] relevant to the magnitude
of flow,” whereas the observed (Q̃m,Ψm) are biased flow
estimators. In the limit 1/n → 0 the two cases would
be identical. A requirement is therefore placed on min-
imum event-wise multiplicity in a “rapidity slice.” To
solve the finite-number problem the paper proposes to
use the event frequency distribution on Q̃2m from event-
wise Fourier analysis to measure flow.
If correlations are zero then
Q̃2m → Q̃2ref =
r2i = n〈r2〉 ≡ σ̃2, (37)
with 〈r2〉 ≡ σ20 . By fitting the frequency distribution on
Q̃2m with a model function it is proposed to obtain Ṽm as
the unbiased flow estimator. The fitting procedure is said
to require sufficiently large event multiplicities to obtain
Ṽm unambiguously.
The distribution on Q̃2m is derived as follows (cf.
Fig. 2 – left panel). The magnitude of statistical ref-
erence ~̃Qref (a random walker) has probability distri-
bution ∝ exp(−Q̃2ref/2Q2ref), with Q2ref = n〈r2〉. But
~̃Qref = ~̃Qm − ~̃Vm, therefore (cf. Fig. 2 – left panel)
Q̃2ref = Q̃
m + Ṽ
m − 2Q̃mṼm cos(m[Ψm −Ψr]), (38)
exp(−Q̃2ref/2Q2ref) → ρ(Q̃m,Ψm; Ṽm,Ψr). (39)
When integrated over cos(m[Ψm −Ψr]) there results the
required probability distribution on Q̃2m, with fit param-
eter Vm. The distribution on Q̃
m is said to show a ‘non-
statistical’ shape change from which Vm can be inferred
by a model fit “free from uncertainties in event-wise de-
termination of the reaction plane.” It is also proposed
to use ρ(Q̃m,Ψm; Ṽm,Ψr) to determine the EP resolu-
tion cos(Ψm − Ψr) by integrating over Q̃2m and using the
resulting projection on cos(Ψm − Ψr) to determine the
ensemble mean (Fig. 3 of [21]).
While one could extract ensemble-average V 2m at some
level of accuracy by fitting the frequency distribution on
Q̃2m with a model function, we ask why go to that trou-
ble when V 2m is easily obtained as a variance difference?
Instead of Eq. (38) we simply write
Q̃2m = Q̃
ref + Ṽ
m, (40)
where the cross term is zero on average and Q̃2ref = n〈r2〉
represents the power-spectrum white noise, which is the
same for all m (i.e., ‘white’). If the vector mean values
are zero that is a relation among variances. Ensemble
mean V 2m = Q
m −Q2ref is therefore simply determined.
For the EP resolution we factor Ṽ 2m
Ṽ 2m = n(n− 1)〈r2 cos2(m[φ−Ψr])〉 (41)
= n(n− 1)〈r2 cos2(m[φ−Ψm])〉 cos2(m[Ψm −Ψr]).
We use Q̃m = n〈r cos(m[φ − Ψm])〉 and the assumption
that random variable Ψm − Ψr is uncorrelated with φ−
Ψm to obtain
cos2(m[Ψm −Ψr]) =
n Ṽ 2m
(n− 1) Q̃2m
, (42)
which defines the EP resolution of the full n-particle
event in terms of power-spectrum elements (cf. App. D).
The square root of that expression is plotted in Fig. 2
(right panel) as the solid curves for several values of n̄.
The solid and dotted curves are nearly identical to those
in Fig. 1.
=       -
vm/σ = √(2nvm
2 ) ~ √(2Vm
2 /n)
n = 5
n = 10, ... , 50
√(n-1)
0 1 2 3
FIG. 2: Left panel: Distribution of event-wise elements of ~̃Qm
components determined by the gaussian-distributed random
walker ~̃Qref and possible correlation component ~̃Vm. Right
panel: Reaction-plane resolution estimator 〈cos(mδΨmr)〉,
with δΨmr = Ψm−Ψr, determined from fits to a distribution
on Q̃2m as in the left panel (dashed curve), and from Eq. 42
for several values of n (solid and dotted curves).
As in other flow papers there is much emphasis on
insuring adequate multiplicities to reduce bias to a man-
agable level, because an easily-determined statistical ref-
erence is not properly subtracted to reveal the contribu-
tion from true two-particle correlations in isolation. For√
2nvm > 1 in Fig. 2 (right panel) the EP is meaningful;
bias of event-wise quantities relative to the EP is man-
ageable. For
2nvm < 1/2 the EP is poorly defined,
and ensemble-averaged two-particle correlations are the
only reliable measure of azimuth correlations. In either
case EP estimation is only justified when a non-flow phe-
nomenon is to be studied relative to the reaction plane.
D. Fourier series II
A more elaborate review of flow analysis methods
based on Fourier series is presented in [22]. The ap-
proach is said to be general. The event plane is ob-
tained for each event. The event-wise Fourier ampli-
tude(s) q̃m ≡ Q̃m/Q̃0 relative to the EP are corrected for
the EP resolution as obtained from subevents. We mod-
ify the notation of the paper to vobsm → q̃m and wi → ri to
maintain consistency within this paper. We distinguish
between the unbiased ṽm ≡ Ṽm/Ṽ0 and the biased q̃m.
According to [22], in the 1/n→ 0 limit (Q̃m → Vm, no
tildes, no random variables) the dependence of the single-
particle density on azimuth angle φ integrated over some
(pt, y) acceptance can be expressed as a Fourier series of
the form
ρ(φ) =
1 + 2
vm cos [m(φ−Ψr)]
, (43)
with reaction-plane angle Ψr. As we have seen, the factor
2 comes from the symmetry on indexm for a real-number
density, not an arbitrary choice as suggested in [22]. In
this definition V0/2π has been factored from the Fourier
series in [21]. The Fourier “coefficients” vm in this form
(actually coefficient ratios) are not easily related to the
power spectrum. In the 1/n → 0 limit the coefficients
are vm = 〈cos(m[φ−Ψr])〉.
In the analysis of finite-multiplicity events, reaction-
plane angle Ψr defined by the collision (beam) axis and
the collision impact parameter is estimated by event
plane (EP) angle Ψm, with Ψm derived from event-wise
azimuth vector ~̃Qm (conventional flow vector)
~̃Qm =
ri~u(mφi) ≡ Q̃m~u(mΨm). (44)
The finite-multiplicity event-wise FT is
ρ̃(φ) =
1 + 2
q̃m cos [m(φ−Ψm)]
. (45)
with q̃m = 〈cos(m[φ−Ψm])〉; e.g., q̃2 ≃ α̃2 (cf. Eq. (34)).
According to the conventional description the EP an-
gle is biased by the presence of self pairs, unfortunately
termed the “autocorrelation effect” or simply “autocor-
relation” in conventional flow analysis [19, 22], whereas
autocorrelations and cross-correlations are distributions
on difference variables used for decades in statistical anal-
ysis to measure correlation structure on time and space.
As in [19], to eliminate “autocorrelations” EP angle Ψmi
is estimated for each particle i from complementary flow
vector ~Qmi =
j 6=i rj~u(2φj) = Qmi~u(mΨmi), a form of
subevent analysis with one particle vs n− 1 particles (cf.
App. D).
In the conventional description the event-plane res-
olution results from fluctuations δΨr ≡ Ψm − Ψr of
the event-plane angle Ψm (or Ψmi) relative to the true
reaction-plane angle Ψr (e.g., due to finite particle num-
ber). The EP resolution is said to reduce the observed
q̃m relative to the true value ṽm:
q̃m =
ṽ2m (46)
= 〈cos (m[Ψr −Ψm])〉 · ṽm.
The EP resolution (first factor, second line) is obtained in
a conventional flow analysis in two ways: the frequency
distribution on Ψm − Ψr discussed in Sec. III C and the
subevent method discussed in App. D. A parameteri-
zation from the frequency-distribution method reported
in [22] is plotted as the dashed curve in Fig. 2 (right
panel). There is good agreement with the simple expres-
n/(n− 1)Vm/Qm obtained in Eq. (42), which also
follows from Eq. (46).
Two methods are described for obtaining vm without
an event-plane estimate, with the proviso that large event
multiplicities in the acceptance are required. The first, in
terms of the conventional flow vector, is expressed (with
ri → 1) as (Eq. (26) of [22])
Q2m = n̄+ n
2 v̄2m. (47)
That expression is constrasted with the exact event-wise
treatment, where for each event we can write
Q̃2m = n+
cos [m(φi − φj)] (48)
= n+ n(n− 1)〈cos(mφ∆)〉
≡ n+ n(n− 1)ṽ2m
→ Q2m = n̄+ n(n− 1)ṽ2m = n̄+ V 2m
Note the similarity with fluctuations measured by num-
ber variance σ2n = n̄+∆σ
n, where the second term on the
RHS is an integral over two-particle number correlations
and the first is the uncorrelated Poisson reference (again,
the self-pair term in the autocorrelation).
There are substantial differences between the two Q2m
formulations above, especially for smaller multiplicities.
Eq. (48) is unbaised for all n and provides a simple way
to obtain V 2m = n(n− 1)ṽ2m = Q2m− n̄. The conventional
method uses a complex fit to the frequency distribution
on Q̃2m to estimate Vm as in [21]. Why do that when such
a simple alternative is available?
E. v2 centrality dependence
In [23] the expected trend of v2 with A-A centrality for
different collision systems is discussed in a hydro context.
It is stated that the v2 centrality trend should reveal the
degree of equilibration in A-A collsions. The centrality
dependence of v2/ǫ should be sensitive to the “physics of
the collision”—the nature of the constituents (hadrons
or partons) and their degree of thermalization or collec-
tivity. “It is understood that such a state requires (at
least local) thermalization of the system brought about
by many rescatterings per particle during the system
evolution...v2 is an indicator of the degree of equilibra-
tion.”
Thermalization is related to the number of rescatter-
ings, which also strongly affects elliptic flow according to
this hydro interpretation. In the full hydro limit, corre-
sponding to full thermalization where the mean free path
λ is much smaller than the flowing system, relation v2 ∝ ǫ
is predicted, with ǫ the space eccentricity of the initial
A-A overlap region [20]. Conversely, in the low-density
limit (LDL), where λ is comparable to or larger than
the system size, a different model predicts the relation
v2 ∝ ǫA1/3/λ, where A1/3/λ estimates the mean num-
ber of collisions per “particle” during system evolution
to kinetic decoupling [24]. In the LDL case v2 ∝ ǫ 1S
where S = πRxRy is the (weighted) cross-section area of
the collision and ǫ =
R2y−R
R2y+R
is the spatial eccentricity.
Those trends are further discussed in App. E.
According to the combined scenario, comparison of
the centrality dependence of v2 at energies from AGS
to RHIC may reveal a transition from hadronic to ther-
malized partonic matter. The key expectation is that at
some combination(s) of energy and centrality v2/ǫ tran-
sitions from an LDL trend (monotonic increase) to hydro
(saturation), indicating (partonic) equilibration.
However, it is important to note two things: 1) That
overall description is contingent on the strict hydro sce-
nario. If the quadrupole component of azimuth correla-
tions arises from some other mechanism then the descrip-
tions in [20] and [24] are invalid, and v2 does not reveal
the degree of thermalization. 2) v2 is a model-dependent
and statistically-biased quantity motivated by the hydro
scenario itself. The model-independent measure of az-
imuth quadrupole structure is V 22 /n̄ ≡ n̄ v22 (defining an
unbiased v2). It is important then to reconsider the az-
imuth quadrupole centrality and energy trends revealed
by that measure to determine whether a hydro interpre-
tation is a) required by or even b) permitted by data.
IV. IS THE EVENT PLANE NECESSARY?
A key element of conventional flow analysis is estima-
tion of the reaction plane and the resolution of the esti-
mate. We stated above that determination of the event
plane is irrelevant if averaged quantities are extracted
from an ensemble of event-wise estimates. The reaction-
plane angle is relevant only for study of nonflow (minijet)
structure relative to the reaction plane on φΣ, the sum
(pair mean azimuth) axis of (φ1, φ2). In contrast to con-
cerns about low multiplicities in conventional flow analy-
sis, proper autocorrelation techniques accurately reveal
“flow” correlations (sinusoids) and any other azimuth
correlations, even in p-p collisions and even within small
kinematic bins. In this section we examine the necessity
of the event plane in more detail.
The reaction plane (RP), nominally defined by the
beam axis and the impact parameter between centers of
colliding nuclei, is determined statistically in each event
by the distribution of participants. The RP is estimated
by the event plane (EP), defined statistically in each
event by the azimuth distribution of final-state particles
in some acceptance. We now consider how to extract
vm relative to a reaction plane estimated by an event
plane in each event. Several different flow measures are
implicitely defined in conventional flow analysis
vm ≡ cos(m[φ−Ψr]) ideal case, 1/n→ 0 (49)
ṽ2m ≡ 〈cos2(m[φ−Ψr])〉 unbiased estimate
q̃m ≡ 〈cos(m[φ−Ψm])〉 self-pair bias
ṽ′m ≡ 〈cos(m[φi −Ψmi])〉 reduced bias
ṽ′m is the event-wise result of a conventional flow analysis.
The ensemble average ṽ′m must be corrected for the “EP
resolution” which we now determine.
The basic event-wise quantities, starting with an inte-
gral over two-particle azimuth space, are
Ṽ 2m ≡
n,n−1
~u(mφi) · ~u(mφj) (50)
= n(n− 1)〈cos(mφ∆)〉
= n(n− 1)〈cos2(m[φ−Ψr])〉
≡ n(n− 1)ṽ2m.
For the limiting case of subevents A and B with A a single
particle the subevent azimuth vector complementary to
particle i is
~̃Qmi ≡
j 6=i
~u(mφj) = Q̃mi~u(mΨmi). (51)
We make the following rearrangment
Ṽ 2m =
~u(mφi) ·
j 6=i
~u(mφj) (52)
Q̃mi cos(m[φi −Ψmi])
≡ nQ̃′m〈cos(m[(φ −Ψ′m])〉
n(n− 1)ṽ2m = nQ̃′m ṽ′m
ṽm = ṽ
and, since Q̃′m ≃
n− 1 + (n− 1)(n− 2)〈cos(mφ∆)〉,
we identify the EP resolution as
〈cos(m[Ψ′m −Ψr])〉 =
where the primes refer to a subevent with n−1 particles.
That expression, for full events with multiplicity n, is
plotted in Fig. 2 (right panel) for several choices of n. An
n-independent universal curve on V 2m/n̄ is multiplied by
n-dependent factor
n/(n− 1), where n is the number
of samples in the event or subevent. The dashed curve is
the parameterization from [22].
The right panel indicates that for large nṽ2m single-
particle reaction-plane estimates can provide a “flow”
measurement with manageable bias. For small nṽ2m the
EP resolution averaged over many events is itself a “flow”
measurement, even though the reaction plane is inacces-
sible in any one event. Ṽ 2m is determined from the same
underlying two-particle correlations by other means—the
only difference is how pairs are grouped across events. In
App. D the EP resolution is determined for the case of
equal subevents A and B
From this exercise we conclude that event-plane deter-
mination is irrelevant for the measurement of cylindrical
multipole moments (“flows”). Following a sequence of
analysis steps in the conventional approach based on de-
termination of and correction for the EP estimate, the
event plane cancels out of the flow measurement. What
results from the conventional method is approximations
to signal components of power-spectrum elements which
can be determined directly in the form V 2m/n̄ = n̄ v
obtained with the autocorrelation method, which defines
an unbiased version of vm. Event-plane determination
can be useful for study of other event-wise phenomena in
relation to azimuth multipoles.
V. 2D (JOINT) AUTOCORRELATIONS
We now return to the more general problem of angular
correlations on (η1, η2, φ1, φ2). We consider the analysis
of azimuth correlations based on autocorrelations, power
spectra and cylindrical multipoles without respect to an
event plane in the context of the general Fourier trans-
form algebra presented in Sec. II. We seek a comprehen-
sive method which treats η and φ equivalently.
In conventional flow analysis there are two concerns
beyond the measurement of flow in a fixed angular accep-
tance: a) study flow phenomena in multiple narrow ra-
pidity bins to characterize the overall “three-dimensional
event shape,” analogous to the sphericity ellipsoid but
admitting of more complex shapes over some rapidity in-
terval, and b) remove nonflow contributions to flow mea-
surements as a systematic error. Maintaining adequate
bin multiplicities to avoid bias is strongly emphasized in
connection with a). The contrast between such individ-
ual Fourier decompositions on single-particle azimuth in
single-particle rapidity bins and a comprehensive analy-
sis in terms of two-particle joint angular autocorrelations
is the subject of this section
A. Stationarity condition
In Fig. 3 we show two-particle pair-density ratios r̂ ≡
ρ/ρref on (η1, η2) (left panel) and (φ1, φ2) (right panel)
for mid-central 130 GeV Au-Au collisions [12]. The hat
on r̂ indicates that the number of mixed pairs in ρref
has been normalized to the number of sibling pairs in
ρ. In each case we observe approximate invariance along
sum axes ηΣ = η1 + η2 and φΣ = φ1 + φ2. In time-series
analysis the equivalent invariance of correlation structure
on the mean time is referred to as stationarity, implying
that averaging pair densities along the sum axes loses no
information. The resulting averages are autocorrelation
distributions on difference axes η∆ and φ∆.
0.9996
0.9998
1.0002
1.0004
1.0006
1.0008
1.001
0.999
0.9992
0.9995
0.9997
1.0002
1.0005
1.0007
1.001
1.0012
1.0015
FIG. 3: Normalized like-sign pair-number ratios r̂ = ρ/ρref
from central Au-Au collisions at 130 GeV for (η1, η2) (left
panel) and (φ1, φ2) (right panel) showing stationarity—
approximate invariance along sum diagonal x1 + x2.
In Fig. 3 (right panel) one can clearly see the cos(2φ∆)
structure conventionally associated with elliptic flow
(quadrupole component). However, there are other con-
tributions to the angular correlations which should be
distinguished from multipole components, accomplished
accurately by combining the two angular correlations into
one joint angular autocorrelation.
B. Joint autocorrelation definition
With the autocorrelation technique the dimensional-
ity of pair density ρ(η1, η2, φ1, φ2) can be reduced from
4D to 2D without information loss provided the dis-
tribution exhibits stationarity. Expressing pair density
ρ(η1, η2, φ2, φ2) → ρ(ηΣ, η∆, φΣ, φ∆) in differential form
d4n/dx4 we define the joint autocorrelation on (η∆, φ∆)
ρA(η∆, φ∆) ≡
dη∆dφ∆
d2n(ηΣ, η∆, φΣ, φ∆)
dηΣdφΣ
ηΣ φΣ
by averaging the 4D density over ηΣ and φΣ within a
detector acceptance. The autocorrelation averaging on
ηΣ is equivalent to 〈dn/dη〉η ≈ n(∆η)/∆η at η = 0. The
magnitude is still the 4D density d4n/dη1dη2dφ1dφ2, but
it varies only on the two difference axes (η∆, φ∆).
In Fig. 4 we illustrate two averaging schemes [18]. In
the left panel we show the averaging procedure applied
to histograms on (x1, x2) as in Fig. 3. Index k denotes
the position of an averaging diagonal on the difference
axis. In the right panel we show a definition involving
pair cuts applied on the difference axes. Pairs are his-
togrammed directly onto (η∆, φ∆). Periodicity on φ im-
plies that the averaging interval on φΣ/2 is 2π indepen-
dent of φ∆. However, η is not periodic and the averaging
interval on ηΣ/2 is ∆η − |η∆|, where ∆η is the single-
particle η acceptance [18]. The autocorrelation value for
average
FIG. 4: Autocorrelation averaging schemes on xΣ = x1 + x2
for a prebinned single-particle space (left panel) and for pairs
accumulated directly into bins on the two-particle difference
axis (right panel).
a given x∆ is the bin sum along that diagonal divided by
the averaging interval.
C. Autocorrelation examples
In Fig. 5 we show pt joint angular autocorrelations for
Hijing central Au-Au collisions at 200 GeV [25]. The left
panel shows quench-on collisions. The right panel shows
quench-off collisions. Aside from an amplitude change
Hijing shows little change in correlation structure from
N-N to central Au-Au collsions. Those results can be
contrasted with examples from an analysis of real RHIC
data [25] shown in Fig. 6.
0.002
0.004
0.006
0.008
0.012
0.002
0.004
0.006
0.008
0.012
FIG. 5: 2D pt angular autocorrelations from Hijing central
Au-Au collisions at 200 GeV for quench on (left panel) and
quench off (right panel) simulations.
In N-N or p-p collisions there is no obvious quadrupole
component. However, that possibility is not quanti-
tatively excluded by the data and requires careful fit-
ting techniques. For correlation structure which is not
sinusoidal there is no point to invoking the Wiener-
Khintchine theorem on φ∆ and transforming to a power
spectrum. Instead, model functions specifically suited to
such structure (e.g., 1D and 2D gaussians) are more ap-
propriate. The power-spectrum model is appropriate for
some parts of ρA, depending on the collision system. We
consider hybrid decompositions in Sec. VII.
D. Comparison with conventional methods
We can compare the power of the autocorrelation tech-
nique with the conventional approach based on EP es-
timation in single-particle η bins. There are concerns
in the single-particle approach about bias or distortion
from low bin multiplicities. In contrast, the autocorrela-
tion method is applicable to any collision system with any
multiplicity, as long as some pairs appear in some events.
The method is minimally biased, and the statistical error
in ∆ρA/
ρref is determined only by the number of bins
and the total number of events in an ensemble.
There is interest within the conventional context in
“correlating” flow results in different η bins. “A three-
dimensional event shape can be obtained by correlat-
ing and combining the Fourier coefficients in different
longitudinal windows” [21]. That goal is automatically
achieved with the joint autocorrelation technique but has
not been implemented with the conventional method.
The content of an autocorrelation bin at some η∆ rep-
resents a covariance averaged over all bin pairs at ηa and
ηb satisfying η∆ = ηa − ηb. If a flow component is iso-
lated (e.g., V 22 /n̄), a given autocorrelation element rep-
resents normalized covariance na nbṽ
n̄an̄b, where
ṽ2mab would be the event-wise result of a ‘subevent’ or
‘scalar-product’ analysis between bins a and b. Averaged
over many events one obtains good resolution on rapid-
ity and azimuth for arbitrary structures, without model
dependence or bias.
VI. FLOW AND MINIJETS (NONFLOW)
The relation between azimuth multipoles and minijets
is a critically important issue in heavy ion physics which
deserves precise study. We should carefully compare
multipole structures conventionally attributed to hydro-
dynamic flows and parton fragmentation dominated by
minijets in the same analysis context. The best arena for
that comparison is the 2D (joint) angular autocorrelation
and corresponding power-spectrum elements. Before pro-
ceeding to autocorrelation structure we consider nonflow
in the conventional flow context
A. Nonflow and conventional flow analysis
In conventional flow analysis azimuth correlation struc-
ture is simply divided into ‘flow’ and ‘nonflow,’ where
the latter is conceived of as non-sinusoidal structure of
indeterminant origin. The premise is that all sinusoidal
structure represents flows of hydrodynamic origin. It is
speculated that nonflow is due to resonances, HBT and
jets, including minijets. Various properties are assigned
to nonflow which are said to distinguish it from flow [26].
Nonflow is by definition non-sinusoidal and is not corre-
lated with the RP, thus it can appear perpendicular to
the RP.
The multiplicity dependence of nonflow is said to be
quite different from flow, where “multiplicity depen-
dence” can sometimes be read as centrality dependence.
For instance, ~̃Qma · ~̃Qmb = Ṽma Ṽmb cos(Ψma−Ψmb) if A,
B are disjoint, since there are no self pairs. The ensemble
average then measures covariance V 2mab. For m = 2 it is
claimed that the nonflow component of V 2
2ab is ∝ n̄ c [22].
Therefore, V 22 /n̄ ∝ c, a constant for nonflow—no cen-
trality dependence. But V 2m/n̄ ∝ ∆ρA/
ρref which, for
the minijets dominating nonflow, is very strongly depen-
dent on centrality [9, 12]; the conventional assumption
is incorrect. The above notation is inadequate because
minijets (nonflow) should not be included in the Fourier
power spectrum. They should be modeled by different
functional forms which we consider in the next section.
B. Cumulants
Another strategy for isolating flow from nonflow is to
use higher cumulants [27]. The basic assumption (a phys-
ical correlation model) is that flow sinusoids are collec-
tive phenomenon characteristic of almost all particles,
whereas nonflow is a property only of pairs, termed “clus-
ters.” That scenario is said to imply that vm should be
the same no matter what the multiplicity, whereas non-
flow should fall off as some inverse power of n.
For instance, by subtracting v2[4] (four-particle cumu-
lant) from v2[2] (two-particle cumulant) one should ob-
tain “nonflow” as the difference (cf. Eq. (10) of [26]). In
Fig. 31 of [26] we find a plot of g2 = Npart (v
2 [2]−v22 [4]) ∝
Npart/nch ×∆ρ/
ρref . Multiplying g2 by n
part we
obtain a measure of minijet correlations per participant
pair. That ‘nonflow’ component increases rapidly with
centrality (and therefore n), consistent with actual mea-
surements of minijet centrality trends. The incorrect fac-
tor in the definition of g2 removes a factor 2x increase
from peripheral to central in the minijet centrality trend,
thus suppressing the centrality dependence of ‘nonflow.’
C. Counterarguments
The conventional flow analysis method requires a com-
plex strategy to distinguish flow from nonflow in projec-
tions onto 1D azimuth difference φ∆. An intricate and
fragile system results, with multiple constraints and as-
sumptions. The assumptions are not a priori justified,
and must be tested. ‘Flow’ isolated with those assump-
tions can and does contain substantial systematic errors.
Claims about the multiplicity (centrality) dependence
of ‘nonflow’ (independent of centrality or slowly varying)
are unsupported speculations without basis in experiment.
In fact, detailed measurements of minijet centrality de-
pendence [9, 12] are quite inconsistent with typical as-
sumptions about nonflow. Multiplicity (centrality) de-
pendence of flow measurements is further compromised
by biases resulting from improper statistical methods,
especially true for small multiplicities or peripheral colli-
sions. Such biases can masquerade as physical phenom-
Finally, it is assumed that nonflow has no correlation
with the RP, thus implying the ability of and need for
the EP to distinguish flow from nonflow. But nonflow
(minijets) should be strongly correlated with the EP (jet
quenching), and such correlations should be measured.
That is the main subject of paper II in this two-part se-
ries. Precise decomposition of angular correlations into
‘flow’ sinusoids and minijet structure is realized with 2D
joint angular autocorrelations combined with proper sta-
tistical techniques.
Non-flow is a suite of physical phenomena, each wor-
thy of detailed study. In the conventional approach this
physics is seen in limited ways by various projections and
poorly-designed measures and described mainly by spec-
ulation. With more powerful analysis methods it is pos-
sibl to separate flow from the various sources of ‘nonflow’
reliably and identify those sources as interesting physical
phenomena.
VII. STRUCTURE OF THE JOINT ANGULAR
AUTOCORRELATION IN A-A COLLISIONS
We now return to the 2D angular autocorrelation. By
separating its structure into a few well-defined compo-
nents we obtain an accurate separation of multipoles,
minijets and other phenomena. Minijets and “flows” can
be compared quantitatively within the same analysis con-
text.
Each bin of an autocorrelation is a comparison of two
“subevents.” The notional term “subevent” represents a
partition element in conventional math terminology (e.g.,
topology, cf. Borel measure theory). An “event” is a dis-
tribution in a bounded region of momentum space (de-
tector acceptance), and a subevent is a partition element
thereof. A distribution can be partitioned in many ways:
by random selection, by binning the momentum space,
by particle type, etc. A uniform partition is a binning,
and the set of bin entries is a histogram.
A bin in an angular autocorrelation represents an av-
erage over all bin pairs in single-particle space separated
by certain angular differences (η∆, φ∆). The bin contents
represent normalized covariances averaged over all such
pairs of bins. The notional “scalar product method,” re-
lating two subevents in conventional flow analysis, is al-
ready incorporated in conventional mathematical meth-
ods developed over the past century as covariances in
bins of an angular autocorrelation. Using 2D angular au-
tocorrelations we easily and accurately separate nonflow
from flow. “Nonflow” so isolated has revealed the physics
of minijets—hadron fragments from the low-momentum
partons which dominate RHIC collisions.
A. Minijet angular correlations
Minijet correlations are equal partners with multi-
pole correlations on difference-variable space (η∆, φ∆).
Minimum-bias jet angular correlations (dominated by
minijets) have been studied extensively for p-p and Au-
Au collisions at 130 and 200 GeV [9, 10, 11, 12]. Those
structures dominate the “nonflow” of conventional flow
analysis. In p-p collisions minijet structure—a same-side
peak (jet cone) and away-side ridge uniform on η∆—
are evident for hadron pairs down to 0.35 GeV/c for
each hadron. Parton fragmentation down to such low
hadron momenta is fully consistent with fragmentation
studies over a broad range of parton energies (e.g., LEP,
HERA) [28].
80-90%
0.001
0.002
0.003
0.004
0.005
45-55%
-0.005
0.005
0.015
80-90%
0.001
0.002
0.003
0.004
45-55%
-0.002
0.002
0.004
0.006
0.008
0.012
FIG. 6: 2D pt angular autocorrelations from Au-Au collisions
at 200 GeV for 80-90% central collisions (left panels) and 45-
55% central collisions (right panels). In the lower panels si-
nusoids cos(φ∆) and cos(2φ∆) have been subtracted to reveal
“nonflow” structure.
In Fig. 6 we show autocorrelations obtained by inver-
sion of pt fluctuation scale (bin-size) dependence [9]. The
upper-left panel is 80-90% central and the upper-right
panel is 45-55% central Au-Au collisions. Correlation
structure is dominated by a same-side peak and mul-
tipole structures (sinusoids). Subtracting the sinusoids
reveals the minijet structure in the bottom panels and
illustrates the precision with which flow and nonflow can
be distinguished. The negative structure surrounding the
same-side peak at lower right is an interesting and unan-
ticipated new feature [9].
B. Decomposing 2D angular autocorrelations: a
controlled comparison
Based on extensive analysis [9, 10, 11, 25] we find
three main contributions to angular correlations in RHIC
nuclear collisions: 1) transverse fragmentation (mainly
minijets), 2) longitudinal fragmentation (modeled as
“string” fragmentation), 3) azimuth multipoles (flows).
Longitudinal fragmentation plays a reduced role in heavy
ion collisions. In this study we focus on the interplay be-
tween 1) and 3), transverse parton fragmentation and
azimuth multipoles, as the critical analysis issue for az-
imuth correlations in A-A collisions.
The 2D joint autocorrelation ρA(η∆, φ∆) is the basis
for decomposition. The criteria for distinguishing az-
imuth multipoles from minijet structure are η∆ depen-
dence and sinusoidal φ∆ dependence. Structure with si-
nusoidal φ∆ dependence and η∆ invariance is assigned to
azimuth multipoles. Other structure, varying generally
on (η∆, φ∆), is assigned in this exercise to minijets. We
adopt the decomposition
ρA(η∆, φ∆) = ρj(η∆, φ∆) + ρm(φ∆), (55)
where j represents (mini)jets and m represents multi-
poles. That decomposition is reasonable within a limited
pseudorapidity acceptance, e.g., the STAR TPC accep-
tance [29]. Over a larger acceptance other separation
criteria must be added.
To illustrate the separation process we construct artifi-
cial autocorrelations combining flow sinusoids and mini-
jet structure with centrality dependence taken from mea-
surements. We add statistical noise appropriate to a typ-
ical event ensemble of a few million Au-Au collisions. We
then fit the autocorrelations with model functions and χ2
 0.1994E-03/    23
P1  0.1456E-01
P2  0.4811E-02
-0.02
-0.01
0.005
0.015
0.025
0 2 4
FIG. 7: Simulated three-component 2D angular autocorrela-
tion for 80-90% central Au-Au collisions at 200 GeV (upper
left), model of data distribution from fitting (upper-right),
autocorrelation with eta-acceptance triangle imposed (lower
left) and sinusoid fit to 1D projection of lower-left panel
(lower-right).
minimization. We compare the resulting fit parameters
with the input parameters. We then project the 2D au-
tocorrelations onto φ∆ and fit the results with a sinusoid.
The result of the 1D fit represents the product of a con-
ventional flow analysis if the EP resolution is perfectly
corrected and there is no statistical bias in the method.
We then compare the resulting sinusoid amplitudes.
In Fig. 7 we show an analysis for 80-90% central Au-Au
collisions, which are nearly N-N collisions. The upper-
right panel is an accurate model of data from p-p col-
lisions [11]. The upper-left panel is the simpler repre-
sentation for this exercise with added statistical noise.
The difference in constant offsets is not relevant to the
exercise. At lower left is the distribution “seen” by a
conventional flow analysis, which simply integrates over
(projects) the η dependence. The projection includes a
triangular acceptance factor on η∆ imposed on the joint
autocorrelation, resulting in distortions that affect the si-
nusoid fit. The lower-right panel is the projection onto
φ∆ with χ
2 fit. The 1D fit gives ∆ρ[2]/
ρref = 0.0052,
compared to 0.0013 from the 2D fit
In Fig. 8 we show the same analysis applied to 5-10%
Au-Au collisions. The model distribution, derived from
observed data trends, is dominated by a dipole term
∝ cos(φ∆) and an elongated same-side jet peak, although
the combination closely mimics a quadrupole ∝ cos(2φ∆)
(elliptic flow). The lower-left panel shows the effect of the
η-acceptance triangle on η∆ applied to the upper-right
panel which is implicit in any projection onto φ∆ (obvious
in this 2D plot). The resulting projection on φ∆ is shown
 0.2331E-01/    23
P1  0.8826E-01
P2  0.1026
0 2 4
FIG. 8: Same as the previous figure but for 5-10% central Au-
Au collisions at 200 GeV. Note the pronounced effect of the
η-acceptance triangle in the lower-left panel (relative to the
upper-right panel) resulting from projection of space (η1, η2)
onto its difference axis.
at lower right. The 1D fit gives ∆ρ[2]/
ρref = 0.103,
compared to 0.060 from the 2D fit. The differences be-
tween 1D and 2D fits are much larger than the differences
between input and fitted parameters in the previous ex-
ercise. They reveal the limitations of conventional flow
analysis.
In Fig. 9 we give a summary of results for twelve cen-
trality classes, nine 10% bins with mean values 95% · · ·
15%, plus two 5% bins with mean values 7.5% and 2.5%.
The twelfth class is b = 0 (0%), constructed by extrapo-
lating the parameterizations from data. The solid curves
and points represent the parameters inferred from 2D fits
to the joint angular autocorrelations. The dashed curves
represent the model parameters used to construct the
simulated distributions. There is good agreement at the
percent level.
m = 1
projection
onto φ∆
peak amplitude
2 4 6
2 4 6
FIG. 9: A parameter summary for the previous two figures.
The solid curves are input model parameters for m = 1 and
m = 2 sinusoids cos(mφ∆) and a same-side 2D gaussian with
two widths and peak amplitude which approximate 200 GeV
Au-Au collisions. The dashed curves are the results of fits to
the model 2D autocorrelations exhibiting excellent accuracy.
The dash-dot curve represents 1D fits to projections on φ∆
(previous lower-right panels) corresponding to conventional
flow analysis.
The fits to 1D projections on φ∆ (dash-dot curve) how-
ever differ markedly from the 2D fit results and the in-
put parameters. The differences are very similar to the
changes of conventional flow measurements with different
strategies to eliminate “nonflow.” This exercise demon-
strates that with the 2D autocorrelation there is no guess-
work. We can distinguish the multipole contributions
from the minijet contributions. The 2D angular autocor-
relation provides precise control of the separation.
C. v2 in various contexts
In Fig. 10 we contrast the results of 1D conven-
tional flow analysis (dashed curves) and extraction of the
quadrupole amplitude from the 2D angular autocorrela-
tion (solid curves). We make the correspondence
∆ρA[2]√
2π n̄
n̄v22
among variables, with n̄/2π modeling ρ0 = d2n/dηdφ.
That expression defines a minimally-biased v2. We make
the comparison in four plotting formats. The upper-right
panel is most familiar from conventional flow analysis,
where v2 is plotted vs participant nucleon number. The
trends can be compared with Fig. 13 of [30] (open circles
vs solid stars). Also included in that panel is the trend
v2 ∼ 0.22 ǫ predicted by hydro for thermalized A-A col-
lisions at 200 GeV (dotted curve, [31]).
npart
0.22 ε
projection
onto φ∆
2D fit
2 4 6
0 100 200 300 400
2 4 6 0
2 4 6
FIG. 10: The quadrupole component in various plotting con-
texts derived from the previous model exercise. Conventional
flow measure v2 is shown in the upper panels. ∆ρ[2]/
based on Pearson’s normalized covariance is shown in the
lower panels. The upper-right panel shows v2 vs npart in the
conventional plotting format. The lower-left panel repeats
the left panel of the previous figure, with conventional 1D fit
results shown as the dashed curve. Dashed and solid curves
correspond in all panels. ν estimates the mean N-N encoun-
ters per participant pair. The dotted curve at lower right is
the error function used to generate the model quadrupole am-
plitudes. It’s correspondent for v2 is at upper left. The dotted
curve at upper right is eccentricity ǫ from a parameterization.
The remaining panels are plotted on parameter ν =
2nbin/npart, the ratio of N-N binary encounters to par-
ticipant pairs which estimates the mean participant path
length in number of encountered nucleons, a geometri-
cal measure. Comparing the upper panels we see that ν
treats peripheral and central collisions equitably, whereas
npart or ncharge compresses the important peripheral re-
gion into a small interval.
In the lower-left panel we plot per-particle density ra-
tio ∆ρA[2]/
ρref vs ν. That quantity, when extracted
from 2D fits, rises from near zero for peripheral collisions
to a maximum for mid-central collisions, falling toward
zero again for b = 0. In contrast to v2, which is the
square root of a per-pair correlation measure, the per-
particle density ratio reflects the trend for “flow” in the
sense of a current density. “Flow” is small for periph-
eral collisions and grows rapidly with increasing nucleon
path length. The trend with centrality is intuitive. The
values obtained from the 1D projection per conventional
flow analysis (dashed curve) are consistently high, espe-
cially for central collisions, exhibiting a strong systematic
bias. The 1D fit procedure (identical to the “standard”
and two-particle flow methods) confuses minijet structure
with quadrupole structure.
In the lower-right panel we show the density ratio di-
vided by initial spatial eccentricity ǫ defined by a pa-
rameterization derived from a Glauber simulation and
plotted as the dotted curve in the upper-right panel [31].
The trend from 2D fits (solid curve) is closely approx-
imated by a simple error function (dotted curve) with
half-maximum point at the center of the ν range. In
fact, the dotted curve is the basis for generating the input
quadrupole amplitudes for our model, and the small devi-
ation of the solid curve from the dotted curves in upper-
left and lower-right panels reveals the systematic error or
bias in the 2D fitting procedure (∼ 20% at ν ∼ 1, < 5%
at ν = 5.8). The dashed curves from the conventional 1D
fits show large relative deviations from the input trend,
especially for peripheral collisions where the emergence of
collectivity is of interest and for central collisions where
the issue of thermalization is most important.
Comparing the solid curve to existing v2 data in the
format of the upper-right panel shown in [30] (Fig. 13)
indicates that our simple formulation at lower right
(dotted curve) is roughly consistent with analysis of
flow data based on four-particle cumulants. The re-
sult from data generated with the same model and ana-
lyzed with the conventional flow analysis method (dashed
curve in upper-right panel) also agrees with the conven-
tional method applied to real RHIC data. Our model
may therefore indicate an underlying simplicity to the
quadrupole mechanism which is not hydrodynamic in ori-
gin. The model centrality trend for v2 is certainly incon-
sistent with the hydro expectation v2 ∝ ǫ [20, 23], as
demonstrated in the upper-right panel (cf. App. E).
VIII. DISCUSSION
A. Conventional flow analysis
The overarching premise of conventional flow analy-
sis is that in the azimuth distribution of each collision
event lies evidence of collective phenomena which must
be discovered to establish event-wise thermalization. The
1/n → 0 hydro limit shapes the analysis strategy, and
flow manifestations are the principal goal. Finite event
multiplicities are seen as a major source of systematic
error, as are correlation structures other than flow. Mul-
tiple strategies are constructed to deal with non-flow and
finite multiplicities. A stated advantage of the conven-
tional method is that the Fourier coefficients can be cor-
rected. The great disadvantage is theymust be corrected.
From the perspective of two-particle correlation anal-
ysis, especially in the context of autocorrelations and
power spectra, the conventional program leaves much to
be desired. The conventional analysis is essentially an
attempt to measure two-particle correlations with single-
particle methods combined with RP estimation, similar
to the use of trigger particles in high-pt jet analysis. But,
by analysis of the algebraic structure of conventional flow
analysis we have demonstrated that RP estimation does
not matter to the end result. Without a proper statistical
reference conventional analysis results contain extrane-
ous contributions from the statistical reference which are
partially ‘corrected’ in a number of ways. The improper
treatment of random variables incorporates sources of
multiplicity-dependent bias in measurements, and the fi-
nal results are questionable.
Flow measure vm is nominally the square root of per-
pair correlation measure V 2m/n(n− 1). vm centrality
trends are thus nonintuitive and misleading (e.g., “el-
liptic flow” decreases with increasing A-A centrality).
The situation is similar to per-pair fluctuation measure
Σpt , which provides a dramatically misleading picture of
pt fluctuation dependence on collision energy [10]. In
contrast, per-particle correlation measures provide intu-
itively clear results and often make dynamical correla-
tion mechanisms immediately obvious. In particular, the
mechanisms behind “nonflow” in the form of minijets are
clearly apparent when correlations are measured by per-
particle normalized covariance density ∆ρ/
ρref in a 2D
autocorrelation.
B. Autocorrelations and nonflow
In conventional flow analysis it is proposed to measure
flow in narrow rapidity bins (“strips”) so as to develop
a three-dimensional picture of event structure. However,
there has been little implementation of that proposal. In
contrast, the joint angular autocorrelation by construc-
tion contains all possible covariances among pseudora-
pidity bins within a detector acceptance. The ideal of
full event characterization is thereby realized.
The angular autocorrelation is the optimum solution to
a geometry problem—how to reduce the six-dimensional
two-particle momentum space to two-dimensional sub-
spaces with minimum distortion or information loss. The
autocorrelation is the unique solution to that problem,
involving no model assumptions. The ensemble-averaged
angular autocorrelation contains all the correlation infor-
mation obtainable from a conventional flow analysis, but
with negligible bias and no sensitivity to individual event
multiplicities.
Because it is a two-dimensional representation the
angular autocorrelation is far superior for separating
“flow” (multipoles) from “nonflow” (minijets), as we have
demonstrated. A simple exercise demonstrates that sep-
aration is complete at the percent level, whereas the con-
ventional method admits crosstalk at the tens of percent
level. Precise separation leads to new physics insights
from the multipoles and minijets so revealed.
C. Collision centrality dependence
Collision centrality dependence is of critical impor-
tance in the comparison of flow and minijets. Parton col-
lisions and hydrodynamic response to early pressure have
very different dependence on impact parameter and col-
lision geometry, especially for peripheral collisions. Pe-
ripheral A-A collisions should approach p-p (N-N) col-
lisions, and correlation structure may change rapidly in
mid-peripheral collisions as collective phenomena develop
there. The possible onset of collective behavior in mid-
peripheral collisions and reduction in more central col-
lisions are of major importance for understanding the
relation of minijets to flow. The conventional flow analy-
sis method is severely limited for peripheral collisions.
In contrast, correlation measure ∆ρ/
ρref , centrality
measure ν and associated centrality techniques described
in [32] are uniquely adapted to cover all centrality regions
down to N-N with excellent accuracy.
D. Physical interpretations
Because similar flow measurement techniques have
been applied at Bevalac and RHIC energies with sim-
ilar motivations it is commonly assumed that azimuth
multipoles have a common source over a broad collision
energy range—hydrodynamic flows, collective response
to early pressure. The hydro mechanism was proposed
as the common element in [20] and persists as the lone
interpretation of azimuth multipoles in HI collisions to
date.
At Bevalac and AGS energies it is indeed likely that
azimuth multipoles result from ‘flow’ of initial-state nu-
cleons in response to early pressure, with consequent
final-state correlations of those nucleons—a true hydro
phenomenon. However, at SPS and RHIC energies the
source of azimuth multipoles inferred from final-state
produced hadrons (mainly pions) may not be hydrody-
namic, in contrast to arguments by analogy with lower
energies. Other sources of multipole structure should be
considered [34, 35]. Multipoles at higher energies could
arise at the partonic or hadronic level, early or late in the
collision, with collective motion or not, and if collective
then implying thermalization or not.
The chain of argument most often associated with
elliptic flow asserts that observation of flow as a col-
lective phenomenon demonstrates that a thermalized
medium (QGP) has been formed which responds hy-
drodynamically to early pressure and converts an ini-
tial configuration-space eccentricity to a corresponding
quadrupole moment in momentum space.
However, nonflow in the form of minijets provides con-
tradictory evidence. Minijet centrality trends indicate
that thermalization is incomplete, and substantial mani-
festations of initial-state parton scattering remain at ki-
netic decoupling [9, 10, 12]. Precision studies of mini-
jet centrality dependence (ν dependence) indicate that a
large fraction of the minijet structure expected from lin-
ear superposition of N-N collisions (no thermalization)
persists in central Au-Au collisions. That contradic-
tion requires more complete experimental characteriza-
tion and careful theoretical study [36].
Arguments based on interpreting the quadrupole com-
ponent as hydrodynamic flow exclude alternative phys-
ical mechanisms. Aside from minijet systematics there
are other hints that a different mechanism might be re-
sponsible for azimuth multipoles. In Fig. 10 we showed
that flow measurements based on four-particle cumulants
(with bias sources and nonflow thereby reduced) are best
described by a trend (solid curves) that is inconsistent
with the hydro expectation v2 ∝ ǫ. The trend is in-
stead simply described in terms of per-particle measure
ρref and two shape parameters relative to ǫ.
We question the theoretical assumption that ǫ should
be simply related to v2 as opposed to some other mea-
sure of the azimuth quadrupole component. We expect
a priori and find experimentally that variance measures,
integrals over two-particle momentum space, more typ-
ically scale linearly with geometry parameters. Thus,
∆ρ[2]/
ρref ∝ n̄v22 may be more closely related to ǫ,
and the relation may or may not be characteristic of a
hydro scenario.
IX. SUMMARY
In conclusion, we have reviewed Fourier transform the-
ory, especially the relation of autocorrelations to power
spectra, essential for analysis of angular correlations in
nuclear collisions. In that context we have reviewed
five papers representative of conventional flow analysis
and have related the methods and results to autocorre-
lation structure and spherical and cylindrical multipole
moments.
We have examined the need for event-plane evaluation
in correlationmeasurements and find that it is extraneous
to measurement of azimuth multipole moments. The EP
estimate drops out of the final ensemble average.
We have introduced the definition of the 2D (joint) an-
gular autocorrelation and considered the distinction be-
tween flows (cylindrical multipoles) and nonflow (domi-
nated by minijet structure) in conventional flow analysis
and criticized the basic assumptions used to distinguish
the two in that context.
Based on measured minijet and flow centrality trends
we have constructed a simulation exercise in which model
autocorrelations of known composition are combined
with statistical noise from a typical event ensemble and
fit with a model function consisting of a few simple com-
ponents, first as a 2D autocorrelation and second as a 1D
projection on azimuth difference axis φ∆. We show that
the 2D fit returns input parameters accurately at the
percent level, whereas the 1D fit, representing conven-
tional flow analysis, deviates systematically and strongly
from the input. Comparisons with published flow data
indicate that the observed bias in the simulation is ex-
actly the difference attributed to “nonflow” in conven-
tional measurements.
By comparing our simple algebraic model of
quadrupole centrality dependence to data we observe
that the trend v2 ∝ ǫ is not met for any collision sys-
tem, nor is there asymptotic approach to such a trend.
That observation raises questions about the relevance of
hydrodynamics to phenomena currently attributed to el-
liptic flow at the SPS and RHIC.
This work was supported in part by the Office of Sci-
ence of the U.S. DoE under grant DE-FG03-97ER41020.
APPENDIX A: BROWNIAN MOTION
There is a close analogy between Brownian motion and
the azimuth structure of nuclear collisions. The long his-
tory of Brownian motion and its mathematical descrip-
tion can thus provide critical guidance for the analysis
of particle distributions. Brownian motion (more gen-
erally, random motion of particles suspended in a fluid)
was modeled by Einstein as a diffusion process (random
walk) [7]. He sought to test the “kinetic-molecular” the-
ory of thermodynamics and provide direct observation of
molecules. Paul Langevin developed a differential equa-
tion to describe such motion, which included a stochastic
term representing random impulses delivered to the sus-
pended particle by molecular collisions. Jean Perrin and
collaborators performed extensive measurements which
confirmed Einstein’s predictions and provided definitive
evidence for the reality of molecules [8].
1. The quasi-random walker
We model a 2D quasi-random walker (including
nonzero correlations) as follows. The walker position
is recorded in equal time intervals δt. After n steps,
with step-wise displacements r sampled randomly from a
bounded distribution, the walker position relative to an
arbitrary starting point is, in the notation of this paper,
i ri~u(φi), where ri is the i
th displacement. The
squared total displacement is then
R2 = n〈r2〉+ n(n− 1)〈r2 cos(φ∆)〉. (A1)
The first term, linear in n (or t), was described by Ein-
stein. The second term could represent “drift” of the
walker due to deterministic response to an external in-
fluence. The composite is then termed “Brownian motion
with drift,” a popular model for stock markets and other
quasi-random processes. Measuring multipole moments
on azimuth in nuclear collisions is formally equivalent
to measuring “drift” terms on time in the quasi-random
walk of a charged particle suspended in a molecular fluid
within a superposition of oscillating electric fields. There
are many other applications for Eq. (A1).
For a true random walk consisting of uncorrelated steps
Einstein expressed 〈r2〉/δt ≡ d ·2D (random walk in d di-
mensions) in terms of diffusion coefficient D. The second
term 〈r2 cos(φ∆)〉 ≡ (δt)2v2x represents a possible deter-
ministic component (correlations), with x̂ the direction of
an applied “force.” In that case successive angles φi are
correlated, and the result is a macroscopic nonstochastic
drift of the walker trajectory.
The fractal dimension of a random walk [first term
in Eq. (A1)] is df = 2. The trajectory is therefore
a “space-filling” curve in 2D configuration space. The
appropriate measure of trajectory size is area, and the
rate of size increase is the diffusion coefficient (rate of
area increase). In contrast, the second term in Eq. (A1)
represents a deterministic trajectory whose nominal di-
mension is 1 (modulo the extent of curvature, which in-
creases the dimension above 1). Therefore, the appropri-
ate measure of trajectory size is length, and speed is the
correct rate measure. For Brownian motion with drift
the trajectory dimension is not well-defined, depending
on the relative magnitudes of the drift and stochastic
terms, and the concept of speed is therefore ambiguous.
Attempts to measure the linear speed of Brownian mo-
tion in the nineteenth century failed because of the frac-
tal structure of random walks. From the structure of
Eq. (A1) the average speed over interval ∆t = nδt is
R2/(∆t)2 ∼ 〈r2/(δt)2〉/n, and the limiting case for
∆t = nδt → 0 is the so-called “infinite speed of diffu-
sion.” That topological oddity is formally equivalent to
the “multiplicity bias” of conventional flow analysis.
2. Brownian motion and nuclear collisions
We now consider the close analogy between Einstein’s
theory of Brownian motion and the measurement of p2x
in a nuclear collision, using directivity as an example.
Just as ~R is the vector total displacement of a quasi-
randomwalker in 2D configuration space, ~Q1 is the vector
total displacement of a quasi-random walker (event-wise
particle ensemble) in 2D momentum space. After n steps
the squared displacements are
R2 = n2 δt2 v′2x = n δt 4D+ n(n− 1) δt2 v2x (A2)
Q21 = n
2 p′2x = n 〈p2t 〉+ n(n− 1)p2x.
4Dδt is the increase in area per step of a random walker
in 2D configuration space. 〈p2t 〉 is the increase in area per
step (per particle) of a random walker in 2D momentum
space, playing the same role as the diffusion coefficient.
The RHS first term in the first line is the subject of Ein-
stein’s 1905 Brownian motion paper. Its measurement by
Perrin confirmed the reality of molecules and the validity
of Boltzmann’s kinetic theory.
As noted, attempts to measure mean speed v′x of a par-
ticle in a fluid failed because speed is the wrong rate mea-
sure for trajectory size increase. Speed measurements
decreased with increasing sample number or observation
time. It was not until Einstein’s formulation and later
mathematical developments that the topology of the ran-
dom walk and its consequences became apparent. Initial
attempts at the Bevalac to measure px in the form p
using directivity failed for the same reason. Corrections
were developed to approximate the unbiased quantity px,
and the failure was attributed to multiplicity bias or ‘au-
tocorrelations.’ Ironically, the autocorrelation distribu-
tion is the ideal method to access the unbiased quantity
in either case.
3. Einstein and autocorrelations
To provide a statistical description of Brownian motion
Einstein introduced the autocorrelation concept with the
following language [7].
Another important consideration can be re-
lated to this method of development. We
have assumed that the single particles are all
referred to the same co-ordinate system. But
this is unnecessary, since the movements of
the single particles are mutually independent.
We will now refer the motion of each parti-
cle to a co-ordinate system whose origin co-
incides at the [arbitrary] time t = 0 with the
[arbitrary] position of the center of gravity of
the particle in question; with this difference,
that [probability distribution] f(x, t)dx now
gives the number of the particles whose x co-
ordinate has increased between the time t = 0
and the time t = t, by a quantity which lies
between x and x+ dx.
Einstein’s function f(ξ, τ) is a 2D autocorrelation
which satisfies the diffusion equation. The solution is
a gaussian on x relative to an arbitrary starting point
(thus defining difference variables ξ = x − xstart and
τ = t − tstart), with 1D variance σ2ξ = 2Dτ . The au-
tocorrelation is sometimes called a two-point correlation
function or two-point autocorrelation. The angular auto-
correlation is a wide-spread and important analysis tool,
e.g., in astrophysics, nuclear collisions and many other
fields.
4. Wiener, Khintchine, Lévy and Kolmogorov
The names Wiener, Lévy, Kolmogorov and Khintchine
figure prominently in the copious mathematics derived
from the Brownian motion problem. Norbert Wiener led
efforts to provide a mathematical description of Brown-
ian motion, abstracted to aWiener process, a special case
of a Lévy process (generalization of a discrete random
walk to a continuous random process) [37]. The Wiener-
Khintchine theorem provides a power-spectrum represen-
tation for stationary stochastic processes such as random
walks, for which a Fourier transform does not exist. We
have acknowledged the theorem with our Eq. (10).
The analysis of azimuth structure in nuclear collisions
in terms of angular autocorrelations is based on power-
ful mathematics developed throughout the past century.
Autocorrelations make it possible to study azimuth struc-
ture for any event multiplicity down to p-p collisions with
as little as two detected particles per event. The effects of
“non-flow” can be eliminated from “flow” measurements
(and vice versa) without model dependence or guesswork.
The Brownian motion problem and Einstein’s fertile so-
lution inform two central issues for studies of the correla-
tion structure of nuclear collisions: analysis methodology
and physics interpretation.
APPENDIX B: RANDOM VARIABLES
A random variable represents a set of samples from
a parent distribution. The outcome of any one sample
is unpredictable (i.e., random), but through statistical
analysis an ensemble of samples can be used to infer prop-
erties (statistics – results of algorithms applied to a set of
samples) of the parent distribution. Sums over particles
and particle pairs of kinematic quantities are the primary
random variables in analysis of nuclear collision data.
1. The algebra of random variables
Products and ratios of random variables behave non-
intuitively because random variables don’t obey the alge-
bra of ordinary variables. E.g., factorization of random
variables results in the spawning of covariances. The
approximation xy ≃ x̄ ȳ common in conventional flow
analysis is a source of systematic error (bias) because
xy = x̄ ȳ + xy − x̄ ȳ. The omitted term is a covariance.
Such covariances play a role in statistics similar to QM
commutators, with 1/n↔ ~. Conventional flow analysis
assumes the 1/n → 0 limit for some random variables,
and the results are undependable for small multiplicities.
Similarly, improper treatment of ratios of random vari-
ables results in infinite series of covariances. E.g.,
x/n =
δx · δn
x̄ n̄
x · (δn)2
x̄ n̄2
+ · · · ), (B1)
with (δn)2/n̄ ≡ σ2n/n̄ ∼ 1 − 2. Thus, the common ap-
proximation x/n ≃ x̄/n̄ can result in significant n- and
physics-dependent (x-n covariances) bias for small n.
In this paper we distinguish between event-wise and
ensemble-averaged quantities and do not employ en-
semble averages of ratios of random variables. We in-
clude event-wise factorizations and ratios only to sug-
gest qualitative connections with conventional flow anal-
ysis. E.g., we consider Ṽ 2m ≡ n(n − 1)〈cos(mφ∆)〉 with
〈cos(mφ∆)〉 = 〈cos2(m[φ − Ψr])〉 ≡ ṽ2m. But, ṽ2m 6=
V 2m/n(n− 1) 6= v̄2m. vm as typically invoked in conven-
tional flow analysis is not a well-defined statistic.
2. Statistical references
The concept of a statistical reference is largely absent
from conventional flow analysis. By ‘statistical reference’
we mean a quantity or distribution which represents an
uncorrelated system, a system consistent with indepen-
dent samples from a fixed parent distribution (central
limit conditions [5]). Concerns about ‘bias’ from low mul-
tiplicities [19, 21, 22] typically relate to the presence of an
unsubtracted and unacknowledged statistical reference
in the final result. Finite multiplicity fluctuations are
then said to produce systematic errors, false azimuthal
anisotropies, a problem masking true collective effects.
In the limit 1/n → 0 the statistical reference may in-
deed become negligible compared to the true correlation
structure. However, its presence for nonzero 1/n is a po-
tential source of systematic error which may block access
to important small-multiplicity systems (peripheral col-
lisions and/or small kinematic bins). In general, if the
statistical reference is not correctly subtracted the result
is increasingly biased with smaller multiplicities. Identi-
fication and subtraction of the proper reference is one of
the most important tasks in statistical analysis.
Use of the term ‘statistical’ to mean ‘uncorrelated’ is
misleading (e.g., ‘statistical’ vs ‘dynamical’). All ran-
dom variables and their fluctuations about the mean are
‘statistical.’ Some random variables and their statistics
are reference quantities, representing systems that are by
construction uncorrelated (independent sampling from a
fixed parent). We therefore label statistical reference
quantities ‘ref,’ not ‘stat.’
3. Random variables and Fourier analysis
In the context of Fourier analysis the basic finite-
number (Poisson) statistical reference is manifested as
the delta-function component in the autocorrelation den-
sity Eq. (13) and the white-noise constant term n〈r2〉 in
the event-wise power spectrum. Other reference compo-
nents may arise from two-particle correlations which are
not of interest to the analysis (e.g., detector effects) and
which may be revealed in mixed-pair distributions. A
clear distinction should always be maintained between
the reference and the sought-after correlation signal.
Careful attention to random-variable algebra is es-
pecially important in a Fourier analysis. The power-
spectrum elements and autocorrelation density must sat-
isfy the transform equations both for each event and after
ensemble averaging. In conventional flow analysis that
condition is often not satisfied. For instance, Ṽ 2m and V
satisfy the FT transforms and Wiener-Khintchine the-
orem before and after ensemble averaging respectively,
whereas the vm do not.
4. Minimally-biased random variables
It is frequently stated in the conventional flow liter-
ature that flow analysis must insure sufficiently large
multiplicities. The operating assumption in the design
of conventional flow methods is the continuum limit
1/n → 0, with inevitable bias for smaller multiplicities.
However, careful reference design and algebraic manipu-
lation of random variables makes possible precise treat-
ment of event-wise multiplicities down to n = 1. Some
statistical measures perform consistently no matter what
the sample number. The full multiplicity range is essen-
tial to measure azimuth multipole evolution with central-
ity down to N-N and p-p collisions, so that A-A “flow”
phenomena may be connected to phenomena observed in
elementary collisions and understood in a QCD context.
Since multiplicity necessarily varies strongly with cen-
trality, multiplicity-dependent bias in flow measurements
is unacceptable, and every means should be used to in-
sure minimally-biased statistics. To achieve that end
analysis methods must carefully transition from safe
event-wise factorizations (as featured in this paper) to
ensemble averages minimally biased for all n. Linear
combinations of powers of random variables, e.g., vari-
ances and covariances, satisfy a linear algebra. Such in-
tegrals of two-particle momentum space are nominally
free of bias.
APPENDIX C: MULTIPOLES AND SPHERICITY
The 1D Fourier transform on azimuth is part of a larger
representation of angular structure. The encompassing
context is a 2D multipole decomposition on (θ, φ) repre-
sented by the sphericity tensor, with the spherical har-
monics Y m2 as elements. In limiting cases submatrices
of the sphericity tensor reduce to “cylindrical harmon-
ics” cos(mφ), part of the 1D Fourier representation on
azimuth.
The central premise of a multipole representation
is that the final-state particle angular distribution on
[θ(yz), φ] is efficiently represented by a few low-order
spherical harmonics (SH) Y ml (θ, φ). At the Bevalac,
sphericity tensor S containing spherical harmonics Y m2 as
elements was introduced. Directivity ~Q1, simply related
to Y 12 , was employed to represent a rotated quadrupole
as a dipole pair antisymmetric about the collision mid-
point. At lower energies (Bevalac, AGS) the quadrupole
principal axis may be rotated to a large angle with re-
spect to the collision axis and Y 12 dominates. At higher
energies and near midrapidity (θ ∼ π/2) the dominant
SH is Y 22 .
1. Spherical harmonics
The spherical harmonics are defined as
Ylm(Ω) =
2l + 1
· (l −m)!
(l +m)!
Pml (cos θ) e
imφ, (C1)
where Pml (θ) is an associated Legendre function [38]. An
event-wise density on the unit sphere can be expanded
ρ̃(Ω) =
Q̃lm Ylm(Ω) (C2)
Q̃lm =
dΩY ∗lm(Ω)ρ̃(Ω)
Y ∗lm(Ωi)
= n〈Y ∗lm(Ω)〉,
where Ω → (θ, φ) and dΩ ≡ d cos(θ)dφ. The FTs on φ
form a special case of those relations when ρ̃(Ω) is peaked
near θ ∼ π/2. The Ylm are orthonormal and complete:
dΩYlm(θ, φ)Yl′m′(θ, φ) = δll′δmm′ (C3)
Ylm(θ, φ)Y
′, φ′) = δ(Ω− Ω′).
2. Multipoles
The spherical harmonics are model functions for single-
particle densities on (θ, φ). The coefficients of the mul-
tipole expansion of a distribution are complex spherical
multipole moments describing 2l poles and defined as en-
semble averages of the spherical harmonics over the unit
sphere weighted by an angular density.
The following relation is defined by analogy with the
expansion of an electric potential in spherical harmonics,
in this case on momentum space ~p rather than configu-
ration space ~r [38]
ρ̃(p′,Ω′)
|~p− ~p ′|
2l+ 1
Ylm(Ω)
. (C4)
The coefficients are the event-wise spherical multipole
moments
q̃lm ≡
p2dp dΩ plρ̃(p,Ω)Y ∗lm(Ω) (C5)
pli Y
lm(Ωi)
= n〈pl Y ∗lm(Ω)〉.
Eq. (C2) is the special case for p restricted to unity (i.e.,
distribution on the unit sphere).
In general, ℜY mm ∝ sinm(θ) cos(mφ), and the cos(mφ)
are by analogy “cylindrical harmonics” [42]. The ensem-
ble average of a cylindrical harmonic over the unit circle
weighted by 1D density ρ(φ) results in complex cylindri-
cal multipole moments Qm. The Fourier coefficients Qm
obtained from analysis of SPS and RHIC data are there-
fore cylindrical multipole moments describing 2m poles.
E.g., m = 2 denotes a quadrupole moment and m = 4
denotes an octupole moment,.
If nonflow contributions (i.e., structure rapidly vary-
ing on η or y) are present, a multipole decomposition of
ρ(θ, φ) is no longer efficient, and the inferred multipole
moments are difficult to interpret physically (e.g., flow in-
ferences per se are biased). In Sec. V we describe a more
differential method for representing angular structure us-
ing two-particle joint angular autocorrelations on differ-
ence axes (η∆, φ∆). Given a decomposition of ρ(θ, φ)
based on variations on η∆ we can distinguish cylindri-
cal multipoles accurately from “nonflow” structure (cf.
Sec. VII).
3. Sphericity
The sphericity tensor has been employed in both jet
physics and flow studies. A normalized 3D sphericity
tensor was defined in [14] to search for initial evidence
of jets in e+-e− collisions. A decade later sphericity was
introduced to the search for collective nucleon flow in
heavy ion collisions [15]. The close connection between
flow and jets continues at RHIC, where we seek the rela-
tion between minijets and “elliptic flow.”
Event-wise sphericity S̃ is a measure of structure in
single-particle density ρ(θ, φ) on the unit sphere. We
use dyadic notation to reduce index complexity, analo-
gous to vector notation ~̃Qm =
i ri~u(mφi). S̃ (with
ri → pi) describes a 3D quadrupole with arbitrary orien-
tation. Given ~p = p [sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)] ≡
p ~u(θ, φ) we have
2S̃ ≡ 2
~pi~pi = 2
p2i ~u(θi, φi) ~u(θi, φi) (C6)
p2i Ũ(θi, φi)
= n〈p2 Ũ(θ, φ)〉,
the last being an event-wise average, where
U(θ, φ) = sin2(θ) I + Y(θ, φ) (C7)
Y(θ, φ) = (C8)
sin2(θ) cos(2φ) sin2(θ) sin(2φ) sin(2θ) cos(φ)
sin2(θ) sin(2φ) − sin2(θ) cos(2φ) sin(2θ) sinφ
sin(2θ) cos(φ) sin(2θ) sinφ 3 cos2(θ) − 1
In terms of event-wise quadrupole moments q̃2m derived
from the Y2m
2S̃ = n
p2 sin2(θ)
I (C9)
2l+ 1
ℜq̃22 ℑq̃22 −ℜq̃21
ℑq̃22 −ℜq̃22 −ℑq̃21
−ℜq̃21 −ℑq̃21
an event-wise estimator of angular structure on the unit
sphere, its reference defined by 2S̃ref = n〈p2 sin2(θ)〉 I,
and p2 sin2(θ) = p2t . The sphericity tensor of [14] was
normalized to Ŝ = S/n〈p2〉.
Note that
Q̃ = 3S̃ − n
I (C10)
is the traceless Cartesian quadrupole tensor appearing in
the Taylor expansion of the (~p equivalent of the) electro-
static potential [38]. We have defined instead
Q̃′ = 3S̃ −
I, (C11)
an alternative quadrupole tensor wherein each element
is a single spherical quadrupole moment. The difference
lies in the diagonal elements: linear combinations of the
ℜq̃2m in the diagonal elements of Q̃ are simplified to
single moments in Q̃′. The ensemble mean of both ten-
sors for an uncorrelated (spherically symmetric) system
or system with event-wise quadrupole orientations ran-
domly varying is the null tensor (all elements zero).
APPENDIX D: SUBEVENTS
The “subevent” is a notional re-invention of partition-
ing/binning, the latter having a history of more than a
century in mathematics. In conventional flow analysis
subevents are groups of particles in an event segregated
on the basis of random selection, charge, strangeness,
PID or a kinematic variable such as pt, y or η. The scalar-
product method [30] is based on a covariance between two
single-particle bins (nominally equal halves of an event).
The subevent method is thus a restricted reinvention of
a common concept in multiparticle correlation analysis:
determining covariances among all pairs of single-particle
bins at some arbitrary binning scale – a two-particle cor-
relation function. Diagonal averages of such distributions
are the elements of autocorrelations.
In the language of conventional flow analysis one way
to eliminate statistical reference Q2ref from Q
m is to par-
tition events into a pair of disjoint (non-overlapping)
subevents A, B [19]. In that case ~̃Qma· ~̃Qmb = ~̃Vma · ~̃Vmb =
nanbṽ
mab, a covariance. The partition may be asym-
metric (unequal particle numbers) and may be as small
as a pair of particles. In addition to eliminating the
self-pair statistical reference such partitioning is said to
reduce nonflow correlation sources, depending on their
physical origins and the partition definition [30]. We as-
sume for simplicity that there is no nonflow contribution.
Subevent pairs can be used to determine the event-plane
resolution for subevents A, B and full events.
First, we consider the symmetric case, defining equiv-
alent subevents A and B with multiplicities nA = nB =
n/2 from an event with n particles. E.g., subevent A
has azimuth vector ~QmA =
i∈A ~u(mφi). The scalar
product is a covariance
~̃Qma · ~̃Qmb ≡ Q̃a Q̃b〈cos(m[Ψa −Ψb])〉 (D1)
na,nb
i∈A,j∈B
cos(m[φi − φj ])
≡ na nb ṽ2mab = Ṽ 2mab,
with e.g. Q̃2a = na + na(na − 1)ṽ2ma = na + Ṽ 2ma. Then
cos(m[Ψma −Ψmb]) =
Ṽ 2mab
na + Ṽ 2ma
nb + Ṽ
.(D2)
If subevents A and B are physically equivalent (e.g., a
random partition of the total of n particles), then
cos(m[Ψma −Ψmb]) = rab
V 2ma
n̄a + V 2ma
= cos(m[Ψma −Ψr]) cos(m[Ψmb −Ψr]),
where rab = V
V 2ma V
mb is Pearson’s normalized co-
variance between subevents A and B for the mth power-
spectrum elements. If A and B are perfectly correlated
(rab = 1) then
cos(m[Ψma −Ψmb]) = cos2(m[Ψma −Ψr]) (D4)
In general, V 2m/n̄ = (1+ rab)V
ma/n̄a, which provides the
exact relation between the EP resolution for subevents
and for composite events A + B. It is not generally cor-
rect that cos(m[Ψm −Ψr]) =
2 · cos(m[Ψma −Ψr]). In
this case
cos2(m[Ψma −Ψr]) =
n̄a − 1
n̄a + V 2ma
and V 2ma = V
m/4 for perfectly correlated subevents.
Second, we consider the most asymmetric case A =
one particle and B = n− 1 particles.
〈cos(m[φi −Ψr]) cos(m[Ψmi −Ψr])〉 = (D6)
(n− 1)ṽ2mi
n− 1 + (n− 1)(n− 2)ṽ2mi
ṽm ·
Ṽ ′m
where Q̃′2m = n− 1+ Ṽ ′2m describes a subevent with n− 1
particles. In general, the EP resolution for a full event of
n particles is given by
cos2(m[Ψm −Ψr]) ≃
nṼ 2m
(n− 1)Q̃2m
. (D7)
Measurement of the EP resolution is simply a measure-
ment of the corresponding power-spectrum element, since
V 2m/n̄ ≃
cos2(m[Ψm −Ψr])
1− cos2(m[Ψm −Ψr])
. (D8)
In [22] the approximation
〈cos(m[Ψm −Ψr])〉2 ≈
V 2m/n̄ (D9)
is given for V 2m/n̄≪ 1 or Q2m ∼ n̄.
Equal subevents, as the largest possible event parti-
tion, imply an expectation that only global (large-scale)
variables are relevant to collision dynamics (e.g., to de-
scribe thermalized events). The possibility of finer struc-
ture in momentum space is overlooked, whereas autocor-
relation studies with finer binnings and the covariances
among those bins discover detailed event structure highly
relevant to collision dynamics.
APPENDIX E: CENTRALITY ISSUES
Accurate A-A centrality determination and the cen-
trality dependence of azimuth multipoles and related pa-
rameters is critical to understanding heavy ion collisions.
We must locate b = 0 accurately in terms of measured
quantities to test theory expectations relative to hydro-
dynamics and thermalization. And we must obtain accu-
rate measurements for peripheral A-A collisions to pro-
vide a solid connection to elementary collisions.
1. Centrality measures
In [32] is described the power-law method of centrality
determination. Because the minimum-bias distribution
on participant-pair number npart/2 goes almost exactly
as (npart/2)
−3/4 the distribution on (npart/2)
1/4 is al-
most exactly uniform, as is the experimental distribution
ch , dominated by participant scaling. Those sim-
ple forms can greatly improve the accuracy of central-
ity determination, especially for peripheral and central
collisions. The cited paper gives simple expressions for
npart/2, nbin and ν relative to fraction of total cross sec-
tion.
In conventional centrality determination the minimum-
bias distribution on nch is divided into several bins rep-
resenting estimated fractions of the total cross section.
The main source of systematic error is uncertainty in
the fraction of total cross section which passes triggering
and event reconstruction. The total efficiency is typi-
cally 95%, the loss being mainly in the peripheral region,
and the most peripheral 10 or 20% bins therefore have
large systematic errors resulting in abandonment. Flow
measurements with EP estimation are also excluded from
peripheral collisions due to low event multiplicities.
In contrast, with the power-law method running inte-
grals of the Glauber parameters and nch can be brought
into asymptotic coincidence for peripheral collisions re-
gardless of the uncertainty in the total cross section. Pa-
rameter ν measures the centrality and greatly reduces
the cross-section error contribution. Centrality accuracy
< 2% on ν is thereby achievable down to N-N collisions.
That capability is essential to determine the correspon-
dence of A-A quadrupole structure in elementary colli-
sions, to test the LDL hypothesis for instance: is there
“collective” behavior in N-N collisions?
For central collisions the upper half-maximum point
on the power-law minimum-bias distribution provides a
precise determination of b = 0 on nch and therefore ν.
The b = 0 point is critical for evaluation of correlation
measures relative to Glauber parameter ǫ in the context
of hydro expectations for v2/ǫ.
2. Geometry parameters and azimuth structure
We consider the several A-A geometry parameters rel-
evant to azimuth structure. In Fig. 11 (left panel) we
plot npart/2 vs ν using the parameterization in [32].
The relation is very nonlinear. The dashed curve is
npart/2 ≃ 2ν2.57. The most peripheral quarter of the
centrality range is compressed into a small interval on
npart/2. Mean path-length ν is the natural geometry
measure for sensitive tests of departure from linear N-
N superposition, whereas important minijet correlations
(nonflow) ∝ ν are severely distorted on npart.
2 4 6
2 4 6
FIG. 11: Left panel: Participant pair number vs mean path-
length ν for 200 GeV Au-Au collisions. Because of the nonlin-
ear relation the peripheral third of collisions is compressed to
a small interval on npart/2. Right panel: Impact parameter b
vs ν. To good approximation the relation is linear over most
of the centrality range.
In Fig. 11 (right panel) we plot impact parameter b vs
ν, again using the parameterization in [32] with fractional
cross section σ/σ0 = (b/b0)
2. We note the interesting
fact that over most of the centrality range b/b0 ≃ (R −
ν)/(R − 1), with b0 ≡ 2R = 14.7 fm for Au-Au. Thus,
any anticipated trends on b are also accessible on ν with
minimal distortion.
In Fig. 12 (left panel) we show the LDL parameter
1/S dnch/dη [24] vs ν for three collision energies. The
energy dependence derives from the multiplicity factor,
which we parameterize in terms of a two-component
model [39]. Weighted cross-section area S(b/b0) (fm
is an optical Glauber parameterization from [31]. Both ν
and 1/S dnch/dη are pathlength measures. They can be
compared with the inverse Knudson number K−1n intro-
duced in [36] as a measure of collision number. The LDL
measure is based on energy-dependent physical particle
collisions, whereas ν is based on A-A geometry alone.
The relation is monotonic and almost linear. Thus, struc-
ture on one parameter should appear only slightly dis-
torted on the other.
17 GeV
62 GeV
200 GeV
2 4 6 1/S  dnch/dη
hydro
0 10 20
FIG. 12: Left panel: Correspondence between LDL parame-
ter 1/S dnch/dη and centrality measure ν for three energies.
Right panel: Theory expectations for two limiting cases at
200 GeV. The solid curve is derived from the solid curve in
Fig. 10 (upper-right panel) using the relation in the left panel.
The hatched region is typically not measured in a con-
ventional flow analysis, due to a combination of large sys-
tematic uncertainty in the centrality determination and
large biases in flow measurements due to small multi-
plicities. However, peripheral collisions provide critical
tests of flow models: e.g., how does collective behavior
(if present) emerge with increasing centrality? In this
paper we describe analysis methods which, when com-
bined with the centrality methods of [32], make all A-A
collisions accessible for accurate measurements down to
In Fig. 12 (right panel) we show v2/ǫ vs 1/S dnch/dη
for theory expectations (hatched bands) and the simula-
tion in Sec. VII C. The latter is based on a simple error
function on ν and is roughly consistent with four-particle
cumulant results at 200 GeV [26]. We observe that the
solid curve is not consistent with either the LDL trend
for peripheral collisions (the LDL slope is arbitrary) or
the hydro trend for central collisions. That provocative
result suggests that accurate analysis of azimuth corre-
lations over a broad range of energies and centralities
with the methods introduced in this paper and [32] may
produce interesting and unanticipated results.
3. Correlation measures
If the centrality dependence of azimuth structure is
to be accurately determined the correlation measure em-
ployed must have little or no multiplicity bias, includ-
ing statistical biases and irrelevant multiplicity factors
which lead to incorrect physical inferences. The quantity
ρref is the unique solution to a measurement prob-
lem subject to multiple constraints. It is the only portable
measure (density ratio) of two-particle correlations ap-
plicable to collision systems with arbitrary multiplicity.
ρref is invariant under linear superposition. If, ac-
cording to that measure, central Au-Au is different from
N-N the difference certainly indicates a unique physical
aspect of Au-Au collisions relative to N-N, exactly what
we require in a correlation measure. Conventional flow
measures do not satisfy that basic requirement.
Drawing a parallel with measures of 〈pt〉 fluctuations
we compare v2 ↔ Σpt [40]. Both are square roots of
per-pair correlation measures which tend to yield mis-
leading systematic trends (on centrality and energy) [10].
In contrast V 2m ↔ ∆Σ2pt:n [5] (the total variance dif-
ference for pt fluctuations) are integrals of two-particle
correlations (azimuth number correlations vs pt correla-
tions). The first is a measure of total azimuth corre-
lations, the second a measure of total pt variance dif-
ference [41], the integral of a two-particle distribution
relative to its reference. In a minimally-biased con-
text we then have vm =
V 2m/n(n− 1), analogous to
Σpt(CERES) =
∆Σ2pt:n/n(n− 1) p̂
t [10] (the similar-
ity of notations is coincidental) as the square roots of
per-pair measures. Both are physically misleading.
[1] R. C. Hwa, Phys. Rev. D 32, 637 (1985); U. W. Heinz
and P. F. Kolb, Nucl. Phys. A 702, 269 (2002).
[2] F. Becattini, M. Gazdzicki and J. Sollfrank, Eur. Phys.
J. C 5, 143 (1998); P. Braun-Munzinger, I. Heppe and
J. Stachel, Phys. Lett. B 465, 15 (1999).
[3] R. Stock, “Event by event analysis of ultrarelativistic
nuclear collisions: A new method to search for criti-
cal fluctuations,” Prepared for NATO Advanced Study
Workshop on Hot Hadronic Matter: Theory and Exper-
iment, Divonne-les-Bains, France, 27 Jun - 1 Jul 1994;
R. Stock, Nucl. Phys. A661, 282c (1999); H. Heiselberg,
Phys. Rep. 351, 161 (2001).
[4] M. Gaździcki, St. Mrówczyński, Z. Phys. C 54, 127
(1992).
[5] T. A. Trainor, hep-ph/0001148.
[6] N. Borghini, P. M. Dinh and J. Y. Ollitrault, Talk given
at 31st International Symposium on Multiparticle Dy-
namics (ISMD 2001), Datong, China, 1-7 Sep 2001. Pub-
http://arxiv.org/abs/hep-ph/0001148
lished in “Datong 2001, Multiparticle dynamics” 192-197,
hep-ph/0111402.
[7] A. Einstein, Ann. Phys. 17, 549, 1905; Investigations on
the Theory of Brownian Movement, ed. R. Fürth, trans-
lated by A. D. Cowper (1926, Dover reprint 1956); Ein-
stein, Collected Papers, 2, 170-82, 206-22.
[8] J. Perrin, Ann. Phys., 18, 5 - 114 (1909), translated
by Frederick Soddy (London: Taylor and Francis, 1910)
reprinted in David M. Knight, ed., Classical scientific pa-
pers: chemistry (New York: American Elsevier, 1968).
[9] J. Adams et al. (STAR Collaboration), J. Phys. G 32,
L37 (2006).
[10] J. Adams et al. (STAR Collaboration), J. Phys. G 33,
451 (2007).
[11] R. J. Porter and T. A. Trainor (STAR Collaboration), J.
Phys. Conf. Ser. 27, 98 (2005).
[12] J. Adams et al. (STAR Collaboration), Phys. Rev. C 73,
064907 (2006).
[13] A. H. Mueller, Nucl. Phys. B 572, 227 (2000); R. Baier,
A. H. Mueller, D. Schiff and D. T. Son, Phys. Lett. B
502, 51 (2001).
[14] J. D. Bjorken and S. J. Brodsky, Phys. Rev. D 1, 1416
(1970).
[15] P. Danielewicz and M. Gyulassy, Phys. Lett. B 129
(1983) 283.
[16] G. Arfken, Mathematical Methods for Physicists, 3rd ed.
(Orlando, FL) Academic Press, 1985; Blackman, R. B.
and Tukey, J. W. The Measurement of Power Spectra,
From the Point of View of Communications Engineering.
New York, Dover, 1959.
[17] H. Nyquist, Trans. AIEE 47, 617 (1928).
[18] T.A. Trainor, R. J. Porter and D. J. Prindle, J. Phys. G
31 809.
[19] P. Danielewicz and G. Odyniec, Phys. Lett. B 157, 146
(1985).
[20] J. Y. Ollitrault, Phys. Rev. D 46, 229 (1992).
[21] S. Voloshin and Y. Zhang, Z. Phys. C 70, 665 (1996).
[22] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58,
1671 (1998).
[23] S. A. Voloshin and A. M. Poskanzer, Phys. Lett. B 474,
27 (2000).
[24] H. Heiselberg and A.-M. Levy, Phys. Rev. C 59, 2716
(1999).
[25] Q. J. Liu, D. J. Prindle and T.A. Trainor, Phys. Lett. B
632, 197 (2006).
[26] J. Adams et al. (STAR Collaboration), Phys. Rev. C 72,
014904 (2005).
[27] N. Borghini, P. M. Dinh and J. Y. Ollitrault, Phys. Rev.
C 64, 054901 (2001)
[28] T. A. Trainor and D. T. Kettler, Phys. Rev. D 74, 034012
(2006).
[29] K. H. Ackermann et al., Nucl. Instrum. Meth. A 499,
624 (2003).
[30] C. Adler et al. (STAR Collaboration), Phys. Rev. C 66,
034904 (2002).
[31] P. Jacobs and G. Cooper, Remarks on the geometry of
high energy nuclear collisions, STAR note SN0402 (1999).
[32] T. A. Trainor and D. J. Prindle, hep-ph/0411217.
[33] N.Borghini, P.M.Dinh, J. Y.Ollitrault, A.M. Poskanzer
and S.A.Voloshin, Phys. Rev. C 66, 014901 (2002)
[34] A. Leonidov and D. Ostrovsky, Eur. Phys. J. C 16, 683
(2000).
[35] A. Leonidov and D. Ostrovsky, Phys. Rev. C 63, 037901
(2001).
[36] N. Borghini, Eur. Phys. J. A 29, 27 (2006).
[37] H. Stark and J. W. Woods, Probability and Random Pro-
cesses with Applications to Signal Processing, 3rd edition,
Prentice Hall (New Jersey, 2002); J. Bertoin, Lévy Pro-
cesses, Cambridge Univ. Press (New York, 1996).
[38] J. D. Jackson, Classical Electrodynamics, sixth edition,
John Wiley & Sons, Inc. (New York), 1967.
[39] D. Kharzeev and M. Nardi, Phys. Lett. B 507, 121
(2001).
[40] D. Adamová et al. (CERES Collaboration), Nucl. Phys.
A727, 97 (2003).
[41] J. Adams et al. (STAR Collaboration), Phys. Rev. C 71
064906 (2005).
[42] The term “cylindrical harmonic” is conventionally ap-
plied to Bessel functions of the first kind, assuming
cylindrical symmetry. However, the term “spherical har-
monic” does not assume spherical symmetry, and we use
the term “cylindrical harmonic” in that sense by analogy
to denote sinusoids on azimuth.
http://arxiv.org/abs/hep-ph/0111402
http://arxiv.org/abs/hep-ph/0411217
