System size and shape dependences of collective flow fluctuations in relativistic nuclear collisions

We start with a three-parameter parametrization of the nucleon density distribution function inside a nucleus as DeVries:1987atn ,

where $r_{p}$ represents the radial position of a nucleon, $\rho_{0}=0.17$ fm^-3 is the equilibrium density of nuclear matter, $R$ and $a$ are the radius and the surface thickness parameter of the nucleus. The equation above returns to the standard two-parameter Woods-Saxon distribution with $w=0$. The model parameters of ¹⁹⁷Au, ⁶⁴Cu, and ¹⁶O nuclei used in our present study are listed in Tab. 1. According to these distributions, we use the Monte-Carlo (MC) method to sample the positions of nucleons inside the projectile and target nuclei. To prevent two nucleons from being too close to each other inside a nucleus, a minimum distance of $d=0.81$ fm is imposed in our MC sampling.

By assuming high-energy nucleons stream along straight lines without changing their directions when being scattered, we let two nucleons (one from projectile along the $+\hat{z}$ direction and one from target along the $-\hat{z}$ direction) collide when their transverse distance is smaller than $(\sigma_{\mathrm{NN}}^{\mathrm{inel}}/\pi)^{-\frac{1}{2}}$, where the inelastic nucleon-nucleon cross section is taken as $\sigma_{\mathrm{NN}}^{\mathrm{inel}}=42$ mb at $\sqrt{s_{\mathrm{NN}}}=200$ GeV PHOBOS:2007vdf . The positions of both binary collisions (taken as the mid-points of nucleon pairs) and nucleons that participate in these collisions (or participant nucleons) are recorded and contribute to particle production from a nucleus-nucleus collision event. The multiplicity of the final state charged particles then follow the form Qin:2010pf

where $\alpha$ is a parameter controlling the relative contribution from the participant nucleon number ($N_{\mathrm{part}}$) and the binary collision number ($N_{\mathrm{bin}}$), which can be determined by the centrality dependence of charged particle yield later. We note that although the contribution from binary collisions might be negligible when the nuclear collisions are less energetic, e.g. $\sqrt{s_{\mathrm{NN}}}\lesssim 62.4$ GeV Wu:2021fjf , it is crucial for a simultaneous description of hadron observables at different centralities at $\sqrt{s_{\mathrm{NN}}}=200$ GeV in our present work. Although this two-component Glauber model has been shown inaccurate in ultra-central collisions STAR:2015mki , it is still considered a reasonable and convenient model for the initial condition of QGP as long as the centrality is not very small.

The local entropy density $s$ and the normalized baryon number density $n_{0}$ then read:

where the former is contributed by both participant nucleons and binary collisions, while the latter is only contributed by participant nucleons. Here, $x$ and $y$ are coordinates in the transverse plane, and $\eta_{\mathrm{s}}$ denotes the spacetime rapidity. Following Ref. Wu:2021fjf , $s_{\mathrm{part}}$ and $n_{\mathrm{part}}$ above at the initial time ($\tau_{0}$) of hydrodynamic evolution are given by

where $\tilde{s}_{\mathrm{P}}$ and $\tilde{s}_{\mathrm{T}}$ represent entropy density in the transverse plane contributed by the projectile and target respectively, $H_{\mathrm{P/T}}^{s}$ and $H_{\mathrm{P/T}}^{n}$ are the longitudinal envelop functions for entropy density and baryon number density respectively, and $K$ is an overall factor that controls the magnitude of the initial entropy density.

The $\tilde{s}_{\mathrm{P/T}}$ function can be constructed using the distribution of participant nucleons in the projectile/target nucleus as

where $(x_{i},y_{i})$ represents the transverse position of the $i$-th participant nucleon obtained from the MC-Glauber model, $\sigma_{r}$ is the transverse Gaussian smearing width. The longitudinal envelop functions take the forms of Bozek:2010bi ; Bozek:2011if ; Bozek:2013uha ; Denicol:2018wdp ; Wu:2021fjf

in which parameters $\eta_{0}^{\mathrm{s}}$, $\sigma_{s}$, $\eta^{0}_{n;\mathrm{P/T}}$, $\sigma_{n;\mathrm{T/P}}$ will be determined by the charged particle distribution along the longitudinal direction, $N$ is the normalization factor for $n_{0}$, and $\eta_{\mathrm{max}}=y_{\mathrm{beam}}=\mathrm{arctanh}(v_{\mathrm{beam}})$ is the rapidity of a projectile beam nucleon, or the maximum $\eta_{\mathrm{s}}$ it can reach with velocity $v_{\mathrm{beam}}=\sqrt{s_{\mathrm{NN}}}/(2m_{\mathrm{N}})$ with $m_{\mathrm{N}}$ the nucleon mass. Alternative ways of introducing the longitudinal profiles of the QGP are discussed in Refs. Hirano:2005xf ; Okai:2017ofp .

Similarly, one can define the $s_{\mathrm{bin}}$ part in Eq. (4) as

with

where $i$ runs over the binary collision points in a nucleus-nucleus collision event, $\eta_{\mathrm{w}}$ and $\sigma_{\mathrm{w}}$ are model parameters.

We summarize in Tab. 2 the parameters we use for constructing the initial entropy density and the baryon number density through Eq. (4) to Eq. (13). They are fitted from the existing data in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV in the next section. We assume these parameters are solely sensitive to the beam energy, and remain the same across different collision systems. The evolution of nuclear matter prior to the QGP phase is not included in this work, which becomes more important for lower energy collisions Dore:2020fiq . To partly compensate effects of this pre-equilibrium evolution, we choose the initial proper time to be larger than the overlap time between the colliding nuclei, $\tau_{0}>2R/\sinh(y_{\mathrm{beam}})$, in order to allow more time for the system to approach local equilibrium before hydrodynamics starts.

With the initial condition constructed in the previous subsection, we use the (3+1)-D viscous hydrodynamics model CLVisc Pang:2018zzo ; Wu:2021fjf to simulate the evolution of the QGP event-by-event. The hydrodynamic equations are based on the energy-momentum conservation and the net baryon number conservation as

where $T^{\mu\nu}$ and $J^{\mu}$ are the energy-momentum tensor and the net baryon current respectively. They can be further decomposed as:

where $e$ is the energy density, $U^{\mu}$ is the 4-velocity of the fluid cell, $P$ is the pressure, $\Delta^{\mu\nu}=g^{\mu\nu}-U^{\mu}U^{\nu}$ is the projector operator, $\pi^{\mu\nu}$ is the shear-stress tensor, $n$ is the net baryon number density, and $V^{\mu}$ is the baryon diffusion current. We neglect the effects of bulk viscosity in this study for the purpose of accelerating computing speed. Earlier studies indicate the bulk viscosity may reduce the radial flow of the QGP Ryu:2015vwa ; Ryu:2017qzn . This is compensated in our work by adjusting the initial proper time $\tau_{0}$ to reproduce the mean transverse momenta of the final state charged particles.

Two model parameters embedded in the $\pi^{\mu\nu}$ and $V^{\mu}$ terms above are the specific shear viscosity $C_{\eta_{\mathrm{v}}}$ and the baryon diffusion coefficient $\kappa_{\mathrm{B}}$. They are related to the shear viscosity $\eta_{\mathrm{v}}$, the baryon chemical potential $\mu_{\mathrm{B}}$ and the medium temperature $T$ via

We set $C_{\eta_{\mathrm{v}}}=0.08$ and $C_{\mathrm{B}}=0.4$ through our calculations. The relaxation times are given by $\tau_{\pi}=5C_{\eta_{\mathrm{v}}}/T$ and $\tau_{V}=C_{\mathrm{B}}/T$.

We employ the NEOS-BQS equation of state Monnai:2019hkn ; Monnai:2021kgu for hydrodynamic evolution, which is based on the lattice QCD calculation at high temperature and vanishing net baryon density, and extended to finite net baryon density using the Taylor expansion. At lower energy density, it transits into the equation of state of hadron gas via a smooth crossover. Detailed discussions on the CLVisc model can be found in Ref. Wu:2021fjf .

When the QCD matter becomes sufficiently dilute, quark and gluon degrees of freedom would become confined into hadrons again, and the hydrodynamic description of the bulk evolution should be switched to a transport description of hadronic rescatterings.

On the hadronization hypersurface (chosen as energy density $e_{\mathrm{frz}}=0.26$ GeV/fm³ in this work), we use the Cooper-Frye formula to obtain the distributions of different hadrons with respect to transverse momentum ($p_{\mathrm{T}}$), azimuthal angle ($\phi$) and rapidity ($Y$) as:

where $h$ denotes the hadron species, $g_{h}$ is its spin degeneracy factor, and $d\Sigma_{\mu}$ is the 3-D hypersurface element inside the 4-D spacetime determined by the Cornelius routine Nabi:2012edo . The thermal equilibrium distribution $f_{\mathrm{eq}}(x,p)$ and its out-of-equilibrium corrections $\delta f_{\pi}(x,p)$ and $\delta f_{V}(x,p)$ can be calculated using thermal quantities from the hydrodynamic model as McNelis:2021acu :

Here, $T_{\mathrm{f}}$ represents the chemical freeze-out temperature corresponding to the $e_{\mathrm{frz}}$ we use, and $B$ is the baryon number of identified particle.

Based on the Cooper-Frye formalism above, we use the Monte-Carlo method to sample hadrons out of the QGP medium and then feed them into the SMASH SMASH:2016zqf ; Schafer:2019edr ; Mohs:2019iee ; Hammelmann:2019vwd ; Mohs:2020awg model for simulating their subsequent scatterings in the hadronic phase. SMASH solves the relativistic Boltzmann equation that includes processes of elastic collisions, resonance excitations, string excitations, and decays for hadrons with masses up to about 2 GeV. To improve statistical accuracy, for each hydrodynamic event, we repeat the particle sampling and SMASH simulation 2000 times for Au+Au and Cu+Au collisions. For O+O collisions with smaller sizes and therefore stronger statistical fluctuations, this number is increased to 20000 to ensure statistically stable results.

We start with validating our hydrodynamic calculation with the charged hadron yields per unit pseudorapidity ($dN_{\mathrm{ch}}/d\eta$) as functions of pseudorapidity ($\eta$) in Fig. 1. As shown in panel (a), with the model setup and parameter tuning presented in the previous section, the hydrodynamic calculation can quantitatively describe both the centrality dependence and the pseudorapidity dependence of charged particle production measured by the PHOBOS Collaboration PHOBOS:2005zhy in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV.

In panels (b) and (c), we present the same calculation for Cu+Au and O+O collisions respectively. The centrality classes are still set according to the previous PHOBOS data. It is important to note that for the asymmetric Cu+Au collisions, the center-of-mass of collisions in our computational frame (located at $Y=0$) needs to be shifted towards the Au-moving ($+\hat{z}$) direction by STAR:2017ykf

in order to compare to experimental data. Here, $N_{\mathrm{part}}^{\mathrm{Au}}$ and $N_{\mathrm{part}}^{\mathrm{Cu}}$ are the average numbers of participant nucleons from Au and Cu nuclei respectively, as evaluated from the MC-Glauber model. This is why we observe a slight shift of each curve’s saddle point toward the Au-moving direction in panel (b). This shift becomes more prominent when the imbalance between the energy deposition from Au and Cu nuclei becomes stronger, i.e., in more central collisions. Comparing among the three panels, one can clearly see a decrease in the particle yield as the system size becomes smaller.

We further present in Fig. 2 the multiplicities of both charged hadrons and identified particles as functions of centrality in the three collision systems. The values of $dN_{\mathrm{ch}}/d\eta$ (for charged hadrons) and $dN/dY$ (for identified particles) here are taken from their mid-(pseudo)rapidity regions. As shown in panel (a), our model calculation provides a good description of the currently available data from the PHENIX Collaboration on charged hadrons, pions, kaons, and protons in Au+Au collisions PHENIX:2003iij ; PHENIX:2015tbb . Predictions for Cu+Au and O+O collisions are provided in panels (b) and (c) respectively. Same as the observation from Fig. 1, the particle yield decreases from central to peripheral collisions, and also from Au+Au to Cu+Au and then to O+O collisions. Our prediction on the particle yields in O+O collisions based on the Glauber initial condition is comparable to that based on the more advanced IP-Glasma model Schenke:2020mbo .

The mean transverse momenta $\langle p_{\mathrm{T}}\rangle$ of final state hadrons help quantify the thermal properties of the QGP and its radial flow developed from hydrodynamic expansion. They are strongly correlated with the $p_{\mathrm{T}}$ spectra of hadrons Pratt:2015zsa ; Sangaline:2015isa , and meanwhile, own the advantage of a better visualization on a linear scale than the $p_{\mathrm{T}}$ spectra that are conventionally plotted on a logarithmic scale.

Shown in Fig. 3 are the $\langle p_{\mathrm{T}}\rangle$’s of identified particles as functions of centrality in Au+Au, Cu+Au, and O+O collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV. As seen in panel (a), our calculation reasonably describes the PHENIX PHENIX:2003iij and STAR STAR:2008med data in Au+Au collisions. Nevertheless, there is a hint of overestimation of the proton $\langle p_{\mathrm{T}}\rangle$ here, possibly due to the neglect of bulk viscosity in hydrodynamic evolution in our current work. This may lead to an underestimate of the proton $v_{2}$ later. Calculations using the IP-Glasma or Trento initial condition models and taking into account the bulk viscosity can be found in Refs. Schenke:2019ruo ; Schenke:2020mbo ; JETSCAPE:2020mzn . The mass ordering in $\langle p_{\mathrm{T}}\rangle$ can be clearly observed in the figure: being produced from the same QGP medium, heavier hadrons acquire larger $p_{\mathrm{T}}$ from the thermal background than lighter hadrons do. Meanwhile, compared to lighter hadrons, it is easier for heavier hadrons to gain additional $p_{\mathrm{T}}$ from the radial flow of the medium, which decreases as centrality increases. Therefore, the $\langle p_{\mathrm{T}}\rangle$ of heavier hadrons has a stronger dependence on the centrality of heavy-ion collisions. Comparing among the three panels, we also observe a weaker centrality dependence of the hadron $\langle p_{\mathrm{T}}\rangle$ in O+O collisions than in Au+Au and Cu+Au collisions. This results from the weak radial flow developed in the small-size O+O system, even in its central collisions. In peripheral collisions, the same species of hadrons exhibit similar magnitudes of $\langle p_{\mathrm{T}}\rangle$ across different collision systems, since the effect of radial flow is negligible in peripheral collisions and the $\langle p_{\mathrm{T}}\rangle$ is mainly determined by the hadronization temperature of the QCD medium at $\sqrt{s_{\mathrm{NN}}}=200$ GeV. Note that this $\langle p_{\mathrm{T}}\rangle$ may depend on beam energy of heavy-ion collisions Wu:2021fjf ; Shen:2020jwv ; Du:2023efk considering the varying boundary between QGP and hadron gas in the QCD phase diagram as $\sqrt{s_{\mathrm{NN}}}$ changes.

As discussed in the introduction section, to avoid the difficulty in determining the event plane in realistic experimental measurements, the multi-particle correlation method is usually preferred in evaluating the collective flow coefficients in heavy-ion collisions.

Consider $m$ is an positive integer, the $n$-th order $2m$-particle azimuthal correlator is defined as Borghini:2001vi

where the inner layer of angle bracket denotes an average over all possible combinations of particles within an event, and the outer layer denotes an average across different events. The corresponding 2-, 4-, and 6-particle cumulants are then given by

with the 2-, 4-, and 6-particle harmonic coefficients given by

By assuming Gaussian distributions of the collective flow fluctuations, one would obtain Voloshin:2008dg ; Bhalerao:2011yg

where $\langle v_{n}\rangle$ is the magnitude of the average of the $\vec{v}_{n}$ vector (with its direction denoting the event plane angle) in the transverse plane, and $\sigma_{n}$ is the Gaussian width of the flow fluctuations. Therefore, one can use the ratio $v_{n}\{4\}/v_{n}\{2\}$ or $v_{n}\{6\}/v_{n}\{2\}$ to quantify the strength of fluctuation in the collective flow coefficient: larger deviation from one implies stronger fluctuation. Equations (32) to (34) are also valid for non-Gaussian distributions of fluctuations as long as their variances satisfy $\sigma_{n}\ll\langle v_{n}\rangle$.

Since a direct evaluation of Eq. (25) requires 2$m$ loops over the particle list in each event, which is computationally inefficient when $m$ is large, we adopt the $Q_{n}$-vector method developed in Ref. Bilandzic:2010jr to compute these correlators. For an event consisting of $M$ particles within a desired kinematic region, $Q_{n}$ is defined as:

The single-event-averaged 2-, 4-, and 6-particle azimuthal correlators can then be written as:

We then average these results over different events and substitute them into Eqs. (26) to (31) to obtain the collective flow coefficients from the cumulant method.

To study the collective flow fluctuations across systems with different shapes and sizes, we implement event-by-event hydrodynamic simulations of Au+Au, Cu+Au, and O+O collisions in this work. In each centrality bin, 500 events are simulated for Au+Au collisions. To enhance statistical accuracy in smaller systems, this number is increased to 1000 for Cu+Au collisions and 1500 for O+O collisions per centrality bin. Additionally, following the STAR Collaboration work STAR:2022gki , we adjust the smallest centrality bin from 0-5% to 1-5% for each collision system to avoid possible positive values of $c_{n}\{4\}$ in very central (0-1%) collisions, considering that geometric fluctuations in ultra-central collisions can be very strongZhou:2018fxx .

In the upper panels of Fig. 4, we first present the collective flow coefficients evaluated from the 2-particle [panel (a)] and 4-particle [panel (b)] cumulant methods. The kinematic cuts are chosen as $|\eta|<1$ and $p_{\mathrm{T}}\in(0.2,4.0)$ GeV/$c$ according to experiment. In general, our calculation provides a reasonable description of the STAR data STAR:2022gki on $v_{2}\{2\}$ and $v_{2}\{4\}$ in both Au+Au and Cu+Au collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV, except for some deviation in $v_{2}\{4\}$ in peripheral Cu+Au collisions. Comparing between these two systems, we observe a slightly larger $v_{2}$ in Cu+Au than in Au+Au in the most central collisions. This results from stronger initial state fluctuations in a smaller system. On the other hand, at larger centrality where $v_{2}$ is dominated by the average medium geometry, it is larger in Au+Au than in Cu+Au collisions. We have verified that at centrality greater than 10%, the eccentricity of Au+Au collisions is larger than that of Cu+Au collisions. Predictions for O+O collisions are also presented here for comparison. Due to the small system size, $v_{2}$ from O+O collisions is mainly driven by fluctuations, and therefore its value is larger in central collisions while smaller in peripheral collisions compared to those seen in Au+Au and Cu+Au collisions. Since the number of particles produced in peripheral O+O collisions is limited, the statistics of its $v_{2}\{4\}$ constructed from 4-particle correlation becomes poor at large centrality.

Shown in the lower panels of Fig. 4 are $v_{2}\{4\}/v_{2}\{2\}$ [panel (c)] and $v_{2}\{6\}/v_{2}\{4\}$ [panel (d)] of charged particles in different collision systems. For Au+Au and Cu+Au collisions, we observe the values of $v_{2}\{4\}/v_{2}\{2\}$ are significantly smaller than one in central collisions, indicating the dominant effect of fluctuations on forming $v_{2}$ in the corresponding region. As centrality increases, this ratio first increases towards one as the average geometry of the medium starts to dominate, and then slightly decreases again as the system becomes small and effect of fluctuations grows again at peripheral collisions. In the mid-central to semi-peripheral region, the eccentricity of the initial state geometry is smaller in Cu+Au than in Au+Au collisions at the same centrality. Meanwhile, fluctuations in the former is stronger due to its smaller system size. As a result, $v_{2}\{4\}/v_{2}\{2\}$ is smaller in Cu+Au than in Au+Au collisions within this centrality region. Fluctuations in the small O+O system is strong across the entire centrality region and therefore the corresponding $v_{2}\{4\}/v_{2}\{2\}$ keeps significantly below one. Here, we can clearly observe the strong dependence of $v_{2}\{4\}/v_{2}\{2\}$ on the size and geometry of the collision system. This is different from the insensitivity of this ratio to the beam energy within the same collision system, as seen in earlier theoretical calculation Wu:2021fjf and experimental data STAR:2022gki . The $v_{2}\{6\}/v_{2}\{4\}$ ratios agree with one for both Au+Au and Cu+Au collisions in panel (d). Possible slight deviation from one in very central and very peripheral regions may result from non-Gaussian form of strong fluctuations there. Because of very limited statistics in O+O collisions for constructing $v_{2}\{6\}$, result for this system is not presented in panel (d).

In Fig. 5, we present $v_{2}\{2\}$, $v_{2}\{4\}$, and $v_{2}\{4\}/v_{2}\{2\}$ for identified particles in the three collision systems at $\sqrt{s_{\mathrm{NN}}}=200$ GeV. Compared to the available data from the STAR Collaboration STAR:2022gki , we provide a good description of $v_{2}\{2\}$ and $v_{2}\{4\}$ of pions and kaons, while slightly underestimating them for protons in Au+Au collisions. Since our previous result on the proton $\langle p_{\mathrm{T}}\rangle$ is close to the upper edges of the data error bars in Fig. 3 (a), selecting the $p_{\mathrm{T}}\in(0.2,2.0)$ GeV/$c$ range here, as used in measurements, excludes protons with $p_{\mathrm{T}}\gtrsim 2.0$ GeV/$c$ from our model which can contribute a larger $v_{2}$ to its average value. It is interesting to note that the mass hierarchy of the $p_{\mathrm{T}}$-averaged $v_{2}$ in Au+Au and Cu+Au collisions here – heavier hadrons show stronger $v_{2}$ – is opposite to what one usually sees in the $p_{\mathrm{T}}$-dependent $v_{2}$ of identified hadrons below 2 GeV/$c$ PHENIX:2006dpn ; ALICE:2014wao ; ALICE:2022zks ; STAR:2022ncy . This is because of the harder $p_{\mathrm{T}}$ spectra of heavier hadrons which add more weights from the higher $p_{\mathrm{T}}$ region to their average $v_{2}$. In other words, one actually compares the $v_{2}$ of heavier hadrons with higher average $p_{\mathrm{T}}$ to the $v_{2}$ of lighter hadrons with lower average $p_{\mathrm{T}}$ in the $p_{\mathrm{T}}$-averaged $v_{2}$ here. Although different species of hadrons show different magnitudes of $v_{2}\{2\}$ and $v_{2}\{4\}$, they have similar $v_{2}\{4\}/v_{2}\{2\}$ ratios in Au+Au and Cu+Au systems, except in very central and very peripheral collision. This indicates these hadrons keep the memory of the momentum space fluctuations of the QGP, and are not strongly affected by additional fluctuations during the hadronization and hadronic afterburner processes. Contrarily, $v_{2}\{4\}/v_{2}\{2\}$ is no longer independent of hadron species in O+O collisions. Although the statistical errors in our results for O+O collisions are large, the dependence of $v_{2}\{4\}/v_{2}\{2\}$ on hadron species is visible across a wide centrality range. This dependence can result from the fluctuations in hadronization and hadronic afterburner, which become more prominent for a smaller system. These non-initial-state fluctuations may also be important in very central Au+Au and Cu+Au collisions, where the collective flows are dominated by fluctuations. To better compare $v_{2}\{4\}/v_{2}\{2\}$ between different hadron species, we show its ratio between kaons and pions, and between protons and pions in the sub-panels in the last column.

Earlier studies Magdy:2018itt ; STAR:2022gki ; Alba:2017hhe ; Schenke:2019ruo ; Magdy:2020gxf ; Magdy:2021sba ; Rao:2019vgy suggest that the collective flow fluctuations, quantified by $v_{2}\{4\}/v_{2}\{2\}$, predominantly originates from the fluctuations of the initial state geometry in large nuclear collision systems, such as Au+Au and Pb+Pb collisions. To characterize the anisotropy of the initial state geometry, one may evaluate its $n$-th order eccentricity as

where $r$ and $\phi$ are the position and azimuthal angle of an element of nuclear matter with respect to its center of mass, and the angle brackets denote average within an event. In this work, both participant nucleons and binary collision points sampled from the MC-Glauber model contribute to the elements in this average, with $\alpha$ as their relative weight as shown in Eq. (3).

Within the hypothesis of linear hydrodynamic response, $v_{n}\propto\varepsilon_{n}$, one may construct the cumulants of eccentricities according to Eqs. (26) to (28)  as Qiu:2011iv ; Ma:2016hkg :

where the angle brackets denote average over different events. The $n$-th order eccentricities defined by the cumulant method then read:

Analogous to the collective flow coefficients determined through the cumulant method, the cumulant eccentricities here encode information of fluctuations in the initial state geometry of the QGP. Therefore, similar to using $v_{n}\{4\}/v_{n}\{2\}$ to measure the fluctuations in the final state collective flow, one can use $\varepsilon_{n}\{4\}/\varepsilon_{n}\{2\}$ to quantify the strength of the initial state fluctuations.

To study the correlation between the initial state geometric fluctuation and the final state flow fluctuation, we compare $\varepsilon_{2}\{4\}/\varepsilon_{2}\{2\}$ and $v_{2}\{4\}/v_{2}\{2\}$ in the upper panel of Fig. 6, and compare $\varepsilon_{2}\{6\}/\varepsilon_{2}\{4\}$ and $v_{2}\{6\}/v_{2}\{4\}$ in the lower panel, for Au+Au, Cu+Au, and O+O collisions at $\sqrt{s_{\mathrm{NN}}}=200$ GeV. The cumulant eccentricities $\varepsilon_{2}\{k\}$ here are obtained from averaging over 50000 MC-Glauber profiles in each centrality bin of each collision system. In the upper three panels, we observe a monotonic increase of $\varepsilon_{2}\{4\}/\varepsilon_{2}\{2\}$ for the three systems from central to peripheral collisions. This indicates, in central collisions, the geometric anisotropy in the initial state is mainly contributed by fluctuations, while in peripheral collisions, it is mainly determined by the average shape of the overlapping region between the two nuclei. In Au+Au and Cu+Au collisions, $\varepsilon_{2}\{4\}/\varepsilon_{2}\{2\}$ agrees well with $v_{2}\{4\}/v_{2}\{2\}$ except in very central and very peripheral regions. This suggests the initial state fluctuations are the main source of the final state collective flow fluctuations in these large enough collision systems. When the centrality is large, the systems become so small that effects of additional fluctuations, e.g. nonlinear hydrodynamic response, hadronization, and hadronic afterburner, become important, leading to smaller values of $v_{2}\{4\}/v_{2}\{2\}$ than $\varepsilon_{2}\{4\}/\varepsilon_{2}\{2\}$. As revealed in Refs. Noronha-Hostler:2015dbi ; Giacalone:2017uqx , nonlinear hydrodynamic response has a significant contribution to the elliptic flow fluctuations in peripheral collisions. For the small O+O system, sources other than the initial state always have strong contributions to the flow fluctuations, and thus $v_{2}\{4\}/v_{2}\{2\}<\varepsilon_{2}\{4\}/\varepsilon_{2}\{2\}$ for the whole centrality range.

In the lower three panels of Fig. 6, we observe the values of $\varepsilon_{2}\{6\}/\varepsilon_{2}\{4\}$ are close to one in all the three systems. No obvious deviation is seen between $\varepsilon_{2}\{6\}/\varepsilon_{2}\{4\}$ and $v_{2}\{6\}/v_{2}\{4\}$ in Au+Au and Cu+Au collisions, except for some hints of slight deviation in very central and very peripheral regions. The result of $v_{2}\{6\}/v_{2}\{4\}$ in O+O collisions is not available due to the limited statistics of hadrons produced in these small systems.