<figure><img src="image_1.png" /></figure>

HAL
open science

# Elliptic flow fluctuations in central collisions of spherical and deformed nuclei

Giuliano Giacalone

► To cite this version:

Giuliano Giacalone. Elliptic flow fluctuations in central collisions of spherical and deformed nuclei. Phys.Rev.C, 2019, 99 (2), pp.024910. 10.1103/PhysRevC.99.024910. hal-01937770

HAL Id: hal-01937770

https://hal.science/hal-01937770

Submitted on 6 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Elliptic flow fluctuations in central collisions of spherical
and deformed nuclei

Giuliano Giacalone

Institut de physique théorique, Université Paris-Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France

<figure><img src="image_3.png" /></figure>

(Received 16 November 2018; published 19 February 2019)

Elliptic flow ($v_2$) fluctuations in central heavy-ion collisions are direct probes of the fluctuating geometry of the quark-gluon plasma and, as such, are strongly sensitive to any deviation from spherical symmetry in the shape of the colliding nuclei. I investigate the consequences of nuclear deformation for $v_2$ fluctuations, and I assess whether current models of medium geometry are able to predict and capture such effects. Assuming linear hydrodynamic response between $v_2$ and the eccentricity of the medium, $\varepsilon_2$, I perform accurate comparisons between model calculations of $\varepsilon_2$ fluctuations and STAR data on cumulants of elliptic flow in central Au+Au and U+U collisions. From these comparisons, I evince that the most distinct signatures of nuclear deformation appear in the non-Gaussianities of $v_2$ fluctuation, and I show, in particular, that the non-Gaussian $v_2$ fluctuations currently observed in central Au+Au collisions are incompatible with model calculations that implement a quadrupole coefficient of order 12% in the <sup>197</sup>Au nuclei. Finally, I make robust predictions for the behavior of higher order cumulants of $v_2$ in collisions of nonspherical nuclei.

DOI: 10.1103/PhysRevC.99.024910

# I. INTRODUCTION

Elliptic flow is the dynamical response of a fluid to an el-
liptic deformation of its geometry. It is a salient feature of the
hydrodynamic expansion of the quark-gluon plasma created
in relativistic heavy-ion collisions, where elliptic anisotropy is
generated, on the one hand, by the fact that the area of overlap
of two nuclei at finite impact parameter looks like an ellipse
[1], which explains why elliptic flow grows quickly with the
impact parameter of the collision, and on the other hand, by
density fluctuations in the fluid [2,3], that explain the striking
emergence of elliptic flow in collisions at small impact param-
eter [4]. Following Teaney and Yan [5], the elliptic anisotropy
of the medium, dubbed $\varepsilon_2$, can be defined rigorously for a
generic heavy-ion collision, and hydrodynamic simulations
show that elliptic flow ($v_2$) is essentially a linear response to
$\varepsilon_2$, i.e., $v_2 = \kappa\varepsilon_2$, where $\kappa$ is a constant [6–8]. This simple
relation implies that $v_2$ and its event-by-event fluctuations can
be used as direct probes of the fluctuating geometry of the
quark-gluon plasma at the beginning of the hydrodynamic
phase.

In central collisions, where linear hydrodynamic response
is exhibited at its best [9,10], an important source of fluctua-
tions that contributes to $\varepsilon_2$ is given by the random orientation
of the colliding nuclei, if they are nonspherical. This explains
the large magnitude of the rms elliptic flow measured in
central collisions of nuclei that have a pronounced deforma-
tion, i.e., U+U collisions at the BNL Relativistic Heavy Ion
Collider (RHIC) [11] and Xe+Xe collisions at the CERN
Large Hadron Collider (LHC) [12]. In the current model-
ing of initial conditions for hydrodynamics, the nonspheri-
cal shape of the nuclei is obtained by adding a quadrupole
deformation in the wave functions of the colliding bodies,
the implementation of which is performed with the guid-

ance of tabulated data on nuclear ground-state deformations
[13,14]. This kind of modeling allows hydrodynamic simu-
lations to reproduce quantitatively the aforementioned large
rms elliptic flow observed in U+U [15] and Xe+Xe [16,17]
collisions.

In this paper, I argue that more interesting and nontrivial
signatures of nuclear deformation can be observed in the
non-Gaussian fluctuations of $v_2$, and that such effects are
currently visible in the experimental data, in particular, in
the fourth-order cumulant of elliptic flow $v_2\{4\}$, accurately
measured by the STAR Collaboration [11] in central colli-
sions of nonspherical nuclei, i.e., <sup>197</sup>Au + <sup>197</sup>Au collisions
and <sup>238</sup>U + <sup>238</sup>U collisions. Note that the present study fills an
important gap in the literature: Theoretical studies devoted to
the consequences of the prolate shape of <sup>238</sup>U nuclei for rel-
evant observables are numerous in the literature [15,18–28],
but a careful assessment of the impact of nuclear deformation
on the non-Gaussianities of $v_2$ fluctuations, which are central
to the phenomenology of flow in heavy-ion collisions [29–32],
is still missing.

I first investigate how nuclear deformation affects the
fluctuations of eccentricity in Monte Carlo simulations and
whether these effects help explain the experimental data. To
achieve this, I use state-of-the-art Monte Carlo models of
initial conditions (described in Sec. II) to perform extensive
calculations of $\varepsilon_2$ fluctuations, that I rescale and compare
(Sec. III) to STAR data on cumulants of elliptic flow fluc-
tuations, $v_2\{2\}$ and $v_2\{4\}$, in central Au+Au and U+U col-
lisions. The outcome of these comparisons is eventually used
(Sec. IV) to explain the prominent consequences of nuclear
deformation for the non-Gaussian fluctuations of elliptic flow,
and I predict how these effects can be observed in the higher
order cumulants of $v_2$.
## II. MODELS OF INITIAL GEOMETRY

In this section, I aim at exhibiting models of initial conditions that are viable for hydrodynamic simulations of Au+Au and U+U collisions at $\sqrt{s} = 200$ GeV. A model of initial conditions is a prescription that provides, event by event, the energy or entropy density, say $S(x, y)$, deposited in the transverse plane (for simplicity, at midrapidity) by a given collision. This is by far the most crucial ingredient in the hydrodynamic framework, because it provides the medium with its spatial anisotropies,<fnref n="1" /> that get eventually converted into momentum anisotropies (i.e., the Fourier coefficients $v_n$) through the hydrodynamic evolution.

In this paper, I shall employ two models of initial conditions that correspond to two parametrizations of $T_R$ENTo [28], a model for the profile of entropy density deposited at midrapidity in heavy-ion collisions. Let me provide a detailed explanation of these models.

### A. $T_R$ENTo initial conditions

The starting point is the modeling of the colliding bodies, which one performs through a random sampling of nucleons. For a generic nonspherical nucleus, the sampling is done using the following Fermi distribution [33],

$$ \rho(x, y, z) = \frac{\rho_0}{1 + \exp\left\{-\frac{1}{a}r(-R[1 + \beta_2 Y_{20} + \beta_4 Y_{40}])\right\}}, \quad (1) $$

where $r = \sqrt{x^2 + y^2 + z^2}$, $a$ and $R$ are the skin depth and radius of the nucleus, respectively, and the coefficients $\beta_2$ and $\beta_4$ multiplying the spherical harmonics, $Y_{20} = \sqrt{\frac{5}{16\pi}}(3 \cos^2 \theta - 1)$ and $Y_{40} = \frac{3}{16\sqrt{\pi}}(35 \cos^4 \theta - 30 \cos^2 \theta + 3)$, are the coefficients of quadrupole and hexadecapole deformation of the nucleus, respectively. Once the transverse coordinates of the nucleons, $(x, y)$, are sampled, the nuclei are overlapped at a random impact parameter, and nucleon-nucleon collisions take place. One simply assumes that a participant nucleon deposits entropy, $s$, according to a Gaussian distribution,

$$ s_i(x, y) = \frac{\omega_i}{2\pi\sigma^2} \exp\left[-\frac{(x-x_i)^2 + (y-y_i)^2}{2\sigma^2}\right], \quad (2) $$

where I take $\sigma = 0.5$ fm for the width of each participant, and the index $i$ refers to the $i$th participant nucleon. The normalization, $\omega_i$, is randomly distributed according to the following gamma distribution,

$$ P(\omega) = \frac{k^k \omega^{k-1} e^{-k}}{\Gamma(k)}, \quad (3) $$

which has mean value equal to unity and variance proportional to $k^{-1}$. Hence, the total entropy profile of a given nucleus, say $A$, is given by

$$ S_A(x, y) = \sum_i s_i. \quad (4) $$

Eventually, for two colliding nuclei, $A$ and $B$, the total entropy profile of the system is given by a generalized mean, i.e.,

$$ S(x, y) = \left( \frac{S_A^p + S_B^p}{2} \right)^{1/p}, \quad (5) $$

where $p$ is any real number.

In this paper, I shall use $S(x, y)$ from two different $T_R$ENTo parametrizations. I take a geometric mean,

$$ S = \left( \frac{S_A^p + S_B^p}{2} \right)^{1/p} \bigg|_{p=0} = \sqrt{S_A S_B}, \quad (6) $$

and an arithmetic mean,

$$ S = \left( \frac{S_A^p + S_B^p}{2} \right)^{1/p} \bigg|_{p=1} = \frac{S_A + S_B}{2}. \quad (7) $$

Let me provide a bit of motivation for these choices.

The case $p=0$ is very successful in phenomenological applications and is the favored value of $p$ resulting from the extensive Bayesian analyses of Refs. [34,35]. The reason for this success is likely the fact that $S$ in this model is proportional to the product $S_A S_B$. This is reminiscent of a class of models inspired by high-energy QCD [9,36,37]. Taking $S(x, y)$ as a product leads typically to systems whose eccentricity follows closely the almond shape of the nuclear overlap [38], leading eventually to a very good description of elliptic flow data [28,39].

The case with $p=1$ corresponds instead to a wounded nucleon model with participant nucleon scaling [40]. This class of models are variations of the original Monte Carlo Glauber [33] model and are typically the models employed in experimental analyses.<fnref n="2" /> Taking $S(x, y)$ as the sum of two components leads to systems that are more scattered in the transverse plane and that have less enhanced eccentricity [38]. The $p=1$ model is essentially ruled out by elliptic flow fluctuations data in Pb+Pb collisions at the LHC [39]. Nevertheless, such a model has never been compared to RHIC data, and RHIC data was not used in the Bayesian analyses of Refs. [34,35]. I deem, then, that one can not a priori state that the $p=1$ model is ruled out as well at RHIC energies.

The bottom line, and I shall come back to this point later, is that these two parametrizations of $T_R$ENTo capture, arguably, all the basic features of the widest classes of initial condition models for nucleus-nucleus collisions that are on the market.

### B. Multiplicity

I discuss now the implementation of the parameter $k$ that regulates the fluctuations of the entropy produced by each participant nucleon. This feature is important for a correct

<footnote n="1">¹Note that for the definition of the eccentricity of the medium, $\epsilon_2$, it does not matter whether one considers anisotropy in the energy density or in the entropy density of the fluid [6].</footnote>

<footnote n="2">²For instance, the model calculations shown by the STAR Collaboration in Ref. [11], i.e., a Glauber-Monte Carlo model with binary collision scaling, and a constituent quark Glaber model, are essentially mild variations of the $p=1$ model for what concerns the medium geometry and should lead to the same kind of $\epsilon_2$ fluctuations. This will be confirmed in Sec. III.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 1. Symbols: Distribution of multiplicity measured by the STAR Collaboration [11] in Au+Au collisions at $\sqrt{s} = 200$ GeV (a) and U+U collisions at $\sqrt{s} = 193$ GeV (b). Lines are model calculations: IP-Glasma simulations [26] (dotted line) and TRENTo with $p = 1$ (solid line) and $p = 0$ (dashed line). The vertical lines indicate the knees of the histograms.</figcaption></figure>

definition of collision centrality in the models. In experiments, centrality classes are defined from the histogram of the multiplicity of particles (or hits, or energy measured in a calorimeter). In hydrodynamics, one typically assumes that the total entropy at the initial condition is, in each event, proportional to multiplicity of particles in the final state. Therefore, in numerical simulations one uses the histogram of total entropy to define the classes of collision centrality. In order to exhibit meaningful model-to-data comparisons, I need to obtain in the models a distribution of total entropy which is in agreement with the distribution of $dN_{ch}/d\eta$ used by the STAR Collaboration [11] to sort their events into centrality bins. In particular, since I shall deal with central collisions, I want the high-multiplicity tails of the measured multiplicity distributions to be captured by the models. This can be achieved by a proper choice of the fluctuation parameter, $k$.

In Fig. 1, I display as circles the distributions of $dN_{ch}/d\eta$ in Au+Au and U+U collisions measured by the STAR Collaboration.<fnref n="3" /> Note that multiplicity distributions collected in the STAR detector present a high-multiplicity tail that is twice as broad as that measured in detectors at the LHC [41]. This is consistent with the fact that detectors at RHIC have a much smaller acceptance, and therefore they are more sensitive to statistical fluctuations. Since in TRENTo the variance of fluctuations scales essentially as $1/k$, I expect that one needs to implement a lower value for the parameter $k$ at RHIC than

at LHC. The choice of this parameter is, essentially, detector dependent.

Let me first discuss the case of the TRENTo model with $p = 0$. In Ref. [16], this model is found to provide an excellent description of ALICE data when $k = 2.0$. To describe RHIC data, i.e., a large-multiplicity tail which is twice as broad, I use $k = 0.5$, which is smaller by a factor of 4, consistent with the fact that the width of the $\Gamma$ distribution is proportional to $k^{-1/2}$. The probability distribution of the total entropy in this model is reported as a dashed line, for both Au+Au and U+U collisions, in Fig. 1. To show a meaningful comparison between the distribution of entropy provided by TRENTo and the measured multiplicities, I rescale the entropy in TRENTo in order to have the knee of the histogram at the same coordinate as the knee of the histogram of experimental data.<fnref n="4" /> Agreement with the data is reasonable,<fnref n="5" /> especially in the high-multiplicity tail, although the comparison is not as excellent as the one observed with LHC data.<fnref n="6" />

Moving to TRENTo with $p = 1$, I find that also in this case an excellent description of data is achieved using $k = 0.5$. The rescaled distributions of entropy in this model are shown as solid lines in Fig. 1. Note that, for most of the histograms, the description provided by this model is better than that observed with $p = 0$.

<footnote n="3">³These distributions can be obtained from the parametrizations of $dN_{ch}/d\eta$ vs centrality (c) provided at the beginning of Ref. [11]. Knowing $dN/d\eta(c)$, and using the fact that centrality is defined as the cumulative of the multiplicity, i.e., $c(dN/d\eta) = \int_{dN/d\eta}^{\infty} P(dN/d\eta)$, the plots of $P(dN/d\eta)$ shown in Fig. 1 can be simply obtained by plotting $(\frac{d(dN/d\eta)}{dc})^{-1}|_c$ vs $dN/d\eta|_c$, with $c \in [0, 1]$.</footnote>

<footnote n="4">⁴The knee is defined as the mean value of multiplicity at zero impact parameter, and for the experimental data I calculate it through the fitting procedure of Ref. [41]. From the STAR parametrizations, I find $dN_{ch}/d\eta|_{knee} = 663$ in Au+Au collisions, and $dN_{ch}/d\eta|_{knee} = 821$ in U+U collisions.</footnote>

<footnote n="5">⁵I have also checked that the model presents the same centrality of the knee, $c_{knee}$, of data, i.e., the area of the histogram on the right of the knee. From the STAR parametrizations, I find $c_{knee} = 0.81\%$ in Au+Au collisions, and $c_{knee} = 0.75\%$ in U+U collisions.</footnote>

<footnote n="6">⁶I have actually tried several values of $k$, and agreement does not get better.</footnote>
I would like to expand a bit further on this point, which is very striking. The model with $p = 0$ allows one to accurately reproduce LHC multiplicity data [16], but here it yields a poor description of RHIC data at multiplicities that are smaller than the location of the knee of the histogram. At such multiplicities, one cannot simply improve agreement with data via a change in $k$: The model is simply unable to reproduce the data. Another model that is able to reproduce with great accuracy the distributions of multiplicity measured at LHC is the IP-Glasma model [36], which was employed in calculations of Au+Au and U+U collisions in Ref. [26]. Rescaled multiplicity distributions in this model are shown as dotted lines in Fig. 1. Note that they are essentially compatible with the curves of $T_{RENTO}$ $p = 0$, and consequently, this model provides a rather bad description of Au+Au data before the knee of the histogram.<fnref n="7" /> Therefore, either the difference between multiplicity distributions at RHIC and at LHC is due to some unknown issue due to the different detectors used or there might actually be some physical motivation for which entropy production within IP-Glasma-like models (e.g., $T_{RENTO}$ $p = 0$) provides a better description of data at the TeV energy scale. Moreover, let me emphasize that the same issue might jeopardize future Bayesian analyses that, using the $T_{RENTO}$ model, are aimed at a determination of the parameters $k$ and $p$ via a simultaneous fit of RHIC and LHC data. My results here indicate that it is very unlikely that, using the same values of $k$ and $p$, one may be able to obtain a good description of both RHIC and LHC data.

That being said: I have exhibited two $T_{RENTO}$ parametrizations whose distributions of total entropy in Au+Au and U+U collisions provide a good description of the multiplicity distributions measured by the STAR Collaboration. Therefore, these are viable models of initial conditions for hydrodynamic calculations, and I use them to compute the eccentricity of the medium in each event, as discussed below.

### C. Eccentricity fluctuations

The elliptic anisotropy of a smooth profile in two dimensions, for instance, the entropy density $S(x, y)$ given in each collision by the $T_{RENTO}$ model, can be computed as indicated by Teaney and Yan [5] (in polar coordinates)<fnref n="8" />:

$$ \mathcal{E}_2 = \varepsilon_2 e^{2i\psi_2} = -\frac{\int r^2 e^{i2\phi} S(r, \phi) r dr d\phi}{\int r^2 S(r, \phi) r dr d\phi}. \quad (8) $$

This complex quantity fluctuates in both magnitude ($\varepsilon_2$) and orientation ($\psi_2$) in each event. Linear hydrodynamic response implies that $\mathcal{E}_2$ is linearly correlated with the complex elliptic

flow coefficient, $V_2 = v_2 e^{2i\Phi_2}$, that also fluctuates in magnitude ($v_2$) and orientation ($\Phi_2$) event to event. If this is the case, the probability distribution of $V_2$ coincides then with that of $\mathcal{E}_2$ up to a factor, and the statistical properties of $V_2$ fluctuations provide direct information about the fluctuations of the initial $\mathcal{E}_2$ [43]. Let me recall the formulas of the first two cumulants of the $v_2$ distribution,

$$ v_2\{2\} = \sqrt{\langle v_2^2 \rangle}, \quad v_2\{4\} = 4\sqrt{2\langle v_2^2 \rangle^2 - \langle v_2^4 \rangle}, \quad (9) $$

where brackets indicate an average over events in a centrality bin. Now, if $v_2 = \kappa\varepsilon_2$, then,

$$ v_2\{2\} = \kappa\varepsilon_2\{2\}, \quad v_2\{4\} = \kappa\varepsilon_2\{4\}. \quad (10) $$

The ratio $v_2\{4\}/v_2\{2\}$ has been used in many studies to observe, in hydrodynamic simulations, the transition between the linear regime, where $v_2\{4\}/v_2\{2\} = \varepsilon_2\{4\}/\varepsilon_2\{2\}$, and the nonlinear regime where this equality breaks down [39,43–45]. In Ref. [39], in particular, the ratio $\varepsilon_2\{4\}/\varepsilon_2\{2\}$ was computed in extensive Monte Carlo calculations of $\varepsilon_2$ in very central collisions, and directly compared to LHC Pb+Pb data, in order to test the validity of different $T_{RENTO}$ parametrizations.

In what follows, I repeat this game, although without taking any ratio: I compute $\varepsilon_2\{2\}$ and $\varepsilon_2\{4\}$ in the $T_{RENTO}$ models, and then rescale them by an appropriate factor $\kappa$, in order to find the best possible agreement with experimental data on $v_2\{2\}$ and $v_2\{4\}$ in central Au+Au and U+U collisions. As anticipated, eccentricity fluctuations in central collisions are strongly sensitive to the event-by-event fluctuations of the spatial orientation of the colliding $^{197}$Au and $^{238}$U nuclei, i.e., to the choice of the deformation parameters that enter in Eq. (1). The goal of performing Monte Carlo calculations and model-to-data comparisons is essentially twofold. First, I want to understand how $\varepsilon_2$ fluctuations are affected by the presence of deformed nuclei in the models. Second, I want to check whether these effects predicted by the models are in agreement with the features of $v_2$ fluctuations observed in the data.

## III. COMPARISON WITH STAR DATA

### A. Setup

Using the $T_{RENTO}$ parametrizations with $p = 0, k = 0.5$, and $p = 1, k = 0.5$, I simulate Au+Au and U+U collisions at $\sqrt{s} = 200$ GeV. The nuclear shape parameters that enter in Eq. (1) are implemented as follows: I take $^{238}$U nuclei with $R = 6.80$ fm, $a = 0.60$ fm, $\beta_2 = 0.236$, and $\beta_4 = 0.098$. Concerning $^{197}$Au nuclei, I use $R = 6.40$ fm and $a = 0.53$ fm, and I shall run calculations for both spherical nuclei, i.e., $\beta_2 = \beta_4 = 0$, and deformed nuclei with $\beta_2 = -0.125$, $\beta_4 = -0.017$,

A couple of comments are in order. The chosen values of radii and skin depths are rather standard. I do not follow the suggestion of Ref. [27] and do not set $a \sim 0.4$ in both $^{238}$U and $^{197}$Au nuclei. As suggested by the results shown in that same reference, this choice has essentially no impact on eccentricity fluctuations in central collisions. Moving on to the deformation parameters, for $^{238}$U nuclei I implement a value of $\beta_2$ that is smaller than in previous studies, as I

<footnote n="7">⁷Note that the IP-Glasma results shown in Fig. 1 present as well a large-multiplicity tail which is twice as steep as the STAR data, suggesting that they miss the detector-dependent part of the fluctuations of multiplicity. In the IP-Glasma formalism, this issue could be solved via the inclusion of fluctuations of the saturation scale at the level of the colliding nucleons [42].</footnote>

<footnote n="8">⁸The expression for the anisotropy of order 3 (triangularity) is completely analogous [5].</footnote>
take it from the most recent table of nuclear deformations [14]. Moving to the deformation of <sup>197</sup>Au nuclei, I cautiously stress that the parameter $\beta_2$ for this nucleus is not a measured quantity. The data tables resulting from the model calculations of Refs. [13,14] yield $|\beta_2| \sim 0.13$, whereas, in Ref. [46], a value of 0.15 is guessed from the measured deformations of neighbor nuclei. My choice $\beta_2 = -0.125$ is taken again from Ref. [14].

For these $T_{\text{RENTO}}$ setups ($p=0$ and $p=1$), and using both spherical and deformed <sup>197</sup>Au nuclei, I simulate $O(10^7)$ U+U and Au+Au collisions, and I compute the fluctuations of $\epsilon_2$ as function of collision multiplicity that I compare to experimental data.

When rescaling the results for $\epsilon_2$ fluctuations, I shall take the same value of $\kappa$ [Eq. (10)] for both Au+Au and U+U systems: This is a good approximation, as the value of $\kappa$ is reduced solely by viscous corrections. Dimensional analysis indicates that viscous corrections scale like $1/R$, where $R$ is the radius of the system. A good approximation for $R$ is given by $A^{-1/3}$. Therefore, moving from U+U collisions to (smaller) Au+Au collisions, I expect a reduction of $\kappa$ of order $3\sqrt{238}/197$, which is a negligible 5% correction. Let me also stress that, in view of this argument concerning the system size, the value of $\kappa$ is expected to decrease as I move from central to peripheral collisions. This effect should be very small in the centrality range I am interested in, as also indicated by hydrodynamic calculations [10]; therefore, I simply assume that $\kappa$ is a constant.

## B. Model vs data

I show STAR data [11] on cumulants of elliptic flow fluctuations as empty symbols in Fig. 2. Figures 2(a) and 2(b) show STAR data in Au+Au collisions, whereas Figs. 2(c) and 2(d) show U+U collisions. In each panel, the minimum multiplicity corresponds to roughly 20% centrality, whereas the maximum multiplicity is around 0.1% centrality.

### 1. Spherical <sup>197</sup>Au nuclei

I start by showing results from the $T_{\text{RENTO}}$ calculations implementing spherical <sup>197</sup>Au nuclei. These results are reported as full symbols in Fig. 2.

$p=0$: The calculations using $T_{\text{RENTO}}$ with $p=0$ are shown as red symbols in Figs. 2(a) and 2(c). The value of the constant $\kappa = 0.165$ is chosen in order to provide the most accurate possible description of the cumulant $v_2\{4\}$ in Au+Au collisions [Fig. 2(a)]. I choose to tune the model to $v_2\{4\}$ because 4-particle cumulants are not affected by non-flow contributions [32]. The agreement between the model and the data in Au+Au collisions in Fig. 2(a) is excellent. The description of $v_2\{4\}$ provided by this $T_{\text{RENTO}}$ parametrization is perfect. It captures the trend of the data all the way up to 20% centrality, and it correctly reproduces the change of sign of the cumulant $v_2\{4\}^4$ observed around 2.5% centrality in experiment. Agreement is equally impressive for $v_2\{2\}$, all the way up

to 20%. The slight shift of the experimental data toward larger values with respect to the model is easily understandable as due to nonflow. Indeed, the nonflow subtraction performed by the STAR Collaboration is not perfect, as two-particle correlations were calculated implementing a small gap of 0.1 units of pseudorapidity in the analysis [11]. Moving to the red symbols shown in Fig. 2(c), the description of U+U data provided by this model is overall good, but less satisfactory. I see that, although the comparison with data is not as quantitatively as good as in Au+Au collisions, the model captures nicely the qualitative features due to the deformation of <sup>238</sup>U nuclei: $v_2\{2\}$ is larger in U+U than in Au+Au and, for reasons that I shall discuss in detail in Sec. IV, $v_2\{4\}$ is observed to be nonzero all the way to the highest multiplicities.

$p=1$: I look now at the red symbols shown in Figs. 2(b) and 2(d), where I show results for $T_{\text{RENTO}}$ with $p=1$. Again, the value of $\kappa = 0.27$ is chosen such to yield the best description of $v_2\{4\}$ in Au+Au collisions.<sup>9</sup> Starting with Au+Au collisions in Fig. 2(b), I see that the description of experimental data provided by this $T_{\text{RENTO}}$ parametrization is much worse than for the $p=0$ case: $v_2\{2\}$ is overestimated and $v_2\{4\}$ reaches zero at too large centrality with respect to the experimental data. This is consistent with the comparisons between this model and LHC data shown in Ref. [39]. Agreement with data is also bad if I look at U+U collisions in Fig. 2(d), where a much smaller value of $\kappa$ would be needed to match the experimental data. The qualitative features of U+U collisions are, on the other hand, reasonably captured, in particular, the fact that the cumulant $v_2\{4\}$ is nonzero all the way to the highest multiplicity.

### 2. Deformed <sup>197</sup>Au nuclei

I discuss now the results obtained in $T_{\text{RENTO}}$ with deformed <sup>197</sup>Au nuclei presenting $\beta_2 = -0.125$ and $\beta_4 = -0.017$ [14]. These results are reported as lines in Fig. 2.

$p=0$: The results for $T_{\text{RENTO}}$ with $p=0$ are the green lines in Figs. 2(a) and 2(c). The constant $\kappa$ that I need to tune the model to $v_2\{4\}$ in Au+Au collisions turns out to be smaller than previously with spherical <sup>197</sup>Au nuclei, since both cumulants increase if I implement $\beta_2 > 0$. Now, in Fig. 2(a) I see that the description of $v_2\{2\}$ is still very good all the way up to 20% centrality, and so is the description of $v_2\{4\}$ essentially above 5% centrality. But switching on the nuclear deformation has a dramatic effect of the fourth-order cumulant: It prevents $v_2\{4\}^4$ from going negative, and one

<footnote n="9"><sup>9</sup>Note that the coefficient $\kappa$ varies a lot moving from $p=0$ to $p=1$. This is due to the fact that the eccentricity grows much faster with centrality for $p=0$ than for $p=1$.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 2. (a) and (b)] and central U+U collisions at $\sqrt{s} = 193$ GeV [panel (c) and (d)], as function of collision multiplicity. Full symbols: Cumulants of initial eccentricity fluctuations in the $T_{RENTO}$ setup that implements spherical $^{197}$Au nuclei, with $p = 0$ [panels (a) and (c)] and $p = 1$ [panels (b) and (d)]. Lines: Cumulants of initial eccentricity fluctuations in the $T_{RENTO}$ setup that implements deformed $^{197}$Au nuclei, with $p = 0$ [panels (a) and (c)] and $p = 1$ [panels (b) and (d)]. The arrows indicate the correspondence between multiplicity and collision centrality.</figcaption></figure>

observes a nonzero $v_2\{4\}$ all the way to the highest multiplicity, much as in the case of U+U collisions. Nothing notable occurs for U+U collisions in Fig. 2(c).

p = 1: As for the results with $p = 1$ with deformed $^{197}$Au nuclei, i.e., the lines shown in Figs. 2(b) and 2(d), I would simply like to point out is that I find the same striking result observed in Fig. 2(a): $v_2\{4\}$ is always nonzero in Au+Au collisions.

### C. Discussion

Let me draw my conclusions from the comparisons shown in Fig. 2.

First, I confirm that the $p = 1$ model is ruled out by elliptic flow fluctuations data at RHIC. This is not a surprise: First, this model is currently ruled out by LHC data [39]; Second, as anticipated, this model is similar to the models used by the STAR Collaboration in their analysis, and those models provide a bad description of data, in the sense that they present a $\kappa$ that varies with centrality even at the highest multiplicities.

Second, my results imply that a value $\beta_2 \sim -0.12$ in $^{197}$Au nuclei is essentially ruled out by experimental data. The description of data achieved with the $T_{RENTO}$ $p = 0$ parametrization implementing spherical $^{197}$Au nuclei is remarkably good and is spoiled by the inclusion of the quadrupole parameter reported in the literature. Note that this is not a small effect, but on the contrary, it is a very visible change in the fourth-order cumulant that does not reach zero
even in the most central collisions, as clearly observed in STAR data. This occurs for both $p=0$ and $p=1$. I would like to stress, once more, that the analysis reported here is exhaustive, in the sense that one cannot easily argue that negative values of $v_2\{4\}^4$ can be obtained with $\beta_2 \sim -0.12$ in a different model setup. As an example, a slightly different class of models are those models that can be fitted by a negative $p$ in $T_{RENTO}$, for instance, the Monte Carlo KLN model [47] or the Monte Carlo rcBK [48] model. But there the problem would remain, as these models present both a larger eccentricity [28] and a larger $\varepsilon_2\{4\}$ [39] than the $p=0$ case.

The bottom line of this section is that the most notable effects of nuclear deformation are visible in $v_2\{4\}$ rather than in $v_2\{2\}$. The typical statement is that nuclear deformation yields larger fluctuations of $v_2$ and thus a larger $v_2\{2\}$. My analysis clearly indicates that experimental data are now precise enough to show distinct signatures of nuclear deformation in the details of the flow fluctuations: $v_2\{4\}^4$ in positive U+U collisions and negative in Au+Au collisions. Note that LHC data seems to point at the same phenomenon. Indeed, preliminary ATLAS data show that $v_2\{4\}$ in Xe+Xe collisions is positive, all the way down to the most central events [49]. This is in contrast with Pb+Pb data, where a change of sign of the fourth-order cumulant of elliptic flow is currently observed [50], and $^{208}$Pb nuclei are perfectly spherical.<fnref n="10" />

How is it, then, that $v_2\{4\}$ is so sensitive to the deformed nuclear shapes?

## IV. NON-GAUSSIAN FLUCTUATIONS

Flow fluctuations are to a good approximation Gaussian, in the sense that the elliptic flow vector $(v_x, v_y)$ has a two-dimensional Gaussian distribution in a given centrality bin. The breakdown of Gaussian behavior leaves observable consequences in the higher order cumulants of elliptic flow. I refer to Refs. [52–54] for exhaustive calculations explaining how the non-Gaussian nature of fluctuations manifests in higher order cumulants and to the detailed analyses of Refs. [29–32] for experimental confirmations of such results.

In this section, I show that nuclear deformation can yield pronounced non-Gaussianity in the event-by-event distribution of elliptic flow, and I discuss the observable consequences of such an effect.

### A. Deformation as a source of non-Gaussianity

Let me start by looking at the sensitivity of the cumulants $\varepsilon_2\{2\}^2$ and $\varepsilon_2\{4\}^4$ under variations of the parameter $\beta_2$, using the $T_{RENTO}$ model with $p=0$ and $k=0.5$. I simulate U+U collisions at zero impact parameter, $b=0$, and I systematically vary the quadrupole deformation of the nuclei. Note that colliding at $b=0$ implies that $\varepsilon_2\{2\}^2$ is equal to the variance of the distribution of eccentricity, whereas $\varepsilon_2\{4\}^4$ measures the

<figure><img src="image_6.png" /><figcaption>FIG. 3. Cumulants of eccentricity fluctuations in U+U collisions at zero impact parameter. The cumulants are plotted as function of the quadrupole coefficient, $\beta_2$. I simulate collisions using the $T_{RENTO}$ model with $p=0$ and $k=0.5$. (a) Variance of eccentricity fluctuations. (b) Kurtosis of eccentricity fluctuations. Lines are polynomial fits.</figcaption></figure>

kurtosis, i.e., whether the distribution has heavier or lighter tails than a Gaussian. In particular, at $b=0$ one expects $\varepsilon_2\{4\}^4 = 0$ if the distribution of the eccentricity vector is a two-dimensional Gaussian.

In Fig. 3(a), I display $\varepsilon_2\{2\}^2$ as function of $\beta_2$ for collisions at zero impact parameter. By the symmetry properties of the spherical harmonic multiplying $\beta_2$ in Eq. (1), any effect that is linear in $\beta_2$ should cancel when averages over events are taken. Therefore, the value of $\langle \varepsilon_2^2 \rangle$ is expected to grow with the square of $\beta_2$. To check this, I perform a parabolic fit of $\varepsilon_2\{2\}^2$, reported as a dashed line in Fig. 3(a). The fit is of excellent quality, and it shows how the variance of the distribution varies with the quadrupole coefficient.

Figure 3(b) shows instead $\varepsilon_2\{4\}^4$ as function of $\beta_2$. I note a great enhancement of this quantity with increasing $\beta_2$. This implies that nuclear deformation does not simply make the distribution of eccentricity broader, but also less Gaussian, as it makes the kurtosis grow by essentially orders of magnitude. The dashed line is a quartic fit, which again confirms the symmetry argument.

The previous result is very intuitive. The kurtosis enters in $\varepsilon_2\{4\}^4$ with a negative sign [54], which means that nuclear deformation makes the kurtosis of $\varepsilon_2$ fluctuations less negative. Negative kurtosis for an azimuthally symmetric eccentricity

<footnote n="10">¹⁰I should stress that the precise mechanism that leads to the change of sign of this cumulant in Pb+Pb collisions is not understood yet. It seems to be a generic feature of systems presenting small eccentricity driven by impact parameter fluctuations [51].</footnote>
distribution can be simply understood as a consequence of the fact that the eccentricity is bounded by unity [55]. Nuclear deformation, then, causes eccentricity fluctuations to get closer to their bound, and thus to present more negative kurtosis.

I conclude that nuclear deformation yields a slight increase in the variance of elliptic flow fluctuations and a dramatic increase in their kurtosis. This explains intuitively why $v_2\{4\}^4$ is much larger in U+U collisions than in Au+Au collision, provided, as experimental data seem to suggest, that $^{197}$Au nuclei are much more spherical than $^{238}$U nuclei.

In the following, I argue that this phenomenon has robust, observable consequences that go beyond the simple positive value of $v_2\{4\}$ in central collisions.

## B. Higher order cumulants

The bottom line of the previous discussion is that nuclear deformation yields stronger non-Gaussianity in the distribution of $v_2$. Here I discuss the implication of this phenomenon for the splitting between higher order cumulants of elliptic flow.

The third cumulant of elliptic flow is of order 6 and is equal to

$$v_2\{6\} = 6\sqrt{\frac{1}{4}(\langle v_2^6 \rangle - 9\langle v_2^2 \rangle\langle v_2^4 \rangle + 12\langle v_2^2 \rangle^3)} \quad (11)$$

Again, in the regime of linear hydrodynamic response, one can write

$$v_2\{6\} = \kappa\epsilon_2\{6\}, \quad (12)$$

that, using Eq. (10), leads to [39,56]

$$\frac{v_2\{6\}}{v_2\{4\}} = \frac{\epsilon_2\{6\}}{\epsilon_2\{4\}}. \quad (13)$$

For collisions of spherical nuclei, e.g., $^{208}$Pb nuclei, the ratio $v_2\{6\}/v_2\{4\}$ is very close to unity (~0.99) in noncentral collisions [29,30]. This is because the probability distribution of elliptic flow, and consequently of the initial eccentricity, is well approximated by a two-dimensional Gaussian. Gaussian fluctuations imply that all higher order cumulants of $\epsilon_2$ are equal to the mean value of the eccentricity projected along the reaction plane [57], which I dub, in a standard notation, $\epsilon_{RP}^{11}$

$$\epsilon_2\{4\} \approx \epsilon_2\{6\} \approx \epsilon_{RP}. \quad (15)$$

This is the famous degeneracy of cumulants of the eccentricity and of elliptic flow. Departure from Gaussian behavior is expected starting from semicentral collisions, due to the negative skewness of $\epsilon_2$ fluctuations in the reaction plane [52], as recently confirmed by LHC data in Pb+Pb collisions [30,31].

$$\epsilon_{RP} = \left\langle \frac{\int r^2 \cos 2\phi S(r, \phi) r dr d\phi}{\int r^2 S(r, \phi) r dr d\phi} \right\rangle, \quad (14)$$

where the average is over events.

I test the validity of Eq. (15) in the $T_R$ENTo model with $p=0$. I show in Fig. 4(a) the values of $\epsilon_2\{4\}$ and $\epsilon_{RP}$ in central Au+Au and U+U collisions. I choose Au+Au collisions implementing spherical nuclei, in order to highlight how results change between collisions of spherical nuclei and collisions of nonspherical nuclei. For Au+Au collisions, I see, as expected, that $\epsilon_2\{4\}$ is essentially equal to $\epsilon_{RP}$ already at 5% centrality. This is the onset of Gaussian fluctuations. Moving on to U+U collisions, I observe that the splitting between $\epsilon_2\{4\}$ and $\epsilon_{RP}$ is much larger than in Au+Au collisions. This nicely illustrates how this cumulant, in the most central events, becomes fully dominated by the enhancement of the kurtosis of the distribution, which is due to the fluctuations of the orientation of the colliding nuclei. In these collisions, then, I do not observe any onset of Gaussian fluctuations, at least up to 10% centrality. I conclude that nuclear deformation breaks the degeneracy of cumulants, Eq. (15), in central events.

Therefore, the robust, model-independent prediction I can easily draw is the following: In collisions of deformed nuclei, I expect a large splitting between $v_2\{4\}$ and $v_2\{6\}$ between ~5% centrality and semicentral collisions. In Fig. 4(b), I show predictions for the ratio $v_2\{6\}/v_2\{4\}$, up to 20% centrality. Note that for U+U collisions, Gaussian fluctuations, i.e., $v_2\{6\}/v_2\{4\} \sim 0.99$, are not observed below 15% centrality. This confirms my expectations: The ratio is significantly lower in U+U collisions than in Au+Au collisions and the effect is very visible. Experimental verification of this feature would provide additional confirmation of the great robustness of the hydrodynamic modeling.

Note that, from the inset in Fig. 4(b), I also predict that $v_2\{6\}^6$, much as $v_2\{4\}^4$, is positive in U+U collisions all the way to the ultracentral events.

## V. CONCLUSIONS

I have shown that the largest effects of nuclear deformation are hidden into the tails of the distribution of elliptic flow in central collisions. Nuclear deformation yields broader distribution of elliptic flow, and, for central collisions, it engenders the mechanism of enhanced negative kurtosis presented in Sec. IV, which allows me to predict and explain the behavior of the higher order cumulants of $v_2$. I have thus explained in very simple terms the striking observation that $v_2\{4\}^4$ is much larger in U+U collisions than in Au+Au collisions and predicted that the splitting between $v_2\{6\}$ and $v_2\{4\}$ is larger in collisions of deformed nuclei.

My results imply that there is no room in the experimental data for a $\beta_2$ of order 0.12 in the $^{197}$Au nuclei, as it would lead to too large values of the cumulant $v_2\{4\}^4$ in a model-independent way. Clearly, one would rather trust the numbers provided in established nuclear physics literature, instead of the crude models of initial conditions for heavy-ion collisions. But, arguably, this would just amount to rejecting evidence. Signatures of nuclear deformation are observed in elliptic flow data in U+U collisions and, more remarkably, in collisions of $^{129}$Xe nuclei that present a $\beta_2$ of order 16%, which is very close to the reported 12% of $^{197}$Au nuclei. As I have shown in this paper, in the hydrodynamic framework one can understand all these observations as simply due to generic and

<footnote>$^{11}$Following Eq. (8), the mean eccentricity along the direction of the reaction plane, which is customarily taken as the x axis, is given by</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 4. (a) $\varepsilon_2\{4\}$ in central Au+Au (full symbols) and U+U (empty symbols) collisions in the $T_R$ENTO model with $p = 0$. Lines indicate the mean eccentricity in the reaction plane, $\varepsilon_{RP}$, in both Au+Au (dashed line) and U+U (solid line) collisions. (b) Model predictions for the ratio $v_2\{6\}/v_2\{4\}$ as function of centrality percentile, in central Au+Au (empty symbols) and U+U (full symbols) collisions. The inset is an enlargement below 5% centrality.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

model-independent features of the fluctuations of the initial geometry. Therefore, any such effects should show up as well in Au+Au data.

This is an interesting puzzle, that, I think, nicely underlines the close link between flow fluctuations in central heavy-ion collisions, whose origin is purely geometric, and the actual shape and structure of the colliding nuclei. This may lead to interesting developments in the future, aimed at matching these two very different areas of nuclear physics. It would be useful, for instance, to collide nuclei

whose structure and deformation is precisely determined experimentally.

## ACKNOWLEDGMENTS

I thank Jean-Yves Ollitrault for useful input and illuminating discussions. I thank Björn Schenke and Prithwish Tribedy for providing me with the results of IP-Glasma simulations and for useful comments on the manuscript. I also thank Jamie Nagle for useful comments on the first draft of this paper.

[1] J. Y. Ollitrault, Phys. Rev. D 46, 229 (1992).

[2] B. Alver, B. B. Back, M. D. Baker, M. Ballintijn, D. S. Barton, R. R. Betts, R. Bindel, W. Busza, V. Chetluru, E. Garcia et al., Phys. Rev. C 77, 014906 (2008).

[3] J. P. Blaizot, W. Broniowski, and J. Y. Ollitrault, Phys. Lett. B 738, 166 (2014).

[4] B. Alver et al. (PHOBOS Collaboration), Phys. Rev. Lett. 98, 242302 (2007).

[5] D. Teaney and L. Yan, Phys. Rev. C 83, 064904 (2011).

[6] F. G. Gardim, F. Grassi, M. Luzum, and J. Y. Ollitrault, Phys. Rev. C 85, 024908 (2012).

[7] H. Niemi, G. S. Denicol, H. Holopainen, and P. Huovinen, Phys. Rev. C 87, 054901 (2013).

[8] F. G. Gardim, J. Noronha-Hostler, M. Luzum, and F. Grassi, Phys. Rev. C 91, 034902 (2015).

[9] H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016).

[10] J. Noronha-Hostler, L. Yan, F. G. Gardim, and J. Y. Ollitrault, Phys. Rev. C 93, 014909 (2016).

[11] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 115, 222301 (2015).

[12] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B 784, 82 (2018).

[13] P. Moller, J. R. Nix, W. D. Myers, and W. J. Swiatecki, Atom. Data Nucl. Data Tabl. 59, 185 (1995).

[14] P. Möller, A. J. Sierk, T. Ichikawa, and H. Sagawa, Atom. Data Nucl. Data Tabl. 109–110, 1 (2016).

[15] A. Goldschmidt, Z. Qiu, C. Shen, and U. Heinz, Phys. Rev. C 92, 044903 (2015).

[16] G. Giacalone, J. Noronha-Hostler, M. Luzum, and J. Y. Ollitrault, Phys. Rev. C 97, 034904 (2018).

[17] K. J. Eskola, H. Niemi, R. Paatelainen, and K. Tuominen, Phys. Rev. C 97, 034911 (2018).

[18] P. F. Kolb, J. Sollfrank, and U. W. Heinz, Phys. Rev. C 62, 054909 (2000).

[19] U. W. Heinz and A. Kuhlman, Phys. Rev. Lett. 94, 132301 (2005).

[20] A. J. Kuhlman and U. W. Heinz, Phys. Rev. C 72, 037901 (2005).

[21] A. Kuhlman, U. W. Heinz, and Y. V. Kovchegov, Phys. Lett. B 638, 171 (2006).

[22] C. Nepali, G. Fai, and D. Keane, Phys. Rev. C 73, 034911 (2006).
[23] H. Masui, B. Mohanty, and N. Xu, Phys. Lett. B 679, 440 (2009).

[24] T. Hirano, P. Huovinen, and Y. Nara, Phys. Rev. C 83, 021902 (2011).

[25] M. Rybczyński, W. Broniowski, and G. Stefanek, Phys. Rev. C 87, 044908 (2013).

[26] B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. C 89, 064908 (2014).

[27] Q. Y. Shou, Y. G. Ma, P. Sorensen, A. H. Tang, F. Videbæk and H. Wang, Phys. Lett. B 749, 215 (2015).

[28] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 92, 011901 (2015).

[29] G. Aad et al. (ATLAS Collaboration), Eur. Phys. J. C 74, 3157 (2014).

[30] A. M. Sirunyan et al. (CMS Collaboration), Phys. Lett. B 789, 643 (2019).

[31] S. Acharya et al. (ALICE Collaboration), J. High Energy Phys. 1807, 103 (2018).

[32] A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 99, 024903 (2019).

[33] C. Loizides, J. Nagle, and P. Steinberg, SoftwareX 1-2, 13 (2015).

[34] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016).

[35] J. S. Moreland, J. E. Bernhard, and S. A. Bass, arXiv:1808.02106 [nucl-th].

[36] B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. Lett. 108, 252301 (2012).

[37] J. L. Nagle and W. A. Zajc, arXiv:1808.01276 [nucl-th].

[38] L. Yan, J. Y. Ollitrault, and A. M. Poskanzer, Phys. Lett. B 742, 290 (2015).

[39] G. Giacalone, J. Noronha-Hostler, and J. Y. Ollitrault, Phys. Rev. C 95, 054910 (2017).

[40] M. Rybczyński, G. Stefanek, W. Broniowski, and P. Bozek, Comput. Phys. Commun. 185, 1759 (2014).

[41] S. J. Das, G. Giacalone, P. A. Monard, and J. Y. Ollitrault, Phys. Rev. C 97, 014905 (2018).

[42] L. McLerran and P. Tribedy, Nucl. Phys. A 945, 216 (2016).

[43] R. S. Bhalerao, M. Luzum, and J. Y. Ollitrault, Phys. Rev. C 84, 034910 (2011).

[44] P. Alba, V. M. Sarti, J. Noronha, J. Noronha-Hostler, P. Parotto, I. P. Vazquez and C. Ratti, Phys. Rev. C 98, 034909 (2018).

[45] C. Chattopadhyay, R. S. Bhalerao, J. Y. Ollitrault, and S. Pal, Phys. Rev. C 97, 034915 (2018).

[46] C. Nair, M. Erhard, A. R. Junghans, D. Bemmerer, R. Beyer, E. Grosse, J. Klug, K. Kosev, G. Rusev, K. D. Schilling, R. Schwengner, and A. Wagner, Phys. Rev. C 78, 055802 (2008).

[47] H.-J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007).

[48] J. L. Albacete and A. Dumitru, arXiv:1011.5161 [hep-ph].

[49] ATLAS Collaboration, ATLAS-CONF-2018-011, https://cds.cern.ch/record/2318870.

[50] ATLAS Collaboration, ATLAS-CONF-2017-066, https://cds.cern.ch/record/2285570.

[51] M. Zhou and J. Jia, Phys. Rev. C 98, 044903 (2018).

[52] G. Giacalone, L. Yan, J. Noronha-Hostler, and J. Y. Ollitrault, Phys. Rev. C 95, 014913 (2017).

[53] N. Abbasi, D. Allahbakhshi, A. Davody, and S. F. Taghavi, Phys. Rev. C 98, 024906 (2018).

[54] R. S. Bhalerao, G. Giacalone, and J. Y. Ollitrault, Phys. Rev. C 99, 014907 (2019).

[55] L. Yan and J. Y. Ollitrault, Phys. Rev. Lett. 112, 082301 (2014).

[56] L. Ma, G. L. Ma, and Y. G. Ma, Phys. Rev. C 94, 044915 (2016).

[57] S. A. Voloshin, A. M. Poskanzer, A. Tang, and G. Wang, Phys. Lett. B 659, 537 (2008).