# Measuring jet energy loss fluctuations in the quark-gluon plasma via
multiparticle correlations

Abraham Holtermann  Jacquelyn Noronha-Hostler  Anne M. Sickles  Xiaoning Wang
Department of Physics, University of Illinois, Urbana, IL 61801, USA

(August 28, 2024)

###### Abstract

The quark-gluon plasma (QGP) is a high temperature state of matter produced in
the collisions of two nuclei at relativistic energies. The properties of this
matter at short distance scales are probed using jets with high transverse
momentum ($p_{T}$) resulting from quarks and gluons scattered with large
momentum transfer in the earliest stages of the collisions. The Fourier
harmonics for anisotropies in the high transverse momentum particle yield,
$v_{n}(p_{T})$, indicate the path length dependence of jet energy loss within
the QGP. We present a framework to build off of measurements of jet energy
loss using $v_{n}(p_{T})$ by characterizing fluctuations in jet energy loss
that are currently not constrained experimentally. In this paper, we utilize a
set of multivariate moments and cumulants as new experimental observables to
measure event-by-event fluctuations in the azimuthal anisotropies of rare
probes, and compare them to the azimuthal anisotropies of soft particles.
Ultimately, these fluctuations can be used to quantify the magnitude and
fluctuations of event-by-event jet energy loss. We relate these quantities to
existing multivariate cumulant observables, highlight their unique properties,
and validate their sensitivities with a Monte Carlo simulation.

##  I Introduction

Sufficient energy densities are reached in collisions at the Relativistic
Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) to produce the
quark-gluon plasma (QGP), a state of hot nuclear matter characterized by the
deconfinement of the quarks and gluons that comprise the hadronic matter in
atomic nuclei Busza:2018rrf . This matter is describable as a nearly ideal
fluid Heinz:2013th ; Luzum:2013yya ; DerradideSouza:2015kpt  and jets
traveling through the medium are suppressed as they interact and exchange
energy with the medium Cunqueiro:2021wls . Jets arise from the hard scattering
between two partons in the early stages of a collision and propagate through
the QGP. The suppression of jets, jet quenching, has been extensively studied
Cunqueiro:2021wls ; Qin:2015srf  because the high momentum jets are sensitive
to short distance scale processes within the QGP. Understanding how short-
distance scale interactions give rise to the emergent hydrodynamic properties
of the QGP is a major focus of the experimental programs at RHIC and the LHC
Citron:2018lsq ; Achenbach:2023pba ; Arslandok:2023utm .

Fluctuations in the geometry of the individual nuclei at the instant of the
collision give rise to anisotropies in the azimuthal angular distribution of
hadrons in the final state through the hydrodynamic evolution of the QGP
Takahashi:2009na ; Alver:2010gr . The size of these anisotropies is
characterized by the Fourier coefficients, $v_{n}$, which are defined as:

$\frac{dN}{d\phi}\propto
1+2\sum_{n=1}^{\infty}\mbox{$v_{n}$}\cos\left[n\left(\phi-\mbox{$\Psi_{n}$}\right)\right]$
(1)

where $\phi$ is the azimuthal angle of a particle, and $\Psi_{n}$ defines the
$n$-th order event plane. These quantities have been measured extensively at
RHIC and the LHC STAR:2004jwm ; PHENIX:2011yyh ; STAR:2013qio ; ALICE:2016ccg
; CMS:2017xnj ; ATLAS:2018ezv . Measurements of $v_{n}$ have been essential to
extractions of the QGP shear viscosity to entropy density ratio $\eta/s$
Bernhard:2019bmu ; Nijs:2020roc ; JETSCAPE:2020shq . Additionally, in order to
understand how the collision geometry fluctuates on an event-by-event basis,
the distributions of the $v_{n}$ values have also been measured at low $p_{T}$
ATLAS:2013xzf , providing insight into fluctuations in the QGP initial state
produced in heavy ion collisions.

The Fourier coefficients $v_{n}$ are also used to describe the azimuthal
anisotropies of more restricted classes of objects, such as jets. Measurements
indicate that the magnitude of jet quenching is sensitive to the path length
of the jet through the QGP ALICE:2015efi ; ATLAS:2021ktw ; CMS:2022nsv , the
type of parton from which the jet originated ATLAS:2022fgb ; ATLAS:2023iad ,
and the structure of the parton shower ALargeIonColliderExperiment:2021mqf ;
ATLAS:2022vii . Measurements of these quantities provide information about the
magnitude of the dependence of energy loss on each quantity, but do not
provide any information on the role of the fluctuations in each energy loss
process. These fluctuations are expected to include fluctuations related to
the initial state geometry (as in the soft sector) but could also include
fluctuations in the energy loss process itself. Some discussion of observables
sensitive to these fluctuations has appeared in Ref. Betz:2016ayq . Similarly,
studies of the angular correlations of charmed mesons have been completed in
Ref. CMS:2021qqk . Most recently, in Ref. Holtermann:2023vwr , the
multiparticle correlation framework was extended and generalized with the goal
of defining experimental observables sensitive to the fluctuations in the
$v_{n}$ distributions of hard processes.

Jet quenching fluctuation measurements require very large data samples of
events with particles of interest (POI) (e.g. from jets, heavy-flavor, etc)
and particles from the entire event. The planned high luminosity data-taking
over the coming years at RHIC and the LHC Citron:2018lsq ; Achenbach:2023pba
will allow for the first measurements of these quantities. The LHC Runs 3 and
4 and ongoing detector upgrades will provide much larger data samples for such
measurements than has been previously available. Finally, the very large event
rate for data from sPHENIX will greatly enhance the ability to perform these
measurements at RHIC PHENIX:2015siv .

In this paper, we show that the framework outlined in Ref. Holtermann:2023vwr
allows for the construction of observables that can evaluate fluctuations and
correlations between azimuthal anisotropies for jets and reference particles.
We develop a toy model which can be used to calculate the explicit
sensitivities of four distinct observables to fluctuations and correlations in
soft and hard particle azimuthal anisotropies. Using these sensitivities, we
show how measurements of multiple observables can be used in tandem with our
toy model to isolate both fluctuations and correlations between azimuthal
anisotropies for particles in the hard sector and soft sector. We anticipate
that the capacity of our model and observables to evaluate fluctuations in the
azimuthal anisotropies of jets can be instrumental in both theoretical and
experimental settings to constrain and validate energy loss properties of the
QGP.

Our paper is outlined as follows. In Sec. II we illustrate how multiparticle
correlations allow for the approximation of different moments of the event-by-
event azimuthal anisotropy distributions for jets and reference particles. We
then detail how these moments can be used to construct a variety of
observables to measure correlations and fluctuations in jet and reference
azimuthal anisotropies. In Sec. III we introduce a bivariate copula toy model
to describe the joint distributions of azimuthal anisotropies for jet and
reference particles. Then, in Sec. IV we demonstrate the sensitivity of four
different observables to second order fluctuations in jet azimuthal anisotropy
and correlations between jet and reference azimuthal anisotropies under this
toy model. In Sec. V we discuss the applicability of the methods and
observables discussed in this paper to constrain fluctuations and correlations
experimentally through the joint measurement of multiple observables. Finally,
we discuss our outlook and conclusions in Sec. VI.

##  II Multiparticle Correlators

###  II.1 Azimuthal Anisotropy

The fundamental variables in the study of azimuthal anisotropy in heavy ion
collisions are $v_{n}$, the $n$th order Fourier coefficients in the angular
distribution of particles within an event:

$\frac{dN}{d\phi}\propto
1+2\sum_{n=1}^{\infty}v_{n}\cos\left[n\left(\phi-\Psi_{n}\right)\right]$ (2)

where the azimuthal angle $\phi_{i}$ for each particle $i$ in an event with
$N$ particles is measured in relation to the event plane $\Psi_{n}$; an
azimuthal angle at which the distribution of $n\phi_{i}$ approximately
symmetric and reaches a maximum Bilandzic:2020csw . When using the symbol
$v_{n}$, we specifically refer to the azimuthal anisotropy coefficients for
reference particles, i.e. all charged particles. The angular distribution
restricted only to POI can likewise be decomposed into Fourier harmonics
$v^{\prime}_{n}$:

$\frac{dN^{\prime}}{d\psi}\propto
1+2\sum_{n=1}^{\infty}v^{\prime}_{n}\cos\left[n\left(\psi-\Psi^{\prime}_{n}\right)\right]$
(3)

where $N^{\prime}$ indicates the total multiplicity of POI for a single event,
and each POI is at azimuthal angle $\psi_{i}$. The event plane is labeled
$\Psi^{\prime}_{n}$, about which the angles $n\psi_{i}$ are symmetric. A
useful quantity that encodes information about both the event planes and the
azimuthal anisotropy are event flow vectors:

$\displaystyle V_{n}$ $\displaystyle=v_{n}e^{in\Psi_{n}}$ (4) $\displaystyle
V^{\prime}_{n}$ $\displaystyle=v^{\prime}_{n}e^{in\Psi^{\prime}_{n}}$ (5)

complex numbers with magnitude $v_{n}$ and $v^{\prime}_{n}$ oriented in the
direction of the $n$th order event plane angles $\Psi_{n}$ and
$\Psi^{\prime}_{n}$.

Event-by-event fluctuations in $v_{n}$ and $v^{\prime}_{n}$ can be calculated
using multiparticle correlators. In heavy ion collisions with large
multiplicity, these correlators can be used to evaluate statistical moments of
$v_{n}$, and $v^{\prime}_{n}$ Bilandzic:2010jr ; Voloshin:2008dg . In Ref.
Holtermann:2023vwr , a formalism for measuring multiparticle correlations
between POI and reference particles was developed that allows for the
construction of event-by-event raw moments with arbitrary dependence on
$v^{\prime}_{n}$ and $v_{n}$.111In Ref. Holtermann:2023vwr  correlations are
described that can take arbitrary dependence on both POI and reference
particles from multiple harmonics: generally these correlations could evaluate
moments of $v^{\prime}_{n}$, $v_{n}$, $v^{\prime}_{m}$, $v_{m}$,
$v^{\prime}_{l}$, $v_{l},$ etc. These correlations also can be used to
construct more complex fluctuation observables as described here in Sec. II.2.
Using flow vectors, we can write correlations that involve arbitrary numbers
of POI and reference particles, and relate them to the values of
$V^{\prime}_{n}$ and $V_{n}$ as seen below,

$\left\langle
e^{in(\psi_{1}+...+\psi_{k^{\prime}}-\psi_{k^{\prime}+1}-...\psi_{m^{\prime}})}e^{in(\phi_{1}+...+\phi_{k}-\phi_{k+1}-...\phi_{m})}\right\rangle=({V^{\prime}}_{n})^{k^{\prime}}({{V^{\prime}_{n}}^{*}})^{m^{\prime}}({V_{n}})^{k}({V_{n}}^{*})^{m}$
(6)

where $k^{\prime}$ and $k$ indicate the POI and reference particles
respectively for which the angles are added within the correlation, while
$m^{\prime}$ and $m$ indicate the number of POI and reference particles with
angles subtracted in the correlation. Here, the angle brackets
$\left\langle\right\rangle$ indicate an average taken over all
$k^{\prime}+m^{\prime}+k+m$ tuples of unique particle angles within an event
in which $k^{\prime}+m^{\prime}$ particles are POI and $k+m$ particles are
reference particles. Note that we must have $k^{\prime}+k=m^{\prime}+m$,
otherwise the above quantity is isotropic, and will average to zero.

To study fluctuations in $v_{n}$, multiparticle particle correlations are
often used to evaluate moments of the event-by-event $v_{n}$ distribution. To
evaluate moments of both $v_{n}$ and $v^{\prime}_{n}$, a weighted event-by-
event average of Eq. (6) is taken,

$\left\langle\left\langle
e^{in(\psi_{1}+...+\psi_{k^{\prime}}-\psi_{k^{\prime}+1}-...\psi_{m^{\prime}})}e^{in(\phi_{1}+...+\phi_{k}-\phi_{k+1}-...\phi_{m})}\right\rangle\right\rangle=\left\langle({V^{\prime}}_{n})^{k^{\prime}}({{V^{\prime}_{n}}^{*}})^{m^{\prime}}({V_{n}})^{k}({V_{n}}^{*})^{m}\right\rangle$
---
$=\frac{\sum_{j}W_{j}\left[({V^{\prime}}_{n})^{k^{\prime}}({{V^{\prime}_{n}}^{*}})^{m^{\prime}}({V_{n}})^{k}({V_{n}}^{*})^{m}\right]_{j}}{\sum_{j}W_{j}}$
(7)

where a weighted average is expressed by double brackets, and the single set
of brackets indicates an expectation value, evaluated on the distribution of
events indexed by $j$. Here,
$\left[({V^{\prime}}_{n})^{k^{\prime}}({{V^{\prime}_{n}}^{*}})^{m^{\prime}}({V_{n}})^{k}({V_{n}}^{*})^{m}\right]_{j}$
refers to the quantity
$({V^{\prime}}_{n})^{k^{\prime}}({{V^{\prime}_{n}}^{*}})^{m^{\prime}}({V_{n}})^{k}({V_{n}}^{*})^{m}$
evaluated within event $j$. In this paper, we restrict ourselves to
multiparticle correlators that correlate a maximum of six particles total
within an event, of which a maximum of two particles are POI. This choice is
based on current measurements, which indicate a capacity to measure
correlations with two jet-like POI Bilandzic:2010jr ; ALICE:2022smy . The 11
unnormalized correlators with up to two POI are detailed in Table 1. Note that
for two POI in correlations with more than two total particles, there are two
ways to construct the correlation: the angles of the two POI are added, or
their difference is taken. Both instances are described in Table 1.

Table 1: Explicit definition of the multiparticle correlations used in this paper, expressing joint moments of the distribution of $v_{n}$ and $v^{\prime}_{n}$ as correlations between up to six particles, and up to two POI. | Two Particles | Four Particles | Six Particles
---|---|---|---
0 POI | $\left\langle\left\langle e^{in(\phi_{1}-\phi_{2})}\right\rangle\right\rangle=\left\langle v_{n}^{2}\right\rangle$ | $\left\langle\left\langle e^{in(\phi_{1}+\phi_{2}-\phi_{3}-\phi_{4})}\right\rangle\right\rangle=\left\langle v_{n}^{4}\right\rangle$ | $\left\langle\left\langle e^{in(\phi_{1}+\phi_{2}+\phi_{3}-\phi_{4}-\phi_{5}-\phi_{6})}\right\rangle\right\rangle=\left\langle v_{n}^{6}\right\rangle$
1 POI | $\left\langle\left\langle e^{in\left(\psi_{1}-\phi_{2}\right)}\right\rangle\right\rangle=\left\langle V_{n}^{\prime}V_{n}^{*}\right\rangle$ | $\left\langle\left\langle e^{in\left(\psi_{1}+\phi_{2}-\phi_{3}-\phi_{4}\right)}\right\rangle\right\rangle=\left\langle V^{\prime}_{n}V^{*}_{n}v_{n}^{2}\right\rangle$ | $\left\langle\left\langle e^{in(\psi_{1}+\phi_{2}+\phi_{3}-\phi_{4}-\phi_{5}-\phi_{6})}\right\rangle\right\rangle=\left\langle V^{\prime}_{n}V_{n}^{*}v_{n}^{4}\right\rangle$
2 POI | $\left\langle\left\langle e^{in(\psi_{1}-\psi_{2})}\right\rangle\right\rangle=\left\langle{v^{\prime}_{n}}^{2}\right\rangle$ |  | $\left\langle\left\langle e^{in(\psi_{1}+\psi_{2}-\phi_{3}-\phi_{4})}\right\rangle\right\rangle=\left\langle(V^{\prime}_{n}V_{n}^{*})^{2}\right\rangle$
---
$\left\langle\left\langle
e^{in(\psi_{1}-\phi_{2}+\psi_{3}-\phi_{4})}\right\rangle\right\rangle=\left\langle{v^{\prime}}_{n}^{2}v_{n}^{2}\right\rangle$
| $\left\langle\left\langle
e^{in(\psi_{1}+\psi_{2}+\phi_{3}-\phi_{4}-\phi_{5}-\phi_{6})}\right\rangle\right\rangle=\left\langle(V^{\prime}_{n}V_{n}^{*})^{2}v_{n}^{2}\right\rangle$
---
$\left\langle\left\langle
e^{in(\psi_{1}+\phi_{2}-\phi_{3}-\psi_{4}-\phi_{5}-\phi_{6})}\right\rangle\right\rangle=\left\langle{v^{\prime}}_{n}^{2}v_{n}^{4}\right\rangle$

As shown in Eq. (7), each correlation can be interpreted as an expectation
value for powers of $V_{n}$ and $V^{\prime}_{n}$, across an ensemble of
events, meaning that they can be considered to be raw moments of a
distribution. However, since $V_{n}$, $V^{\prime}_{n}$, and their odd powers
cannot be directly measured, we opt to use stochastic variables
$v_{n}^{2},{v^{\prime}}_{n}^{2}$, and the dot product
$V^{\prime}_{n}V^{*}_{n}$. This choice is consistent with the convention for
symmetric cumulants ($SC$), which use stochastic variables $v_{n}^{2}$ and
$v_{m}^{2}$, azimuthal anisotropies at two different harmonics
Bilandzic:2013kga .

Raw moments of a distribution often roughly correspond to the product of the
averages of their stochastic variables:

$\left\langle
X^{\nu_{x}}Y^{\nu_{y}}Z^{\nu_{z}}...\right\rangle\sim\left\langle
X\right\rangle^{\nu_{x}}\left\langle Y\right\rangle^{\nu_{y}}\left\langle
Z\right\rangle^{\nu_{z}}...$ (8)

where correlations and decorrelations between the variables above result in a
deviation of the above expression from equality. To ensure that each of these
moments are comparable, and to highlight correlations between each variable,
we can normalize each correlation by the product of the averages of its
stochastic variables

$n\left\langle
X^{\nu_{x}}Y^{\nu_{y}}Z^{\nu_{z}}...\right\rangle\equiv\frac{\left\langle
X^{\nu_{x}}Y^{\nu_{y}}Z^{\nu_{z}}...\right\rangle}{\left\langle
X\right\rangle^{\nu_{x}}\left\langle Y\right\rangle^{\nu_{y}}\left\langle
Z\right\rangle^{\nu_{z}}...}\sim 1.$ (9)

Here and below, we use the notation $n\left\langle\right\rangle$ to indicate a
normalized moment. If $n\left\langle
X^{\nu_{x}}Y^{\nu_{y}}Z^{\nu_{z}}...\right\rangle>1$, we expect significant
correlations between the stochastic variables $X,Y,$ and $Z$, while
$n\left\langle X^{\nu_{x}}Y^{\nu_{y}}Z^{\nu_{z}}...\right\rangle<1$ indicates
an anticorrelation.

The two-particle correlators are only comprised of one stochastic variable,
and normalize trivially to $n\left\langle v_{n}^{2}\right\rangle=1$. Each
four-particle correlation correlates two stochastic variables. Likewise, each
six-particle correlation correlates three stochastic variables. We demonstrate
the normalization scheme for each of the four- and six-particle correlators in
Table 2. These normalizable multiparticle correlators with dependence up to
${V^{\prime}_{n}}^{2}$ on POI angles can be used to construct a variety of
angular correlations Holtermann:2023vwr , each with their own sensitivities to
moments of the distributions of $v^{\prime}_{n}$ and $v_{n}$.

Table 2: We show normalization schemes for each four- and six-particle correlation. Note that for two POI, there are two possible normalization schemes, both of which we show here. | Four Particles | Six Particles
---|---|---
0 POI | $n\left\langle v_{n}^{4}\right\rangle=\frac{\left\langle v_{n}^{4}\right\rangle}{\left\langle v_{n}^{2}\right\rangle^{2}}$ | $n\left\langle v_{n}^{6}\right\rangle=\frac{\left\langle v_{n}^{6}\right\rangle}{\left\langle v_{n}^{2}\right\rangle^{3}}$
1 POI | $n\left\langle V^{\prime}_{n}V_{n}^{*}v_{n}^{2}\right\rangle=\frac{\left\langle V^{\prime}_{n}V_{n}^{*}v_{n}^{2}\right\rangle}{\left\langle V^{\prime}_{n}V_{n}^{*}\right\rangle\left\langle v_{n}^{2}\right\rangle}$ | $n\left\langle V^{\prime}_{n}V_{n}^{*}v_{n}^{4}\right\rangle=\frac{\left\langle V^{\prime}_{n}V_{n}^{*}v_{n}^{4}\right\rangle}{\left\langle v_{n}^{2}\right\rangle^{2}\left\langle V^{\prime}_{n}V_{n}^{*}\right\rangle}$
2 POI |  | $n\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)^{2}\right\rangle=\frac{\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)^{2}\right\rangle}{\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)\right\rangle^{2}}$
---
$n\left\langle{v^{\prime}}_{n}^{2}v_{n}^{2}\right\rangle=\frac{\left\langle{v^{\prime}}_{n}^{2}v_{n}^{2}\right\rangle}{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle}$
|
$n\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)^{2}v_{n}^{2}\right\rangle=\frac{\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)^{2}\right\rangle}{\left\langle\left(V^{\prime}_{n}V_{n}^{*}\right)\right\rangle^{2}\left\langle
v_{n}^{2}\right\rangle}$
---
$n\left\langle{v^{\prime}}_{n}^{2}v_{n}^{4}\right\rangle=\frac{\left\langle{v^{\prime}}_{n}^{2}v_{n}^{2}\right\rangle}{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle^{2}}$

###  II.2 Observables

Following Ref. Holtermann:2023vwr , we consider three types of observables to
quantify fluctuations and correlations in the stochastic variables
$\left\langle v_{n}^{2}\right\rangle,\left\langle
V^{\prime}_{n}V_{n}^{*}\right\rangle$, and
$\left\langle{v^{\prime}}_{n}^{2}\right\rangle$, using the correlations
detailed in Table 1. In this paper, we examine the sensitivity of these
quantities to fluctuations and correlations in $v^{\prime}_{n}$ and $v_{n}$ to
demonstrate the feasibility for these quantities to be used experimentally for
the same purpose. Each observable has a “key value” that aids interpretation,
which is described in Table 3. Comparing any observable with its key value
indicates whether or not the stochastic variables are correlated (the
observable exceeds the key value), uncorrelated (the observable is equivalent
to the key value), or anticorrelated (the observable’s value is lower than the
key value). The three types of observables are:

  * •

Moments: The observables we evaluate include each of the 11 raw moments
designated in Table 1, which can be combined easily to create normalized
moments, detailed in Table 2. While a raw or normalized raw moment has some
standalone interpretive value for quantifying the correlations between
variables, they are difficult to interpret and tend to have nonlinear
relationships with fluctuation measurements, meaning they are generally not
used alone in flow analyses. However, they are necessary to construct more
sophisticated fluctuation observables, such as cumulants, and have nontrivial
sensitivities to the fluctuations we seek to evaluate.

  * •

Cumulants: We evaluate cumulant expansions for each of the eight “higher
order” moments corresponding to correlators of either four or six particles in
Table 1, again with up to two particles being POI. Cumulants identify the
genuine correlation between a set of stochastic variables
$\\{x_{1},...,x_{n}\\}$ by subtracting away cumulants for every smaller subset
of stochastic variables
$\\{x_{i_{1}},...,x_{i_{m}}\\}\subset\\{x_{1},...,x_{n}\\}$. The derivation
for these cumulant expansions is described in Ref. Holtermann:2023vwr , but
for up to six-particle correlations, they all take the form of either a
symmetric cumulant or their generalization to higher order, an asymmetric
cumulant ($ASC$) as defined in Refs. Holtermann:2023vwr ; Bilandzic:2013kga ;
Bilandzic:2021rgb :

$\displaystyle SC(X,Y)=\left\langle XY\right\rangle-\left\langle
X\right\rangle\left\langle Y\right\rangle$ (10) $\displaystyle
ASC(X,Y,Z)=\left\langle XYZ\right\rangle-\left\langle
X\right\rangle\left\langle YZ\right\rangle-\left\langle
Y\right\rangle\left\langle ZW\right\rangle-\left\langle
Z\right\rangle\left\langle XY\right\rangle+2\left\langle
X\right\rangle\left\langle Y\right\rangle\left\langle Z\right\rangle$ (11)

where each stochastic variable $X,Y,Z$ can be any stochastic variable
$v_{n}^{2},{v^{\prime}}_{n}^{2},$ or $V^{\prime}_{n}V_{n}^{*}$. The cumulants
can be normalized simply as:

$\displaystyle nSC(X,Y)=\frac{SC(X,Y)}{\left\langle X\right\rangle\left\langle
Y\right\rangle}$ (12) $\displaystyle
nASC(X,Y,Z)=\frac{ASC(X,Y,Z)}{\left\langle X\right\rangle\left\langle
Y\right\rangle\left\langle Z\right\rangle}.$ (13)

Finally, we note that since these cumulants involve three or fewer stochastic
variables (since we restrict ourselves to correlations of no more than six
particles), they are equivalent to the central co-moments of these quantities.
The covariance between $X$ and $Y$ is measured by $SC(X,Y)$, and the “co-
skewness,” the correlation between deviations in the mean for $X,Y,$ and $Z$,
is measured by $ASC(X,Y,Z)$.

  * •

Comparisons of Fluctuations: We finally define four quantities ($\Gamma$ and
$\zeta$ for two and three stochastic variables each) that compare the
fluctuations in $v_{n}^{2}$ to fluctuations with ${v^{\prime}}_{n}^{2}$ and
$V^{\prime}_{n}V_{n}^{*}$. They were initially described in Ref.
Holtermann:2023vwr  and are written below as functions of normalized moments:

$\displaystyle\Gamma(X,Y;v_{n}^{2},v_{n}^{2})$ $\displaystyle=$ $\displaystyle
n\left\langle XY\right\rangle-n\left\langle v_{n}^{4}\right\rangle$ (14)
$\displaystyle\Gamma(X,Y,Z;v_{n}^{2},v_{n}^{2},v_{n}^{2})$ $\displaystyle=$
$\displaystyle n\left\langle XYZ\right\rangle-n\left\langle
v_{n}^{6}\right\rangle$ (15) $\displaystyle\zeta(X,Y;v_{n}^{2},v_{n}^{2})$
$\displaystyle=$ $\displaystyle\frac{n\left\langle
XY\right\rangle}{n\left\langle v_{n}^{2},v_{n}^{2}\right\rangle}$ (16)
$\displaystyle\zeta(X,Y,Z;v_{n}^{2},v_{n}^{2},v_{n}^{2})$ $\displaystyle=$
$\displaystyle\frac{n\left\langle XYZ\right\rangle}{n\left\langle
v_{n}^{6}\right\rangle}$ (17)

where $X$, $Y$, and $Z$ can be any stochastic variable
$v_{n}^{2},{v^{\prime}}_{n}^{2},$ or $V^{\prime}_{n}V^{*}_{n}$, so long as the
total dependence on ${v^{\prime}}_{n}^{2}\leq 1$ and the total dependence on
$V^{\prime}_{n}V_{n}^{*}\leq 2$. This constraint allows for 12 unique
observables which compare the normalized raw moments for a set of arbitrary
stochastic variables to normalized raw moments of the same order in
$v_{n}^{2}$. As detailed in Table 2, values larger than 0 displayed by
$\Gamma$ indicate more significant fluctuations displayed by the random
variables $X,Y,$ and $Z$ than are displayed by $v_{n}^{2}$. Negative values
indicate that $v_{n}^{2}$ displays larger fluctuations than are displayed by
$X,Y,$ and $Z$. For $\zeta$ a similar relationship exists: greater
fluctuations are displayed by $X,Y,$ and $Z$ than by $v_{n}^{2}$ when
$\zeta>1$, and $v_{n}^{2}$ displays larger fluctuations than $X,Y$ and $Z$
when $\zeta<1$. Additionally, while $\zeta$ is nonlinear and describes
differences in the overall magnitude of fluctuations, the difference in
central moments between two sets of stochastic variables can be written as a
sum of different $\Gamma$ observables. $\Gamma$ and $\zeta$ illustrate the
differences and ratios between the four- and six-particle contributions to
fluctuations in $v^{\prime}_{n}$ and $v_{n}$. For this paper, we have
exclusively focused on correlations within one harmonic, and no more than six
particles.

Table 3: We display a list of key values and their interpretations for each observable discussed in this section. The stochastic variables $X,Y,$ and $Z$ are arbitrary.  Observable |  Key Value |  Explanation
---|---|---
$n\left\langle XY\right\rangle,n\left\langle XYZ\right\rangle$ |  1 |  Normalized moments with a value greater than unity indicate that $X,Y,$ and $Z$ are positively correlated, while values of less than unity indicate anticorrelation. A value of exactly one indicates that $X,Y,$ and $Z$ are independent.
nSC$(X,Y),$ nASC$(X,Y,Z)$, |  0 |  Both normalized and unnormalized symmetric cumulants or asymmetric cumulants with values greater than 0 indicate genuine positive correlations between $X$ and $Y$, or $X,Y,$ and $Z$ respectively. A value of 0 indicates independence between at least two of $X,Y$ or $Z$, and a negative value for a symmetric cumulant or asymmetric cumulant indicates anticorrelation between $X$ and $Y$, or $X,Y,$ and $Z$.
$\Gamma(X,Y;v_{n}^{2},v_{n}^{2})$ $\Gamma(X,Y,Z;v_{n}^{2},v_{n}^{2},v_{n}^{2})$ |  0 |  Positive values of $\Gamma$ indicate that the normalized raw moments $n\left\langle XY\right\rangle>n\left\langle v_{n}^{4}\right\rangle$ or $n\left\langle XYZ\right\rangle>n\left\langle v_{n}^{6}\right\rangle$, meaning that the fluctuations between $X$ and $Y$, or $X,Y,$ and $Z$ are larger than the second or third order fluctuations between $v_{n}^{2}$. $\Gamma=0$ means that the fluctuations have no difference, and $\Gamma<0$ indicates that the second or third order fluctuations in $v_{n}^{2}$ are greater than the fluctuations displayed by $X$ and $Y$, or $X,Y,$ and $Z$.
$\zeta(X,Y;v_{n}^{2},v_{n}^{2})$, $\zeta(X,Y,Z;v_{n}^{2},v_{n}^{2},v_{n}^{2})$ |  1 |  Values of $\zeta$ greater than unity indicate that the overall magnitude of the fluctuations in $X$ and $Y$, or $X,Y,$ and $Z$ are larger than the second or third order fluctuations displayed by $v_{n}^{2}$. Likewise, for $\zeta<1$ we have that the fluctuations in $v_{n}^{2}$ are larger in magnitude than the fluctuations between $X$ and $Y$, or $X,Y,$ and $Z$

##  III Modeling Azimuthal Anisotropy

###  III.1 Parametrizing Azimuthal Anisotropies

Using multiparticle correlators, we have identified a way to evaluate event-
by-event moments of stochastic variables with varying dependence on
$v^{\prime}_{n}$ and $v_{n}$ Holtermann:2023vwr , and construct fluctuation
observables with them. We will show that these fluctuation observables can
evaluate correlations and fluctuations in the actual values of
$v^{\prime}_{n}$ and $v_{n}$. To do this, we parametrize $v^{\prime}_{n}$ and
$v_{n}$ using a bivariate distribution, and by sampling these bivariate
distributions, demonstrate how various observables can be used to isolate
parameters related to the correlation between $v^{\prime}_{n}$ and $v_{n}$,
and the fluctuations within $v^{\prime}_{n}$.

Our approach is not new for the study of fluctuations in the soft sector; in
Refs. Voloshin:2007pc ; Voloshin:2008dg , the authors have used assumptions
about initial state geometry to create ansatzes for the $v_{n}$ distribution,
and in Refs. ATLAS:2013xzf ; CMS:2017glf  cumulants and unfolding methods are
used to approximate the $v_{n}$ distribution $p(v_{n})$, and then solve for
the parameter values used in parametrizations of $p(v_{n})$. We extend this
approach to the study of fluctuations in $v^{\prime}_{n}$ as well as
correlations between $v^{\prime}_{n}$ and $v_{n}$, using the cumulants and
correlations defined in the previous section.

The parametrizations used in our analysis are validated by their ability to
represent simulated values of eccentricity using MCGlauber Loizides:2017ack ,
a heavy ion initial state model. Eccentricity is defined similarly to $v_{n}$
as the coefficients in the Fourier expansion of the energy density
$\rho(r,\varphi)$ of the initial state of a heavy ion collision:

$\varepsilon_{n}e^{in\psi_{n}}\equiv-\frac{\int
r^{n}e^{in\varphi}\rho(r,\varphi)r\mathrm{\leavevmode\nobreak\
d}r\mathrm{\leavevmode\nobreak\ d}\varphi}{\int
r^{n}\rho(r,\varphi)r\mathrm{\leavevmode\nobreak\
d}r\mathrm{\leavevmode\nobreak\ d}\varphi}$ (18)

where $\varphi$ is the azimuthal angle, and $\psi_{n}$ is an $n$th order event
plane. In Refs. Noronha-Hostler:2015dbi ; Teaney:2010vd ; Gardim:2011xv ;
Niemi:2012aj ; Teaney:2012ke ; Qiu:2011iv ; Gardim:2014tya , it is shown that
$v_{n}\propto\varepsilon_{n}$ is a reasonable approximation for $v_{2}$ and
$v_{3}$ in central to mid-central collisions. Thus, a distribution that
accurately reflects MCGlauber simulations also reflects $p(v_{n})$ for $n=2$
or $3$ in central to midcentral collisions.222As is noted in Refs.
Renk:2014jja ; Denicol:2014ywa , MCGlauber cannot fully describe the
eccentricity and $v_{n}$ distributions measured from the LHC. However, for the
purposes of this analysis, we claim that using MCGlauber as a proxy for
$p(v_{n})$ is adequate for constraining simple quantities like the mean and
standard deviations of $p(v_{n})$.

For our analysis, we consider the Gaussian distribution and the elliptic power
distribution. The Gaussian distribution lends itself naturally to this
analysis because it easily generalizes into a bivariate distribution with an
extra parameter corresponding to the Pearson correlation coefficient $\rho$.
Additionally, its fluctuations are typically parameterized by $\sigma$, which,
when squared, is the exact value of the variance, $\sigma^{2}$, of the
distribution. Under a Gaussian $\varepsilon_{n}$ distribution parametrized by
the mean $\mu$ and the standard deviation $\sigma$,

$p(\varepsilon_{n})=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{1}{2}\left(\frac{\varepsilon_{n}-\mu}{\sigma}\right)^{2}}$
(19)

the following equations represent the second and fourth order moments of
$v_{n}$, using

$\displaystyle\left\langle\varepsilon_{n}^{2}\right\rangle=\mu^{2}+\sigma^{2}$
(20)
$\displaystyle\left\langle\varepsilon_{n}^{4}\right\rangle=\mu^{4}+6\mu^{2}\sigma^{2}+3\sigma^{4}$

$\left\langle(kX)^{n}\right\rangle=k^{n}\left\langle X^{n}\right\rangle$. In
general, a measurement of both of the above moments can be used to solve for
$\mu$ and $\sigma$, which will provide information about event-by-event
fluctuations in the $v_{n}$ distribution, since $\sigma^{2}(v_{n})$, can be
isolated.

The elliptic power distribution was explicitly constructed to represent the
eccentricity distributions generated by MCGlauber Yan:2014afa . The elliptic
power distribution with parameters $\varepsilon_{0}$, related to the mean
$\left\langle\varepsilon_{n}\right\rangle$, and $\alpha$, related to the
spread of $p(\varepsilon_{n})$ is written as:

$p(\varepsilon_{n})=2\varepsilon_{n}\alpha\left(1-\varepsilon_{n}^{2}\right)^{\alpha-1}\left(1-\varepsilon_{n}\varepsilon_{0}\right)^{-1-2\alpha}\left(1-\varepsilon_{0}^{2}\right)^{\alpha+\frac{1}{2}}{}_{2}F_{1}\left(\frac{1}{2},1+2\alpha;1;\frac{2\varepsilon_{n}\varepsilon_{0}}{\varepsilon_{n}\varepsilon_{0}-1}\right)$
(21)

where ${}_{2}F_{1}(a,b;c;z)$ is the Gaussian hypergeometric function. The
moments of this distribution are only calculable numerically, but they can be
approximated, as shown in Ref. Yan:2014afa

$\displaystyle\left\langle\varepsilon_{n}^{2}\right\rangle\approx\varepsilon_{0}^{2}+\frac{\left(1-\varepsilon_{0}^{2}\right)\left(1-\frac{3}{2}\varepsilon_{0}^{2}\right)}{\alpha}+\mathcal{O}\left(\frac{1}{\alpha^{2}}\right)$
(22)
$\displaystyle\left\langle\varepsilon_{n}^{4}\right\rangle\approx\varepsilon_{0}^{4}+\frac{\varepsilon_{0}^{2}\left(1-\varepsilon_{0}^{2}\right)\left(4-5\varepsilon_{0}^{2}\right)}{\alpha}+\mathcal{O}\left(\frac{1}{\alpha^{2}}\right).$

By numerically or analytically estimating the various moments of the elliptic
power distribution, we can solve for $\varepsilon_{0}$ and $\alpha$ in the
same way as $\mu$ and $\sigma$ were determined for the Gaussian distribution.
Since the relationship between $v_{n}$ and $\varepsilon_{n}$ is approximately
linear in central to mid-central collisions, a probability distribution
$p(\varepsilon_{n})$ can also describe the distribution $p(v_{n})$. Moreover,
due to cancellations, each observable using normalized correlators can be
calculated using $p(v_{n})$ just as well as with $p(\varepsilon_{n})$. To
indicate the capacity of both the elliptic power and Gaussian parametrizations
to characterize realistic event-by-event distributions of $v_{n}$, sampled
values for $\varepsilon_{n}$ drawn from both elliptic power and Gaussian
distributions are shown in Fig. 1. The fit between MCGlauber $\varepsilon_{2}$
simulated values obtained from Ref. Loizides:2017ack , and $\varepsilon_{2}$
from the elliptic power and Gaussian distributions is demonstrated for two
centrality classes: 5-10% and 40-45%. In each plot, the correspondence between
the MCGlauber simulation results and the elliptic power and Gaussian
distributions is evident, and both parametrizations are shown to accurately
characterize significant differences in the fluctuations of $\varepsilon_{n}$
between the central and peripheral MCGlauber results. Note that the Gaussian
fit displays nonzero probability for $v^{\prime}_{n}$ to have values less than
zero. We acknowledge that this is unrealistic. However, it poses no major
issues for our analysis, and does not significantly affect our results, as
seen by the correspondences in results obtained for both the elliptic power
distribution and Gaussian distributions.

|
---|---
Figure 1:  Marginal distributions sampled from elliptic power (green) and
Gaussian (red) parametrizations are plotted next to MCGlauber eccentricity
simulation results at 5-10% (left) and 40-45% (right) centrality bins.

###  III.2 Bivariate flow models

Up until now, we have concentrated on the distribution of $v_{n}$ values,
$p(v_{n})$. This distribution is already well studied ATLAS:2013xzf ;
CMS:2017glf , both by approximating the underlying distribution and measuring
cumulants. In this work, we show that this approach can be extended to the
consideration of bivariate flow models that relate $v^{\prime}_{n}$ and
$v_{n}$. The distribution $p(v^{\prime}_{n})$ is not currently known, so we
introduce a method to construct bivariate azimuthal anisotropy distributions,
$p(v_{n},v^{\prime}_{n})$ with the same marginal parametrizations for
$p(v^{\prime}_{n})$ (although allowing different values and fluctuations) as
$p(v_{n})$, displaying correlations of varying magnitudes.

In a previous analysis Yi:2011hs , differential flow $v^{\prime}_{n}$ was
evaluated using a Gaussian distribution for both $p(v_{n})$ and
$p(v^{\prime}_{n})$, and assuming a perfect linear correlation between $v_{n}$
and $v^{\prime}_{n}$ to generate $v^{\prime}_{n}$ samples by sampling
$p(v_{n})$ Yi:2011hs . While a reasonable approximation, recent measurements
from ALICE of $p_{\textrm{T}}$ dependent fluctuations of both the magnitude
and the angle associated with $v_{2}$ and $v^{\prime}_{2}$ in PbPb collisions
ALICE:2022smy  demonstrate that the assumption of perfect correlation is not
accurate. In light of this, allowing the correlation coefficient between
$v^{\prime}_{n}$ and $v_{n}$ to vary away from unity in our model is necessary
to capture realistic decorrelation effects for POI and reference particles.

In this analysis we use the bivariate Gaussian copula Takeuchi_2010 ;
Meyer_2013  to model the bivariate distribution of $v^{\prime}_{n}$ and
$v_{n}$ with a fixed correlation coefficient $\rho$ and marginal distributions
$p(v_{n})$ and $p({v^{\prime}}_{n})$. Specifically, we examine two cases:

  * •

Bivariate Gaussian Distribution: The joint probability
$p(v_{n},v^{\prime}_{n})$ is assumed to have Gaussian marginal distributions,
where the standard deviation and mean of the reference $v_{n}$ distribution
are written as $\sigma(v_{n})$ and $\mu(v_{n})$, and the the same quantities
for the $v^{\prime}_{n}$ distribution are $\sigma(v^{\prime}_{n})$ and
$\mu(v^{\prime}_{n})$.

  * •

Joint Elliptic Power Distribution: The joint probability
$p(v_{n},v^{\prime}_{n})$ is assumed to have elliptic power marginal
distributions. In this case, the $\varepsilon_{0}$ and $\alpha$ parameters for
reference flow are written as $\varepsilon_{0}(v_{n})$ and $\alpha(v_{n})$,
and for differential flow, they are written $\varepsilon_{0}(v^{\prime}_{n})$
and $\alpha(v^{\prime}_{n})$.

Using the cumulative distribution function for the Gaussian distribution
denoted by $\Phi(x)$:

$\Phi(X)=\int_{-\infty}^{X}\frac{1}{\sqrt{2\pi}}\exp(\frac{-x^{2}}{2})dx$ (23)

we can define the bivariate Gaussian copula density. For a given $\rho$, the
bivariate Gaussian copula density $c(v_{n},v^{\prime}_{n},\rho)$ is given by
Meyer_2013 ; Takeuchi_2010 :

$c(v_{n},v^{\prime}_{n},\rho)=\frac{1}{\sqrt{1-\rho^{2}}}\exp\left(\frac{2\rho\Phi^{-1}(F_{1}(v_{n}))\Phi^{-1}(F_{2}(v^{\prime}_{n}))-\rho^{2}\left(\Phi^{-1}(F_{1}(v_{n}))^{2}+\Phi^{-1}(F_{2}(v^{\prime}_{n}))^{2}\right)}{2\left(1-\rho^{2}\right)}\right)$
(24)

where $F_{1}(v_{n})$ and $F_{2}(v^{\prime}_{n})$ are the respective cumulative
distribution functions for the desired marginal distributions of $v_{n}$ and
$v^{\prime}_{n}$. The distribution $c(v_{n},v^{\prime}n,\rho)$ is a bivariate
probability distribution with marginal distributions characterized by
$F_{1}(v_{n})$ and $F_{2}(v_{n})$ and Pearson correlation correlation
coefficient $\rho$, that characterizes the probability density for any pair
$v_{n}$, $v^{\prime}_{n}$ by their quantile within their marginal
distributions, and the selected correlation coefficient. Note that if the
marginal distributions are Gaussians, Eq. (24) reduces to a bivariate Gaussian
distribution with correlation parameter $\rho$. Fig. 2 contains example joint
distributions for $\varepsilon_{n}$ and $\varepsilon^{\prime}_{n}$ with
Gaussian and elliptic power marginal distributions using the same parameter
values illustrated in Fig. 1 for $p(v_{n})$. Both figures have similar mean
values and ranges for $v_{n}$ and $v^{\prime}_{n}$. Additionally they
illustrate the same values for $\rho$. However, it is clear that constraints
imposed by the marginal distributions force qualitative differences between
the two distributions, especially for $v_{n},v^{\prime}_{n}\rightarrow 0$.

Figure 2: Example joint distributions for $\varepsilon_{n}$ and
$\varepsilon^{\prime}_{n}$ for the bivariate Gaussian distribution (left) and
a joint elliptic power distribution linked with a Gaussian copula (right). The
color indicates the magnitude of the probability density for both
distributions.

##  IV Results

###  IV.1 Toy Model

In the previous section we have detailed how to construct bivariate
distributions $p(v_{n},v^{\prime}_{n})$ using a bivariate copula model. In
this section we motivate and describe a toy model using the previously
discussed joint distributions to evaluate the sensitivities of the observables
developed in Sec. II.2 to fluctuations and correlations in $v_{n}$ and
$v^{\prime}_{n}$.

First, given a distribution $p(v_{n},v^{\prime}_{n})$, each correlator can be
measured by evaluating products of sampled points: $\left\langle
v_{n}^{2}V^{\prime}_{n}V_{n}\right\rangle$ can be evaluated by simply sampling
$p(v_{n},v^{\prime}_{n})$ and taking the average value of the product
$v_{n}^{3}$ and $v^{\prime}_{n}$ for each sampled
$(v_{n},v^{\prime}_{n})$.333This method relies on the assumption that there is
no angular decorrelation between the $n$th order event planes for
$\Psi^{\prime}_{n}$ and $\Psi_{n}$. Not only is this a good approximation at
low $p_{T}$ ALICE:2022smy , but we can also allow the decorrelations described
by $\rho$ in the copula models to describe any desired combination of angular
and magnitude decorrelations resulting in a weaker correlation between
measured values of the product $v^{\prime}_{n}v_{n}$. To evaluate the
sensitivity of each observable to fluctuations in $v^{\prime}_{n}$ and
$v_{n}$, we examine the value of each observable for an array of different
parameter values for $p(v_{n},v^{\prime}_{n})$. To determine which parameters
are changed, we consider that both the bivariate Gaussian and joint elliptic
power distributions have five parameters: two parameters define $p(v_{n})$,
two define $p(v^{\prime}_{n})$, and $\rho$ describes the correlation between
$v^{\prime}_{n}$ and $v_{n}$. Cumulant analyses as described in Ref.
Voloshin:2007pc  can constrain both parameters for $p(v_{n})$, and
differential cumulant studies, as well as scalar product methods Luzum:2012da
; Bilandzic:2010jr  with only one POI can approximate $\left\langle
v^{\prime}_{n}\right\rangle$, which is typically described by one of the
parameters for $p(v^{\prime}_{n})$. This leaves the remaining parameter (the
one associated with fluctuations either in $\alpha$ or $\sigma$) for
$p(v^{\prime}_{n})$, and $\rho$ to be varied when evaluating the observables
in Sec. II.2.

Table 4: Parameter values for the bivariate Gaussian and joint elliptic power distributions are displayed. Note that the top three rows are fixed, and the bottom two parameters listed for both distributions are incremented between a range of values. Bivariate Gaussian | Joint Elliptic Power
---|---
$\mu(v_{n})=0.3$ | $\varepsilon_{0}(v_{n})=0.28$
$\sigma(v_{n})=0.11$ | $\alpha(v_{n})=27$
$\mu(v^{\prime}_{n})=0.08$ | $\varepsilon_{0}(v^{\prime}_{n})=0.03$
$\sigma(v^{\prime}_{n})=0.8,\cdots,2.0$ | $\alpha(v^{\prime}_{n})=3,\cdots,500$
$\rho=0.0,\cdots,1.0$ | $\rho=0.0,\cdots,1.0$

Specifically, for the joint elliptic power parametrization, we fix
$\varepsilon(v_{n}),$ $\varepsilon(v^{\prime}_{n})$, and $\alpha(v_{n})$,
using values obtained for the $40-45\%$ centrality class in Fig. 1 before
iterating through $\rho$ and $\alpha(v^{\prime}_{n})$. Likewise, for the
bivariate Gaussian parametrization, we fix $\mu(v_{n})$, $\sigma(v_{n})$, and
$\mu(v^{\prime}_{n})$, and iterate through values for $\rho$ and
$\sigma(v^{\prime}_{n})$. The values for each parameter are found in Table 4.

Our goal is to determine the sensitivity of each observable to fluctuations in
$v^{\prime}_{n}$ and correlations between $v^{\prime}_{n}$ and $v_{n}$, and
show that through the measurement of these observables, correlations between
$v^{\prime}_{n}$ and $v_{n}$, and fluctuations in $v^{\prime}_{n}$ can be
quantified and constrained. To characterize the dependence of each observable
on these fluctuations, we require measurements of fluctuations in
$v^{\prime}_{n}$, and correlations between $v^{\prime}_{n}$ and $v_{n}$. We
select $\rho$ as our proxy for correlation between $v^{\prime}_{n}$ and
$v_{n}$ because it is among the most commonly used and simplest measurements
for correlation between two variables. Additionally, our choice of copula in
Eq. (24) enables for $\rho$ to be selected as a direct input regardless of
parametrization. To evaluate the sensitivity of observables to fluctuations in
$v^{\prime}_{n}$ we use the relative spread RS:

$\mathrm{RS}\equiv\frac{\sigma(v^{\prime}_{n})/\mu(v^{\prime}_{n})}{\sigma(v_{n})/\mu(v_{n})}$
(25)

which quantifies the fluctuations in $v^{\prime}_{n}$ relative to fluctuations
within $v_{n}$ by taking the ratio of their standard deviations divided by
their means. This quantity is useful since it is unitless and normalized, and
evidently describes fluctuations in $v^{\prime}_{n}$ and $v_{n}$ relative to
their mean.

Figure 3: A cartoon illustrating how RS and $\rho$ define a two dimensional
space, over which the sensitivity of observables can be evaluated. Example
distributions of bivariate Gaussian data distributed according to $\rho$ and
RS values are depicted.

RS and $\rho$ are independent quantities, and define a two dimensional space
occupied by all bivariate density parametrizations. By evaluating each
observable for bivariate distributions at a variety of locations within this
space, the sensitivity of each observable to RS and $\rho$ can be identified.
Fig. 3 details this notion, showing example bivariate Gaussian distributions
arranged in RS, $\rho$ space. Along the RS axis, the bivariate Gaussian
distribution displays increasing values of $\sigma(v^{\prime}_{n})$, and the
elliptic power distribution displays increasing values for
$\alpha(v^{\prime}_{n})$. Along the $\rho$ axis, $\rho$ changes for both
parametrizations. Examples of bivariate Gaussian data for different values of
RS and $\rho$ are plotted in the corners, illustrating the qualitative
features of the relationship between the bivariate Gaussian distribution, RS
and $\rho$, as shown in Fig. 3.

Figure 4: The RS values as defined in Eq. (25) for the joint elliptic power
distribution as a function of $\alpha(v^{\prime}_{n})$ are displayed. The
dependence is nonlinear, and is non-monotonic for $\alpha(v^{\prime}_{n})\gg
1$. Our selection of values for $\alpha(v^{\prime}_{n})$ spans the whole range
of RS accessible for this parametrization.

While RS is easily varied for the bivariate Gaussian distribution, we note
that it is more complicated to understand RS in the context of parameters in
the elliptic power distribution. Fig. 4 shows the correspondence between RS
and $\alpha(v^{\prime}_{n})$, where the RS is calculated for the elliptic
power distribution at values of $\alpha(v^{\prime}_{n})$ varying from 2 to
500. The fixed values for parameters $\epsilon_{0}(v_{n})=0.28,$
$\epsilon_{0}(v^{\prime}_{n})=0.03,$ and $\alpha(v_{n})=27$ are selected
because they produce the best fit to MCGlauber simulation results for 40-45%
central collisions. From this plot, it is clear that our sampling of
$\alpha(v^{\prime}_{n})$ covers the entire range of RS that the joint elliptic
power distribution can display with $\epsilon_{0}(v_{n})=0.28,$
$\epsilon_{0}(v^{\prime}_{n})=0.03,$ and $\alpha(v_{n})=27$.

The code for this toy model is in Ref. citation-key . To run it, a user
specifies the desired bivariate flow model, and initial values for the fixed
parameters. Then, the code iterates along a defined set of parameter values,
at each junction constructing a new bivariate Gaussian or joint elliptic power
distribution. The code then samples the newly created distribution ($\sim
10^{6}$ times for the results displayed here), and evaluates $RS$ and $\rho$
by calculating
$\rho(v^{\prime}_{n},v_{n}),\sigma(v^{\prime}_{n}),\sigma(v_{n}),\mu(v_{n})$,
and $\mu(v^{\prime}_{n})$ from sampled data. The code then evaluates each
correlation detailed in Sec. II.2 by averaging products of $v_{n}$ and
$v^{\prime}_{n}$. Finally, it evaluates each of the 31 observables using these
correlations, and records the value of each observable at its place in the two
dimensional space defined by RS and $\rho$.

###  IV.2 Sensitivity of Four Observables

Of the observables initially proposed, we present detailed analysis and
interpretation of four specific observables that showcase the various types of
observable, choices for random variable, and number of particles:
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$,
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$,
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2}),$ and
$nASC({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$.

We present $\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$
because it requires the least statistical precision, and fewest POIs to obtain
a measurement. Additionally, predictions for this quantity are evaluated in
Ref. Betz:2016ayq , which allows extra interpretive context, indicating
different constraints on $\rho$ and RS imposed by various hydrodynamical
observables coupled to jet energy loss models. To complement $\Gamma$, which
provides a difference between fluctuations with a quantity sensitive to
relative fluctuations from four particles, we select
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$.
This observable has a higher order of dependence on $v^{\prime}_{n}$, but
correlates the same stochastic variables ($v_{n}^{2}$ and
$V^{\prime}_{n}V_{n}^{*}$) as were shown for $\Gamma$.

As referenced in Sec. II, two-POI correlations have two distinct
interpretations: they can correlate $V^{\prime}_{n}V_{n}^{*}$ and $v_{n}^{2}$,
or they can correlate ${v^{\prime}}_{n}^{2}$ and $v_{n}^{2}$. To compare with
our study of $V^{\prime}_{n}V_{n}^{*}$ fluctuations using
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ and
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$,
we use nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ to examine differences in
sensitivity to RS and $\rho$ displayed by the stochastic variable
${v^{\prime}}_{n}^{2}$. Finally, to build off our choice for
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$, we consider a higher order cumulant
using the same variables: nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$, to
illustrate the viability of higher order cumulants in this analysis, and
examine the extent to which higher order cumulants can actually isolate $\rho$
and RS under our toy model.

We evaluate each observable by sampling joint elliptic power distributions and
bivariate Gaussian distributions that display varying values of RS,
representing the fluctuations in $v^{\prime}_{n}$ relative to $v_{n}$, and
$\rho$, representing the correlation between $v^{\prime}_{n}$ and $v_{n}$. In
Figs. 5 \- 8, we display the sensitivity of $\Gamma,\zeta,nSC,$ and nASC, to
RS and $\rho$ using the joint elliptic power and bivariate Gaussian
distributions. RS and $\rho$ values are represented for each distribution on
the horizontal and vertical axes respectively, while the value of the
observable is indicated by color, delineated on a third axis. The level curves
(bands within RS, $\rho$ space that are the same color) indicate regions of
RS, $\rho$ space over which the observable maintains approximately the same
value. In the following text we discuss results for each observable in detail.

####  IV.2.1  $\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$

$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ was originally
proposed in Ref. Betz:2016ayq  (referred to there as $\Delta_{n}^{sh}$) with
the intent of evaluating a difference in fluctuations between $v^{\prime}_{n}$
and $v_{n}$. Unfortunately, it does not suffice in this goal on its own, as it
demonstrates sensitivity to both $\rho$ and to RS, which is apparent for both
of the Gaussian and elliptic power parametrizations seen in Fig. 5. However,
it does still have significant standalone importance. The authors of Ref.
Betz:2016ayq  produced predictions for this quantity using a jet quenching
hydrodynamical model, which can be used to constrain both $\rho$ and RS.

Fig. 5 shows explicitly the sensitivity of
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ as defined
below,

$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n};v_{n}^{2},v_{n}^{2})=\frac{\left\langle
v_{n}^{2}V^{\prime}_{n}V_{n}^{*}\right\rangle}{\left\langle
V^{\prime}_{n}V_{n}\right\rangle\left\langle
v_{n}^{2}\right\rangle}-\frac{\left\langle
v_{n}^{4}\right\rangle}{\left\langle
v_{n}^{2}\right\rangle^{2}}=nSC(V^{\prime}_{n}V_{n},v_{n}^{2})-nSC(v_{n}^{2},v_{n}^{2})$
(26)

to fluctuations in $\rho$ and RS. Since
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ can be written
as a difference of symmetric cumulants, we can see that positive and negative
values for $\Gamma$ show whether $cov(V^{\prime}_{n}V_{n}^{*},v_{n}^{2})$ is
larger or smaller than $\sigma^{2}(v_{n}^{2})$. In this way, $\Gamma$ is a
measurement of the magnitude of fluctuations in $V^{\prime}_{n}V_{n}^{*}$
relative to $v_{n}^{2}$. $\Gamma$ also has considerable importance because it
only has a first order dependence on $v^{\prime}_{n}$, making it more feasible
to measure with jets or other rare probes, as it has the same statistical
requirements as the correlations used to create differential cumulants
Borghini:2001vi .

|
---|---
Figure 5: The sensitivity of
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ is shown as a
function of $\rho(v_{n},v^{\prime}_{n})$ and RS
$=\sigma(v^{\prime}_{n})\mu(v_{n})/\sigma(v_{n})\mu(v^{\prime}_{n})$ for
$v^{\prime}_{n}$ and $v_{n}$ parametrized using the joint elliptic power
distribution (left), and a bivariate Gaussian distribution (right). Note that
the range of RS is different in the two plots.

The general behavior of
$\Gamma(v_{n}^{2},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$ is desirable;
it is largely monotonic in both $\rho$ and RS, as shown in Fig. 5, and
demonstrates comparable values for both parametrizations. While there is a
clear difference between the right and left panels, we see that the values and
behavior of the observable are similar when considering that the RS for the
elliptic power distribution varies over a much smaller range of values
compared to RS for the bivariate Gaussian. The level of correspondence between
both $p(v_{n})$ parametrizations is discussed in depth in Sec. V.

####  IV.2.2
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$

The quantity
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$
can be understood as the ratio of the non-constant terms in two normalized
symmetric cumulants. While $\zeta$ can not in general express the differences
between central moments, it describes the relative magnitude of the
fluctuations between $V^{\prime}_{n}V_{n}^{*}$ and $v_{n}^{2}$, as detailed
below,

$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})=\frac{\left\langle(V^{\prime}_{n}V_{n}^{*})^{2}\right\rangle}{\left\langle
V^{\prime}_{n}V_{n}^{*}\right\rangle^{2}}\bigg{/}\frac{\left\langle
v_{n}^{4}\right\rangle}{\left\langle v_{n}^{2}\right\rangle^{2}}$ (27)

where a higher order dependence on $v^{\prime}_{n}$ can be seen for this
quantity, allowing for a different sensitivity to both RS and $\rho$.

|
---|---
Figure 6: The sensitivity of
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*};v_{n}^{2},v_{n}^{2})$
is shown as a function of $\rho(v_{n},v^{\prime}_{n})$ and RS for
$v^{\prime}_{n}$ and $v_{n}$ parametrized using the joint elliptic power
distribution (left), and a bivariate Gaussian distribution (right). Note that
the range of RS is different in the two plots.

Fig. 6 shows the sensitivities displayed by $\zeta$ to both $\rho$ and RS. We
find that $\Gamma$ strongly varies with RS, but is less sensitive to $\rho$.
This behavior is quite different compared to $\Gamma$, which increases more
with $\rho$ and is nearly flat in RS. Additionally, the range of values
covered by $\zeta$ is larger than that of $\Gamma$. This comes from a
comparatively stronger second order dependence on $v^{\prime}_{n}$, as well as
from $\zeta$ expressing a ratio rather than a difference in fluctuations. Fig.
6 demonstrates that $\zeta$, like $\Gamma$, monotonically increases along with
both RS and $\rho$, but with very different sensitivities. The difference in
sensitivities between $\Gamma$ and $\zeta$ ensures that there is little
redundancy in simultaneous measurements, and, as described in Sec. V, may
assist in the discernment of values for both $\rho$ and RS.

####  IV.2.3 nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$

Symmetric cumulants have been used to probe flow coefficient magnitude
correlations Bilandzic:2013kga ; Giacalone:2016afq ; ATLAS:2018ngv , and their
interpretation is very straightforward: larger values indicate a stronger
correlation between ${v^{\prime}}_{n}^{2}$ and $v_{n}^{2}$. Since covariance
is a co-central moment, we can also interpret that the value of
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ indicates the correlation between a
departure from the mean in ${v^{\prime}}_{n}^{2}$ with a departure from the
mean in $v_{n}^{2}$. Fig. 7 shows the sensitivity of a normalized symmetric
cumulant nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ measuring a normalized
covariance between ${v^{\prime}}_{n}^{2}$ and $v_{n}^{2}$:

$nSC({v^{\prime}}_{n}^{2},v_{n}^{2})=\frac{\left\langle{v^{\prime}}_{n}^{2}{v}_{n}^{2}\right\rangle}{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle{v}_{n}^{2}\right\rangle}-1=\frac{\left\langle\left({v^{\prime}}_{n}^{2}-\left\langle{v^{\prime}}_{n}^{2}\right\rangle\right)\left(v_{n}^{2}-\left\langle
v_{n}^{2}\right\rangle\right)\right\rangle}{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle}.$ (28)

|
---|---
Figure 7: The sensitivity of nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ is shown as
a function of $\rho(v_{n},v^{\prime}_{n})$ and RS for $v^{\prime}_{n}$ and
$v_{n}$ parametrized using the joint elliptic power distribution (left), and a
bivariate Gaussian distribution (right). Note that the range of RS is
different in the two plots.

The sensitivities of the observable nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ to
$\rho$ and RS, are shown in Fig. 7. The observable
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ is a function of the stochastic variable
${v^{\prime}}_{n}^{2}$ in contrast to the variable $V^{\prime}_{n}V_{n}^{*}$
which was used for the past two observables. Like the other observables,
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ is monotonically increasing in both RS
and $\rho$. Interestingly, the dependence of
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$ is qualitatively very similar to
$\Gamma$, which requires a smaller sample of POIs to measure accurately.
However, nSC displays a larger range of values, and has a fundamentally
different interpretation, owing to a dependence on different random variables
than $\Gamma$. Finally, when considering the restricted domain for the joint
elliptic power distribution, we can see that the observable yields very
similar values across both $\rho$ and RS.

####  IV.2.4 nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$

The sensitivities of the normalized asymmetric cumulant
nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$ are shown in Fig. 8, and have
been used to describe higher order fluctuations between flow harmonics
Bilandzic:2021rgb ; ALICE:2021klf ; ATLAS:2018ngv . Asymmetric cumulants have
also been generalized to incorporate multi-differential correlations in Ref.
Holtermann:2023vwr , and provide a measurement of the genuine correlations
between larger sets of stochastic variables. The observable is included here
to illustrate a higher-order generalization of
nSC$({v^{\prime}}_{n}^{2},v_{n}^{2})$. The asymmetric cumulant we study,
nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$, measures fluctuations
between ${v^{\prime}_{n}}^{2}$ and $v_{n}^{2}$, displaying second order
dependence on $v_{n}^{2}$:

$nASC({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})=\frac{\left\langle{v^{\prime}}_{n}^{2}v_{n}^{4}\right\rangle-\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{4}\right\rangle-2\left\langle v_{n}^{2}\right\rangle\left\langle
v_{n}^{2}{v^{\prime}}_{n}^{2}\right\rangle+2\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle^{2}}{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle^{2}}=\frac{\left\langle\left({v^{\prime}_{n}}^{2}-\left\langle{v^{\prime}_{n}}^{2}\right\rangle\right)\left({v_{n}}^{2}-\left\langle{v_{n}}^{2}\right\rangle\right)^{2}\right\rangle}{{\left\langle{v^{\prime}}_{n}^{2}\right\rangle\left\langle
v_{n}^{2}\right\rangle^{2}}}.$ (29)

The sensitivities of nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$ to RS
and $\rho$ in our toy model can be seen in Fig. 8. nASC clearly increases
along the horizontal and vertical axes, indicating significant sensitivity to
both $\rho$ and RS. The dependence of nASC on both RS and $\rho$ appears to
increase more non-linearly compared to nSC and $\Gamma$, but overall has a
similar qualitative appearance, meaning that in this context, a measurement of
nASC may not add significant information to that provided by a measurement of
$\Gamma$ or nSC. Recalling that this asymmetric cumulant is a central moment,
we can understand the positive result seen at high values for RS as indicative
of a very strong relationship between deviations of ${v^{\prime}}_{n}^{2}$
from the mean, and the overall variance in the distribution of $v_{n}^{2}$,
$\left\langle(v_{n}^{2}-\left\langle
v_{n}^{2}\right\rangle)^{2}\right\rangle.$

|
---|---
Figure 8: The sensitivity of nASC$({v^{\prime}}_{n}^{2},v_{n}^{2},v_{n}^{2})$
is shown as a function of $\rho(v_{n},v^{\prime}_{n})$ and RS for
$v^{\prime}_{n}$ and $v_{n}$ parametrized using the joint elliptic power
distribution (left), and a bivariate Gaussian distribution (right). Note that
the range of RS is different in the two plots.

##  V Applicability

Now that we have discussed the sensitivities of four different correlation
observables using our toy model, we explain how this phenomenological work can
be used to extract values for $\rho$ and RS from experimental data. We then
address considerations about the choices of parametrizations and correlation
models used throughout this paper.

###  V.1 Constraining RS and $\rho$

Given a measurement of any observable detailed in this paper and sufficient
inputs for the toy model specified in Sec. IV (bivariate azimuthal anisotropy
model $p(v^{\prime}_{n},v_{n})$ and a set of fixed parameters), one would be
able to constrain the values of RS and $\rho$ to a curve in RS, $\rho$ space.
Depending on the observables and their sensitivity, this alone may be a
substantial constraint, especially in the context of theoretical predictions.

Figure 9: Values for
$\Gamma(V^{\prime}_{n}V_{n}^{*},v_{n}^{2};v_{n}^{2},v_{n}^{2})$ are displayed
as a function of $\rho$, for both elliptic power (red) and Gaussian (blue)
marginal distributions, and at RS = 1.1 (dotted lines) and RS = 1.4 (solid
lines), alongside a prediction for
$\Gamma(V^{\prime}_{n}V_{n}^{*},v_{n}^{2};v_{n}^{2},v_{n}^{2})$ found in Ref.
Betz:2016ayq  for PbPb collisions at $\sqrt{s_{NN}}=5.02$ TeV, with
$\eta/s=0.05$, at $40-50\%$ centrality for jets at $p_{T}=$ 100 GeV.

As an example, in Fig. 9, we plot values for
$\Gamma(V^{\prime}_{n}V_{n}^{*},v_{n}^{2};v_{n}^{2},v_{n}^{2})$ both with
elliptic power and Gaussian marginals, as a function of $\rho$ for two
different RS values, RS = $1.1$ and RS = $1.4$. These RS values were selected
because they cover most of the range of RS for the joint elliptic power
distribution, shown in Fig. 4. Predictions for this quantity using a combined
jet-hydrodynamical model to predict energy loss were published in Ref.
Betz:2016ayq , for PbPb collisions at $\sqrt{s_{NN}}=5.02$ TeV, measuring
$\Gamma(V^{\prime}_{n}V_{n}^{*},v_{n}^{2};v_{n}^{2},v_{n}^{2})$ for
$\eta/s=0.05$, using jets at $p_{T}=$ 100 GeV. We show the value obtained in
Ref. Betz:2016ayq  for $40-50\%$ centrality, using the elliptic power
distribution parameters from Fig. 1. The predicted value from Ref.
Betz:2016ayq  intersects computed values of
$\Gamma(V^{\prime}_{n}V_{n}^{*},v_{n}^{2};v_{n}^{2},v_{n}^{2})$ at $\rho\geq
0.6$, indicating that for $v_{n}$ and $v^{\prime}_{n}$ modeled using the joint
elliptic power distribution with fixed values detailed in Table 4, $\rho\geq
0.6$. While this relationship differs at different centralities, it suggests
that experimental models predict relatively large correlations between
$v^{\prime}_{n}$ and $v_{n}$ even for jets at high transverse momentum. It is
worth noting that this relationship persists for RS significantly greater than
1, where despite having mean values and fluctuations of different magnitudes,
$v_{n}$ and $v^{\prime}_{n}$ still remain strongly correlated.

Owing to the unique sensitivities displayed by each observable to RS and
$\rho$ as shown in Figs. 5-8, for any two or more observables that display
monotonicity for the inputs, we can obtain RS and $\rho$ values by determining
where the level curves (curves in RS, $\rho$ space where an observable has a
constant value) for each observable intersect. Using plots such as Fig. 8,
alongside experimental measurements of a given observable, we can identify the
level curves in RS, $\rho$ space occupied by each observable. Combining
measurements of multiple observables with different sensitivities allows for
the confinement of the possible values for RS and $\rho$ to separate curves in
RS,$\rho$ space, each corresponding to the level curve of a measured
observable. This means that the actual values for RS and $\rho$ can be
determined by locating where all of the level curves intersect. For a given
parametrization, only two level curves are necessary, but to check for
consistency, or to distinguish a larger number of fluctuation quantities (e.g.
RS, $\rho$, and $\sigma(v_{n})$), any number of simultaneous measurements can
be used. The uncertainties measured for any observable detailed in this paper
will likewise constrain $RS$ and $\rho$ to bands in RS,$\rho$ space containing
all level curves within the uncertainties measured for the observable. The
intersection of these bands will provide a two dimensional overlap region in
RS,$\rho$ space covering all possible points for which both observables take
values within their respective uncertainties. Measuring three or more
observables simultaneously will force the overlap region of uncertainty bands
for each observable to be further reduced, decreasing the region in which all
observables remain within their uncertainties, and providing a stronger
constraint on RS and $\rho$.

A hypothetical example of the intersection of level curves for nASC and
$\Gamma$ is shown in Fig. 10. In this figure, the level curves corresponding
to hypothetical measurements of
$\Gamma(v_{n}^{2},V^{\prime}_{n}V^{*}_{n};v_{n}^{2},v_{n}^{2})=0.02$ and
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*},v_{n}^{2},v_{n}^{2})=0.3$
are shown. Their intersection indicates that in this toy model, RS $\approx
1.85$ and $\rho$ $\approx 0.6$ is the only point in RS, $\rho$ space
satisfying
$\Gamma(v_{n}^{2},V^{\prime}_{n}V^{*}_{n};v_{n}^{2},v_{n}^{2})=0.02$ and
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*},v_{n}^{2},v_{n}^{2})=0.3$.

Figure 10: The level curves for
$\Gamma(v_{n}^{2},V^{\prime}_{n}V^{*}_{n};v_{n}^{2},v_{n}^{2})=0.02$ (orange)
and
$\zeta(V^{\prime}_{n}V_{n}^{*},V^{\prime}_{n}V_{n}^{*},v_{n}^{2},v_{n}^{2})=0.3$
(purple) for varying values of $\rho$ and RS for the bivariate Gaussian
distribution are shown. The values for RS and $\rho$ where the level curves
intersect are also labeled (green).

In general, simply measuring a four-particle correlation with two POI enables
the measurement of enough unique observables, (with different normalizations,
stochastic variables, construction, and sensitivity) such as those detailed in
Sec. II.2, to make the analysis of two or more intersecting level curves
possible.

The toy model detailed in this paper uses joint distributions,
$p(v_{n},v^{\prime}_{n}),$ parametrized with five variables, as seen in Table
4. For each parametrization, only two parameters are varied to evaluate
sensitivities to RS and $\rho$. These two variables were selected due to their
ability to characterize fluctuations in $p(v^{\prime}_{n})$, and correlations
between $v_{n}$ and $v^{\prime}_{n}$. Additionally, most existing azimuthal
anisotropy measurement techniques and observables (such as the scalar product,
or reference multiparticle cumulants) are not sensitive to either RS or
$\rho$. However, in each parametrization, the remaining fixed parameters can
also be important to the determination of an observable’s sensitivity to both
RS and $\rho$. In Appendix B, we detail a number of nonlinear relationships
between the parameter values and multivariate moments for a bivariate Gaussian
distribution $p(v_{n},v^{\prime}_{n})$, illustrating how even the fixed
parameters $\mu(v_{n}),\mu(v^{\prime}_{n}),$ and $\sigma(v_{n})$, can have
significant effects on the observables described in this paper.

The fixed parameters in this model, for both the joint elliptic power and
bivariate Gaussian distributions, characterize $p(v_{n})$ and $\left\langle
v^{\prime}_{n}\right\rangle$. However, there may exist quantities that drive
changes in the parameters defining $p(v_{n})$ and $\left\langle
v^{\prime}_{n}\right\rangle$. An example of such a quantity is centrality, as
$p(v_{n}),$ and parameters $\sigma(v_{n})$ and $\mu(v_{n})$ take on
significantly different values for events with differing centrality
ATLAS:2018ngv . Another example of such a quantity is the harmonic $n$. The
shape of the probability distribution $p(v_{n})$ changes significantly between
$v_{2},v_{3},$ and $v_{4}$ as documented in Ref. ATLAS:2013xzf . To accurately
apply this toy model in any analysis, the “fixed parameters” must be
reevaluated for variations in every quantity that drives changes in $p(v_{n})$
or $\left\langle v^{\prime}_{n}\right\rangle$. As an example, parameters
characterizing $p(v_{n})$ should be evaluated separately for each centrality
bin and for each harmonic before reapplying the toy model.

###  V.2 Additional Considerations

A bivariate parametrization $p(v_{n},v^{\prime}_{n})$ is necessary to generate
data for the toy model defined in this paper. In this section, we address the
effects of different choices for both the marginal distributions, $p(v_{n})$
and $p(v^{\prime}_{n})$, and the correlation model used to generalize
$p(v_{n}),p(v^{\prime}_{n})$ into a bivariate distribution.

While the Gaussian distribution represents Glauber simulation results less
accurately compared to the elliptic power distribution as seen in Fig. 1, we
have found the Gaussian distribution to be more comprehensive as a model for
both $p(v_{n})$ and $p(v^{\prime}_{n})$. Using the Gaussian distribution,
fluctuations in $p(v_{n})$ and $p(v^{\prime}_{n})$ are easily characterized
and varied, and the marginal distributions $p(v_{n})$ and $p(v^{\prime}_{n})$
generalize naturally with the Gaussian copula model into a bivariate Gaussian
distribution. The elliptic power distribution, more accurately represents the
eccentricity distributions generated by the Glauber model Yan:2014afa  as
shown in Fig. 1. However, it demonstrates a limited range of RS for realistic
parameter values, and demonstrates qualitative consistency with sensitivities
obtained from the bivariate Gaussian distributions. Figures 5 \- 8 indicate
that for each observable, the values from the bivariate Gaussian and joint
elliptic power distributions are qualitatively similar, and display similar
sensitivities to both RS and $\rho$.

|
---|---
|

Figure 11: The values for $\Gamma$ (top left), $\zeta$ (top right),nSC (bottom
left), and nASC (bottom right) are plotted as functions of the relative spread
for selected values of $\rho$, for both elliptic power marginal and Gaussian
marginal distributions.

To better assess the sensitivity of each observable to the choice of
parametrization, we provide a direct comparison. In Fig. 11, we display the
values of each observable discussed in Sec. IV.2 for both the joint elliptic
power distribution and the bivariate Gaussian distribution as a function of RS
for four selected values of $\rho$. In each plot, the correspondence between
the values of each observable derived from the joint elliptic power
distribution and the bivariate Gaussian distribution can be assessed.

The observables nSC, $\zeta$, and $\Gamma$ demonstrate differences in value
between both parametrizations, but this difference is far greater for nASC.
This is unsurprising because nASC is constructed using six-particle
correlations, and thus has sensitivity to a joint moment of order six. This
sixth order moment is more sensitive to the differing tail behavior between
the elliptic power and Gaussian parametrizations than any moment used in the
construction of nSC$,\zeta$ or $\Gamma$, since they only rely on moments of
order at most four. Additionally, Fig. 11 shows that the clearest differences
between our two parameterizations show up at high $\rho$, whereas our two
distributions provide nearly identical results at low $\rho$ values. This is
because highly correlated bivariate distributions $p(v_{n},v^{\prime}_{n})$
are constrained more narrowly onto a line within their domain. As a result,
their probability densities become more compact, and typically display more
extreme values due to normalization. This process amplifies differences in the
parametrization shapes, which are reflected by more significant differences
between nSC, $\zeta$, $\Gamma$ and nASC for the bivariate and elliptic power
distributions at higher values of $\rho$.

Experimentally, while there are constraints for the distribution of reference
particle $v_{n}$ ATLAS:2013xzf , there is currently little theoretical insight
or experimental data to guide the development of parametrizations for
$p(v^{\prime}_{n})$. For our toy model, we relied on the assumption that
$v^{\prime}_{n}$ and $v_{n}$ have the same parametrization, but different
parameter values. There is no way to validate the quality of this assumption,
and in reality, $p(v^{\prime}_{n})$ could be fundamentally different than what
our model describes. However, when considering distributions that specifically
are parametrized to reflect fluctuations in $v^{\prime}_{n}$, the use of a
Gaussian distribution for $p(v^{\prime}_{n})$ with a fixed
$\sigma(v^{\prime}_{n})$ is a reasonable approximation, because even if the
shape of the distribution $p(v^{\prime}_{n})$ is vastly different, it will
still have a value for $\sigma(v^{\prime}_{n})$, and will display identical
second order fluctuations to a Gaussian distribution parametrized by the same
$\sigma(v^{\prime}_{n}).$

The correlation model between $v^{\prime}_{n}$ and $v_{n}$ dictates precisely
how the joint fluctuations between $v^{\prime}_{n}$ and $v_{n}$ are produced
and evaluated by observables, and different choices for this correlation model
will likely have a significant effect on the results of the analysis. The
Gaussian copula Meyer_2013 ; Takeuchi_2010  was selected for this analysis
because it allows for the direct specification of $\rho$, and simply reduces
to the equation for a bivariate Gaussian distribution when each marginal is a
Gaussian distribution Takeuchi_2010 . Generally, there is no reason that the
correlations between $v^{\prime}_{n}$ and $v_{n}$ should be described using a
Gaussian copula model over any other type of correlation. However, regardless
of the underlying distribution function of the correlation, the Pearson
correlation coefficient $\rho(v^{\prime}_{n},v_{n})$ representing second order
fluctuations in $v_{n}$ and $v^{\prime}_{n}$ will always exist for a given
joint probability density $p(v^{\prime}_{n},v_{n})$, as long as up to second
order moments for each variable can be evaluated. While there are plenty of
cases in which the correlation coefficient fails to accurately represent
dependency between $v^{\prime}_{n}$ and $v_{n}$, we can always produce a model
that displays a given $\rho(v^{\prime}_{n},v_{n})$ using the Gaussian copula,
and the correlation coefficient is extremely common as a linear correlation
measure for joint probability distributions.

##  VI Discussion and Conclusion

Multiparticle correlations are one of the most powerful and fundamental tools
for studying the QGP produced in heavy ion collisions. In this study, for the
first time, we describe how these multiparticle correlations can be harnessed
to measure fluctuations in jet energy loss, among other phenomena.

In this paper, we have developed a toy model to construct and evaluate
distributions of $v_{n}$ and $v^{\prime}_{n}$ with varying fluctuations in
$v^{\prime}_{n}$, and varying correlations between $v_{n}$ and
$v^{\prime}_{n}$. Using this toy model we characterize the sensitivity of four
unique observables to fluctuations and correlations in $v^{\prime}_{n}$ and
$v_{n}$. While they depend on parametrizations and correlation models, we show
how the sensitivities of each observable can be used to constrain values for
fluctuations in hard probe azimuthal anisotropies, and how they relate to the
fluctuations of the soft sector azimuthal anisotropies within heavy ion
collisions. Greater accuracy can be obtained when constraining the
correlations and fluctuations in $v^{\prime}_{n}$ and $v_{n}$ using a
simultaneous measurement of multiple observables. The method described in this
paper presents a unique opportunity to constrain different theoretical models
for jet energy loss processes and fluctuations by evaluating observables
detailed in Sec. II.2 in both simulation and experimental settings.
Theoretical calculations of the observables in this paper can create baselines
for both the fluctuation and correlations between $v^{\prime}_{n}$ and
$v_{n}$. Futhermore, they can relate these fluctuations to hydrodynamical
variables and jet energy loss models. Experimental measurements of these
observables will provide concrete information about the fluctuations and
correlations in $v^{\prime}_{n}$ and $v_{n}$, as well as constraints on the
associated physics processes characterized by theoretical work.

Recent and ongoing high luminosity LHC runs, coupled with sPHENIX and STAR
data from the upcoming AuAu RHIC run will provide unprecedented precision for
jet measurements. This abundance of data will provide an opportunity to apply
the observables and analysis method discussed in this paper, and isolate
estimates for the fluctuations in jet energy loss within the QGP medium.
Moreover, while we have focused primarily on jet POIs, our approach and
observables can also be used with other types of POI to probe fluctuations in
rare probe production and interactions with the QGP.

##  VII Acknowledgements

A.H., A.M.S., and X. W. acknowledge support from National Science Foundation
Award Number 2111046. J.N.H. acknowledges the support from the US-DOE Nuclear
Science Grant No. DE-SC0023861 and within the framework of the Saturated Glue
(SURGE) Topical Theory Collaboration

## References

  * (1) Wit Busza, Krishna Rajagopal, and Wilke van der Schee.  Heavy Ion Collisions: The Big Picture, and the Big Questions.  Ann. Rev. Nucl. Part. Sci., 68:339–376, 2018.  arXiv:1802.04801, doi:10.1146/annurev-nucl-101917-020852.
  * (2) Ulrich Heinz and Raimond Snellings.  Collective flow and viscosity in relativistic heavy-ion collisions.  Ann. Rev. Nucl. Part. Sci., 63:123–151, 2013.  arXiv:1301.2826, doi:10.1146/annurev-nucl-102212-170540.
  * (3) Matthew Luzum and Hannah Petersen.  Initial State Fluctuations and Final State Correlations in Relativistic Heavy-Ion Collisions.  J. Phys. G, 41:063102, 2014.  arXiv:1312.5503, doi:10.1088/0954-3899/41/6/063102.
  * (4) R. Derradi de Souza, Tomoi Koide, and Takeshi Kodama.  Hydrodynamic Approaches in Relativistic Heavy Ion Reactions.  Prog. Part. Nucl. Phys., 86:35–85, 2016.  arXiv:1506.03863, doi:10.1016/j.ppnp.2015.09.002.
  * (5) Leticia Cunqueiro and Anne M. Sickles.  Studying the QGP with Jets at the LHC and RHIC.  Prog. Part. Nucl. Phys., 124:103940, 2022.  arXiv:2110.14490, doi:10.1016/j.ppnp.2022.103940.
  * (6) Guang-You Qin and Xin-Nian Wang.  Jet quenching in high-energy heavy-ion collisions.  Int. J. Mod. Phys. E, 24(11):1530014, 2015.  arXiv:1511.00790, doi:10.1142/S0218301315300143.
  * (7) Z. Citron et al.  Report from Working Group 5: Future physics opportunities for high-density QCD at the LHC with heavy-ion and proton beams.  CERN Yellow Rep. Monogr., 7:1159–1410, 2019.  arXiv:1812.06772, doi:10.23731/CYRM-2019-007.1159.
  * (8) P. Achenbach et al.  The Present and Future of QCD.  3 2023.  arXiv:2303.02579.
  * (9) M. Arslandok et al.  Hot QCD White Paper.  3 2023.  arXiv:2303.17254.
  * (10) J. Takahashi, B. M. Tavares, W. L. Qian, R. Andrade, F. Grassi, Yogiro Hama, T. Kodama, and N. Xu.  Topology studies of hydrodynamics using two particle correlation analysis.  Phys. Rev. Lett., 103:242301, 2009.  arXiv:0902.4870, doi:10.1103/PhysRevLett.103.242301.
  * (11) B. Alver and G. Roland.  Collision geometry fluctuations and triangular flow in heavy-ion collisions.  Phys. Rev. C, 81:054905, 2010.  [Erratum: Phys.Rev.C 82, 039903 (2010)].  arXiv:1003.0194, doi:10.1103/PhysRevC.82.039903.
  * (12) J. Adams et al.  Azimuthal anisotropy in Au+Au collisions at s(NN)**(1/2) = 200-GeV.  Phys. Rev. C, 72:014904, 2005.  arXiv:nucl-ex/0409033, doi:10.1103/PhysRevC.72.014904.
  * (13) A. Adare et al.  Measurements of Higher-Order Flow Harmonics in Au+Au Collisions at $\sqrt{s_{NN}}=200$ GeV.  Phys. Rev. Lett., 107:252301, 2011.  arXiv:1105.3928, doi:10.1103/PhysRevLett.107.252301.
  * (14) L. Adamczyk et al.  Third Harmonic Flow of Charged Particles in Au+Au Collisions at sqrtsNN = 200 GeV.  Phys. Rev. C, 88(1):014904, 2013.  arXiv:1301.2187, doi:10.1103/PhysRevC.88.014904.
  * (15) Jaroslav Adam et al.  Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV.  Phys. Rev. Lett., 116(13):132302, 2016.  arXiv:1602.01119, doi:10.1103/PhysRevLett.116.132302.
  * (16) Albert M Sirunyan et al.  Pseudorapidity and transverse momentum dependence of flow harmonics in pPb and PbPb collisions.  Phys. Rev. C, 98(4):044902, 2018.  arXiv:1710.07864, doi:10.1103/PhysRevC.98.044902.
  * (17) Morad Aaboud et al.  Measurement of the azimuthal anisotropy of charged particles produced in $\sqrt{s_{{}_{\text{NN}}}}$ = 5.02 TeV Pb+Pb collisions with the ATLAS detector.  Eur. Phys. J. C, 78(12):997, 2018.  arXiv:1808.03951, doi:10.1140/epjc/s10052-018-6468-7.
  * (18) Jonah E. Bernhard, J. Scott Moreland, and Steffen A. Bass.  Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma.  Nature Phys., 15(11):1113–1117, 2019.  doi:10.1038/s41567-019-0611-8.
  * (19) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings.  Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum.  Phys. Rev. C, 103(5):054909, 2021.  arXiv:2010.15134, doi:10.1103/PhysRevC.103.054909.
  * (20) D. Everett et al.  Phenomenological constraints on the transport properties of QCD matter with data-driven model averaging.  Phys. Rev. Lett., 126(24):242301, 2021.  arXiv:2010.03928, doi:10.1103/PhysRevLett.126.242301.
  * (21) Georges Aad et al.  Measurement of the distributions of event-by-event flow harmonics in lead-lead collisions at = 2.76 TeV with the ATLAS detector at the LHC.  JHEP, 11:183, 2013.  arXiv:1305.2942, doi:10.1007/JHEP11(2013)183.
  * (22) Jaroslav Adam et al.  Azimuthal anisotropy of charged jet production in $\sqrt{s_{\rm NN}}$ = 2.76 TeV Pb-Pb collisions.  Phys. Lett. B, 753:511–525, 2016.  arXiv:1509.07334, doi:10.1016/j.physletb.2015.12.047.
  * (23) Georges Aad et al.  Measurements of azimuthal anisotropies of jet production in Pb+Pb collisions at $\sqrt{s_{NN}}=$ 5.02 TeV with the ATLAS detector.  Phys. Rev. C, 105(6):064903, 2022.  arXiv:2111.06606, doi:10.1103/PhysRevC.105.064903.
  * (24) Armen Tumasyan et al.  Azimuthal anisotropy of dijet events in PbPb collisions at $\sqrt{s_{\textrm{NN}}}$ = 5.02 TeV.  JHEP, 07:139, 2023.  arXiv:2210.08325, doi:10.1007/JHEP07(2023)139.
  * (25) Measurement of the nuclear modification factor of $b$-jets in 5.02 TeV Pb+Pb collisions with the ATLAS detector.  4 2022.  arXiv:2204.13530.
  * (26) Georges Aad et al.  Comparison of inclusive and photon-tagged jet suppression in 5.02 TeV Pb+Pb collisions with ATLAS.  3 2023.  arXiv:2303.10090.
  * (27) Shreyasi Acharya et al.  Measurement of the groomed jet radius and momentum splitting fraction in pp and Pb$-$Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV.  Phys. Rev. Lett., 128(10):102001, 2022.  arXiv:2107.12984, doi:10.1103/PhysRevLett.128.102001.
  * (28) Measurement of substructure-dependent jet suppression in Pb+Pb collisions at 5.02 TeV with the ATLAS detector.  11 2022.  arXiv:2211.11470.
  * (29) Barbara Betz, Miklos Gyulassy, Matthew Luzum, Jorge Noronha, Jacquelyn Noronha-Hostler, Israel Portillo, and Claudia Ratti.  Cumulants and nonlinear response of high $p_{T}$ harmonic flow at $\sqrt{s_{NN}}=5.02$ TeV.  Phys. Rev. C, 95(4):044901, 2017.  arXiv:1609.05171, doi:10.1103/PhysRevC.95.044901.
  * (30) Armen Tumasyan et al.  Probing Charm Quark Dynamics via Multiparticle Correlations in Pb-Pb Collisions at $\sqrt{s_{\mathrm{NN}}}$ =5.02 TeV.  Phys. Rev. Lett., 129(2):022001, 2022.  arXiv:2112.12236, doi:10.1103/PhysRevLett.129.022001.
  * (31) Abraham Holtermann, Jacquelyn Noronha-Hostler, Anne M. Sickles, and Xiaoning Wang.  Multiparticle correlations, cumulants, and moments sensitive to fluctuations in rare-probe azimuthal anisotropy in heavy ion collisions.  Phys. Rev. C, 108(6):064901, 2023.  arXiv:2307.16796, doi:10.1103/PhysRevC.108.064901.
  * (32) A. Adare et al.  An Upgrade Proposal from the PHENIX Collaboration.  1 2015.  arXiv:1501.06197.
  * (33) Ante Bilandzic, Marcel Lesch, and Seyed Farid Taghavi.  New estimator for symmetry plane correlations in anisotropic flow analyses.  Phys. Rev. C, 102(2):024910, 2020.  arXiv:2004.01066, doi:10.1103/PhysRevC.102.024910.
  * (34) Ante Bilandzic, Raimond Snellings, and Sergei Voloshin.  Flow analysis with cumulants: Direct calculations.  Phys. Rev. C, 83:044913, 2011.  arXiv:1010.0233, doi:10.1103/PhysRevC.83.044913.
  * (35) Sergei A. Voloshin, Arthur M. Poskanzer, and Raimond Snellings.  Collective phenomena in non-central nuclear collisions.  Landolt-Bornstein, 23:293–333, 2010.  arXiv:0809.2949, doi:10.1007/978-3-642-01539-7_10.
  * (36) Observation of flow angle and flow magnitude fluctuations in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV at the LHC.  6 2022.  arXiv:2206.04574.
  * (37) Ante Bilandzic, Christian Holm Christensen, Kristjan Gulbrandsen, Alexander Hansen, and You Zhou.  Generic framework for anisotropic flow analyses with multiparticle azimuthal correlations.  Phys. Rev. C, 89(6):064904, 2014.  arXiv:1312.3572, doi:10.1103/PhysRevC.89.064904.
  * (38) Ante Bilandzic, Marcel Lesch, Cindy Mordasini, and Seyed Farid Taghavi.  Multivariate cumulants in flow analyses: The next generation.  Phys. Rev. C, 105(2):024912, 2022.  arXiv:2101.05619, doi:10.1103/PhysRevC.105.024912.
  * (39) Sergei A. Voloshin, Arthur M. Poskanzer, Aihong Tang, and Gang Wang.  Elliptic flow in the Gaussian model of eccentricity fluctuations.  Phys. Lett. B, 659:537–541, 2008.  arXiv:0708.0800, doi:10.1016/j.physletb.2007.11.043.
  * (40) Albert M Sirunyan et al.  Non-Gaussian elliptic-flow fluctuations in PbPb collisions at $\sqrt{\smash[b]{s_{{}_{\text{NN}}}}}=5.02$ TeV.  Phys. Lett. B, 789:643–665, 2019.  arXiv:1711.05594, doi:10.1016/j.physletb.2018.11.063.
  * (41) Constantin Loizides, Jason Kamin, and David d’Enterria.  Improved Monte Carlo Glauber predictions at present and future nuclear colliders.  Phys. Rev. C, 97(5):054910, 2018.  [Erratum: Phys.Rev.C 99, 019901 (2019)].  arXiv:1710.07098, doi:10.1103/PhysRevC.97.054910.
  * (42) Jacquelyn Noronha-Hostler, Li Yan, Fernando G. Gardim, and Jean-Yves Ollitrault.  Linear and cubic response to the initial eccentricity in heavy-ion collisions.  Phys. Rev. C, 93(1):014909, 2016.  arXiv:1511.03896, doi:10.1103/PhysRevC.93.014909.
  * (43) Derek Teaney and Li Yan.  Triangularity and Dipole Asymmetry in Heavy Ion Collisions.  Phys. Rev. C, 83:064904, 2011.  arXiv:1010.1876, doi:10.1103/PhysRevC.83.064904.
  * (44) Fernando G. Gardim, Frederique Grassi, Matthew Luzum, and Jean-Yves Ollitrault.  Mapping the hydrodynamic response to the initial geometry in heavy-ion collisions.  Phys. Rev. C, 85:024908, 2012.  arXiv:1111.6538, doi:10.1103/PhysRevC.85.024908.
  * (45) H. Niemi, G. S. Denicol, H. Holopainen, and P. Huovinen.  Event-by-event distributions of azimuthal asymmetries in ultrarelativistic heavy-ion collisions.  Phys. Rev. C, 87(5):054901, 2013.  arXiv:1212.1008, doi:10.1103/PhysRevC.87.054901.
  * (46) Derek Teaney and Li Yan.  Non linearities in the harmonic spectrum of heavy ion collisions with ideal and viscous hydrodynamics.  Phys. Rev. C, 86:044908, 2012.  arXiv:1206.1905, doi:10.1103/PhysRevC.86.044908.
  * (47) Zhi Qiu and Ulrich W. Heinz.  Event-by-event shape and flow fluctuations of relativistic heavy-ion collision fireballs.  Phys. Rev. C, 84:024911, 2011.  arXiv:1104.0650, doi:10.1103/PhysRevC.84.024911.
  * (48) Fernando G. Gardim, Jacquelyn Noronha-Hostler, Matthew Luzum, and Frédérique Grassi.  Effects of viscosity on the mapping of initial to final state in heavy ion collisions.  Phys. Rev. C, 91(3):034902, 2015.  arXiv:1411.2574, doi:10.1103/PhysRevC.91.034902.
  * (49) Thorsten Renk and Harri Niemi.  Constraints from $v_{2}$ fluctuations for the initial-state geometry of heavy-ion collisions.  Phys. Rev. C, 89(6):064907, 2014.  arXiv:1401.2069, doi:10.1103/PhysRevC.89.064907.
  * (50) G. S. Denicol, C. Gale, S. Jeon, J. F. Paquet, and B. Schenke.  Effect of initial-state nucleon-nucleon correlations on collective flow in ultra-central heavy-ion collisions.  6 2014.  arXiv:1406.7792.
  * (51) Li Yan, Jean-Yves Ollitrault, and Arthur M. Poskanzer.  Eccentricity distributions in nucleus-nucleus collisions.  Phys. Rev. C, 90(2):024903, 2014.  arXiv:1405.6595, doi:10.1103/PhysRevC.90.024903.
  * (52) Li Yi, Fuqiang Wang, and Aihong Tang.  Possibility to disentangle anisotropic flow, flow fluctuation, and nonflow assuming Gaussian fluctuations.  1 2011.  arXiv:1101.4646.
  * (53) Tsutomu T. Takeuchi.  Constructing a bivariate distribution function with given marginals and correlation: application to the galaxy luminosity function.  Monthly Notices of the Royal Astronomical Society, pages no–no, may 2010.  URL: https://doi.org/10.1111%2Fj.1365-2966.2010.16778.x, doi:10.1111/j.1365-2966.2010.16778.x.
  * (54) Christian Meyer.  The bivariate normal copula.  Communications in Statistics - Theory and Methods, 42(13):2402–2422, jul 2013.  URL: https://doi.org/10.1080%2F03610926.2011.611316, doi:10.1080/03610926.2011.611316.
  * (55) Matthew Luzum and Jean-Yves Ollitrault.  Eliminating experimental bias in anisotropic-flow measurements of high-energy nuclear collisions.  Phys. Rev. C, 87(4):044907, 2013.  arXiv:1209.2323, doi:10.1103/PhysRevC.87.044907.
  * (56) Abraham Holtermann.  2023\.  URL: https://github.com/AHoltermann/multiparticle_probe_flutuation.
  * (57) Nicolas Borghini, Phuong Mai Dinh, and Jean-Yves Ollitrault.  Flow analysis from multiparticle azimuthal correlations.  Phys. Rev. C, 64:054901, 2001.  arXiv:nucl-th/0105040, doi:10.1103/PhysRevC.64.054901.
  * (58) Giuliano Giacalone, Li Yan, Jacquelyn Noronha-Hostler, and Jean-Yves Ollitrault.  Symmetric cumulants and event-plane correlations in Pb + Pb collisions.  Phys. Rev. C, 94(1):014906, 2016.  arXiv:1605.08303, doi:10.1103/PhysRevC.94.014906.
  * (59) Morad Aaboud et al.  Correlated long-range mixed-harmonic fluctuations measured in $pp$, $p$+Pb and low-multiplicity Pb+Pb collisions with the ATLAS detector.  Phys. Lett. B, 789:444–471, 2019.  arXiv:1807.02012, doi:10.1016/j.physletb.2018.11.065.
  * (60) Shreyasi Acharya et al.  Multiharmonic Correlations of Different Flow Amplitudes in Pb-Pb Collisions at $\sqrt{s_{{}_{NN}}}=2.76$ TeV.  Phys. Rev. Lett., 127(9):092302, 2021.  arXiv:2101.02579, doi:10.1103/PhysRevLett.127.092302.
  * (61) Ante Bilandzic.  Anisotropic flow measurements in ALICE at the large hadron collider.  PhD thesis, Utrecht U., 2012.

##  Appendix A Evaluation of Multiparticle correlations

The $Q$-Cumulant method Bilandzic:2010jr ; Bilandzic:2012wva  is a more
computationally efficient method used to measure angular particle correlations
without double counting particles, and ignoring contributions from higher
harmonics. The $Q$-Cumulant method writes the even moments of $v_{n}$,
$\langle v_{n}^{2k}\rangle$, as a function of event $Q_{n}$ vectors with
weights multiplied by different powers $r$, a sum over all $N$ reference
particles in an event:

$Q_{n,r}=\sum_{i}w_{i}^{r}e^{in\phi_{i}}$ (30)

where $\phi_{i}$ indicates the azimuthal angle of each particle. We include a
weighting $w_{i}$ for each particle, which is raised to the $r$ power. To
incorporate the azimuthal anisotropies displayed by jets, or any other
differential particle species, we use the “reduced” vector:

$q_{n,r}=\sum_{j}w_{j}^{r}e^{in\psi_{j}}.$ (31)

where $j$ is an index for for each differential particle, featuring azimuthal
angle $\psi.$

For this analysis we also assume there is no overlap between the POI and the
reference particles, as specifically jets and soft particles will have
significantly different selection criteria.

###  A.1 Two Particle Correlations

Following the formalism from Bilandzic:2013kga ; Holtermann:2023vwr , we
obtain expressions for the two-particle and four-particle correlations used to
construct the observables discussed in this paper. The expressions are arrived
at by evaluating the overlap between particles, as illustrated by a derivation
of the expression for $\left\langle{v^{\prime}}_{n}^{2}\right\rangle$,

$\displaystyle\left\langle{v^{\prime}}_{n}^{2}\right\rangle$
$\displaystyle=\left\langle e^{in(\psi_{1}-\psi_{2})}\right\rangle$
$\displaystyle=\frac{\sum_{k=1}^{N}\sum_{j=1,j\neq
i}^{N}w_{k}w_{j}e^{in(\phi_{k}-\phi_{j})}}{\sum_{k=1}^{N}\sum_{j=1,j\neq
i}^{N}w_{k}w_{j}}.$ (32)

Then, we note that the iterated sum over unique pairs of particles can be
written as a subtraction of all duplicates ($\psi_{a}$,$\psi_{a}$) from every
possible combination of angle pairs (note that double counting of angle pairs
$(\psi_{a},\psi_{b}),(\psi_{b},\psi_{a})$ is accounted for in the denominator
of the fraction):

$\displaystyle\frac{\sum_{k=1}^{N}\sum_{j=1,j\neq
i}^{N}w_{k}w_{j}e^{in(\psi_{k}-\psi_{j})}}{\sum_{k=1}^{N}\sum_{j=1,j\neq
i}^{N}w_{k}w_{j}}=\frac{\sum_{k=1}^{N}\sum_{j=1}^{N}w_{j}w_{k}e^{in(\psi_{k}-\psi_{j})}-\sum_{l=1}^{N}w_{l}^{2}e^{in(\psi_{l}-\psi_{l})}}{\sum_{k=1}^{N}\sum_{j=1}^{N}w_{j}w_{k}-\sum_{l=1}^{N}w_{l}^{2}}.$
(33)

Using the definition of $q_{n,r}$ given in Eq. (31), we obtain:

$\displaystyle\frac{\sum_{k=1}^{N}\sum_{j=1}^{N}w_{j}w_{k}e^{in(\psi_{k}-\psi_{j})}-\sum_{l=1}^{N}w_{l}^{2}e^{in(\psi_{l}-\psi_{l})}}{\sum_{k=1}^{N}\sum_{j=1}^{N}w_{j}w_{k}-\sum_{l=1}^{N}w_{l}^{2}}=\frac{q_{n,1}q_{n,-1}-q_{0,2}}{q_{0,1}^{2}-q_{0,2}}.$
(34)

Following this procedure for each two-particle and four-particle correlation
in Table 1, we present their form as a function of the $q_{n,r}$ and $Q_{n,r}$
vectors. Note that all averages here, indicated by angle brackets, indicate an
average taken over an ensemble of events.

$\left\langle
v_{n}^{2}\right\rangle=\left\langle\frac{|Q_{n,1}|^{2}-Q_{0,2}}{|Q_{0,1}|^{2}-Q_{0,2}}\right\rangle$
(35)

$\left\langle
V^{\prime}_{n}V_{n}^{*}\right\rangle=\left\langle\frac{q_{n,1}Q_{n,1}^{*}}{q_{0,1}Q_{0,1}^{*}}\right\rangle$
(36)

$\left\langle{v^{\prime}}_{n}^{2}\right\rangle=\left\langle\frac{|q_{n,1}|^{2}-q_{0,2}}{|q_{0,1}|^{2}-q_{0,2}}\right\rangle$
(37)

###  A.2 Four-Particle Correlations

Following a similar process as for the two particle correlations, and
described in detail in Ref. Holtermann:2023vwr , we obtain the following
expression for the various four-particle correlations using $Q_{n}$ vectors:

$\begin{split}\left\langle
v_{n}^{4}\right\rangle&=\left\langle\left(\left|Q_{n,1}\right|^{4}-Q_{2n,2}Q_{n,1}^{*2}-4Q_{0,2}\left|Q_{n,1}\right|^{2}-Q_{2n,2}^{*}Q_{n,1}^{2}+|Q_{2n,2}|^{2}+2\left|Q_{0,2}\right|^{2}+\right.\right.\\\
&\left.\left.+4Q_{n,3}Q_{n,1}^{*}+4Q_{n,3}^{*}Q_{n,1}-6Q_{0,4}\right)\right.\\\
&\times\left.\left(Q_{0,1}^{4}-6Q_{0,1}^{2}Q_{0,2}+3Q_{0,2}^{2}+8Q_{0,1}Q_{0,3}-6Q_{0,4}\right)^{-1}\right\rangle\\\
\end{split}$ (38)

$\begin{split}\left\langle
V^{\prime}_{n}V_{n}^{*}v_{n}^{2}\right\rangle&=\left\langle\left(q_{n,1}Q_{n,1}^{*}|Q_{n,1}|^{2}-2q_{n,1}Q_{n,1}^{*}Q_{0,2}-q_{n,1}Q_{n,1}Q_{2n,2}^{*}+2q_{n,1}Q_{n,3}^{*}\right)\right.\\\
&\times\left.\left(q_{0,1}Q_{0,1}^{3}-3q_{0,1}Q_{0,1}Q_{0,2}+2q_{0,1}Q_{0,3}\right)^{-1}\right\rangle\end{split}$
(39)

$\begin{split}\left\langle(v^{\prime})_{n}^{2}v_{n}^{2}\right\rangle&=\left\langle\left(|q_{n,1}|^{2}|Q_{n,1}|^{2}-q_{0,2}|Q_{n,1}|^{2}-|q_{n,1}|^{2}Q_{0,2}+q_{0,2}Q_{0,2}\right)\right.\\\
&\times\left.\left(q_{0,1}^{2}Q_{0,1}^{2}-q_{0,2}Q_{0,1}^{2}-q_{0,1}^{2}Q_{0,2}+q_{0,2}Q_{0,2}\right)^{-1}\right\rangle\end{split}$
(40)

$\begin{split}\left\langle(V^{\prime}_{n}V_{n}^{*})^{2}\right\rangle&=\left\langle\left(q_{n,1}^{2}Q_{n,1}^{*2}-q_{2n,2}Q_{n,1}^{*2}-q_{n,1}^{2}Q_{2n,2}^{*}+q_{2n,2}Q_{2n,2}^{*}\right)\right.\\\
&\times\left.\left(q_{0,1}^{2}Q_{0,1}^{*2}-q_{0,2}Q_{0,1}^{*2}-q_{0,1}^{2}Q_{0,2}^{*}+q_{0,2}Q_{0,2}^{*}\right)^{-1}\right\rangle.\end{split}$
(41)

##  Appendix B  Analytic formulas for joint moments of a bivariate Gaussian
Distribution

The probability density function for two arbitrary stochastic variables $X$
and $Y$ parametrized with a bivariate Gaussian distribution with parameters
($\mu_{x}$, $\mu_{y}$, $\sigma_{x}$, $\sigma_{y}$, $\rho$) is given by:

$p(X,Y)=\frac{1}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rho^{2}}}\exp\left(\frac{1}{2(1-\rho)^{2}}\left(\frac{(X-\mu_{x})^{2}}{\sigma_{x}^{2}}+\frac{(Y-\mu_{y})^{2}}{\sigma_{y}^{2}}-\frac{2\rho(X-\mu_{x})(Y-\mu_{y})}{\sigma(v_{n})\sigma(v^{\prime}_{n})}\right)\right)$
(42)

The moment generating function is used to evaluate the moments of a
distribution. Generally for a two variable probability distribution, the
moment generating function is a function of two dummy variables, $t_{x}$ and
$t_{y}$, and is written as follows

$M_{X,Y}(t_{x},t_{y})=\left\langle e^{t_{x}X+t_{y}Y}\right\rangle=\iint
e^{t_{x}X+t_{y}Y}p(X,Y)dXdY$ (43)

where in this case, an average indicated by angle brackets is an integral over
all values of $X$ and $Y$. For a bivariate Gaussian distribution the moment
generating function has the form:

$M_{X,Y}(t_{x},t_{y})=\exp{\left(\mu_{x}t_{x}+\mu_{y}t_{y}+\frac{1}{2}\sigma_{x}^{2}t_{x}^{2}+\frac{1}{2}\sigma_{y}^{2}t_{y}^{2}+\rho\sigma_{x}\sigma_{y}t_{x}t_{y}\right)}$
(44)

Finally, the moments are extracted from the moment generating function by
evaluating terms in the Taylor expansion of $M_{X,Y}(t_{x},t_{y})$ around the
origin as follows:

$\left\langle
X^{n}Y^{m}\right\rangle=\frac{\partial^{n+m}M_{XY}(t_{x},t_{y})}{\partial
X^{n}\partial Y^{n}}\Bigg{|}_{t_{x}=t_{y}=0}.$ (45)

For a concise example, we evaluate $\left\langle XY\right\rangle$ using this
method.

$\left\langle XY\right\rangle=\frac{\partial^{2}M_{XY}(t_{x},t_{y})}{\partial
X\partial Y}\Bigg{|}_{t_{x}=t_{y}=0}=\frac{\partial}{\partial Y}\frac{\partial
M_{XY}(t_{x},t_{y})}{\partial
x}\Bigg{|}_{t_{x}=t_{y}=0}=\frac{\partial}{\partial
Y}\left(M_{XY}(t_{x},t_{y})\left[\mu_{x}+\sigma_{x}^{2}t_{x}+\rho\sigma_{x}\sigma_{y}t_{y}\right]\right)\Bigg{|}_{t_{x}=t_{y}=0}$
(46)
$=M_{XY}(t_{x},t_{y})\left[\mu_{x}+\sigma_{x}^{2}t_{x}+\rho\sigma_{x}\sigma_{y}t_{y}\right]\left[\mu_{y}+\sigma_{y}^{2}t_{y}+\rho\sigma_{x}\sigma_{y}t_{x}\right]+M_{XY}(t_{x},t_{y})\left[\rho\sigma_{x}\sigma_{y}\right]\Bigg{|}_{t_{x}=t_{y}=0}$
(47)

Then, we evaluate at $t_{x}=t_{y}=0$:

$\left\langle XY\right\rangle=\frac{\partial^{2}M_{X,Y}(0,0)}{\partial
x\partial y}=\mu_{x}\mu_{y}+\rho\sigma_{x}\sigma_{y}$ (48)

Using this procedure, we derive formulas for the moments used in this paper up
to four particles (detailed in Table 1).

$\displaystyle\left\langle X^{2}\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{x}^{2}+\sigma_{x}^{2}$ (49) $\displaystyle\left\langle
Y^{2}\right\rangle$ $\displaystyle=$ $\displaystyle\mu_{y}^{2}+\sigma_{y}^{2}$
(50) $\displaystyle\left\langle XY\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{x}\mu_{y}+\rho\sigma_{x}\sigma_{y}$ (51)
$\displaystyle\left\langle X^{4}\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{x}^{4}+6\sigma_{x}^{2}\mu_{x}^{2}+3\sigma_{x}^{4}$ (52)
$\displaystyle\left\langle X^{3}Y\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{x}^{3}\mu_{y}+3\mu_{x}^{2}\rho\sigma_{x}\sigma_{y}+3\sigma_{x}^{2}\mu_{x}\mu_{y}+3\sigma_{x}^{2}\rho\sigma_{x}\sigma_{y}$
(53) $\displaystyle\left\langle X^{2}Y^{2}\right\rangle$ $\displaystyle=$
$\displaystyle\mu_{x}^{2}\mu_{y}^{2}+4\rho\sigma_{x}\sigma_{y}\mu_{x}\mu_{y}+\mu_{x}^{2}\sigma_{y}^{2}+\sigma_{y}^{2}\mu_{x}^{2}+\sigma_{x}^{2}\sigma_{y}^{2}+2\rho^{2}\sigma_{x}^{2}\sigma_{y}^{2}$
(54)