# Constraining the quadrupole deformation of atomic nuclei with relativistic nuclear collisions

Giuliano Giacalone

► To cite this version:

Giuliano Giacalone. Constraining the quadrupole deformation of atomic nuclei with relativistic nuclear collisions. Phys.Rev.C, 2020, 102 (2), pp.024901. 10.1103/PhysRevC.102.024901. hal-02613260

HAL Id: hal-02613260

https://hal.science/hal-02613260

Submitted on 23 Jun 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the quadrupole deformation of atomic nuclei with relativistic nuclear collisions

Giuliano Giacalone*

Université Paris Saclay, CNRS, CEA, Institut de Physique Théorique, 91191 Gif-sur-Yvette, France

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

(Received 1 May 2020; accepted 21 July 2020; published 3 August 2020)

Preliminary data by the STAR collaboration at the BNL Relativistic Heavy Ion Collider shows that the elliptic flow, $v_2$, and the average transverse momentum, $\langle p_t \rangle$, of final-state hadrons produced in high-multiplicity $^{238}$U + $^{238}$U collisions are negatively correlated. This observation brings experimental evidence of a significant prolate deformation, $\beta \approx 0.3$, in the colliding $^{238}$U nuclei. I show that a quantitative description of this new phenomenon can be achieved within the hydrodynamic framework of heavy-ion collisions, and that thus such kind of data in the context of high-energy nuclear experiments can help constrain the quadrupole deformation of the colliding species.

DOI: 10.1103/PhysRevC.102.024901

## I. INTRODUCTION

Deformation is a fundamental property of atomic nuclei, reflecting the greatly collective and correlated nature of the dynamics of nucleons within the quantum many-body system. The majority of atomic nuclei possess in fact an intrinsic deformation, most notably, an axial quadrupole, or ellipsoidal, deformation. Such deformation exists when a nucleus, described by a charge (or mass) density $\rho(\mathbf{r})$, has a nonvanishing electric (or mass) quadrupole moment:

$$ \langle Y_2^0(\Theta, \Phi)r^2\rho(\mathbf{r}) \rangle \neq 0, \quad (1) $$

where $\mathbf{r} = (\Theta, \Phi, r)$, $Y_2^0 \propto 3 \cos^2 \Theta - 1$, and angular brackets denote an expectation value with respect to the nuclear wave function. The magnitude of the quadrupole deformation, i.e., the eccentricity of the nuclear body, can be roughly obtained by dividing Eq. (1) by the mean-squared nuclear radius, and is typically quantified by a dimensionless deformation parameter, $\beta$ [1].

In the context of low-energy nuclear experiments, the value of $\beta$ is usually inferred from measurements of the transition probability of the electric quadrupole operator from the $0^+$ ground state to the first excited $2^+$ state [2–4]. Simple formulas relating this measurable quantity to the value of $\beta$ are derived under model assumptions, and a comprehensive collection of such experimentally-inferred deformation parameters is provided in Ref. [3]. Alternatively, values of $\beta$ for essentially the whole spectrum of known nuclei have been tabulated in extensive predictions, or extrapolations of theoretical models [5,6]. These results do not necessarily agree with the figures reported in Ref. [3].

In this article, I argue that high-energy nuclear experiments, or relativistic heavy-ion collisions, can help place new constraints on the values of $\beta$. In these experiments, substantial evidence of a quadrupole deformation in the colliding

species has emerged only recently. At the BNL Relativistic Heavy Ion Collider (RHIC), the STAR collaboration observed unambiguous signatures of nuclear deformation by means of accurate comparisons between data in $^{197}$Au + $^{197}$Au collisions and data in $^{238}$U + $^{238}$U collisions [7], the latter nuclei being much more deformed in their ground state. The same kind of observations were later made as well at the CERN Large Hadron Collider (LHC), by comparing data in collisions of spherical $^{208}$Pb nuclei with collisions of deformed $^{129}$Xe nuclei [8–10].

It is natural that the observables measured in heavy-ion collisions enable us to probe the geometry of the colliding nuclei. Two nuclei colliding at relativistic energy produce a quark-gluon plasma [11], the high-temperature fluidlike state of strong-interaction matter. This system is produced at rest, and is set in motion by pressure gradients which are determined by its geometry at the onset of the hydrodynamic behavior. A head-on collision between deformed nuclei can yield significant spatial asymmetry in the shape of the created medium [12–18], which in turn triggers an anisotropic flow of matter towards the detectors. This anisotropy is then observed in experiments as an enhanced quadrupole component in the angular distribution of final-state hadrons, dubbed elliptic flow, $v_2$ [19].

However, it was recently realized [20] that new observables exhibiting a great sensitivity to the deformation of the colliding nuclei can be constructed by looking at the magnitude of the isotropic flow of particles, which is given by the average momentum of hadrons, $\langle p_t \rangle$, in the plane orthogonal to the collision axis. This quantity does not carry information about the spatial anisotropy of the quark-gluon plasma, but is sensitive to its size and temperature (or energy). The realization of Ref. [20] is that systems emitting the same number of particles can be classified according to their temperature by looking at their $\langle p_t \rangle$. By doing so, one can get an experimental handle on the orientation of the nuclear ellipsoids at the time of interaction. This idea is reviewed in detail in Sec. III. The outcome is that in nearly head-on collisions at fixed final-state

<footnote>*giuliano.giacalone@universite-paris-saclay.fr</footnote>
multiplicity, there should exist a negative correlation between elliptic flow, $v_2$, and the average transverse momentum of hadrons, $\langle p_t \rangle$, an effect which is entirely engendered by the nonzero value of $\beta$ in the colliding species. Remarkably enough, this prediction has been already investigated by the STAR collaboration at RHIC, and preliminary experimental data confirms [21] the existence of a negative correlation between $v_2$ and $\langle p_t \rangle$ in central $^{238}$U + $^{238}$U collisions.

Here, I show that this new measurement represents a crucial leap towards a quantitative phenomenology of nuclear structure at high-energy nuclear colliders. I take the idea of Ref. [20], and I improve its phenomenological ingredients to show that, within the hydrodynamic framework of heavy-ion collisions, it is indeed possible to achieve a quantitative description of the new phenomenon observed by the STAR collaboration. The most important feature introduced in this work is the implementation of a more realistic initial-state predictor for $\langle p_t \rangle$, as I shall explain through Sec. II. Secondly, I implement in a realistic way the statistical fluctuations that plague the determination of $\langle p_t \rangle$ in a single nucleus-nucleus collision event, due to the finite hadron multiplicity, as I explain in Sec. III D. This improved setup allows me to obtain results that can be considered as quantitative predictions, as shown in Sec. IV, and that can be compared with future data. Further assumptions made in Ref. [20], e.g., about the relation between initial and final anisotropies in hydrodynamics, are investigated by means of state-of-the-art hydrodynamic simulations, presented in Appendix B.

## II. AVERAGE TRANSVERSE MOMENTUM AND ITS FLUCTUATIONS

Two nuclei colliding at relativistic energy produce a quark-gluon plasma, the hot state of strong-interaction matter. Despite being as large as an atomic nucleus, the dynamics of this entity is driven by macroscopic physics, i.e., pressure gradient forces and velocity fields. Due to an immense Lorentz contraction of the colliding nuclei in the laboratory frame, the created medium can be considered as invariant under longitudinal boosts [22], so that for most practical purposes it is enough to look at the dynamics of the system in a rapidity slice around the interaction point (midrapidity).

When the plasma decouples to particles, a spectrum of hadrons is produced at midrapidity, $\frac{dN}{d^2\mathbf{p}_t}$, where $\mathbf{p}_t$ is the particle momentum in the plane transverse to the collision axis, or transverse plane. The quantity I shall focus on is the average transverse momentum of the hadrons emitted in a collision event:<fnref n="1" />

$$ \langle p_t \rangle = \frac{1}{N} \int_{\mathbf{p}_t} \frac{dN}{d^2\mathbf{p}_t} p_t, \qquad (2) $$

where $p_t \equiv |\mathbf{p}_t|$, and $N = \int_{\mathbf{p}_t} dN/d^2\mathbf{p}_t$ is the total number of particles (or multiplicity) detected in the event. The average transverse momentum can be evaluated on an event-by-event

<figure><img src="image_3.png" /><figcaption>FIG. 1. Relation between system size, $R$, and average transverse momentum, $\langle p_t \rangle$, in a system with given entropy, $S$. Systems A and B have the same total entropy, but B has a larger volume. The temperature in A is therefore larger. Assuming an ideal classical gas of massless particles, this implies that the particles in A carry more momentum, $p_t$.</figcaption></figure>

basis, and its probability distribution can be measured in experiments.

This preliminary section is devoted to elucidating the physical origin of $\langle p_t \rangle$, and of its fluctuations in nucleus-nucleus collisions. This knowledge will be crucial in the subsequent discussion of the phenomenology of nuclear structure in heavy-ion collisions.

### A. Relation between system size and $\langle p_t \rangle$

Consider a uniform ideal gas of massless particles with Boltzmann statistics. The energy of a particle, $E$, in such a system is equal to $3T$, where $T$ is the temperature, while the total entropy, $S$, is proportional to the number of particles [23]. If the system is relativistic, the energy of a particle coincides with its momentum, $p$, and one has

$$ p \approx E = 3T. \qquad (3) $$

The quark-gluon plasma created in ultrarelativistic heavy-ion collisions is essentially a relativistic ideal gas of massless particles, and the existence of similar relations in the context of heavy-ion collisions, where one replaces $p$ with the final $\langle p_t \rangle$, has been recently established by means of modern hydrodynamic simulations [24]:

$$ \langle p_t \rangle \approx 3T, \qquad (4) $$

where $T$ is an effective temperature of the system.

Suppose one is given two uniform quark-gluon fluids that have the same entropy, but contained within different volumes. The system with a smaller volume is denser, it has larger temperature, and thus, according to Eq. (4), it yields larger $\langle p_t \rangle$ in the final state. This discussion is summarized in Fig. 1, for two systems with the same entropy but different volumes. For uniform systems with fixed number of particles, then, there is a one-to-one correspondence between the system size, which I dub $R$, and $\langle p_t \rangle$, and in fact all other thermodynamic quantities (such as the energy). It is natural, then, that the origin of the fluctuations of $\langle p_t \rangle$ in the literature was identified with the fluctuations of $R$ [25–28].

In the limit of small fluctuations, one can write an explicit relation between the fluctuations of these quantities.

<footnote n="1">¹In an actual experiment, the average transverse momentum comes from a discrete sum of a finite number of terms (particles), as discussed in Sec. III D.</footnote>
First write the following thermodynamic identity [29]:

$$c_s^2 = \frac{dP}{d\epsilon} = \frac{d \ln T}{d \ln s}, \qquad (5)$$

where P, ε, T, and s are, respectively, the pressure, energy density, temperature, and entropy density of the fluid. Then, dimensional analysis imply that $s \propto R^{-3}$, while $\langle p_t \rangle \propto T$ from Eq. (4), so that

$$\frac{d\langle p_t \rangle}{\langle p_t \rangle} = -3c_s^2 \frac{dR}{\langle R \rangle}, \qquad (6)$$

where the minus sign comes from the fact that $\langle p_t \rangle$ is an-ticorrelated with R. In natural units ($c = 1$), the factor $3c_s^2$ is typically of order unity. This is the relation employed in Ref. [20] to transform the event-by-event distribution of R into an event-by-event distribution of $\langle p_t \rangle$. However, this picture is modified when the system is not uniform. The first important result of the present paper, described in the next subsection, is an improvement of Eq. (6) via the definition of a better initial-state predictor for $\langle p_t \rangle$.

## B. Initial energy as a predictor for $\langle p_t \rangle$

The density profile of the quark-gluon plasma formed in relativistic nuclear collisions is neither uniform nor smooth. It is a rough landscape, that fluctuates significantly on an event-by-event basis. In presence of density fluctuations, Eq. (3) remains valid, and $\langle p_t \rangle$ remains a measure of the temperature, or of the energy, of the system; But the system size, R, is no longer in a one-to-one correspondence with the other thermodynamic quantities, and can not be considered as a natural predictor of $\langle p_t \rangle$.

I give an explicit illustration of this phenomenon. I define the size of the quark-gluon plasma as follows:

$$R^2 = \frac{\int_{\mathbf{x}} |\mathbf{x}|^2 s(\mathbf{x}, \tau_0)}{\int_{\mathbf{x}} s(\mathbf{x}, \tau_0)}, \qquad (7)$$

where $\mathbf{x}$ is a coordinate in the transverse plane, and $s(\mathbf{x}, \tau_0)$ is the transverse profile of entropy density at the time $\tau_0$ when hydrodynamics becomes applicable. I also define the initial energy of the medium by

$$E_0 \equiv \tau_0 \int_{\mathbf{x}} e(\mathbf{x}, \tau_0), \qquad (8)$$

where $e(\mathbf{x}, \tau_0)$ is the initial energy-density profile, obtained from $s(\mathbf{x}, \tau_0)$ via an equation of state. In Refs. [30,31], it was indeed shown that, at fixed total entropy (per unit rapidity), $S = \tau_0 \int_{\mathbf{x}} s(\mathbf{x}, \tau_0)$, the value of $\langle p_t \rangle$ at the end of hydrodynamics is in an almost perfect one-to-one correspondence with the value of $E_0$. This result was obtained in ideal hydrodynam-ics simulations. Here, in Appendix B, in Fig. 13, I check explicitly that this is true as well in viscous hydrodynamics. The correlation between $\langle p_t \rangle$ and R is instead much weaker. The reason is precisely that there is a significant dispersion between R and $E_0$.

I show this dispersion in simulations of Au + Au collisions at top RHIC energy, corresponding to nucleon-nucleon interactions at a center-of-mass energy $\sqrt{s} = 200$ GeV. I do this by means of the TRENTo model of initial condi-tions, a parametrization of the initial entropy density, $s(\mathbf{x}, \tau_0)$,

<figure><img src="image_4.png" /><figcaption>FIG. 2. Joint distribution of initial energy, $E_0$, and system size, R, in ultracentral Au-Au collisions at top RHIC energy simulated with the TRENTo model (described in Appendix A). Darker pixels correspond to larger numbers of counts.</figcaption></figure>

whose working principles and parameters are reported in Appendix A. The entropy density returned by the TRENTo calculation is converted into an energy density by means, for simplicity, of the equation of state of conformal high-temperature QCD ($c_s^2(T) = 1/3$):

$$e(\mathbf{x}, \tau_0) = s(\mathbf{x}, \tau_0)^{4/3} \left(\frac{3}{4}\right)^{4/3} \left(\nu_{\text{QCD}} \frac{\pi^2}{30}\right)^{-1/3} \qquad (9)$$

with a number of degrees of freedom $\nu_{QCD} = 40$. $E_0$ is obtained upon integration of the previous equation following Eq. (8). Figure 2 shows the joint distribution of $E_0$ and R in a sample of ultracentral Au + Au collisions, as defined in Appendix A, where the value of S is essentially fixed. The energy and the size are clearly negatively correlated, but they are not in a one-to-one correspondence. Therefore, if $E_0$ is the same thing as $\langle p_t \rangle$, then $\langle p_t \rangle$ is not the same thing as R.

In this paper, I use $E_0$ as a predictor for $\langle p_t \rangle$ to improve Eq. (6), and thus the results of Ref. [20]. Note that in Fig. 2, and in general throughout this manuscript, the value of S is not strictly fixed, but the considered events belong to a narrow interval of S (as explained in detail in Appendix A). This is a better representation of the centrality selection performed in experiments. To take the small fluctuations of S into account, one can simply replace $E_0$ with $E_0/S$.<fnref n="2" /> The improved Eq. (6) is finally given by

$$\frac{d\langle p_t \rangle}{\langle p_t \rangle} = \kappa_0 \frac{d(E_0/S)}{\langle E_0/S \rangle}, \qquad (10)$$

where $\kappa_0$ is a phenomenological parameter, which can be inferred from experimental data on the relative fluctuation of $\langle p_t \rangle$. The latter was determined recently by the STAR

<footnote n="2">²Note that this choice is not unique. The authors of Ref. [32] show, for instance, that $E_0/A$, where A is the area of the system, works just as well, or better, than $E_0/S$.</footnote>
collaboration [33]. They measured the relative fluctuation of $\langle p_t \rangle$ originating from genuine collective effects, such as those predicted by the TRENTo model. They found that the left-hand side of Eq. (10) is equal to

$$ \frac{\sigma_{\text{dynamical}}(\langle p_t \rangle)}{\langle\langle p_t \rangle\rangle} = 0.012, \qquad (11) $$

in central collisions at top RHIC energy. The relative fluctuation of $E_0/S$, i.e., the right-hand side of Eq. (10), in my TRENTo calculation is instead of order 0.03, so that in practice one has $\kappa_0 \approx 0.4$.

The conclusion of this section is that a natural predictor of $\langle p_t \rangle$ in presence of initial density fluctuations is provided by the initial energy, and not by the initial size. Throughout this manuscript, I assume that the distribution of $\langle p_t \rangle / \langle\langle p_t \rangle\rangle - 1$ at a given centrality is equal to the distribution of $(E_0/S)/\langle E_0/S \rangle - 1$, multiplied by a factor $\kappa_0$ to match the width of the dynamical $\langle p_t \rangle$ distribution observed in data.

### III. NUCLEAR DEFORMATION IN HEAVY-ION COLLISIONS

The goal of this section is to reproduce the results of Ref. [20] with improved theoretical ingredients, in particular, implementing the predictor for the fluctuations of $\langle p_t \rangle$ introduced in the previous section. I establish a phenomenology of nuclear structure in relativistic heavy-ion collisions centered around the fact that $\langle p_t \rangle$ can be used as a tool to access the shape of the colliding nuclear species. This provides the platform for the quantitative results and model-to-data comparisons shown in Sec. IV.

#### A. Modeling the colliding bodies

The modeling of the structure of the colliding ions is simple, and follows typically the Bohr-Mottelson model [1]. The average mass density in the nucleus is defined by

$$ \rho(\mathbf{x}', z') = \frac{\rho_0}{1 + \exp\left\{\frac{1}{a}\left[\sqrt{|\mathbf{x}'|^2 + z'^2} - R\left(1 + \beta Y_2^0\right)\right]\right\}}, \quad (12) $$

where $z'$ is the direction of the nuclear axis, $\mathbf{x}'$ is the plane orthogonal to it. Spherical symmetry is broken by the spherical harmonic $Y_2^0$. R is the nuclear radius, a is the diffusiveness, and $\beta$ is the quadrupole deformation parameter.<fnref n="3" /> In particular, one has [34,35]

$$ \beta \approx \frac{4\pi}{5} \frac{\int d^3\mathbf{r} Y_2^0(\Theta, \Phi) r^2 \rho(\mathbf{r})}{\int d^3\mathbf{r} r^2 \rho(\mathbf{r})}, \quad (13) $$

where $\rho(\mathbf{r})$ corresponds to Eq. (12). If $\beta > 0$, the nucleus is prolate, while it is oblate for $\beta < 0$.

The fluctuations of the mass density are simply obtained by a Monte Carlo sampling of nucleon coordinates [36],

<figure><img src="image_5.png" /><figcaption>FIG. 3. A deformed nucleus with $\beta > 0$ randomly oriented in the laboratory frame. With respect to the intrinsic frame of the nucleus, the ellipsoid is rotated by a polar tilt, $\theta$, and by an azimuthal spin, $\phi$. In a heavy-ion collision experiment, $z$ represents the beam axis, while $(x, y)$ is the transverse plane.</figcaption></figure>

distributed according to Eq. (12), performed on an event-by-event basis. Note that the quadrupole deformation is static, in the sense that the sampling of nucleons is performed from an average density which has always the same value of $\beta$. Doing so one neglects the fact that the $\beta$ can have a significant degree of softness [37,38]. In this paper, I will be concerned mostly with $^{238}$U nuclei, for which this is a good approximation [38], but it would be interesting to include fluctuating values of $\beta$, as well as of a potential triaxial deformation parameter, in collisions of smaller nuclei, such as $^{96}$Ru, $^{96}$Zr, or $^{129}$Xe, for which experimental data are available.

The deformed nucleus is now injected in the beam pipe of a particle accelerator. The orientation of the ellipsoid is random, so that in general the intrinsic nuclear frame $(z', \mathbf{x}')$ differ from the laboratory frame $(z, \mathbf{x})$ by a polar tilt, $\theta$, and by an azimuthal spin, $\phi$. This is illustrated in Fig. 3. Consequently, the collision geometry for two colliding nuclei, A and B, in the laboratory frame is characterized by two polar tilts, $\theta_A$ and $\theta_B$, and two azimuthal rotations, $\phi_A$ and $\phi_B$.

#### B. $\langle p_t \rangle$ and collision geometry: tip-tip vs. body-body collisions

I explain now the argument introduced in Ref. [20] that the average transverse momentum of hadrons produced in a heavy-ion collision can be used as a tool to somehow freeze the orientation of the colliding nuclei. The focus is always on ultracentral collisions, i.e., on collisions where almost all the nucleons from the colliding ions participate in the interaction, and thus the overlap of the two nuclei is maximal. Two kinds of fully overlapping collisions are particularly interesting, as illustrated in Fig. 4.

(i) There are body-body configurations, shown in the left panel of Fig. 4, in which the axes of the two nuclei are both orthogonal to the beam direction, i.e., $\theta_A = \theta_B = \pi/2$, and both nuclei are rotated by the same azimuthal angle, $\phi_A = \phi_B$. As shown in the Fig. 4, in such configurations the area of overlap in the transverse plane has an enhanced elliptical deformation, which originates from the shape of the colliding nuclei.

(ii) There are tip-tip configurations, shown in the right panel of Fig. 4, in which the axis of both nuclei is

<footnote n="3">³Note that it is more customary to truncate the multipole expansion at the hexadecapole deformation, $Y_4^0$, which introduces an additional deformation parameter, $\beta_4$. Here, I neglect the effects of this deformation, which should be unimportant in studies of ultracentral collisions [4].</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Left: illustration of a fully-overlapping body-body collision (top). The nuclei collide along the beam axis, z, with their axes aligned orthogonally to z ($\theta_A = \theta_B = \pi/2$), and with the same azimuthal orientation, $\phi$. The transverse area of nuclear overlap in a body-body collision (bottom) possesses an enhanced quadrupolar asymmetry ($R_y > R_x$), as it follows closely the deformed shape of the colliding bodies. Right: illustration of a fully overlapping tip-tip collision (top). The nuclei collide with $\theta_A = \theta_B = 0$. The area of overlap (bottom) has the symmetry of a disk, $R_x = R_y$.</figcaption></figure>

aligned with the beam axis, $z' = z$, or $\theta_A = \theta_B = 0$. The resulting area of overlap is circular.

The realization of Ref. [20] can be formulated as follows. If we look at fully overlapping events at fixed final-state multiplicity (i.e., fixed initial entropy), then tip-tip collisions correspond to the configurations where the largest densities are achieved. Compared to body-body collisions, a tip-tip event deposits the same amount of entropy in a smaller volume. Consequently, the temperature (and the energy) of the tip-tip event is larger, and thus, if the number of particles is fixed, one naturally expects $\langle p_t \rangle$ to be also larger [Eq. (4)]. The idea is that, in a given sample of events at fixed multiplicity, large values of $\langle p_t \rangle$ are mainly produced by configurations that are tip-tip-like, while low values of $\langle p_t \rangle$ come from body-body configurations.<fnref n="4" />

Following Ref. [20], I provide now an explicit confirmation that this idea works in the TRENTo model. I collide $^{238}$U nuclei at top RHIC energy. In Ref. [3], it is reported that these nuclei have a significant prolate deformation, $\beta = 0.286$, which I round up to 0.3 for simplicity. I focus on ultracentral collisions, as described in Appendix A. I compute the distribution of $(E_0/S)/\langle E_0/S \rangle - 1$, in my sample, and I

turn it into the distribution of $\langle p_t \rangle / \langle \langle p_t \rangle \rangle - 1$ by means of Eq. (10).

In Fig. 5(a), I show the sine of the polar tilts, $\sin \theta_{A,B}$, as a function of the relative variation of $\langle p_t \rangle$. The results are shown as circles. We see that the value of the sine quickly grows to unity as soon as we move towards low values of $\langle p_t \rangle$, which confirms our expectation that low $\langle p_t \rangle$ selects body-body configurations, i.e., $\theta_A = \theta_B = \pi/2$. Note that $\sin \theta$ does not go to zero at the largest values of $\langle p_t \rangle$, meaning that, on average, tip-tip collisions are not achieved.

In Fig. 5(b), I look instead at the alignment of the colliding nuclei in the azimuthal plane, quantified by the correlation $\cos 2(\phi_A - \phi_B)$. This result is shown as circles. Moving towards the low-$\langle p_t \rangle$ tail, I observe that the correlation becomes positive, reaching a value of 0.5. This implies that a selection of fully-overlapping body-body collisions is effectively taking place at low $\langle p_t \rangle$.

I conclude that one can indeed have an handle on the orientation of the colliding bodies by sorting ultracentral events according to their value of $\langle p_t \rangle$. Note that the curves in both Figs. 5(a) and 5(b) are significantly flatter than those shown in Fig. 3 of Ref. [20]. This is due to the fact that I am now consistently taking into account the fact that $R$ and $\langle p_t \rangle$ are not in a one-to-one correspondence.

### C. Anticorrelation between $\varepsilon_2$ and $\langle p_t \rangle$

The second point made in Ref. [20] is that the capability of discerning between body-body and quasi-tip-tip collisions has a very important phenomenological consequence. As illustrated in Fig. 4, body-body collisions have a different geometry than tip-tip collisions, as they are characterized by an elliptical asymmetry, depending on the value of $\beta$. Let me introduce the initial eccentricity of the quark-gluon plasma, defined as the second Fourier harmonic of the entropy density profile [39]:

$$ \varepsilon_2 = \frac{|\int_{\mathbf{x}} \mathbf{x}^2 s(\mathbf{x}, \tau_0)|}{\int_{\mathbf{x}} |\mathbf{x}|^2 s(\mathbf{x}, \tau_0)}, \qquad (14) $$

where $\mathbf{x}^2$ in the numerator should be read in complex notation, $\mathbf{x}^2 = (x + iy)^2$.

The results of Fig. 5, combined with the illustration in Fig. 4, imply that collisions at low $\langle p_t \rangle$, which probe body-body geometries, should have a larger $\varepsilon_2$ than collisions at high $\langle p_t \rangle$. Using my sample of ultracentral U + U collisions from the TRENTo calculation, I plot the rms value of $\varepsilon_2$ as a function of $\langle p_t \rangle$ in Fig. 6. The result is shown as black circles. One sees that $\varepsilon_2$ decreases quickly as a function of $\langle p_t \rangle$, confirming the intuitive picture.

The result in Fig. 6 is one of the most important in this paper, because one can easily relate the initial spatial anisotropy, $\varepsilon_2$, to the measured final-state anisotropy in momentum space, $v_2$, to be defined in Sec. IV, and thus obtain a result expressed solely in terms of measurable observables. That is eventually my goal, but before doing so, an additional ingredient needs to be included in the present analysis.

<footnote n="4">⁴In Appendix B, I show the initial energy-density profile of a tip-tip collision and of a body-body collision as an illustration of the above-mentioned phenomenon for which tip-tip events achieve larger densities of energy (see Fig. 12). Additionally, I carry out hydrodynamic simulations to show explicitly that tip-tip collisions yield larger $\langle p_t \rangle$ than body-body collisions (see Fig. 13).</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Sine of $\theta_A$ and $\theta_B$ as function of $\langle p_t \rangle$ in ultracentral collisions of deformed $^{238}$U nuclei ($\beta = 0.3$) at top RHIC energy. Circles represent the results obtained without the inclusion of statistical fluctuations, as explained in Sec. III D (full circles: $\theta_A$, empty circles: $\theta_B$). Such fluctuations are instead included in the results shown as diamonds (same color coding for A and B). Symbols are shifted horizontally for readability. (b) Azimuthal alignment, $\cos 2(\phi_A - \phi_B)$, of the colliding $^{238}$U nuclei. Diamonds include the effect of statistical fluctuations. See Appendix A for details about the TRENTo implementation leading to these results.</figcaption></figure>

### D. Effect of statistical fluctuations

Evaluating quantities as function of $\langle p_t \rangle$ requires the knowledge of the mean transverse momentum in each event. While this is not an issue in a hydrodynamic simulation, where the output of one event is a continuous spectrum in momentum space, the situation is different in an experiment, where the equivalent of Eq. (2) is given by a discrete average:

$$ \langle p_t \rangle = \frac{1}{N} \sum_{i=1}^{N} p_{t,i}, \qquad (15) $$

where $N$ is the number of particles detected in the event, and $p_{t,i}$ is the transverse momentum of the $i$th particle. The average momentum is evaluated from the sum of a finite number of terms, typically of order 1000 in central

<figure><img src="image_8.png" /><figcaption>FIG. 6. rms eccentricity as function of $\langle p_t \rangle$ in ultracentral U + U collisions at top RHIC energy. Diamonds include the effect of statistical fluctuations explained in Sec. III D.</figcaption></figure>

U + U collisions, and its determination is therefore affected by a statistical error. This implies in particular that, in an experiment, the event-by-event correspondence between $\langle p_t \rangle$ and $E_0$ is smeared by the finite number of detected particles. Statistical fluctuations are naturally proportional to $1/\sqrt{N}$, and, relative to the mean value of $N$, are typically of order 1% for $N \approx 1000$. However, this is as large as the relative dynamical fluctuation of $\langle p_t \rangle$, which is about 1.2% in central collisions at RHIC [33] [see Eq. (11)].

The magnitude of statistical fluctuations can be evaluated from the law of large numbers. In a heavy-ion collisions, particles are emitted independently from the decoupling surface, with a random value of transverse momentum, $p_t$, taken from an underlying probability distribution, i.e., the $p_t$ spectrum, which can be measured at a given collision centrality. Therefore, the fluctuation of the average transverse momentum due to the finite number of particles is given by

$$ \sigma_{\text{stat}} = \frac{1}{\sqrt{N}} \sqrt{\langle p_t^2 \rangle - \langle p_t \rangle^2}, \qquad (16) $$

where the average is weighted with respect to the measured spectrum, $\langle \cdots \rangle = \frac{1}{N} [\int \cdots \frac{dN}{d^2p_t}]$. Considering that at top RHIC energy in the full acceptance of the STAR detector, $|\eta| < 1$, one detects $N \approx 1000$ particles, and using the spectrum measured in central Au + Au collisions [40], I obtain

$$ \sigma_{\text{stat}} = 0.01 \text{ GeV}. \qquad (17) $$

The relative statistical fluctuation is, hence,

$$ \sigma_{\text{stat}} / \langle\langle p_t \rangle\rangle = 0.18, \qquad (18) $$

where I used $\langle\langle p_t \rangle\rangle = 0.57$ GeV [41]. This is the relative fluctuation of $\langle p_t \rangle$ coming from the simple fact that $N$ is not infinite.

This result has to be compared with the magnitude of the relative dynamical fluctuation of $\langle p_t \rangle$. Comparing Eq. (18)
with Eq. (11), one finds

$$ \frac{\sigma_{\text{stat}}(\langle p_t \rangle)}{\langle p_t \rangle} = 1.5 \times \frac{\sigma_{\text{dynamical}}(\langle p_t \rangle)}{\langle p_t \rangle}. \quad (19) $$

Relative statistical fluctuations are, hence, somewhat larger than relative dynamical fluctuations. Any theoretical prediction requiring the evaluation of $\langle p_t \rangle$ on an event-by-event basis is meaningless unless it takes this statistical smearing into account.

Including the effect of these statistical fluctuations in a theoretical calculation is in fact trivial. Since particles are emitted independently at freezeout, their number follows a Poisson distribution, or equivalently, a Gaussian distribution, because $N \gg 1$. Therefore, I can readily correct my TRENTo results for statistical fluctuations. First, I use Eq. (10) to compute the distribution of $\langle p_t \rangle / \langle p_t \rangle$, whose width, say $\sigma$, is given in Eq. (11). Then, I multiply each entry of this distribution by a random number drawn from a Gaussian distribution of unit mean and width given by $1.5 \times \sigma$, following Eq. (19). This produces a distribution of fictitious values of $\langle p_t \rangle / \langle p_t \rangle$, and if I now evaluate the final observables as function of this corrected quantity, the decorrelation between $E_0/S$ and $\langle p_t \rangle$ due the finite-$N$ effect will be properly included.

I implement this correction, and I show its impact on the observables analyzed in this section. In Fig. 5(a), the new results are shown as diamonds. The effect of the statistical smearing is very visible. As expected, it washes out much of the correlation between $\sin \theta_{A,B}$ and $\langle p_t \rangle$, yielding a curve with a smaller slope. An even larger effect is observed in Fig. 5(b) (diamonds), where the alignment between azimuthal angles is almost entirely washed out by the statistical fluctuations. The slope of $\epsilon_n$ as a function of $\langle p_t \rangle$, now shown as diamonds in Fig. 6, is also significantly reduced. Note that statistical fluctuations make the $\langle p_t \rangle$ distribution broader.

However, despite this large correction, the conclusion of Ref. [20] remains valid. The TRENTo model clearly predicts that $\langle p_t \rangle$ is anticorrelated with $\epsilon_2$ in ultracentral collisions of $^{238}$U nuclei with $\beta = 0.3$. Note that one can also define a measure of the correlation between anisotropy and $\langle p_t \rangle$ which is by construction insensitive to the statistical fluctuations analyzed here. This will be discussed in Sec. IV B.

Let me move on, then, to a quantitative phenomenological study.

## IV. CORRELATION BETWEEN $\langle p_t \rangle$ AND $v_n$: QUANTITATIVE ANALYSIS

To exhibit results that can be compared to experimental data, I need to convert $\epsilon_2$ into a quantity that can be measured. This can be done by means of the elliptic flow of the final-state spectrum. In polar coordinates,

$$ \frac{dN}{d^2\mathbf{p}_t} = \frac{dN}{p_t d\mathbf{p}_t d\phi_p}, \quad (20) $$

and elliptic flow is defined as the second (complex) harmonics of the angular part:

$$ V_2 = \frac{1}{N} \int_{\mathbf{p}_t} \frac{1}{2\pi} \frac{dN}{p_t d\mathbf{p}_t} e^{-i2\phi_p}. \quad (21) $$

Elliptic flow originates from a quadrupole-like imbalance of pressure gradient forces, $\vec{F} = -\nabla P$, within the quark-gluon medium [42]. This occurs when the spatial distribution of deposited matter has an elliptical asymmetry. Following Fig. 4, consider a medium that is elongated along a certain transverse direction, e.g., $R_y > R_x$. Pressure gradient forces in the fluid scale with the inverse of the length, so that one has

$$ R_x < R_y \Rightarrow F_x > F_y. \quad (22) $$

Under this condition, then, the hydrodynamic flow builds more momentum along $x$ than along $y$, so that the particles emitted at the end of hydrodynamics have more momentum along $x$ than along the direction orthogonal to it. This phenomenon defines elliptic flow. Hence, body-body collisions in Fig. 4 naturally produce a larger elliptic flow than tip-tip collisions.

Hydrodynamic simulations [43] show that, much as $E_0$ can be used as a predictor of $\langle p_t \rangle$, the value of $v_2 \equiv |V_2|$ in each event can be traced back to the value of $\epsilon_2$, given by Eq. (14), at the beginning of the hydrodynamic expansion. At a given collision centrality, their rms values can be related through a simple linear scaling, $v_2 = \kappa_2 \epsilon_2$. This provides a convenient and powerful method to estimate elliptic flow at a given centrality without running numerically-expensive hydrodynamic simulations, provided that the value of $\kappa_2$ is known. Note that linear scaling is especially valid when $\epsilon_2$ is not large [44], which is typically the case in ultracentral collisions.<fnref n="5" />

The same procedure works as well for the third harmonic. A nonzero initial triangular deformation, $\epsilon_3$, yields a triangular flow of quark-gluon matter, and thus a nonzero third harmonic in the final hadron spectrum, $v_3$ [45]. The two quantities are in a linear relation, $v_3 = \kappa_3 \epsilon_3$. Note that unlike $\epsilon_2$, the triangular deformation is not generated by a symmetry of the problem (unless one collides head-on nuclei that have an intrinsic octupole moment), but solely from density fluctuations.

The remainder of this paper is devoted to a phenomenological study of the correlation between $\langle p_t \rangle$ and $v_n$ in U + U and Au + Au collisions. I evaluate $v_2$ and $v_3$ through linear scaling,

$$ v_n = \kappa_n \epsilon_n, \quad (23) $$

where, depending on the context, $v_n$ indicates either the elliptic flow computed in one event or its rms value at a given $\langle p_t \rangle$. I use

$$ \begin{aligned} \kappa_2[\text{U} + \text{U}] &= 0.165, & \kappa_2[\text{Au} + \text{Au}] &= 0.155, \\ \kappa_3[\text{U} + \text{U}] &= 0.110, & \kappa_3[\text{Au} + \text{Au}] &= 0.100. \end{aligned} \quad (24) $$

The choice of these values is motivated by the results of full hydrodynamic simulations that I carry out in Appendix B. Note that in Ref. [20] the response coefficients were only guessed on the basis of earlier comparisons between the TRENTo model and experimental data. Here, I am explicitly

<footnote n="5">⁵This is not necessarily true in body-body collisions, as I discuss in Appendix B.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 7. rms elliptic flow coefficient, $v_2$, as a function of $\langle p_t \rangle$ in $\sqrt{s} = 200$ GeV U + U and Au + Au collisions. (a) $v_2$ in U + U collisions. (b) $v_2$ in Au + Au collisions. The response coefficient $\kappa_2 = v_2/\epsilon_2$, given by Eq. (24), is specified in each panel. The implementation of the TRENTo model used to obtain these results is detailed in Appendix A.</figcaption></figure>

checking consistency between the employed values of $\kappa_n$ and the hydrodynamic framework of nucleus-nucleus collisions.

### A. $v_n$ as a function of $\langle p_t \rangle$

Figure 7 shows quantitative predictions for the rms elliptic flow in ultracentral U + U and Au + Au collisions at top RHIC energy, as function of the relative variation of $\langle p_t \rangle$. The prediction in Fig. 7(a) corresponds to the results of Fig. 6, rescaled by a factor $\kappa_2 = 0.165$, and including the effect of statistical fluctuations. Figure 7(b) shows instead predictions for Au + Au collisions. The correlation is in this case positive, and this is not surprising. The correlation between mean transverse momentum and anisotropic flow has been recently studied in collisions of spherical $^{208}$Pb nuclei at LHC energy by the ATLAS collaboration [46]. They have indeed observed that $v_2$ and $\langle p_t \rangle$ are positively correlated in ultracentral collisions.

The results in Fig. 7 should, hence, be compared to preliminary STAR data [21]. One finds that my predictions provide a good description of data. The magnitude of elliptic flow is captured within an accuracy of 10%, and the slopes of $v_2$ versus $\langle p_t \rangle$ are also nicely reproduced. In particular, STAR preliminary data present a nontrivial negative slope in U + U collisions, confirming the overall picture of this paper, and so the idea of Ref. [20] about the selection of body-body-like events at low $\langle p_t \rangle$. I stress that my predictions do not have any free parameters. The features of the model are constrained by other sets of data, except for the value of $\beta$, which is taken from nuclear data literature. Note that this calculation implements mildly oblate $^{197}$Au nuclei with $\beta = -0.13$. The argument explained in Sec. III B suggests that the oblate deformation yields a positive slope, enhancing the correlation between $\langle p_t \rangle$ and $v_2$. In Appendix C, I assess the impact of $\beta$ on this result, by studying collisions of spherical $^{197}$Au nuclei.

Figure 8 shows my predictions for the $\langle p_t \rangle$ dependence of triangular flow, $v_3$. As expected from the results of Pb + Pb collisions at LHC energy [46], the correlation between $v_3$ and $\langle p_t \rangle$ is positive in ultracentral collisions. Note that $v_3$ is

slightly larger in U + U collisions, in agreement with STAR data [7]. Note that the comparison between this result and preliminary STAR data [21] is actually excellent. This result is very important. Since $v_3$ is not affected by the deformation parameter, obtaining a good description of data without any additional tuning provides a nontrivial confirmation of the goodness of the model implementation. I emphasize that the correlation between $v_3$ and $\langle p_t \rangle$ predicted by my calculation would be perfectly flat, and so in stark disagreement with data, if I used $R$ as a predictor of $\langle p_t \rangle$, instead of $E_0$, as shown explicitly in Ref. [31].

### B. Removing statistical fluctuations: Božek correlation coefficient

In analyses of heavy-ion collisions, one typically does not evaluate final-state observables, such as $\langle p_t \rangle$ or $v_n$, on

<figure><img src="image_10.png" /><figcaption>FIG. 8. Predictions for the rms $v_3$ as a function of $\langle p_t \rangle$ in ultracentral U + U collisions (diamonds) and Au + Au collisions (circles) at top RHIC energy. The TRENTo implementation that yields these results is detailed in Appendix A.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Predictions for $\rho_n(v_n^2, \langle p_t \rangle)$, as defined by Eq. (27), as a function of collision centrality for $^{197}$Au + $^{197}$Au and $^{238}$U + $^{238}$U collisions at RHIC energy, and $^{208}$Pb + $^{208}$Pb and $^{129}$Xe + $^{129}$Xe at LHC energy. (a) $\rho_2(v_2^2, \langle p_t \rangle)$. (b) $\rho_3(v_3^2, \langle p_t \rangle)$. These results are obtained with the TRENTo model (see Appendix A). For Pb + Pb collisions, experimental data by the ATLAS collaboration are available in Ref. [46], and were compared to the results of the same TRENTo calculation shown here in Ref. [31].</figcaption></figure>

an event-by-event basis, because of the large statistical fluctuations associated with the finite number of particles, as explained in Sec. III D. Observables are typically obtained by statistical averages over many events, constructed in such a way to remove the trivial (and overwhelming) finite-N effects.

For the correlation between $v_n$ and $\langle p_t \rangle$, the corresponding observable was first studied in Ref. [27] in the context of a principal component analysis, and later by Božek [47], who formulated it as a simple correlation coefficient:

$$ \rho_n(v_n^2, \langle p_t \rangle) = \frac{\langle \langle p_t \rangle v_n^2 \rangle - \langle \langle p_t \rangle \rangle v_n^2}{\sigma_{p_t} \sigma_{v_n^2}}, \quad (25) $$

where outer angular brackets denote an average over events in a given multiplicity (centrality) window. The standard deviations, $\sigma_{p_t}$ and $\sigma_{v_n}$, are given by

$$ \sigma_{p_t} = \sqrt{\langle \langle p_t \rangle \rangle^2 - \langle \langle p_t \rangle \rangle^2}, \quad (26) $$

$$ \sigma_{v_n^2} = \sqrt{\langle v_n^4 \rangle - \langle v_n^2 \rangle^2}. $$

In Ref. [31], it has been shown that if one replaces $v_n$ with $\epsilon_n$, and $\langle p_t \rangle$ with $E_0/S$, then the corresponding correlator

$$ \frac{\langle E_0/S \epsilon_n^2 \rangle - \langle E_0/S \rangle \langle \epsilon_n^2 \rangle}{\sigma_{E_0/S} \sigma_{\epsilon_n^2}}, \quad (27) $$

with an appropriate redefinition of the standard deviations in Eq. (26), provides a good description of experimental data on the correlation (25) measured in Pb + Pb collisions by the ATLAS collaboration [46].

The Božek coefficient is extremely robust. Since it does not depend on the value of $\kappa_n$, it is essentially insensitive to the medium properties of the quark-gluon plasma, as confirmed to a great extent by both the results of Ref. [31], and the results shown in Fig. 9 of Ref. [32]. Furthermore, as anticipated, the correlation coefficient in Eq. (25) is by construction insensitive to trivial statistical fluctuations, and isolates the correlations coming from genuine collective effects. This implies that

when I evaluate Eq. (27), I use directly the distribution of $E_0/S$ coming from the TRENTo model. I do not have to rescale the distribution by the factor $\kappa_0$ in Eq. (10), nor to perform the correction described in Sec. III D.

I use Eq. (27) to make quantitative predictions for the Božek correlation coefficient in U + U and Au + Au collisions at top RHIC energy, as function of the collision centrality. The results are presented in Fig. 9. In panel (a), I observe that $\rho_2(v_2^2, \langle p_t \rangle)$ is negative in central U + U collisions, while it is positive in central Au + Au collisions. This effect is compatible with Fig. 7, and is caused by the prolate deformation of $^{238}$U nuclei. Note that the value of $\rho_2(v_2^2, \langle p_t \rangle)$ in U + U collisions is lower than in Au + Au collisions across essentially the full range of centrality. In panel (b) I show instead my predictions for $\rho_3(v_3^2, \langle p_t \rangle)$. This quantity should not be sensitive to the value of $\beta_2$, and I observe, accordingly, only a minor difference between Au + Au and U + U.

It is insightful to look as well at the results for LHC systems, in particular, at the difference between collisions of spherical $^{208}$Pb nuclei and collisions of prolate $^{129}$Xe nuclei. In Fig. 9(a), we note that the comparison between Pb + Pb and Xe+Xe is similar to that found for the RHIC systems. The reason is clearly that $^{129}$Xe nuclei are prolate ($\beta = 0.18$). Equally striking, though, is the result shown in Fig. 9(b), where we observe now four curves that fall essentially into two distinct categories. RHIC systems present $\rho_3(v_3^2, \langle p_t \rangle)$ of order 0.1 and mildly dependent on centrality, while LHC systems present a steep decrease with the centrality percentile. This is a remarkable byproduct of this analysis. The feature of the TRENTo model that drives this result is the fluctuation parameter, $k$, discussed in Appendix A. This quantity determines the magnitude of initial-state fluctuation, and is here tuned in such a way that fluctuations at RHIC energy are essentially twice as large as fluctuations at LHC energy (as explained in Appendix A). Therefore, my results suggest that $\rho_3(v_3^2, \langle p_t \rangle)$ could serve as a fine probe of the evolution of fluctuations from RHIC energy to LHC energy.
<figure><img src="image_12.png" /><figcaption>FIG. 10. rms anisotropic flow coefficients, $v_n$, as a function of $\langle p_t \rangle$ in ultracentral U + U collisions at $\sqrt{s} = 200$ GeV, for different values of the deformation parameter: $\beta = 0.30$ [diamonds, corresponding to Fig. 7(a)], and $\beta = 0.15$ (circles). (a) $v_2$. (b) $v_3$.</figcaption></figure>

### C. Role of the quadrupole deformation

My claim in this paper is that future STAR data will allow us to place a constraint on the value of $\beta$. To check this, I assess the role of $\beta$ by repeating the previous calculations with a quadrupole deformation reduced by a factor 2, i.e., $\beta = 0.15$, in the colliding $^{238}$U nuclei. All the other model parameters are kept fixed.<fnref n="6" />

The resulting flow coefficients as function of $\langle p_t \rangle$ are shown as full circles in Fig. 10. Panel (a) may be considered as the most important result of this paper. I find that, with $\beta = 0.15$, the correlation between $\langle p_t \rangle$ and $v_2$ is positive, and the effect of the deformation is essentially not visible. Panel (b) shows, on the other hand, that modifying $\beta$ does not affect $v_3$. In Fig. 11, I show as well the correlation coefficient of Bożek for two values of $\beta$. The result for $\rho_2(v_2^2, \langle p_t \rangle)$ is in panel (a). The sensitivity of this observable to the value of the quadrupole parameter is impressive. Differences between the curves corresponding to $\beta = 0.15$ and $\beta = 0.3$ are visible up to 40% centrality, whereas usually the effects of deformation on quantities such as the cumulants of elliptic flow are typically inferrable only from central collisions [8]. In panel (b) I show results for $\rho_3(v_3^2, \langle p_t \rangle)$. As expected, the two curves overlap.

Let me stress once more that this calculation has no free parameters. The most uncertain ingredient in Fig. 10(a) is the value of $\kappa_2$, but this quantity changes only the normalization of the results, and not their shape. The results in Fig. 10 suggest then that it would not be possible for state-of-the-art simulations of heavy-ion collision to get a reasonable description

of preliminary STAR data [21] without implementing a value of $\beta$ close to 0.3. My conclusion is that heavy-ion collision data are in agreement with the result of Raman et. al [3], and that $^{238}$U nuclei have indeed a quadrupole deformation parameter of order 0.3. Future investigations of $\rho_2(v_2^2, \langle p_t \rangle)$ at RHIC will further support (or refute) this conclusion.

### V. CONCLUSION AND OUTLOOK

The observation of an anticorrelation between $\langle p_t \rangle$ and $v_2$ in ultracentral U + U collisions provides, arguably, the most striking signature of nuclear deformation so far observed in high-energy nuclear experiments, and should be regarded as compelling evidence of a significant prolate deformation, $\beta \approx 0.3$, in $^{238}$U nuclei.

The ingredient missing in the qualitative analysis of Ref. [20] is a realistic predictor for the event-by-event $\langle p_t \rangle$. The present works fixes this problem by implementing a predictor based on the initial energy, which in turn allows me to carry out quantitative predictions for the phenomenon observed in preliminary STAR data. For large nuclei with a significant prolate deformation, such as $^{238}$U, spectacular results can be obtained by measuring $v_2$ as a function of $\langle p_t \rangle$ [Fig. 7(a)], although one has to carefully consider the impact of statistical fluctuations. The Bożek correlation coefficient [Fig. 9(a)] is instead insensitive to trivial finite-$N$ effects, and also largely independent of the details of the hydrodynamical modeling. At the same time, it displays an outstanding sensitivity to the value of $\beta$ [Fig. 11(a)], and, therefore, it should be investigated at RHIC, and in Xe + Xe collisions at LHC.

My conclusion is that, if we look at observables that have a great sensitivity to the quadrupole deformation of the colliding nuclei, then we can use relativistic nuclear collisions to constrain the value of $\beta$. The correlation between $v_2$ and $\langle p_t \rangle$ provides such an observable.

The realization of Ref. [20] concerning the role of $\langle p_t \rangle$ for the phenomenology of nuclear deformation opens a new direction for studies of nuclear structure at high-energy nu-

<footnote n="6">⁶Note that, while varying $\beta$, one should make sure that the resulting parametrization of the nuclear density remains consistent with the observed transition probability of the electric quadrupole operator, as explained in Ref. [4], and otherwise retune the other shape parameters, *a* and *R*, to meet this condition. For the parametrizations used in this paper, I have not performed any such checks. However, this issue will have to be taken into proper account in future quantitative extractions of the value of $\beta$ from experimental data.</footnote>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Predictions for the Božek correlation coefficient, $\rho_n(v_n^2, \langle p_t \rangle)$, as a function of collision centrality for U + U collisions at top RHIC energy, utilizing two distinct values of $\beta$. (a) $n = 2$. (b) $n = 3$. Solid lines: $\beta = 0.3$. Dashed lines: $\beta = 0.15$.</figcaption></figure>

clear colliders. I expect more observables based on $\langle p_t \rangle$, and sensitive to the deformation of the colliding nuclei, to be invented in the near future. One should also consider the possibility of performing new experiments with new species. Heavy-ion collisions could be used to provide independent confirmations of the results of low-energy experiments, as well as to put constraints on the value of $\beta$ for nuclei that have not yet been investigated experimentally. Such studies would be feasible thanks to the great versatility of the RHIC machine.

## ACKNOWLEDGMENTS

I thank Shengli Huang and Prithwish Tribedy for their continued interest, and Jean-Yves Ollitrault for the careful reading of the manuscript. Useful discussions with Thomas Duguet, Fernando Gardim, Jiangyong Jia, Matt Luzum, Aleksas Mazeliauskas, Jean-Yves Ollitrault, Björn Schenke, Vittorio Somà, and Derek Teaney are kindly acknowledged.

## APPENDIX A: THE TRENTO MODEL

The TRENTO model [48] is a flexible generalization of the Glauber Monte Carlo model [36]. The modeling of the collision process starts with a random sampling of nucleons within the nuclear volume. The coordinates of the nucleons are generated according to the two-parameter Fermi density shown in Eq. (12):

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

where $r = \sqrt{x^2 + y^2 + z^2}$, $a$ and $R_0$ are the diffusiveness and radius of the nucleus, respectively, the spherical harmonic $Y_{20} = \sqrt{\frac{5}{16\pi}}(3 \cos^2 \Theta - 1)$ induces a dependence of the density on the solid angle $\Theta$ in the intrinsic nuclear frame, while the parameter $\beta$ quantifies the magnitude of the ellipsoidal deformation. The parameters used to model the nuclei considered in this paper are reported in Table I.

Once the transverse coordinates, $(x, y)$, of the nucleons are known, the nuclei are overlapped at a random impact parameter, and nucleon-nucleon collisions take place. A participant nucleon carries a Gaussian profile of matter density, $\rho_s(\mathbf{x})$, in the transverse plane:

$$ \rho_{s,i}(\mathbf{x}) = \frac{\omega_i}{2\pi\sigma^2} \exp\left[-\frac{(\mathbf{x} - \mathbf{x}_i)^2}{2\sigma^2}\right], \quad (A2) $$

where the index $i$ labels the $i$th participant nucleon, and I use $\sigma = 0.5$ fm. The weight of each participant, $\omega_i$, is random and follows a $\Gamma$ distribution of unit mean:

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

Note that the variance of this distribution is equal to $k^{-1}$.

The total density in a given nucleus, say A, is given by

$$ s_A(\mathbf{x}) = \sum_i \rho_{s,i}(\mathbf{x}), \quad (A4) $$

while the total entropy density profile for two colliding nuclei is proportional to a generalized mean:

$$ s(\mathbf{x}) \propto \left( \frac{s_A(\mathbf{x})^p + s_B(\mathbf{x})^p}{2} \right)^{1/p}, \quad (A5) $$

where $p$ is any real number. I use $p = 0$, which is the value favored by the comprehensive Bayesian analyses of

<table><caption>TABLE I. Parameters used in Eq. (12) for different species.</caption><thead><tr><th>species</th><th><i>a</i> [fm]</th><th><i>R</i> [fm]</th><th><i>β</i></th></tr></thead><tbody><tr><td><sup>238</sup>U [56]</td><td>0.60</td><td>6.80</td><td>0.30</td></tr><tr><td><sup>208</sup>Pb [56]</td><td>0.55</td><td>6.62</td><td>0</td></tr><tr><td><sup>197</sup>Au [56]</td><td>0.53</td><td>6.40</td><td>-0.13 [5,6]</td></tr><tr><td><sup>129</sup>Xe [57]</td><td>0.59</td><td>5.40</td><td>0.18 [57]</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>FIG. 12. Left: initial energy density profile for a body-body collision, with $S = 4040$, $E_0 = 1294$ GeV, $\epsilon_2 = 0.478$, $\epsilon_3 = 0.191$, and average temperature $\langle T \rangle = 0.433$ GeV. At the end of hydrodynamics, this event yields $dN_{ch}/d\eta_{|\eta|<1} = 1296$, $\langle p_t \rangle = 0.587$ GeV, $v_2 = 0.083$, $v_3 = 0.016$, where the flow coefficients are obtained with the kinematic cuts of the STAR collaboration: $|\eta| < 1$, and $0.2 < p_t < 2$ GeV. Right: a tip-tip collision, with $S = 4072$, $E_0 = 1429$ GeV, $\epsilon_2 = 0.096$, $\epsilon_3 = 0.089$, and average temperature $\langle T \rangle = 0.475$ GeV. This event yields $dN_{ch}/d\eta_{|\eta|<1} = 1280$, $\langle p_t \rangle = 0.651$ GeV, $v_2 = 0.027$, $v_3 = 0.009$.</figcaption></figure>

Refs. [49,50], corresponding to a geometric mean:

$$ s(\mathbf{x}, \tau_0) = \frac{N_0}{\tau_0} \left( \frac{s_A(\mathbf{x})^p + s_B(\mathbf{x})^p}{2} \right)^{1/p} \Big|_{p=0} \\ = \frac{N_0}{\tau_0} \sqrt{s_A(\mathbf{x})s_B(\mathbf{x})}, \qquad (A6) $$

where $\tau_0$ is the initial time, while the constant $N_0$ is added in order for the model to yield the right final-state multiplicity at the end of hydrodynamics. Hydrodynamic simulations are presented in Appendix B, where I implement $\tau_0 = 0.2$ fm/c, $N_0 = 21.6$ to reproduce the charged-particle multiplicity measured at top RHIC energy. Note that the $p=0$ Ansatz for the entropy density yields a correlation between anisotropy and multiplicity at fixed number of participant nucleons that is consistent with that reported by the STAR collaboration in ultracentral U + U collisions [7,48], thus making it a naturally viable prescription for modeling collisions of deformed nuclei.

I discuss now the choice of the fluctuation parameter, $k$. Substantial experimental evidence, coming from analyses of anisotropic flow fluctuations, multiplicity fluctuations, and dynamical $\langle p_t \rangle$ fluctuations, points to the fact that initial-state fluctuations are larger at RHIC energy than at LHC energy (see, e.g., Ref. [51] for a short review). I follow, accordingly, the phenomenological applications of Refs. [52,53], and implement $k = 0.5$ for collisions at RHIC energy, and $k = 2$ for collisions at LHC energy. This implies that, in my model, initial-state fluctuations at RHIC energy are essentially larger than at LHC energy by a factor 2.

Finally, the TRENTTo model is used to sort events into centrality classes. I generate about $10^7$ events for each collision system, and evaluate the corresponding distribution of the total entropy per unit rapidity, $S = \tau_0 \int_x s(\mathbf{x}, \tau_0)$. Events are then classified according to their value of $S$ to mimic the centrality selection performed in experiment, where events are classified according to their final-state multiplicity [54,55]. Ultracentral events at top RHIC energy have typically

$S \sim 4000$ at midrapidity. The class defining ultracentral events in this paper, 0.4–0.8 %, corresponds to $4072 < S < 4231$ for U + U collisions, and to $3419 < S < 3559$ for Au + Au collisions. The hydrodynamic simulations of U + U collisions performed in Appendix B evolve events taken from the 0.78–0.96 % class, which corresponds to $4023 < S < 4078$.

## APPENDIX B: HYDRODYNAMIC SIMULATIONS

In this Appendix I perform hydrodynamic simulations of central U + U collisions with a threefold purpose:

(i) Confirm that $\langle p_t \rangle$ is strongly correlated with $E_0$ in viscous hydrodynamics.

(ii) Confirm that body-body collisions yield lower $\langle p_t \rangle$ than tip-tip collisions.

(iii) Confirm that the response coefficients $\kappa_n$ used in the phenomenological applications of this paper [see Eq. (24)] are in agreement with the results of full hydrodynamic simulations.

To do this, I evolve hydrodynamically 60 profiles of entropy density coming from the TRENTTo simulations described in Appendix A. The calculation is carried out for both body-body and tip-tip configurations. The initial profiles are selected in a narrow bin of entropy corresponding to the 0.78–0.96 % centrality class, as explained in Appendix A.

It is instructive to have a look at an example of initial density profile for both a body-body and a tip-tip event. In the left panel of Fig. 12, I show the initial energy density profile, $e(\mathbf{x}, \tau_0)$, for a body-body collision. On the right panel, I show instead the energy density profile of a tip-tip collision. One immediately observes the different global geometry of these profiles, with an enhanced elliptical asymmetry in the body-body event. Secondly, we can infer essentially by eye that the tip-tip collision covers a smaller area in the transverse plane, and presents on average a larger value of energy density
<figure><img src="image_15.png" /><figcaption>FIG. 13. Average transverse momentum, $\langle p_t \rangle$ as a function of the initial energy, $E_0$, in hydrodynamic simulations of central U + U collisions at top RHIC energy. The charged-particle multiplicity in the $|\eta| < 1$ window corresponds approximately to 1% centrality in the STAR analysis. Full symbols: body-body collisions. Empty symbols: tip-tip collisions.</figcaption></figure>

(or temperature). This is the reason why the tip-tip event yields a larger value of $\langle p_t \rangle$ at the end of the hydrodynamic expansion.

The boost-invariant hydrodynamic evolution of the initial profiles is carried out by means of the MUSIC hydrodynamic code [58–60]. In my setup, I neglect the pre-equilibrium phase [61,62], and I start hydrodynamics at initial time $\tau_0 = 0.2$ fm/c. The hydrodynamic expansion is viscous. I implement a constant value of the specific shear viscosity, $\eta/s = 0.16$, and a temperature-dependent specific bulk viscosity, $\zeta/s$. The current status of the implementation of $\zeta/s$ in heavy-ion collisions is fairly uncertain [63]. On the one hand, calculations starting with IP-Glasma [64] initial conditions

implement a bulk viscosity that peaks around $\zeta/s \approx 0.3$, to yield values of $\langle p_t \rangle$ in agreement with experimental data, at both RHIC and LHC energy [65–67]. On the other hand, calculations starting with TRENTo initial condition are able to fit data on $\langle p_t \rangle$ with a $\zeta/s$ that is smaller by essentially one order of magnitude ($\zeta/s \approx 0.03$ at the peak [49,50]), although results are available only at LHC energy. Here, I work at RHIC energy starting with TRENTo initial conditions, so I set up a sort of hybrid scenario: I implement a bulk viscosity having the same functional form as in Refs. [66,68], but I reduce the value of $\zeta/s$ at its peak by a factor 10. The shear and bulk corrections to the momentum distribution functions at freezeout are the same as in Ref. [66].

The medium has the equation of state of lattice QCD [69], and fluid cells convert into hadrons when their temperature falls below $T = 0.15$ GeV. At freezeout, all hadronic resonances can be formed [70,71], and their decays [72] to stable hadrons are taken into account. The outcome is a final boost-invariant spectrum of charged hadrons, $\frac{dN}{d^2p_t}$, which is then integrated to calculate the charged-particle multiplicity, the average transverse momentum, and the flow coefficients, see, e.g., Eqs. (4) and (21).

I start by showing the values of $\langle p_t \rangle$ as a function of $E_0$ in Fig. 13. We observe the anticipated nearly one-to-one correspondence between these quantities. The correlation is essentially as strong as that observed in Refs. [30,31], meaning that the simple picture of $\langle p_t \rangle$ as a measure of the energy per particle in the initial state is only mildly disrupted by the viscous corrections. The second remarkable feature displayed by the results of Fig. 13 is that body-body and tip-tip events, though falling essentially on the same curve, cover distinct intervals in $E_0$ and $\langle p_t \rangle$, which confirms the overall picture of this paper, i.e., that $\langle p_t \rangle$ provides a powerful handle on the orientation of the colliding ellipsoids. Note that the correlation in Fig. 13 would be even stronger if one replaced $E_0$ with $E_0/S$, which is in fact the predictor used in this paper.

<figure><img src="image_16.png" /><figcaption>FIG. 14. Flow coefficients, $v_n$, as a function of the initial spatial anisotropy, $\epsilon_n$, in high-multiplicity U + U collisions. The flow coefficients are evaluated within kinematic cuts $|\eta| < 1$ and $0.2 < p_t < 2$ GeV, following the STAR collaboration. Shaded bands are meant to provide the range of viable values for the response coefficients $\kappa_n = v_n/\epsilon_n$. Dashed lines correspond to the values implemented in Figs. 7 and 8. Left panel: elliptic flow, $v_2$. Right panel: triangular flow, $v_3$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIG. 15. Left: $v_2$ as a function of the relative variation of $\langle p_t \rangle$ in ultracentral Au + Au collisions at top RHIC energy. Full circles are obtained with oblate $^{197}$Au nuclei ($\beta = -0.13$), while triangles are obtained with spherical nuclei. Right: Božek coefficient, $\rho_2(v_n^2, \langle p_t \rangle)$. Solid line: $\beta = -0.13$. Dashed line: $\beta = 0$.</figcaption></figure>

Finally, I calculate flow coefficients, $v_n$, that I plot as functions of $\epsilon_n$ in Fig. 14.

Figure 14(a) shows results for elliptic flow. The first notable feature is that, as expected, body-body and tip-tip collision cover essentially two distinct intervals in $\epsilon_2$ and $v_2$, as discussed in Sec. IIIC. The shaded band shows a range of viable values for the response coefficient $\kappa_2 = v_2/\epsilon_2$. The dashed line corresponds to $\kappa_2 = 0.165$, which is the value suggested in Ref. [53], and used here in Figs. 7 and 10. One sees that $\kappa_2 = 0.165$ captures a significant fraction of the body-body points, but only a minor fraction of tip-tip events, meaning that the prediction shown in Fig. 7 is smaller than the full hydrodynamic result by about 10%. Note also that, contrary to Fig. 13, body-body and tip-tip collisions do not fall on the same curve. The response coefficient is slightly larger, by a few percents, in the tip-tip case. This is consistent with previous studies showing that the linear response coefficient decreases when the initial-state anisotropy is large [44,73], as it occurs, e.g., moving from central to peripheral nucleus-nucleus collisions. The coefficient $\kappa_2$ presents, as a consequence, a small dependence on $\langle p_t \rangle$, which should however be negligible in Fig. 7 thanks to the strong smearing due to statistical fluctuations.

Figure 14(b) shows, on the other hand, results for triangular flow. Body-body and tip-tip collisions overlap to a good extent, although both $\epsilon_3$ and $v_3$ are larger in tip-tip events. The value of the response coefficient chosen in Figs. 8 and 10, $\kappa_3 = 0.110$, provides a good description of the relation between $\epsilon_3$ and $v_3$ observed in Fig. 14.

A final word about the response coefficients, $\kappa_2$ and $\kappa_3$, employed in Figs. 7 and 8 for Au + Au collisions. I do not estimate them by means of hydrodynamic simulations. The value $\kappa_2 = 0.155$ is suggested by Ref. [53], where it

is used to match the second-order cumulant of the eccentricity, $\epsilon_2\{2\}$, computed in the same TRENTo model used here, to STAR data on the second-order cumulant of elliptic flow, $v_2\{2\}$, in ultracentral Au + Au collisions. I play the same game for the third harmonic. In STAR ultracentral data [7], $v_3\{2\} = 0.0138$, while my TRENTo calculation yields $\epsilon_3\{2\} = 0.1330$, hence $\kappa_3 \simeq 0.100$. Note that these coefficients are both smaller than those implemented in U + U collisions. This is consistent with the fact that $\kappa_n$ is damped by viscous corrections, which are larger in Au + Au collisions due to the smaller system size.

## APPENDIX C: EFFECT OF OBLATE DEFORMATION IN $^{197}$Au

I assess the effect of $\beta$ on the correlation between $v_2$ and $\langle p_t \rangle$ in Au + Au collisions. This analysis is motivated by Ref. [53], where it was found that a deformation parameter $\beta = -0.13$ in the colliding $^{197}$Au nuclei is at variance with the experimental result that the fourth-order cumulant of elliptic flow, $v_2\{4\}^4$, is negative in central Au + Au collisions [7]. It is relevant, hence, to figure out whether $\beta$ plays any role for the observable studied in this article.

I look at the correlation between $\langle p_t \rangle$ and $v_2$, implementing both oblate and spherical $^{197}$Au nuclei. Figure 15(a) shows $v_2 = 0.155\epsilon_2$ as a function of $\langle p_t \rangle$. As expected, $v_2$ is smaller in magnitude in the calculation implementing spherical nuclei. However, the slope of the curve is just as large as with deformed nuclei. The effect of $\beta$ can be simply accounted, then, by a rescaling of the value of $\kappa_2$ by few percents. This is also confirmed in Fig. 15(b), where the Božek coefficient, independent of $\kappa_2$, turns out to be insensitive to the value of $\beta$.

<footnote>[1] A. Bohr and B. Mottelson, *Nuclear Structure*, Vol. II (Benjamin, Reading, MA, 1975).</footnote>

<footnote>[2] D. Cline, Annu. Rev. Nucl. Part. Sci. 36, 683 (1986).</footnote>
[3] S. Raman, C. W. G. Nestor, Jr., and P. Tikkanen, At. Data Nucl. Data Tables 78, 1 (2001).

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

[5] S. Hilaire and M. Girod, Eur. Phys. J. A 33, 237 (2007).

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

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

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

[9] A. M. Sirunyan et al. (CMS Collaboration), Phys. Rev. C 100, 044902 (2019).

[10] G. Aad et al. (ATLAS Collaboration), Phys. Rev. C 101, 024906 (2020).

[11] E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017).

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

[13] P. Filip, R. Lednicky, H. Masui, and N. Xu, Phys. Rev. C 80, 054903 (2009).

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

[15] S. A. Voloshin, Phys. Rev. Lett. 105, 172301 (2010).

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

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

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

[19] U. Heinz and R. Snellings, Annu. Rev. Nucl. Part. Sci. 63, 123 (2013).

[20] G. Giacalone, Phys. Rev. Lett. 124, 202301 (2020).

[21] S. Huang (STAR Collaboration), Contribution to the 36th Winter Workshop on Nuclear Dynamics, Puerto Vallarta, Mexico, March 2020, https://indico.cern.ch/event/841247/contributions/3740391/.

[22] J. D. Bjorken, Phys. Rev. D 27, 140 (1983).

[23] J.-Y. Ollitrault, Eur. J. Phys. 29, 275 (2008).

[24] F. G. Gardim, G. Giacalone, M. Luzum, and J.-Y. Ollitrault, Nat. Phys. 16, 615 (2020).

[25] W. Broniowski, M. Chojnacki, and Ł. Obara, Phys. Rev. C 80, 051902(R) (2009).

[26] P. Bożek and W. Broniowski, Phys. Rev. C 85, 044910 (2012).

[27] A. Mazeliauskas and D. Teaney, Phys. Rev. C 93, 024913 (2016).

[28] P. Bożek and W. Broniowski, Phys. Rev. C 96, 014904 (2017).

[29] F. G. Gardim, G. Giacalone, and J. Y. Ollitrault, arXiv:1909.11609 [nucl-th].

[30] F. G. Gardim, G. Giacalone, M. Luzum, and J. Y. Ollitrault, arXiv:2002.07008 [nucl-th].

[31] G. Giacalone, F. G. Gardim, J. Noronha-Hostler, and J. Y. Ollitrault, arXiv:2004.01765 [nucl-th].

[32] B. Schenke, C. Shen, and D. Teaney, arXiv:2004.00690 [nucl-th].

[33] J. Adam et al. (STAR Collaboration), Phys. Rev. C 99, 044918 (2019).

[34] K. Kumar, Phys. Rev. Lett. 28, 249 (1972).

[35] T. Schmidt, K. L. G. Heyde, A. Blazhev, and J. Jolie, Phys. Rev. C 96, 014302 (2017).

[36] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Annu. Rev. Nucl. Part. Sci. 57, 205 (2007).

[37] M. T. Mustonen, C. N. Gilbreth, Y. Alhassid, and G. F. Bertsch, Phys. Rev. C 98, 034317 (2018).

[38] A. Poves, F. Nowacki, and Y. Alhassid, Phys. Rev. C 101, 054307 (2020).

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

[40] B. B. Back et al. (PHOBOS Collaboration), Phys. Lett. B 578, 297 (2004).

[41] J. Adams et al. (STAR Collaboration), Phys. Rev. C 72, 044902 (2005).

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

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

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

[45] B. Alver and G. Roland, Phys. Rev. C 81, 054905 (2010); 82, 039903(E) (2010).

[46] G. Aad et al. (ATLAS Collaboration), Eur. Phys. J. C 79, 985 (2019).

[47] P. Bożek, Phys. Rev. C 93, 044908 (2016).

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

[49] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 101, 024911 (2020).

[50] J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nat. Phys. 15, 1113 (2019).

[51] G. Giacalone, F. Gelis, P. Guerrero-Rodríguez, M. Luzum, C. Marquet, and J. Y. Ollitrault, arXiv:1911.04720 [nucl-th].

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

[53] G. Giacalone, Phys. Rev. C 99, 024910 (2019).

[54] B. Abelev et al. (ALICE Collaboration), Phys. Rev. C 88, 044909 (2013).

[55] M. Aaboud et al. (ATLAS Collaboration), J. High Energy Phys. 01 (2020) 051.

[56] H. De Vries, C. W. De Jager, and C. De Vries, At. Data Nucl. Data Tables 36, 495 (1987).

[57] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B 790, 35 (2019).

[58] B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010).

[59] B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 106, 042301 (2011).

[60] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 93, 044906 (2016).

[61] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney, Phys. Rev. Lett. 122, 122302 (2019).

[62] S. Schlichting and D. Teaney, Annu. Rev. Nucl. Part. Sci. 69, 447 (2019).

[63] M. Byres, S. H. Lim, C. McGinn, J. Ouellette, and J. L. Nagle, Phys. Rev. C 101, 044902 (2020).

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

[65] S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Phys. Rev. C 95, 064913 (2017).

[66] S. Ryu, J.-F. Paquet, C. Shen, G. Denicol, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 97, 034910 (2018).

[67] B. Schenke, C. Shen, and P. Tribedy, Phys. Rev. C 99, 044908 (2019).
[68] S. Ryu, J.-F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 115, 132301 (2015).

[69] S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabó, Phys. Lett. B 730, 99 (2014).

[70] W. Broniowski and W. Florkowski, Phys. Rev. Lett. 87, 272302 (2001).

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

[72] A. Mazeliauskas, S. Floerchinger, E. Grossi, and D. Teaney, Eur. Phys. J. C 79, 284 (2019).

[73] M. Sievert and J. Noronha-Hostler, Phys. Rev. C 100, 024904 (2019).