# Shape of atomic nuclei in heavy ion collisions

Jiangyong Jia<EMAIL_ADDRESS>Department of Chemistry, Stony
Brook University, Stony Brook, NY 11794, USA  Physics Department, Brookhaven
National Laboratory, Upton, NY 11976, USA

###### Abstract

In the hydrodynamic model description of heavy ion collisions, the final-state
anisotropic flow $v_{n}$ are linearly related to the strength of the multi-
pole shape of the distribution of nucleons in the transverse plane
$\varepsilon_{n}$, $v_{n}\propto\varepsilon_{n}$. The $\varepsilon_{n}$, for
$n=1,2,3,4$, are sensitive to the shape of the colliding ions, characterized
by the quadrupole $\beta_{2}$, octupole $\beta_{3}$ and hexadecapole
$\beta_{4}$ deformations. This sensitivity is investigated analytically and
also in a Monte Carlo Glauber model. One observes a robust linear relation,
$\langle\varepsilon_{n}^{2}\rangle=a_{n}^{\prime}+b_{n}^{\prime}\beta_{n}^{2}$,
for events in a fixed centrality. The $\langle\varepsilon_{1}^{2}\rangle$ has
a contribution from $\beta_{3}$ and $\beta_{4}$, and
$\langle\varepsilon_{3}^{2}\rangle$ from $\beta_{4}$. In the ultra-central
collisions, there are little cross contributions between $\beta_{2}$ and
$\varepsilon_{3}$ and between $\beta_{3}$ and $\varepsilon_{2}$, but clear
cross contributions are present in non-central collisions. Additionally,
$\langle\varepsilon_{n}^{2}\rangle$ are insensitive to non-axial shape
parameters such as the triaxiality. This is good news because the measurements
of $v_{2}$, $v_{3}$ and $v_{4}$ can be used to constrain simultaneously the
$\beta_{2}$, $\beta_{3}$, and $\beta_{4}$ values. This is best done by
comparing two colliding ions with similar mass numbers and therefore nearly
identical $a_{n}^{\prime}$, to obtain simple equation that relates the
$\beta_{n}$ of the two species. This opens up the possibility to map the shape
of the atomic nuclei at a timescale ($<10^{-24}$s) much shorter than probed by
low-energy nuclear structure physics ($<10^{-21}$s), which ultimately may
provide information complementary to those obtained in the nuclear structure
experiments.

###### pacs:

25.75.Gz, 25.75.Ld, 25.75.-1

###### Contents

  1. I Introduction
  2. II Analytical estimate in ultra-central collisions and Glauber model setup
  3. III Results
    1. III.1 Influence of axial-symmetric multipole deformation
    2. III.2 Influence of non-axial deformations
    3. III.3 Simultaneous presence of quadrupole-octupole-hexadecapole deformations
    4. III.4 Constraining quadrupole, octupole, and hexadecapole deformations using collective flow data
  4. IV Discussion and summary
  5. A An analytical estimate in head-on collisions
  6. B More detailed results

##  I Introduction

Most of the atomic nuclei in their ground state are deformed from a well
defined spherical shape. The deformation has non-trivial dependence on the
proton and neutron number, especially in the vicinity of full shell or
subshell, reflecting collective motion induced by interaction between valence
nucleons and shell structure Bohr and Mottelson (1998). The collective motion
leads to characteristic rotational spectra of nuclear excited state, where the
electric multi-pole transition probability $B(En)$ between low-lying
rotational states with $n\hbar$ difference in angular momentum can be used to
infer the shape parameters. Past efforts have led to the discovery of a rich
variety of phenomena, such as quadrupole deformation, shape evolution,
triaxiality/shape coexistence, octupole deformation, hexadecapole deformation
and other exotic shapes Heyde and Wood (2011); Togashi  _et al._ (2016); Heyde
and Wood (2016); Frauendorf (2018); Zhou (2016).

No-one has directly observed the deformed nucleus, however. This is because
the nucleus is deformed in the so-called intrinsic (body-fixed) frame, and its
wave function in the laboratory frame actually does not pick a particular
direction. The typical scattering experiments probe the nuclear form factors
averaged over all orientations, and the static deformation appears mostly as
an increased surface thickness De Vries  _et al._ (1987). On the other hand,
high-energy heavy ion collisions at RHIC and the LHC, as illustrated in Fig.
1, can image the shape of the nucleus by colliding them together and looking
at the collective expansion of the produced system responding to the geometry
of the overlap. In these collisions, two Lorentz-contracted nuclei, by a
factor of 100 at RHIC and more than a factor of 1000 at the LHC, cross each
other over a time scale $\tau<0.1$fm/$c\approx 3\times 10^{-24}$s, forming a
hot and dense quark-gluon plasma (QGP) Busza  _et al._ (2018) in the overlap
region, whose initial shape is correlated with the deformed shape of the
nuclei. Driven by the large pressure gradient forces, the QGP expands
hydrodynamically, converting the spatial anisotropies into azimuthal
anisotropies of final-state particles in the momentum space Heinz (2013).
Nuclear shape imaging is possible because each collision probes simultaneously
the entire mass distribution of the nuclei, and one can use particle
correlations among thousands of produced particles to infer the two-point and
multi-point correlations of this mass distribution and hence its spatial
shape. Since the time scales involved in these collisions are much shorter
($<10^{-24}$s) than the typical timescale of the rotational bands ($10^{-21}$s
Nakatsukasa  _et al._ (2016)), this raises an important question of whether
the manifestation of nuclear deformation – a collective feature of the nuclear
many-body system – is the same across energy scales Giacalone  _et al._
(2021a).

Figure 1:  The cartoon of collision of nuclei with quadrupole (left), octupole
(middle) and hexadecapole (right) deformations including either the $Y_{n,0}$
mode (top row) or the $Y_{n,n}$ mode (middle row) and with $\beta_{n}=0.25$.
The Lorentz contraction in the $z$-direction, by factor of 100 at RHIC and
more than a factor of 1000 at the LHC, are not shown. The bottom row shows how
the initial condition of the QGP formed after the collision looks like in the
transverse plane. The hallow arrows indicate the direction of maximum pressure
gradients along which the medium expand with largest velocity, leading to
final state harmonic flow $v_{n}$ with $n$-fold symmetry.

The shape of the nucleus in nuclear physics is often modeled though a nucleon
density profile of the Woods-Saxon form,

$\displaystyle\rho(r,\theta,\phi)=\frac{\rho_{0}}{1+e^{\left[r-R(\theta,\phi)/a_{0}\right]}},\;R(\theta,\phi)=R_{0}\left(1+\beta_{2}[\cos\gamma
Y_{2,0}+\sin\gamma
Y_{2,2}]+\beta_{3}\sum_{m=-3}^{3}\alpha_{3,m}Y_{3,m}+\beta_{4}\sum_{m=-4}^{4}\alpha_{4,m}Y_{4,m}\right),$
(1)

where the nuclear surface $R(\theta,\phi)$ is expanded in spherical harmonics
$Y_{n,m}$, but keeping only the terms that are most relevant in nuclear
structure physics, i.e. quadrupole $n=2$, octupole $n=3$ and hexadecapole
$n=4$. Note that the $Y_{2,-1},Y_{2,1}$ and $Y_{2,-2}$ are used to define the
intrinsic frame, leaving $Y_{2,0}$ and $Y_{2,2}$ as the only relevant
quadrupole components (higher order deformations defined in this frame should
have all components as relevant degrees of freedom). The positive number
$\beta_{2}$ describes the overall quadrupole deformation, and the triaxiality
parameter $\gamma$ controls the relative order of the three radii
$R_{a},R_{b},R_{c}$ of the nucleus in the intrinsic frame. It has the range
$0\leq\gamma\leq\pi/3$, with $\gamma=0$, $\gamma=\pi/3$, and $\gamma=\pi/6$
corresponding, respectively, to prolate ($R_{a}=R_{b}<R_{c})$, oblate
($R_{a}<R_{b}=R_{c}$) or maximum triaxiality ($R_{a}<R_{b}<R_{c}$ and
$2R_{b}=R_{a}+R_{c}$). Similarly, $\beta_{3}$ and $\beta_{4}$ control the
overall octupole and hexadecapole deformations, respectively. The
$\alpha_{3,m}$ and $\alpha_{4,m}$, in analogue to $\gamma$, are internal
“angular” parameters describing deviation from axial and/or reflection
symmetry, and they satisfy the normalization condition
$\sum_{m=-3}^{3}\alpha_{3m}^{2}=1$ and $\sum_{m=-4}^{4}\alpha_{4m}^{2}=1$.

In heavy ion collisions, the initial condition and dynamics of QGP are
naturally formulated in a cylindrical coordinate system with the $z$-axis
coincide with the beam-line. The initial condition is determined by the
distribution of nucleons in the transverse plane $\rho(r_{\perp},\phi)$, which
drives the collective flow of final-state particles, reflected by the momentum
spectra $N(p_{\mathrm{T}},\phi)$. The $N(p_{\mathrm{T}},\phi)$ is often
analyzed in terms of a Fourier expansion $dN/d\phi\propto
1+2\sum_{n}v_{n}(p_{\mathrm{T}})\cos\,n(\phi-\Psi_{n}(p_{\mathrm{T}}))$. The
$\rho(r_{\perp},\phi)$ is fully characterized via a 2D multi-pole expansion,
whose leading radial modes have the following expression C.G.Joslin and Gray
(1983); Teaney and Yan (2011) in the center of mass frame,

$\displaystyle\varepsilon_{1}e^{i\Phi_{1}}=-\int
d^{2}r_{\perp}r_{\perp}^{3}e^{i\phi}\rho(\vec{r}_{\perp})/\int
d^{2}r_{\perp}r_{\perp}^{3}\rho(\vec{r}_{\perp})\;,\;\varepsilon_{n}e^{in\Phi_{n}}|_{n>1}=-\int
d^{2}r_{\perp}r_{\perp}^{n}e^{in\phi}\rho(\vec{r}_{\perp})/\int
d^{2}r_{\perp}r_{\perp}^{n}\rho(\vec{r}_{\perp}),$ (2)

The 2D eccentricity vectors $\varepsilon_{n}e^{in\Phi_{n}}$ are a close
analogue of $Y_{n,m}$ in the 3D. In fact, the eccentricity vectors for $n>1$
are directly related to the multi-pole moments of the mass distribution,
$\varepsilon_{n}e^{in\Phi_{n}}\propto-\left\langle Y_{n}^{n}\right\rangle$,
which I will show later lead to a simple relation between $\varepsilon_{n}$
and $\beta_{n}$. Note that the radial weight of $\varepsilon_{1}$ is
$r_{\perp}^{3}$ instead of the naively expected $r_{\perp}$, because the
latter contribution vanishes in the center of mass frame and $r_{\perp}^{3}$
weighting gives the next radial mode.

Study of the relation between the initial-state $\varepsilon_{n}$ and final-
state $v_{n}$, within the relativistic viscous hydrodynamics or transport
model framework, has always been one central focus of the heavy ion community.
Comprehensive model and data comparisons Niemi  _et al._ (2016) show very good
linear relations, $v_{n}=k_{n}\varepsilon_{n}$ not only on average but also
for each event 111The linear relation is very good for $n=2$ and 3 in general
and for $n=1$ and $n=4$ in the case of central collisions Teaney and Yan
(2012)., where the response coefficients $k_{n}$ capture the transport
properties of the QGP produced in the collision. Thanks to the precision
measurements of $v_{n}$ and its event-by-event fluctuations $p(v_{n})$ Jia
(2014); Busza  _et al._ (2018), and detailed understanding of the properties
of $k_{n}$ Teaney and Yan (2012); Bernhard  _et al._ (2016, 2019); Nijs  _et
al._ (2021), hydrodynamic models can now determine the $\varepsilon_{n}$ and
$p(\varepsilon_{n})$ with enough precision to constrain the deformation
parameters.

Influence of nuclear deformation on dynamics of heavy ion collisions has been
considered early on Rosenhauer  _et al._ (1986); Li (2000); Gupta and Gale
(2000). More recent studies focused on the relation between $\beta_{2}$ and
$v_{2}$ Heinz and Kuhlman (2005); Filip  _et al._ (2009); Shou  _et al._
(2015); Goldschmidt  _et al._ (2015); Giacalone  _et al._ (2018); Giacalone
(2019). Experimental evidences for quadrupole deformation appear as large
differences of $v_{2}$ between ultra-central collisions (UCC) of different
systems, in particular between 197Au+197Au and 238U+238U collisions at RHIC
Adamczyk  _et al._ (2015) and between 129Xe+129Xe and 208Pb+208Pb collisions
at the LHC ALICE Collaboration (2018); Sirunyan  _et al._ (2019); Aad  _et
al._ (2020). Ref. Giacalone  _et al._ (2021a) explored the parametric
dependence of various flow observables on $\beta_{2}$, and find that both
$\epsilon_{2}^{2}$ and $v_{2}^{2}$ depend linearly on $\beta_{2}^{2}$; a
simple formula is derived relating the $\beta_{2}$ in the two collision
systems to the ratio of $v_{2}$. The influence of octupole deformation is
considered recently in Pb+Pb collisions to explain the order of $v_{2}$ and
$v_{3}$ in the UCC region Carzon  _et al._ (2020).

Another observable showing a strong sensitivity to the nuclear deformation is
the Pearson correlation coefficient, $\rho(v_{2}^{2},[p_{\mathrm{T}}])$,
between $v_{2}$ and the mean transverse momentum, $[p_{\mathrm{T}}]$, which
probes both the $\beta_{2}$ Giacalone (2020a, b) and its triaxiality $\gamma$
Jia  _et al._ (2021) of the colliding ions. Recent measurement from the STAR
collaboration Jia  established unambiguously the large and dominating
influence of the nuclear quadrupole deformation of 238U. The large prolate
deformation of 238U yields a strong negative contribution to the
$\rho(v_{2}^{2},[p_{\mathrm{T}}])$, enough to make it change sign. Large
influence of deformation is also observed in the fluctuations of
$[p_{\mathrm{T}}]$ Jia .

Continuing this line of work, I explore the parametric dependence of
$\epsilon_{1}$, $\epsilon_{2}$, $\epsilon_{3}$ and $\epsilon_{4}$ on various
deformation parameters $\beta_{n}$ and deviations from axial and reflection
symmetries ($\gamma$ and combinations of $\alpha_{n,m}$ in Eq. (1)). In the
UCC region, the mean square (ms) of $\epsilon_{n}$,
$\left\langle\epsilon_{n}^{2}\right\rangle$, are found to be driven primarily
by $\beta_{n}$ of the same order and the triaxiality parameter $\gamma$ only
has very modest impact on $\epsilon_{2}$. Away from the UCC region,
$\left\langle\epsilon_{n}^{2}\right\rangle$ are insensitive to $\gamma$, but
they receive contributions from $\beta_{m}$ of a different order $m\neq n$. In
other words, I establish the following empirical relation,

$\displaystyle\left\langle\epsilon_{n}^{2}\right\rangle=a^{\prime}_{n}+b^{\prime}_{n}\beta_{n}^{2}+\sum_{m\neq
n}b^{\prime}_{n,m}\beta_{m}^{2}\;,\;\;n=1,2,3,4,\;m=2,3,4\;\mathrm{and}\;b^{\prime}_{1}=0\;,$
(3)

with significant values of $b^{\prime}_{n,m}$ observed for $b^{\prime}_{1,3}$,
$b^{\prime}_{1,4}$ and $b^{\prime}_{3,4}$. Since $v_{n}\propto\epsilon_{n}$,
one expects similar parametric dependencies to hold between $\left\langle
v_{n}^{2}\right\rangle$ and $\beta_{m}$. This simple scaling relation provides
a strong motivation for a collision system scan to map out the shape of atomic
nuclei in most interesting region of nuclear chart and compare with the
knowledge from nuclear structure physics.

##  II Analytical estimate in ultra-central collisions and Glauber model setup

To gain some intuitive insight on Eq. (3) it is useful to demonstrate its
validity using a simpler version of nuclear surface. Here I consider a density
distribution of liquid-drop model with a sharp surface:
$\rho(r,\theta,\phi)=\rho_{0}$ when $r<R(\theta,\phi)$ and zero otherwise, and
I assume that the energy density distribution is given by the distribution of
participating nucleons. I limit the discussion to head-on collisions with
nearly maximum overlap i.e. the two nuclei not only have zero impact
parameter, but are also required to align in a way to ensure the overlap
region contains all the nucleons $N_{\mathrm{part}}=2A$. In reality, the
selection of UCC events naturally encompasses a wider range of rotation angles
and also a finite range of $N_{\mathrm{part}}$, therefore I also study a
second case which requires zero impact parameter but independent rotations for
the two nuclei. The details of the calculation can be found in Appendix A.

As illustrated in Fig. 1, the maximum $\varepsilon_{n}$ for $Y_{n,0}$ is
reached in a “body-body” configuration, when the symmetry-axis of the nuclei
is perpendicular to the beam. In the case of $Y_{n,n}$, the maximum
$\varepsilon_{n}$ is reached in a “tip-tip” configuration, when the $z$-axis
of the nuclei is aligned with the beam. For these two configurations, it is
straightforward to calculate eccentricities, they are listed in the first two
rows of Table 1. However, one is more interested in the eccentricity values
averaged over random orientations. The first non-trivial and the most
important moment is $\left\langle\varepsilon_{n}^{2}\right\rangle$, which
relates directly to the $\left\langle v_{n}^{2}\right\rangle$ measured by the
two-particle correlation method. The results obtained by requiring same random
rotations for the two nuclei are listed in the third row of Table 1, and those
obtained by requiring independent random rotations the two nuclei are listed
in the last row of Table 1. The two cases have the same dependencies but the
coefficients are a factor of two smaller in the second case.

| $n=1$ | $n=2$ | $n=3$ | $n=4$
---|---|---|---|---
| $\frac{16}{5\pi\sqrt{7\pi}}\beta_{3}+$ | $\sqrt{\frac{45}{16\pi}}\beta_{2}+\frac{15}{112\pi}\beta_{2}^{2}+$ | $\frac{16}{\pi\sqrt{7\pi}}\beta_{3}+$ | $\frac{35}{16\sqrt{\pi}}\beta_{4}+\frac{45}{16\pi}\beta_{2}^{2}+$
$1+\sum_{m=2}^{4}\beta_{m}Y_{m,0}$
| $\frac{302}{3\sqrt{35}\pi^{2}}\beta_{2}\beta_{3}\\!+\\!\frac{493967\sqrt{7}}{73920\pi^{2}}\beta_{3}\beta_{4}$ | ${\frac{9\sqrt{5}}{7\pi}\\!\beta_{2}\beta_{4}+\\!\frac{75}{77\pi}\\!\beta_{4}^{2}\\!+\\!\frac{1}{\pi}\beta_{3}^{2}}$ | $\frac{2\sqrt{35}}{3\pi^{2}}\\!\beta_{2}\beta_{3}\\!+\\!\frac{2067\sqrt{7}}{704\pi^{2}}\\!\beta_{3}\beta_{4}$ | $\frac{215\sqrt{5}}{352\pi}\beta_{2}\beta_{4}+\frac{315}{176\pi}\beta_{3}^{2}+\frac{43305}{36608\pi}\beta_{4}^{2}$
$\varepsilon_{n}$ (body-body)
$=0.22\beta_{3}+$ |  $\\!=\\!0.95\beta_{2}\\!+\\!0.043\beta_{2}^{2}+$ |  $=1.08\beta_{3}+$ | $=1.23\beta_{4}+0.90\beta_{2}^{2}+$
| $1.7\beta_{2}\beta_{3}+1.8\beta_{3}\beta_{4}$ | $0.92\beta_{2}\beta_{4}+\\!0.31\beta_{4}^{2}\\!+\\!0.32\beta_{3}^{2}\\!$ | $0.40\beta_{2}\beta_{3}+0.79\beta_{3}\beta_{4}$ | $0.43\beta_{2}\beta_{4}+0.57\beta_{3}^{2}+0.38\beta_{4}^{2}$
| $\frac{320}{\sqrt{378}\pi^{2}}\beta_{2}\beta_{3}+\frac{800}{33\pi^{2}\sqrt{2}}\beta_{3}\beta_{4}$ | $\\!\sqrt{\frac{15}{4\pi}}\beta_{2}\\!+\\!\frac{15}{\sqrt{21}\pi}\beta_{2}\beta_{4}\\!-\\!\frac{15}{8\pi}\beta_{2}^{2}\\!$ | $\frac{64}{\pi\sqrt{70\pi}}\beta_{3}-$ | $\sqrt{\frac{35}{4\pi}}\beta_{4}+\frac{15}{4\pi}\beta_{2}^{2}-$
$1+\sum_{m=2}^{4}\beta_{m}Y_{m,m}$
| $=1.67\beta_{2}\beta_{3}+1.74\beta_{3}\beta_{4}$ | $\\!=\\!1.09\beta_{2}+1.0\beta_{2}\beta_{4}-0.6\beta_{2}^{2}\\!$ | $\frac{12\sqrt{6}}{\sqrt{7}\pi^{2}}\beta_{2}\beta_{3}-\frac{81}{16\sqrt{2}\pi^{2}}\beta_{3}\beta_{4}$ | $\frac{\sqrt{525}}{4\pi}\beta_{2}\beta_{4}-\frac{105}{32\pi}\beta_{4}^{2}$
$\varepsilon_{n}$ (tip-tip)
|  | $\\!=1.37\beta_{3}-\\!$ | $=1.67\beta_{4}+1.19\beta_{2}^{2}-$
|  |  | $1.13\beta_{2}\beta_{3}-0.36\beta_{3}\beta_{4}$ | $-1.82\beta_{2}\beta_{4}-1.04\beta_{4}^{2}$
$\left\langle\varepsilon_{n}^{2}\right\rangle$ | $\frac{4096}{3675\pi^{3}}\beta_{3}^{2}=0.036\beta_{3}^{2}$ | $\frac{3}{2\pi}\beta_{2}^{2}=0.477\beta_{2}^{2}$ | $\frac{4096}{245\pi^{3}}\beta_{3}^{2}=0.539\beta_{3}^{2}$ | $\frac{35}{18\pi}\beta_{4}^{2}+\frac{45}{14\pi^{2}}\beta_{2}^{4}$
(same rotation)
|  |  | $=0.62\beta_{4}^{2}+0.32\beta_{2}^{4}$
$\left\langle\varepsilon_{n}^{2}\right\rangle$ | $\frac{2048}{3675\pi^{3}}\beta_{3}^{2}=0.018\beta_{3}^{2}$ | $\frac{3}{4\pi}\beta_{2}^{2}=0.239\beta_{2}^{2}$ | $\frac{2048}{245\pi^{3}}\beta_{3}^{2}=0.270\beta_{3}^{2}$ | $\frac{35}{36\pi}\beta_{4}^{2}+\frac{45}{28\pi^{2}}\beta_{2}^{4}$
(indep. rotation)
|  |  | $=0.31\beta_{4}^{2}+0.16\beta_{2}^{4}$
Table 1:  The value of eccentricity generated by the deformation component
$Y_{n,0}$ (first row) and $Y_{n,n}$ (second row) for the special alignment of
two colliding nuclei similar to those shown in Fig. 1 which maximizes the
eccentricity values, obtained within an optical Glauber model with sharp
surface by setting $a_{0}=0$ in Eq. (1). Here only the leading and subleading
order contributions are included. The leading-order mean square values
$\left\langle\varepsilon_{n}^{2}\right\rangle$ obtained by averaging over
common random orientations for the two nuclei and independent random
orientations for the two nuclei are given in the third row and the last row,
respectively. The values in the latter case are a factor of two smaller, but
in both cases they are independent of $\gamma$, $\alpha_{3,m}$ or
$\alpha_{4,m}$.

A few remarks are in order. The maximum possible $\varepsilon_{n}$ values are
different between $Y_{n,n}$ and $Y_{n,0}$, and they are generally comparable
to the corresponding $\beta_{n}$ value. However, the ms values after averaging
over random orientations have exactly the same quadratic dependence on
$\beta_{n}$, showing no explicit dependence on the internal angular variables
$\gamma$ and $\alpha_{n,m}$ in Eq. 1 to the leading order 222The actual
probability density distribution $p(\varepsilon_{n})$ is different between
deformation described by $Y_{n,n}$ and by $Y_{n,0}$. This difference can be
captured by the fourth- and higher-order cumulants.. Remarkably, the octupole
deformation also gives rise to a dipolar eccentricity, following the same
quadratic dependence on $\beta_{3}$ but with a coefficient that is a factor of
ten smaller. Furthermore, the quadrupole deformation gives rise to a quartic
contribution to $\varepsilon_{4}$, and in an analogy to the well-known non-
linear contribution of $\varepsilon_{2}$ to $\varepsilon_{4}$ Teaney and Yan
(2011), scales as
$\left\langle\varepsilon_{4}^{2}\right\rangle\approx(0.7-1.4)\left\langle\varepsilon_{2}^{2}\right\rangle^{2}$.
Lastly, the coefficients $b_{n,m}^{\prime}$ listed in the table are derived
under a simplified scenario. In a more realistic Monte-Carlo Glauber model
calculation based on the Woods-Saxon nuclear profile, the coefficients in the
UCC region as shown in bottom row of Fig. 2, are comparable or slightly
smaller than those obtained by requiring zero impact parameter and independent
rotations.

For a more realistic estimation of influence of nuclear deformation, a Monte-
Carlo Glauber model Miller  _et al._ (2007) is used to simulate collisions of
238U and 96Zr systems and calculate $\varepsilon_{n}$ in each event. These
systems are chosen because the experimental collision data exist already. The
nucleons are assumed to have a hard-core of 0.4 fm in radii, with a density
described by Eq. (1). The nuclear radius $R_{0}$ and the surface thickness
$a_{0}$ are chosen to be $R_{0}=6.81$ fm and $a_{0}=0.55$ fm for 238U and
$R_{0}=5.09$ fm and $a_{0}=0.52$ fm for 96Zr, respectively. The nucleon-
nucleon inelastic cross-section are chosen to be $\sigma_{\mathrm{nn}}=42$ mb
at $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV. In each collision event,
nucleons are generated in each nuclei at a random impact parameter from each
other. Each nucleus is then rotated by randomly generated three Euler angles
before they are set on a straight line trajectory towards each other along the
$z$ direction. From this, the nucleons in the overlap region, known as
participants, are identified. The $\varepsilon_{n}$ are calculated from
nucleon participants according to Eq. (2), and the results are studied as a
function of $N_{\mathrm{part}}$. For a systematic study of the influence of
different shapes, one deformation component or particular combination of
components of the same $n$ is enabled at a time, the latter is useful to
understand the influence of departure from axial and/or reflection symmetry. A
special study is performed to also investigate the presence of shapes of
different $n$, where two or three non-zero values for $\beta_{2}$, $\beta_{3}$
and $\beta_{4}$ are enabled simultanously.

It is well known that the particle production in nucleus-nucleus collisions
only scales approximately with $N_{\mathrm{part}}$. A better scaling can be
achieved by considering the constituent quarks as effective degree-of-freedom
for particle production Adler  _et al._ (2014); Lacey  _et al._ (2018);
Loizides (2016); Bożek  _et al._ (2016); Acharya  _et al._ (2019), which would
naturally give rise to different $\varepsilon_{n}$ in each event. Defining
centrality with constituent quarks is also expected to change the fluctuations
of eccentricity Zhou and Jia (2018), and provide a way to quantify the
centrality smearing effects. For this purpose, a quark Glauber model from Ref.
Loizides (2016) is used. Three quark constituents are generated for each
nucleon according to the “mod” configuration Mitchell  _et al._ (2016), which
ensures that the radial distribution of the three constituents after re-
centering follows the proton form factor
$\rho_{\mathrm{proton}}(r)=e^{-r/r_{0}}$ with $r_{0}=0.234$ fm De Forest and
Walecka (1966). The value of quark-quark cross-section is chosen to be
$\sigma_{\mathrm{qq}}=8.2$ mb in order to match the $\sigma_{\mathrm{nn}}$.
The $\varepsilon_{n}$ are then calculated from the list of quark participants
in the overlap region, and the number of quark participants
$N_{\mathrm{quark}}$ is used as an alternative centrality estimator. In the
quark Glauber model, I also keep track explicitly the participant nucleons,
i.e. a nucleon is counted as participant as long as one of its quark
participate in the collision. This paper presents and compares results
obtained from both nucleon participants and quark participants.

In the presence of large deformation, the total volume of the nucleus
increases slightly for fixed $R_{0}$ Myers and Schmidt (1983). Considering the
quadrupole deformation only, for the largest value considered,
$\beta_{2}=0.34$, the ratio to the original volume is approximately
$1+\frac{3}{4\pi}\beta_{2}^{2}+\frac{\sqrt{5}}{28\pi^{3/2}}\cos(3\gamma)\beta_{2}^{3}=1.021+0.0004\cos(3\gamma)$.
In order to keep the overall volume fixed, it would require a small less than
1% decrease of the $R_{0}$, which is safely ignored in the present study.

##  III Results

The goal of this paper is to explore the relation between
$\left\langle\varepsilon_{n}^{2}\right\rangle$ and various deformation
parameters in Eq. (1), and to provide insights on the deformation dependence
of experimentally measured $\left\langle v_{n}^{2}\right\rangle$. The
influence of nuclear deformation on higher-order cumulants of
$\varepsilon_{n}$ will be explored in a separate study. Section III.1
establishes the quadratic relation Eq. (3) by considering the axial-symmetric
deformation $Y_{n,0}$, $n=2,3$ and 4. The influences of non-axial deformation,
$Y_{n,m\neq 0}$, characterized by the triaxiality parameter $\gamma$ and
$\alpha_{n,m}$ parameters are considered in Section III.2. One finds that the
slope parameters $b_{n,m}^{\prime}$ have a very weak dependence on the
$\gamma$ and $\alpha_{n,m}$. Section III.3 considers the presence of multiple
shape components $\beta_{2}$, $\beta_{3}$, and $\beta_{4}$, which is generally
expected in nuclear structure physics. For moderate deformation values, one
finds that the non-linear contributions, terms like $\beta_{n}\beta_{m},m\neq
n$ is subdominant and Eq. (3) still holds well. Section III.4 discusses ways
to constrain these deformations simultaneously using flow measurements.

###  III.1 Influence of axial-symmetric multipole deformation

The top row of Fig. 2 shows the $N_{\mathrm{part}}$ dependence of
$\left\langle\varepsilon_{n}^{2}\right\rangle$, $n=2$, 3 and 4 for various
$\beta_{n}$ values in U+U collisions, calculated from the participating
nucleons according to Eq. (2). One observes that the larger $\beta_{2}$ values
increase $\varepsilon_{2}$ mostly in the central region, while larger
$\beta_{3}$ values increase $\varepsilon_{3}$ over the full centrality range.
The reason is that the contribution associated with the average elliptic
geometry to $\varepsilon_{2}$ dominates over the deformation effects in the
mid-central and peripheral collisions. On the other hand, the
$\varepsilon_{3}$ without nuclear deformation arises solely from random
fluctuations of nucleon positions and has a much smaller value, therefore the
$\varepsilon_{3}$ is more sensitive to $\beta_{3}$. The increase of
$\varepsilon_{4}$ with $\beta_{4}$ is observed only in the central region. In
fact, the peripheral region shows a slight decrease of $\varepsilon_{4}$ with
$\beta_{4}$. Overall, the influence of deformation on $\varepsilon_{n}$ is
largest in the most central region for all harmonics.

Figure 2:  Top row: The $N_{\mathrm{part}}$ dependence of mean squre
eccentricity $\left\langle\varepsilon_{n}^{2}\right\rangle$ for several
$\beta_{n}$ values considering only the $Y_{n,0}$ component, with markers and
lines correspond $d_{\perp}$ obtained with nucleons and quarks, respectively.
Middle row: The $\beta_{n}^{2}$ dependence of
$\left\langle\varepsilon_{n}^{2}\right\rangle$ in several centrality ranges
based on $N_{\mathrm{part}}$, which can be nicely described by a linear
function
$\left\langle\varepsilon_{n}^{2}\right\rangle-\left\langle\varepsilon_{n}^{2}\right\rangle_{\beta_{n}=0}=b^{\prime}_{n}\beta_{n}^{2}$.
Bottom row: the centrality dependence of the slope parameter $b^{\prime}_{n}$
in U+U (black) and Zr+Zr (red) systems for $\varepsilon_{n}$ calculated based
on nucleons (markers) or quarks (lines). The results are shown for $n=2$ (left
column), $n=3$ (middle column) and $n=4$ (right column). The first three
points in the bottom panels correspond to 0–0.2%, 0.2–0.5% and 0.5–1%,
respectively.

In the same plots, I also show the
$\left\langle\varepsilon_{n}^{2}\right\rangle$ calculated from quark
participants as solid lines, with the same color as those calculated from
nucleon participants. Small differences are observed in the UCC region when
$\beta_{n}$ are small, or in the peripheral region for $n=2$, implying that
the influences of deformation are insensitive to nucleon substructures.

To quantify these dependencies, $\left\langle\varepsilon_{n}^{2}\right\rangle$
values obtained for fixed $N_{\mathrm{part}}$ are averaged in narrow
centrality ranges, and plotted as a function of $\beta_{n}^{2}$ in the middle
row of Fig. 2. An linear dependence is observed in all cases, confirming the
first part of the relation in Eq. (3) involving $a^{\prime}_{n}$ and
$b^{\prime}_{n}$. Note that $a^{\prime}_{n}$ correspond to eccentricities in
the absence of deformation
$a^{\prime}_{n}=\left\langle\varepsilon_{n}^{2}\right\rangle_{|\beta_{n}=0}$
shown by the black solid circles in the top row, while $b^{\prime}_{n}$
describes the slope of the $\beta_{n}^{2}$ dependencies in the middle row .

The bottom row of Fig. 2 shows the centrality dependence of $b^{\prime}_{n}$
for U+U and Zr+Zr collisions. The values of $b^{\prime}_{n}$ are largest in
the UCC region and decrease toward mid-central and peripheral region. It is
quite remarkable that the value of $b^{\prime}_{n}$ starts at around 0.2–0.3
for all harmonics in both collision systems. This value of $b^{\prime}_{n}$
reflects an effect that is purely geometrical. If the two nuclei were to
collide head-on in the direction perpendicular to maximum deformation as shown
in Fig. 1, $b^{\prime}_{n}$ should be on the order of unity, see Table 1. In
reality, after averaging over all possible random orientations, the
$b^{\prime}_{n}$ values are reduced to about 0.2–0.3. In a Monte-Carlo Glauber
with finite number of nucleons, the random fluctuation of nucleon positions
smear the correlation between the shape of the overlap region and the
$N_{\mathrm{part}}$. This smearing is expected to be larger for smaller
system, leading to a slightly smaller $b^{\prime}_{n}$ in the Zr+Zr collisions
than in the U+U collisions. Note that the centrality and collision system
dependencies of $b_{n}^{\prime}$ are just the opposite of $a^{\prime}_{n}$
(the latter corresponds to the $\left\langle\varepsilon_{n}^{2}\right\rangle$
without deformation shown in the top row). The values for $a^{\prime}_{n}$ are
smallest in the UCC region and increase towards more peripheral region and
exhibit a much larger difference between Zr+Zr and U+U.

The bottom row of Fig. 2 also compares the $b_{n}^{\prime}$ calculated from
nucleon participants with that calculated from quark participants. In the
central and mid-central collisions, the differences are negligible for $n=2$,
but for $n=3$ and $n=4$ the results based on quark participants are
systematically smaller. At this point, one may wonder if the $\varepsilon_{n}$
are also affected by the $\beta_{m}$ of different order, $m\neq n$. I have
performed such calculations. In most cases, the influences are small. But one
identifies three cases for which the influences are quite large. In
particular, one finds that the octupole deformation contributes strongly to
the dipolar eccentricity in all centrality, and the hexadecapole deformation
contributes to both dipolar and triangular eccentricities in the mid-central
collisions. They are presented in Fig. 3 with a similar layout as Fig. 2. As
these contributions are a global geometry effect with a radial distribution
different from contributions arising from random fluctuation in nucleon
positions in each event, they probably will be damped differently by viscous
effects in comparison to $a_{n}^{\prime}$. Similar effects are known in the
context of hydrodynamic model studies as leading and sub-leading
eccentricities, which characterize different length scale in radial direction
for $\varepsilon_{n}$, and sub-leading $\varepsilon_{n}$ with higher frequency
in radial direction is more damped than the leading $\varepsilon_{n}$
Mazeliauskas and Teaney (2015, 2016). Results for other $b_{n,m}^{\prime}$ can
be found in Figs. 10–13 in the Appendix B, including the quartic dependence of
$\left\langle\varepsilon_{4}^{2}\right\rangle$ on $\beta_{2}$ predicted in
Tab. 1.

Figure 3:  Characterization of the influence of $\beta_{3}$ to
$\left\langle\varepsilon_{1}^{2}\right\rangle$ (left column), $\beta_{4}$ to
$\left\langle\varepsilon_{1}^{2}\right\rangle$ (middle column) and $\beta_{4}$
to $\left\langle\varepsilon_{3}^{2}\right\rangle$ (right column). The top row
shows the centrality dependence of
$\left\langle\varepsilon_{n}^{2}\right\rangle$ for several values of
$\beta_{3}$ or $\beta_{4}$ as indicated in the panels, with markers and lines
correspond $d_{\perp}$ obtained with nucleons and quarks, respectively. The
middle row shows the $\beta_{3}^{2}$ or $\beta_{4}^{2}$ dependence of
$\left\langle\varepsilon_{1}^{2}\right\rangle$ or
$\left\langle\varepsilon_{1}^{2}\right\rangle$ in several centrality ranges.
The bottom row shows the centrality dependence of the extracted slope
parameter $b^{\prime}_{1,3}$ (left), $b^{\prime}_{1,4}$ (middle) or
$b^{\prime}_{3,4}$ (right) for U+U (black) and Zr+Zr (red) systems for
$\varepsilon_{n}$ calculated based on nucleons (markers) or quarks (lines).
The first three points in the bottom panels correspond to 0–0.2%, 0.2–0.5% and
0.5–1%, respectively.

The middle row of Fig. 3 shows that contributions between different orders
also follow a quadratic dependence, confirming the second part of the Eq. (3).
The slopes, $b^{\prime}_{1,3}$, $b^{\prime}_{1,4}$ and $b^{\prime}_{2,4}$, are
summarized in the bottom row of Fig. 3. One should not be tricked by the
apparent small value of $b_{1,3}^{\prime}$, though. Since the value of
$a^{\prime}_{1}$, the $\left\langle\varepsilon_{1}^{2}\right\rangle$ without
deformation, is very small in the UCC region, even a value of
$b^{\prime}_{1,3}=0.015$ together with a modest octupole deformation of
$\beta_{3}=0.1$ could increase the
$\left\langle\varepsilon_{1}^{2}\right\rangle$ by about 15%. In mid-central
collisions, due to a much larger $a^{\prime}_{1}$, the combined contributions
from $\beta_{3}=0.1$ and $\beta_{4}=0.1$ are less than 10%. This result
suggests that the dipolar flow in the UCC region could in principle be used to
probe the octupole deformation. On the other hand, the influence of
$\beta_{4}$ on $\varepsilon_{3}$ is significant in the mid-central collisions,
and is negligible in the UCC region. The bottom row of Fig. 3 also compares
the slope parameters between U+U and Zr+Zr collisions, they are very similar
in the UCC region, but values in Zr+Zr are about 20% smaller in the mid-
central collisions.

To summarize the main message of Fig. 3, the $\varepsilon_{n}$ for $n=2$,3,
and 4 in the UCC region is not affected by deformation of different order
$\beta_{m}$, $m\neq n$, leading to a particularly simple expression,
$\left\langle\epsilon_{n}^{2}\right\rangle_{\mathrm{UCC}}=a^{\prime}_{n}+b^{\prime}_{n}\beta_{n}^{2}$.
Exploiting this relation in the UCC collisions from experimentally measured
$v_{n}$ values provides a clean way to constrain the $\beta_{n}$ parameters as
will be discussed in Section III.4.

###  III.2 Influence of non-axial deformations

Let’s first consider the influence of triaxiality parameter $\gamma$, which
mixes the contribution from $Y_{20}$ and $Y_{22}$ components, while keeping
the overall magnitude of quadrupole deformation $\beta_{2}$ fixed. The top-
left panel of Fig. 4 shows the $N_{\mathrm{part}}$ dependence of
$\left\langle\varepsilon_{2}^{2}\right\rangle$ for $\beta_{2}=0.28$ but
different $\gamma$ values in the U+U collisions. They are contrasted to the
case for spherical nuclei $\beta_{2}=0$. It is clear that over most of the
centrality range, $\left\langle\varepsilon_{2}^{2}\right\rangle$ have very
little sensitivity to $\gamma$. The
$\left\langle\varepsilon_{2}^{2}\right\rangle$ calculated with quark
participants, shown as solid lines in the same panel, also give very similar
results.

Figure 4:  Characterization of the influence of triaxiality $\gamma$ on
$\left\langle\varepsilon_{2}^{2}\right\rangle$. Top row: The
$N_{\mathrm{part}}$ (left) and $N_{\mathrm{quark}}$ (right) dependencies of
$\left\langle\varepsilon_{2}^{2}\right\rangle$ for different values of
$\gamma$, where markers and lines correspond $d_{\perp}$ obtained with
nucleons and quarks, respectively. Middle row: The
$\left\langle\varepsilon_{2}^{2}\right\rangle-\left\langle\varepsilon_{2}^{2}(\gamma=\pi/6)\right\rangle$
as a function of $\cos(3\gamma)$ in several centrality ranges based on
$N_{\mathrm{part}}$ (left) and $N_{\mathrm{quark}}$ (right), which follows a
linear function of $\cos(3\gamma)$, whose slope is parametrized as
$c_{2}^{\prime}\beta_{2}^{2}$ (see Eq. (4)). Bottom row: the extracted
$c^{\prime}_{2}$ as a function of centrality for U+U (black) and Zr+Zr (red)
systems for $\varepsilon_{n}$ calculated based on nucleons (markers) or quarks
(lines). The first three points in the bottom panels correspond to 0–0.2%,
0.2–0.5% and 0.5–1%, respectively.

To test the influence of volume/centrality fluctuations, results obtained
using $N_{\mathrm{quark}}$ as centrality is shown in the top-right panel.
Large splittings between different $\gamma$ cases are observed in the UCC
region of 0–1%. Namely, the $\left\langle\varepsilon_{2}^{2}\right\rangle$ for
oblate deformation $\gamma=\pi/3$ shows a stronger increase as a function of
$N_{\mathrm{quark}}$ before they all start to decrease slightly at the largest
$N_{\mathrm{quark}}$ values. This behavior suggests that the events selection
based on $N_{\mathrm{part}}$ or $N_{\mathrm{quark}}$ have different
correlation with, and therefore different sensitivity to, the triaxiality of
the nucleus. The largest difference is reached between the prolate deformation
and the oblate deformation, consistent with a previous study based on the AMPT
model Giacalone  _et al._ (2021a); Jia  _et al._ (2021). In that study, a
similar dependence on $\gamma$ is observed for the final-state $v_{2}$.
Interestingly, the U+U $v_{2}$ data from the STAR Collaboration show a
decreasing behavior as a function of $N_{\mathrm{ch}}$ in the UCC region,
while the Au+Au $v_{2}$ data show a slight rising trend Adamczyk  _et al._
(2015), compatible with a prolate deformation of 238U nucleus and a oblate
deformation of 197Au. In summary, our results suggest that the two-particle
correlators $\left\langle v_{2}^{2}\right\rangle$ are sensitive to the
triaxiality only in the UCC region, and the level of sensitivity depends on
the choice of centrality estimator.

The middle row of Fig. 4 quantifies the $\gamma$ dependence of
$\left\langle\varepsilon_{2}^{2}\right\rangle$ in several centrality ranges in
$N_{\mathrm{part}}$ on the left and $N_{\mathrm{quark}}$ on the right in the
U+U collisions. The $\gamma$ dependence is well described by a linear function
of $\cos(3\gamma)$, reflecting the expected three-fold periodicity. Similar
observation is also made in the Zr+Zr collisions, although the sensitivity to
$\gamma$ is observed over a larger centrality range. Based on this finding, I
arrive the following empirical formula that accounts for the dependence on
both $\beta_{2}$ and $\gamma$,

$\displaystyle\left\langle\varepsilon_{2}^{2}\right\rangle=a^{\prime}_{2}+(b^{\prime}_{2}+c^{\prime}_{2}\cos(3\gamma))\beta_{2}^{2}\;.$
(4)

Note that the $c^{\prime}_{2}\beta_{2}^{2}$ is the slope of the
$\cos(3\gamma)$ dependence in the middle row of Fig. 4. The bottom row of Fig.
4 shows the centrality dependence of $c^{\prime}_{2}$. The value of
$c^{\prime}_{2}$ is generally much smaller than $b^{\prime}_{2}$,
$c^{\prime}_{2}\ll b^{\prime}_{2}$, and approaches zero in the peripheral
collisions. However, its value in the UCC region changes sign and could gain a
sizable magnitude depending on the centrality estimator. Lastly, for other
eccentricities $\varepsilon_{n},n\neq 2$, only very small dependencies on the
triaxiality are observed, typically less than 5%. However, these dependencies
to a good extent can also be described by a $\cos(3\gamma)$ function (see the
left column of Figs. 10-13).

This result begs the question of whether the finding about triaxiality also
applies for the octupole and hexadecapole deformations. These higher-order
deformations have many more shape parameters. In the intrinsic frame defined
by the quadrupole deformation, after taking out the $\beta_{n}$ that describe
the overall strength of the deformation, there are still six and eight
independent shape variables for octupole and hexadecapole deformations,
respectively. Besides, there are also large redundancies in the parameter
space due to spatial symmetry of spherical harmonics. For example, if the
underlying quadrupole deformation is axial-symmetric, terms like
$\cos(\delta)Y_{n,m}+\sin(\delta)Y_{n,-m}$ can be absorbed by an azimuthal
rotation of $Y_{n,m}$ without real physical consequence. There have been
several attempts to find efficient parameterizations to reduce this
redundancy, see Refs.Hamamoto  _et al._ (1991); Rohozinski (1997). Our paper
follows a more relaxed approach, where I just test special cases of the
octupole and hexadecapole shapes. It is reassuring that
$\left\langle\varepsilon_{n}^{2}\right\rangle$ has very small sensitivity to
these internal angular parameters, as will be described below.

For this study, I consider all real valued spherical harmonics $Y_{3,m}$ and
$Y_{4,m}$. They are introduced one at a time in the Glauber model and the
resulting $\varepsilon_{n}$ are calculated. The results are summarized in Fig.
5. I found that $\varepsilon_{3}$ values are the same for all $Y_{3,m}$
components, except for small differences in the UCC region. I also tried
several combinations, such as $\cos(\delta)Y_{3,0}+\sin(\delta)Y_{3,1}$ with
$\delta$ a free mixing angle, and the conclusion remains the same. I suspect
that this is true for general mixing of all components $\alpha_{3,m}$ in Eq.
(1), as long as $\sum_{m=-3}^{3}\alpha_{3m}^{2}=1$ is satisfied. This
independence is also observed for $\beta_{3}$ contribution to the dipolar
eccentricity $\varepsilon_{1}$ and probably is a property for all odd-order
deformations.

On the other hand, the $\left\langle\varepsilon_{4}^{2}\right\rangle$ values
show a modest, at a level of 15%, differences among different $Y_{4,m}$
components for $\beta_{4}=0.2$. Results for non-axial components
$1+\beta_{4}Y_{4,m},m\neq 0$ lie exactly between $1+\beta_{4}Y_{4,0}$ and
$1-\beta_{4}Y_{4,0}$. The differences are largest in central collisions but
are present throughout the entire centrality range. This is different from
$\left\langle\varepsilon_{2}^{2}\right\rangle$, for which the dependence on
the triaxiality is observed only in the UCC region. I initially thought that
the $\left\langle\varepsilon_{4}^{2}\right\rangle$ for the most general
hexadecapole shape should be in between the results for $Y_{40}$ and
$-Y_{40}$, which turns out is not the case. In fact, the extrema of
$\left\langle\varepsilon_{4}^{2}\right\rangle$ are reached for deformation
described by $\sqrt{7/12}Y_{4,0}+\sqrt{5/12}Y_{4,4}$ and
$\sqrt{5/12}Y_{4,0}-\sqrt{7/12}Y_{4,4}$, with the maximum for
$\beta_{4}=|\beta_{4}|$ and minimum for $\beta_{4}=-|\beta_{4}|$ Rohozinski
(1997). Identifying heavy ion observables that are sensitive to the sign of
$\beta_{4}$ will be particularly useful for understanding the nuclear fission
data Lemmon  _et al._ (1993). A more detailed investigation of this topic is
given in the Appendix B.

Figure 5:  Left: The $N_{\mathrm{part}}$ dependence of
$\left\langle\varepsilon_{3}^{2}\right\rangle$ for different component of
octupole deformation $Y_{3,m}$ with a deformation value of $\beta_{3}=0.2$.
Right: The $N_{\mathrm{part}}$ dependence of
$\left\langle\varepsilon_{4}^{2}\right\rangle$ for different component of
octupole deformation $Y_{4,m}$ with a deformation value of $\beta_{4}=0.2$.
The $\varepsilon_{n}$ are calculated either from nucleons (symbols) or quarks
(lines). For all cases $\beta_{2}=0$.

###  III.3 Simultaneous presence of quadrupole-octupole-hexadecapole
deformations

Although the axial quadrupole distortion is the nuclear deformation of primary
importance, secondary contributions from octupole and hexadecapole components
often coexist and can be important in some regions of nuclear chart Butler
(2016). One example is the pear-shaped 224Ra P  _et al._ (2013) with a
$(\beta_{2},\beta_{3},\beta_{4})$ value of $(0.1545,0.097,0.080)$ Nazarewicz
_et al._ (1984). A summary of the deformation parameters for the large systems
collided at RHIC and the LHC are listed in Tab. 2, highlighting the importance
of possible higher-order deformations. It would be interesting to study how
the eccentricities depend on the simultaneous presence of these different
deformations, in particular, whether the contribution from each component to
$\varepsilon_{n}$ is independent of each other.

| $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$
---|---|---|---
238U | 0.286 Raman  _et al._ (2001) | 0.078 Agbemava  _et al._ (2016) | 0.07–0.09 Libert  _et al._ (1982); Möller  _et al._ (2016)
208Pb | 0.05 Raman  _et al._ (2001) | 0.04Robledo and Bertsch (2011) | ?
197Au | -(0.13-0.16) Möller  _et al._ (2016); Hilaire and Girod (2007) | ? | -0.03 Möller  _et al._ (2016)
129Xe | 0.16 Möller  _et al._ (2016) | ? | ?
96Ru | 0.05-0.16 Möller  _et al._ (2016); Raman  _et al._ (2001) | ? | ?
96Zr | 0.08 Raman  _et al._ (2001) | 0.20-0.27 KIBÉDI and SPEAR (2002) | 0.06 Möller  _et al._ (2016)
Table 2:  Some estimates of the deformation values $\beta_{2},\beta_{3}$, and
$\beta_{4}$ for the large nuclei collided at RHIC and the LHC with references
given, mostly based on global analysis of the $B(En)$ transition data.

For this exploratory study, only combinations of axial-symmetric components
$Y_{n,0},n=2,3,4$ are considered, from which the $\varepsilon_{1}$,
$\varepsilon_{2}$, $\varepsilon_{3}$ and $\varepsilon_{4}$ are calculated. The
analysis is carried out for different combination of
$(\beta_{2},\beta_{3},\beta_{4})$ from the values $\beta_{2}=\pm 0.1,0$,
$\beta_{3}=0.1,0$ and $\beta_{4}=0.1,0$. The results for U+U and Zr+Zr, in
terms of ratios to the spherical nuclei, are shown in Fig. 6. The
contributions to eccentricities from different deformation components are
almost independent of each other, i.e. following Eq. (3). Modest deviations
are observed in a few cases, however. In particular, the $\varepsilon_{2}$ is
observed to increase with $\beta_{3}$ in non-central region, and the
difference of $\varepsilon_{n}$ between $\beta_{2}=0.1$ and $-0.1$ is also
larger when $\beta_{3}$ and/or $\beta_{4}$ are non-zero. Remarkably, such non-
linear effects are very small in the UCC region, where $\varepsilon_{n}$ is
only sensitive to $\beta_{n}$ except for $n=1$. For the odd harmonics
$\varepsilon_{1}$ and $\varepsilon_{3}$, both $\beta_{3}$ and $\beta_{4}$ can
have large contribution in non-central collisions.

Figure 6:  Relative change of $\left\langle\varepsilon_{1}^{2}\right\rangle$
(left column), $\left\langle\varepsilon_{2}^{2}\right\rangle$ (second column),
$\left\langle\varepsilon_{3}^{2}\right\rangle$ (third column)
$\left\langle\varepsilon_{4}^{2}\right\rangle$ (right most column) for U+U
(top row) and Zr+Zr (bottom row) collisions, relative undeformed case, for
different combinations of $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$ as
indicated in the left panels. Only the axial component of the deformation,
$Y_{2,0}$,$Y_{3,0}$, and $Y_{4,0}$ are considered.

Figure 7 considers a different scenario where the quadrupole component is much
larger than the octupole and hexadecapole. For this case, I increase the
$\beta_{2}$ to the value of 0.28. Most trends remain qualitatively the same as
Fig. 6. In particular, the $\varepsilon_{1}$ and $\varepsilon_{3}$ over most
the centrality range, as well as $\varepsilon_{4}$ in the UCC region, are
still dominated by the $\beta_{3}$ and $\beta_{4}$. The behaviors for
$(\beta_{2},\beta_{3},\beta_{4})=(0.28,0.1,0.1)$ for 238U, comparable to the
values obtained from nuclear structure calculations in Table 2, are
particularly interesting. A significant enhancement of
$\left\langle\varepsilon_{3}^{2}\right\rangle$ of about 40% is expected in the
central collisions relative to the case of no deformation. Since
$\left\langle\varepsilon_{3}^{2}\right\rangle\propto 1/A$ without deformation
in large system, the $\left\langle\varepsilon_{3}^{2}\right\rangle$ in the UCC
Au+Au collisions should be 238/197-1= 21% larger than those in the UCC U+U
collisions. Therefore in the presence of non-zero $\beta_{3}$ and $\beta_{4}$,
the ordering is expected to be flipped: the
$\left\langle\varepsilon_{3}^{2}\right\rangle_{\rm U}$ is expected to be 20%
larger than $\left\langle\varepsilon_{3}^{2}\right\rangle_{\rm Au}$, and
consequently $\left\langle v_{3}^{2}\right\rangle_{\rm U}$ is expected to be
larger than $\left\langle v_{3}^{2}\right\rangle_{\rm Au}$. The reverse
ordering of $v_{3}$ between U+U and Au+Au collisions, if observed, would be a
strong indication for the presence of octupole deformation in 238U nucleus.

Another useful example is the 96Zr+ 96Zr and 96Ru+ 96Ru isobar collisions
taken by the STAR Collaboration in 2018. The ratio $\left\langle
v_{n}^{2}\right\rangle_{\rm Zr}/\left\langle v_{n}^{2}\right\rangle_{\rm Ru}$
will directly constrain the relative ordering of $\beta_{n,\rm Zr}$ and
$\beta_{n,\rm Ru}$ Giacalone  _et al._ (2021b), especially in the UCC region,
where other effects associated with the radial distribution of nucleons, such
as neutron skin, are less important Hammelmann  _et al._ (2020); Xu  _et al._
(2021). Therefore observation of significant deviation of the ratio
$\left\langle v_{3}^{2}\right\rangle_{\rm Zr}/\left\langle
v_{3}^{2}\right\rangle_{\rm Ru}$ from unity with the characteristic centrality
dependence similar to those shown in bottom panels of Fig. 6 would be a strong
evidence for the presence of octupole correlations in these isobar systems.

In summary, the contributions of deformation to $\varepsilon_{n}$ arise mainly
from $\beta_{n}$ for $n=2$,3 and 4. In particular, there are small cross-
contributions and non-linear effects between $\beta_{2}$ and $\varepsilon_{3}$
and between $\beta_{3}$ and $\varepsilon_{2}$, especially in the UCC region.
This should be contrasted to the well-known anti-correlation between
$a_{2}^{\prime}=\left\langle\varepsilon_{2}^{2}\right\rangle_{|\beta_{n}=0}$
and
$a_{3}^{\prime}=\left\langle\varepsilon_{3}^{2}\right\rangle_{|\beta_{n}=0}$
in the absence of nuclear deformation Huo  _et al._ (2014), where
$\varepsilon_{3}$ arises from random position fluctuations of nucleons, the
latter have opposite effects for $\varepsilon_{3}$ and $\varepsilon_{2}$.
Therefore, one can constrain the value of $\beta_{2}$, $\beta_{3}$, and
$\beta_{4}$ by combining the information from $v_{1}$, $v_{2}$, $v_{3}$ and
$v_{4}$ in the ultra-central collisions as will be discussed next.

Figure 7:  Relative change of $\left\langle\varepsilon_{1}^{2}\right\rangle$
(left column), $\left\langle\varepsilon_{2}^{2}\right\rangle$ (second column),
$\left\langle\varepsilon_{3}^{2}\right\rangle$ (third column)
$\left\langle\varepsilon_{4}^{2}\right\rangle$ (right most column) for U+U
(top row) and Zr+Zr (bottom row) collisions, relative to undeformed case, for
different combinations of $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$ as
indicated in the left panels. Only the axial component of the deformation,
$Y_{2,0}$, $Y_{3,0}$, and $Y_{4,0}$ are considered. The only difference from
Fig. 6 is that a larger $\beta_{2}$ value of 0.28 is considered.

###  III.4 Constraining quadrupole, octupole, and hexadecapole deformations
using collective flow data

Exploiting the linear dependence of between $\varepsilon_{2}^{2}$ and
$\beta_{2}^{2}$, I have previously proposed a method Giacalone  _et al._
(2021a) to relate the quadrupole deformation between two collision systems of
similar sizes. This method can be straightforwardly generalized to octupole
and hexadecapole deformations.

Recall that in the UCC region, $v_{n}$ for $n=2-4$ depends linearly on
$\varepsilon_{n}$, and $\varepsilon_{n}$ is sensitive only to $\beta_{n}$
according to Eq. (3). Therefore, one expects that the $v_{n}$ in the UCC
region to also follow a similar dependence on $\beta_{n}$,

$\left\langle v_{n}^{2}\right\rangle=a_{n}+b_{n}\beta_{n}^{2},\;n=2,3,4,$ (5)

where averages are performed over events in a narrow centrality class and
$a_{n}=\left\langle v_{n}^{2}\right\rangle_{|\beta_{n}=0}$ . Following the
argument of Ref. Giacalone  _et al._ (2021a), I write down a simple equation
relating the deformation and flow in two collision systems X+X and Y+Y that
are close in mass number, with subscript X(Y) indicating a quantity evaluated
in X+X(Y+Y) collisions:

$\displaystyle\beta^{2}_{n,\mathrm{Y}}=\left(\frac{r_{v_{n}^{2}}r_{a,n}-1}{r_{n,\mathrm{Y}}}\right)+\left(r_{v_{n}^{2}}r_{b,n}\right)\beta^{2}_{n,\mathrm{X}}\;,\;\;r_{v_{n}^{2}}\equiv\frac{\left\langle
v_{n}^{2}\right\rangle_{\rm Y}}{\left\langle v_{n}^{2}\right\rangle_{\rm
X}}\;,\;r_{b,n}=\frac{b_{n,\rm X}}{b_{n,\rm Y}}\;,\;r_{a,n}=\frac{a_{n,\rm
X}}{a_{n,\rm Y}},\;r_{n,\rm Y}=\frac{b_{n,\rm Y}}{a_{n,\rm Y}}.$ (6)

As shown in the bottom panels of Fig. 2, $b_{n}^{\prime}$ has very weak
dependence on system size, therefore one expects it is true also for $b_{n}$
and therefore $r_{b,n}\approx 1$. In the absence of deformation $\beta_{n,\rm
X}=\beta_{n,\rm Y}=0$, using the linear response relation $\left\langle
v_{n}^{2}\right\rangle=k_{n}^{2}\left\langle\varepsilon_{n}^{2}\right\rangle$,
the relative difference of harmonic flow between X+X and Y+Y collisions,
$\Delta\left\langle v_{n}^{2}\right\rangle/\left\langle
v_{n}^{2}\right\rangle=\left(\left\langle v_{n}^{2}\right\rangle_{\rm
X}-\left\langle v_{n}^{2}\right\rangle_{\rm Y}\right)/\left\langle
v_{2}^{2}\right\rangle_{\rm Y}$, can be decomposed as

$\displaystyle\frac{\Delta\left\langle
v_{n}^{2}\right\rangle_{|\beta_{n}=0}}{\left\langle
v_{n}^{2}\right\rangle_{|\beta_{n}=0}}=\frac{\Delta
k_{n}^{2}}{k_{n}^{2}}+\frac{\Delta\left\langle\varepsilon_{n}^{2}\right\rangle_{|\beta_{n}=0}}{\left\langle\varepsilon_{n}^{2}\right\rangle_{|\beta_{n}=0}}\rightarrow\frac{\Delta
a_{n}}{a_{n}}=\frac{\Delta k_{n}^{2}}{k_{n}^{2}}+\frac{\Delta
a^{\prime}_{n}}{a^{\prime}_{n}}\;.$ (7)

In the UCC region, the eccentricities are dominated by the random fluctuations
of nucleon positions and to a good extent can be approximated by
$\left\langle\varepsilon_{n}^{2}\right\rangle\propto 1/A$ Alver  _et al._
(2010); Bhalerao  _et al._ (2011), and therefore
$\Delta\left\langle\varepsilon_{n}^{2}\right\rangle/\left\langle\varepsilon_{n}^{2}\right\rangle\approx\Delta\frac{1}{A}/\frac{1}{A}$.
The response coefficient $k_{n}$ is damped with the respect to the ideal
hydrodynamic value, $k_{n,{\rm ih}}$. In the simplified acoustic scaling
scenario of Ref. Gubser and Yarom (2011); Staig and Shuryak (2011); Teaney and
Yan (2012); Lacey  _et al._ (2016), one has $k_{n}/k_{n,\mathrm{ih}}\approx
1-Kn^{2}$, where $K$ encodes the viscous correction. This leads to $\Delta
k_{n}/k_{n}\approx-\Delta Kn^{2}k_{n,\mathrm{ih}}/k_{n}$. For central
collisions of large systems Alver  _et al._ (2010), $k_{n}$ is close to the
ideal hydro limit, and $k_{n,\mathrm{ih}}/k_{n}$ are nearly independent of
$n$, therefore one obtains

$\frac{\Delta k_{n}^{2}}{k_{n}^{2}}=\frac{n^{2}}{m^{2}}\frac{\Delta
k_{m}^{2}}{k_{m}^{2}}\;.$ (8)

Combining Eqs. (7) and (8) yield two coupled equations,

$r_{a,2}^{\prime}-r_{a,2}=x_{3}(r_{a,3}^{\prime}-r_{a,3})=x_{4}(r_{a,4}^{\prime}-r_{a,4})\;,\;\;x_{3}\approx\frac{4}{9},x_{4}\approx\frac{4}{16}\;.$
(9)

These equations involve only ratios of quantities between two systems close in
size. All these ratios are close to unity and can be reliably estimated from
the hydrodynamic model. I have verified these relations explicitly in the AMPT
model simulation of Au+Au and U+U collisions without deformations in a
previous study Giacalone  _et al._ (2021a); Jia  _et al._ (2021). The
centrality dependence of these ratios are shown in Fig. 8 with ${\rm
X}=^{197}$Au, ${\rm Y}=^{238}$U. It is immediately clear that $r_{a,n}$
follows closely the centrality dependence trends of $r^{\prime}_{a,n}$, but
has smaller values due to viscous damping. The difference grows with $n$,
reflecting the stronger viscous damping for higher-order flow harmonics. In
the 0–1% most central collisions (the rightmost point for each dataset), one
has $r^{\prime}_{a,n}=1.23$ independent of the harmonic number. This number is
very close to the expected ratio of atomic numbers
$(1/A_{\rm{Au}})/(1/A_{\rm{U}})=238/197=1.21$.

Figure 8:  The ratio of mean square eccentricity and mean square flow
harmonics without deformation
$a_{n}^{\prime}=\left\langle\varepsilon_{n}^{2}\right\rangle_{\beta_{n}=0}$
and $a_{n}=\left\langle v_{n}^{2}\right\rangle_{\beta_{n}=0}$ between Au+Au
and U+U collisions for $n=2$ (top-left), $n=3$ (top-right), $n=4$ (bottom
left) as a function of centrality based on $N_{\mathrm{ch}}$ obtained from
AMPT model in Ref. Giacalone  _et al._ (2021a). The arrows indicate the
expected ratios of ms eccentricity in the large system limit:
$(1/A_{Au})/(1/A_{U})=238/197=1.21$. The bottom right panel shows the
difference between $r_{a,n}^{\prime}=a_{n,Au}^{\prime}/a_{n,U}^{\prime}$ and
$r_{a,n}=a_{n,Au}/a_{n,U}$, scaled by the expected viscous damping factor
according to Eq. (9).

The bottom right panel demonstrates the robustness of Eq. (9). Using the
viscous damping relation Eq.(8), the predicted values in the 0–1% most central
collisions converge remarkably within 0.005 between different harmonics. But
deviations from this scaling appear away from the most central region. For
$v_{3}$, one finds a value of $x_{3}=5/9\sim 0.55$ achieves best agreements
over the 1%–25% centrality range as shown. This value is close to
$x_{3}=0.57\sim 5/9$ from a recent state-of-the-art hydrodynamic simulation
Gardim and Ollitrault (2021). For the $v_{4}$, I found $x_{4}=5/16\sim 0.31$
has the best agreement in the 1%–25% centrality range, although the
interpretation may be complicated by the mode-mixing contribution from
elliptic flow that scales like $v_{4}\sim v_{2}^{2}$, for which a smaller
damping $x_{4}=8/16$ is expected Teaney and Yan (2012). For the remaining
discussion, I shall focus simply on the 0–1% most central bin.

First, I use the approximation
$\left\langle\varepsilon_{n}^{2}\right\rangle_{|\beta_{n}=0}\propto 1/A$ and
rewrite Eq. (9) and Eq. (6) as:

$\displaystyle\beta^{2}_{n,\mathrm{Y}}=\frac{r_{v_{n}^{2}}r_{a,n}-1}{r_{n,\mathrm{Y}}}+r_{v_{n}^{2}}\beta^{2}_{n,\mathrm{X}}$
$\displaystyle
r_{a,2}=(1-x_{3})\frac{A_{Y}}{A_{X}}+x_{3}r_{a,3}=(1-x_{4})\frac{A_{Y}}{A_{X}}+x_{4}r_{a,4},\;\;x_{3}\approx\frac{4}{9},x_{4}\approx\frac{4}{16}\;.$
(10)

Presumably, if one deformation e.g $n=3$ is absent, the $r_{a,3}=\left\langle
v_{3,X}^{2}\right\rangle_{|\beta_{3}=0}/\left\langle
v_{3,Y}^{2}\right\rangle_{|\beta_{3}=0}$ can be obtained directly from
experiments, which allow us to fix $r_{a,2}$ and $r_{a,4}$ values.
Alternatively, $r_{a,n}$ can be cross-calibrated by picking nuclei with
similar mass number, therefore all of them are very close to unity. One such
example is the Zr+Zr and Ru+Ru isobar datasets for which both $r_{b,n}$ and
$r_{a,n}$ should be very close to unity 333A few percent difference in
$\varepsilon_{2}$ might arise because the difference in neutron skin effects
between 96Zr and 96Ru Hammelmann  _et al._ (2020); Xu  _et al._ (2021), but
these effects are much smaller than the influence of nuclear deformation in
the UCC region.. The only variable that needs to be evaluated numerically in
hydrodynamic model is $r_{n,Y}=b_{n,Y}/a_{n,Y}$, which is the property of a
single collision system.

One such numerical analysis has been performed in Ref. Giacalone  _et al._
(2021a), here I offer a bit more discussion on the expected behavior. Defining
two response coefficients, $k_{b,n}=\sqrt{b_{n}/b_{n}^{\prime}}$ and
$k_{a,n}=\sqrt{a_{n}/a_{n}^{\prime}}$, $r_{n}$ can be rewritten as

$\displaystyle
r_{n}=\frac{b_{n}}{a_{n}}=\frac{k_{b,n}^{2}}{k_{a,n}^{2}}\frac{b^{\prime}_{n}}{a^{\prime}_{n}}$
(11)

The $k_{a,n}^{2}$ is the usual viscous damping coefficient for $\left\langle
v_{n}^{2}\right\rangle$ in the absence of deformation, while $k_{b,n}^{2}$
describes the damping of the $\beta_{n}^{2}$-dependent part of the
$\left\langle v_{n}^{2}\right\rangle$ in Eq. (5). From the AMPT model study, I
have found that the ratio of the two damping coefficients for $n=2$,
$k_{b,2}/k_{a,2}\approx 0.75$ in the UCC region in U+U collisions, and only
has a very weak dependence on $N_{\mathrm{part}}$, suggesting that
$k_{b,2}/k_{a,2}$ is not very sensitive to viscosity. If this is the case, the
model dependence lies in the ratio
$r_{n}^{\prime}=b^{\prime}_{n}/a^{\prime}_{n}$, whose uncertainty arises
mainly from centrality smearing effects, e.g. the relative smearing of
centrality based on the final state charged particle multiplicity
$N_{\mathrm{ch}}$ and the $N_{\mathrm{part}}$.

##  IV Discussion and summary

The main finding of the paper is the simple parametric relation between
$\varepsilon_{n}$ and multi-pole deformation of nuclei $\beta_{n}$,
$\left\langle\varepsilon_{n}^{2}\right\rangle=a^{\prime}_{n}+b^{\prime}_{n}\beta_{n}^{2}$
for $n=2$,3 and 4, valid in all centrality and different collision systems.
The $a_{n}^{\prime}$ reflects the eccentricities for spherical nuclei, i.e.
$a_{n}^{\prime}=\left\langle\varepsilon_{2}^{2}\right\rangle_{|\beta_{n}=0}$.
The $a_{2}^{\prime}$ is dominated by elliptic shape of the overlap region,
which starts at a small value in central collisions and grows rapidly toward
mid-central and peripheral collisions. Other $a_{n}^{\prime},n\neq 2$ are
generated by random fluctuations of participating nucleons and typically
scales as $1/N_{\mathrm{part}}$. On the other hand, $\beta_{n}$ influences the
global shape of the overlap region on an event-by-event bases, and its
contribution to eccentricity $b^{\prime}_{n}\beta_{n}^{2}$ plays a similar
role as the so called reaction plane ellipticity $\varepsilon_{2,\mathrm{RP}}$
associated with the average elliptic shape of the overlap. Due to linear
response $v_{n}=k_{n}\varepsilon_{n}$ predicted by hydrodynamic models, I
expect a similar dependence for $v_{n}$, $\left\langle
v_{n}^{2}\right\rangle=a_{n}+b_{n}\beta_{n}^{2}$. From these, I define a
deformation-dependent and deformation-independent hydrodynamic response
coefficients $k_{b,n}=\sqrt{b_{n}/b_{n}^{\prime}}$ and
$k_{a,n}=\sqrt{a_{n}/a_{n}^{\prime}}$. It would be insightful to investigate
and compare $k_{a,n}$ and $k_{b,n}$, which will provide new kind of test on
the hydrodynamic models.

The best place to reveal nuclear deformation is the ultra-central collisions
(UCC) of large systems, where the deformation-driven components become
comparable or even larger than the values without deformation. For this
purpose, I propose a collision-system scan of a few species of similar size at
RHIC to systematically establish the influence of deformation, see the sketch
in Fig. 9.

First, it would be useful to scan two nuclei, e.g. 208Pb and another species,
in the vicinity of 197Au, to improve the modeling of Au+Au collisions, an
information which is crucial for the precision interpretation of high-
statistics flow data. Comparison between Pb+Pb at RHIC and the LHC will
constrain any possible energy dependence of the initial state effects and pre-
equilibrium dynamics. Since 208Pb is nearly spherical, a comparison of Pb+Pb
with Au+Au collisions at the same energy will also allow us to better
understand the impact of the moderate deformation of 197Au in Au+Au
collisions. The collisions of another species e.g 198Hg+198Hg
($\beta_{2}=-0.11$) would then probe more deeply the nature of the deformation
of 197Au, which, being an odd-mass nucleus, hasn’t been directly measured in
low-energy experiments. Having additional systems also provides an independent
cross-check on the initial state, for example one can setup three relations
like Eq.(10) to “triangulate” the consistency of the three deformation values.

In the second step, one can use flow measurements in conjunction with
hydrodynamics to map out the evolution of the quadrupole deformation along the
chain of stable samarium isotopes. As proposed in Ref. Giacalone  _et al._
(2021a), it would be useful to collide three isotopes: 144Sm
($\beta_{2}=0.08$, as spherical as 208Pb), 148Sm ($\beta_{2}=0.18$, triaxial
much as 129Xe and 197Au), and 154Sm ($\beta_{2}=0.32$ well-deformed like
238U). The evolution of the quadrupole deformation can be mapped precisely at
RHIC, thus offering a valuable test of nuclear structure knowledge. This scan
also enables a search for enhanced octupole correlations, i.e., $\beta_{3}$
values, which are predicted to be present in the region $Z\sim 56$/$N\sim 88$
Butler (2016) including the samarium isotopes. The influence of octupole
correlations would manifest in high-energy collisions as enhanced $v_{3}$, as
well as modified $\rho(v_{3}^{2},[p_{\mathrm{T}}])$ correlator. Evidence of
static octupole moments at low energies is rather sparse, and heavy ion
collisions might be a more sensitive approach.

Figure 9:  Correlation of mean square elliptic eccentricity and $\beta_{2}$ in
0–1% centrality,
$\left\langle\varepsilon_{2}^{2}\right\rangle=a^{\prime}_{2}+b^{\prime}_{2}\beta_{n}$,
estimated from Glauber model for various large collision species, where is
$a^{\prime}_{2}$ and $b^{\prime}_{2}$ are approximately $a^{\prime}_{2}=1.3/A$
and $b_{2}^{\prime}=0.23$. For even-even nuclei, the values of deformation
come from low-energy data, as in Table 2. For 197Au, whose $\beta_{2}$ is not
directly measured, I also include recent estimate based on the $v_{2}$ data
from Au+Au and U+U collisions, $|\beta_{2,\rm Au}|\approx 0.18$ Giacalone  _et
al._ (2021a). Offsets between different species are due to differences in mass
numbers $A$, therefore each group of nuclei with very similar masses, (197Au,
198Hg,208Pb) or (144Sm, 148Sm,154Sm), almost fall on the same curve.

In summary, I have studied the parametric dependence of eccentricity
$\varepsilon_{n}$ on the quadrupole $\beta_{2}$, octupole $\beta_{3}$ and
hexadecapole $\beta_{4}$ deformation of the nucleus in heavy ion collisions.
The mean square eccentricity $\left\langle\varepsilon_{n}^{2}\right\rangle$
are found to depend linearly on $\beta_{n}^{2}$ for $n$=2, 3 and 4. I also
find that $\beta_{3}$ contributes significantly to $\varepsilon_{1}$, and in
non-central collisions, the $\beta_{4}$ also contributes to $\varepsilon_{1}$
and $\varepsilon_{3}$. In central collisions, there are very little cross
contribution between $\beta_{2}$ and $\varepsilon_{3}$ and between $\beta_{3}$
and $\varepsilon_{2}$, although in non-central collisions, $\beta_{3}$
($\beta_{2}$) contributes modestly to $\varepsilon_{2}$ ($\varepsilon_{3}$).
Since harmonic flow $v_{n}$ is directly driven by the corresponding
eccentricity, $v_{n}\propto\varepsilon_{n}$, one expects very similar
parametric dependencies between $v_{n}$ and $\beta_{m}$ for both $m=n$ and
$m\neq n$. These findings provide a strong motivation to use heavy ion
collisions as a precision tool to scan and map out the ground state nuclear
deformations and compare with low energy nuclear measurements, not only for
the quadrupole deformation, but also for the octupole deformation, whose
evidence is quite sparse in nuclear structure physics. The procedure for doing
this is presented. One finds that the non-axial deformations, in particular
the triaxiality of quadrupole deformation, do not influence
$\left\langle\varepsilon_{n}^{2}\right\rangle$, but they can be probed by
other observables such as $v_{n}-[p_{\mathrm{T}}]$ correlation, and possibly
higher-order cumulants of $v_{n}$ and $[p_{\mathrm{T}}]$. The information
about the shape of atomic nuclei obtained in heavy-ion collisions is fully
complementary to that obtained in nuclear structure experiments. A carefully
planned system scan of stable species in the nuclear chart at RHIC, the LHC,
and other collider facilities could open new direction of research in nuclear
physics.

Acknowledgements: I thank Giuliano Giacalone for stimulating discussions and
recent collaborations that motivated this work. I thank Shengli Huang and
Chunjian Zhang for valuable discussions. This work is supported by DOE
DEFG0287ER40331 and NSF PHY-1913138.

##  Appendix A An analytical estimate in head-on collisions

I consider liquid-drop nuclear potential with a sharp surface,
$\rho(r,\theta,\phi)=\rho_{0}$ when $r<R(\theta,\phi)$ and zero otherwise, and
focus on head-on collisions with maximum overlap, i.e the two nuclei not only
have zero impact parameter, but also need to be aligned in a way to ensure the
overlap region contains all the nucleons $N_{\mathrm{part}}=2A$. However, at
the end of this section, I also relax the requirement to consider only zero
impact parameter, which corresponds more closely to the realistic scenario
when the effects of centrality resolution is taken into account. Our goal is
to establish the simple parametric dependence between
$\left\langle\varepsilon_{n}^{2}\right\rangle$ and $\beta_{m}$ in Eq. (3).

For this discussion, I switch back to using complex spherical harmonics
$Y_{n}^{m}$ with the normalization condition on the coefficients
$\sum_{m=-n}^{n}|\alpha_{n,m}|^{2}=1$ and $\alpha_{n,m}=\alpha_{n,-m}^{*}$.
Let’s first re-express Eq. (2) as an integral in the 3D, using the relation
$r_{\perp}=r\sin\theta$:

$\displaystyle\varepsilon_{n}e^{in\Phi_{n}}=-\frac{\int r^{n}\sin^{n}\\!\theta
e^{in\phi}\rho(\vec{r})d^{3}\vec{r}}{\int
r^{n}\sin^{n}\\!\theta\rho(\vec{r})d^{3}\vec{r}}=-\sqrt{\frac{4\pi(2n)!!}{(2n+1)!!}}\frac{\int(1+\sum_{l,m}\beta_{l}\alpha_{l,m}Y_{l}^{m})^{n+3}Y_{n}^{n}\sin\\!\theta
d\theta
d\phi}{\int(1+\sum_{l,m}\beta_{l}\alpha_{l,m}Y_{l}^{m})^{n+3}\sin^{n+1}\\!\theta
d\theta d\phi},$ (12)

where I have used
$Y_{n}^{n}=\sqrt{\frac{(2n+1)!!}{4\pi(2n)!!}}\sin^{n}\\!\theta e^{in\phi}$ and
the fact that the range of integration along radial direction is
$r\in[0,R_{0}(1+\sum_{l,m}\beta_{l}\alpha_{l,m}Y_{l}^{m})]$. Keeping the
integration to leading order in $\beta_{n}$, one has

$\displaystyle\varepsilon_{n}e^{in\Phi_{n}}=-A_{n}\int(\sum_{l,m}\beta_{l}\alpha_{l,m}Y_{l}^{m}Y_{n}^{n})\sin\\!\theta
d\theta
d\phi=-A_{n}\beta_{n}\alpha_{n,n},\;A_{n}\equiv\frac{(n+3)\Gamma(1+1/2+n/2)}{\pi\Gamma(1+n/2)}\sqrt{\frac{(2n)!!}{(2n+1)!!}}$
(13)

This result is easy to understand, for tip-tip collision where the $z$-axis is
aligned with beam direction (middle row of Fig. 1), only the $Y_{n}^{n}$
component can contribute to the eccentricity. The results,
$\varepsilon_{n}=A_{n}\beta_{n}/\sqrt{2}$, for nuclear surface containing only
this mode
$R=R_{0}(1+\beta_{n}Y_{n,n})=R_{0}(1+\beta_{n}/\sqrt{2}(Y_{n}^{-n}+(-1)^{n}Y_{n}^{n}))$,
are listed as the first term of each entry in the second row of Table 1.

In order to calculate the body-body collision for $Y_{n,0}=Y_{n}^{0}$ shown in
top row of Fig. 1, the nuclear surface, or equivalently the direction of
projection $Y_{n}^{n}$, need to be rotated by Euler angles
$(\alpha_{e},\beta_{e},\gamma_{e})=(0,\pi/2,0)$, i.e.
$Y_{n}^{n}=\sum_{m^{\prime}}D_{n,m^{\prime}}^{n}(0,\pi/2,0)Y_{n}^{m^{\prime}}$,
where $D_{n,m^{\prime}}^{n}$ is the Wigner rotational matrix. Plugging this
into Eq. (13) and considering only axial deformation
$R=R_{0}(1+\beta_{n}Y_{n}^{0})$, give,
$\varepsilon_{n}=A_{n}D_{n,0}^{n}(0,\pi/2,0)\beta_{n}=\sqrt{(2n)!}/(n!2^{n})A_{n}\beta_{n}$.
The values are provided as the first term of each entry in the top row of
Table 1.

The calculation of Eq. (12) including higher-order terms in $\beta_{n}$ is
straightforward; and the full expression up to the second order expansion are
listed in the table. Interestingly, It is found that if one only considers the
expansion for nuclear surface described by axial deformations
$R=R_{0}(1+\beta_{2}Y_{2,0}+\beta_{3}Y_{3,0}+\beta_{4}Y_{4,0})$ in the
numerator of an equation similar to Eq. (12), the ratios of the coefficients
of the high-order terms to that of the leading order in the first row of Table
1 are exactly the same as the Eq. (13) in Ref. Ryssens  _et al._ (2019).
However, including also the influence of the denominator of Eq. (12) modifies
the coefficients of some of the second-order terms.

Next, I consider random orientation of the nucleus around its center of mass
described by Euler angles $\Omega=(\alpha_{e},\beta_{e},\gamma_{e})$, and the
rotations of the two nuclei are required to be the same. For the purpose to
calculating the $\varepsilon_{n}$, it is equivalent to rotate the transverse
plane by applying the substitution,
$Y_{n}^{n}\rightarrow\sum_{m^{\prime}}D_{n,m^{\prime}}^{n}(\alpha,\beta,\gamma)Y_{n}^{m^{\prime}}$,
in Eq. (13),

$\displaystyle\varepsilon_{n}e^{in\Phi_{n}}=-A_{n}\int(\sum_{l,m,m^{\prime}}\beta_{l}\alpha_{l,m}D_{n,m^{\prime}}^{n}Y_{l}^{m}Y_{n}^{m^{\prime}})\sin\\!\theta
d\theta
d\phi=-A_{n}(\sum_{l,m,m^{\prime}}\beta_{l}\alpha_{l,m}D_{n,m^{\prime}}^{n}\delta_{ln}\delta_{m,-m^{\prime}})=-A_{n}\beta_{n}\sum_{m}\alpha_{n,m}D_{n,m}^{n}$
(14)

Using the orthogonality relation for D-matrix, the mean square average over
Euler angle is

$\displaystyle\left\langle\varepsilon_{n}^{2}\right\rangle=A_{n}^{2}\beta_{n}^{2}\int(\sum_{m}\alpha_{n,m}D_{n,m}^{n})(\sum_{m^{\prime}}\alpha_{n,m^{\prime}}D_{n,m^{\prime}}^{n})^{*}\frac{d\Omega}{8\pi^{2}}=\frac{A_{n}^{2}}{2n+1}\beta_{n}^{2}.$
(15)

The result apparently is independent of the mixture of different shape
component $Y_{n}^{m}$, as long as the overall magnitude $\beta_{n}$ remains
the same. This numerical values are listed in the third row of Table 1.

The dipolar eccentricity can be calculated in a similar way. The angular
weight involved is decomposed into two spherical harmonics, $\sin^{3}\\!\theta
e^{i\phi}=8/5\sqrt{\pi/21}(Y_{3}^{1}-\sqrt{14}Y_{1}^{1})$, which implies that
a octupole nuclear shape can give rise to a dipole eccentricity. As the
original nuclear surface has no dipole component, $Y_{1}^{1}$ drops out, and
Ione obtains,

$\displaystyle\varepsilon_{1}e^{i\Phi_{1}}=-\frac{\int r^{3}\sin^{3}\\!\theta
e^{i\phi}\rho(\vec{r})d^{3}\vec{r}}{\int
r^{3}\sin^{3}\\!\theta\rho(\vec{r})d^{3}\vec{r}}\approx-\frac{64}{5\sqrt{21\pi^{3}}}\beta_{3}\int\sum_{m}\alpha_{3,m}Y_{3}^{m}(\sum_{m^{\prime}}D_{1,m^{\prime}}^{3}Y_{3}^{m^{\prime}})=-\frac{64}{5\sqrt{21\pi^{3}}}\beta_{3}\sum_{m}\alpha_{3,m}D_{1,m}^{3}.$
(16)

Therefore the ms average of dipolar eccentricity over all orientations gives

$\displaystyle\left\langle\varepsilon_{1}^{2}\right\rangle=\frac{64^{2}}{525\pi^{3}}\beta_{3}^{2}\frac{1}{8\pi^{2}}\int\sum_{m,m^{\prime}}\alpha_{3,m}\alpha_{3,m^{\prime}}^{*}D_{1,m}^{3}D_{1,m}^{3*}d\Omega=\frac{4096}{3625\pi^{3}}\beta_{3}^{2},$
(17)

again independent of the mixture of different shape component $Y_{3}^{m}$.
Note that, this remarkable contribution is present entirely because the
$r_{\perp}^{3}$ weight in the definition of $\varepsilon_{1}$ in Eq. (2),
which is naturally required by the cumulant framework Teaney and Yan (2011).

Next, let us consider the possible contribution of quadrupole deformation to
the $\varepsilon_{4}$. For this purpose, one expands the nuclear shape in the
numerator of Eq. (12) and keeping terms that are proportional to
$\beta_{2}^{2}$,

$\displaystyle\varepsilon_{4}$ $\displaystyle\approx
3A_{4}\beta_{2}^{2}\int\sum_{m,m^{\prime}}\alpha_{2,m}\alpha_{2,m^{\prime}}D_{4,-m-m^{\prime}}^{4}Y_{2}^{m}Y_{2}^{m^{\prime}}Y_{4}^{-m-m^{\prime}}\sin\\!\theta
d\theta d\phi$
$\displaystyle=3A_{4}\beta_{2}^{2}\frac{15}{\sqrt{4\pi}}\sum_{m,m^{\prime}}\alpha_{2,m}\alpha_{2,m^{\prime}}D_{4,-m-m^{\prime}}^{4}\left(\begin{array}[]{ccc}2&2&4\\\
0&0&0\end{array}\right)\left(\begin{array}[]{ccc}2&2&4\\\
m&m^{\prime}&-m-m^{\prime}\end{array}\right)$ (22)
$\displaystyle=\frac{45}{2\pi}\sqrt{\frac{2}{35}}\beta_{2}^{2}[\alpha_{4,0}^{\prime}D_{4,0}^{4}+\frac{\alpha_{4,2}^{\prime}}{\sqrt{2}}(D_{4,2}^{4}+D_{4,-2}^{4})+\frac{\alpha_{4,4}^{\prime}}{\sqrt{2}}(D_{4,4}^{4}+D_{4,-4}^{4})]$
(23)

where $\alpha_{4,0}^{\prime}=(7+5\cos(2\gamma))/12$,
$\alpha_{4,2}^{\prime}=\sqrt{5/12}\sin(2\gamma)$ and
$\alpha_{4,4}^{\prime}=(1-\cos(2\gamma))\sqrt{35}/12$, satisfying
$(\alpha_{4,0}^{\prime})^{2}+(\alpha_{4,2}^{\prime})^{2}+(\alpha_{4,4}^{\prime})^{2}=1$.
The average over Euler angles gives

$\displaystyle\left\langle\varepsilon_{4}^{2}\right\rangle=\frac{45^{2}}{4\pi^{2}}\frac{2}{35}\beta_{2}^{4}\int(\sum_{m}\alpha_{4,m}^{\prime}D_{4,m}^{4})(\sum_{m^{\prime}}\alpha_{4,m^{\prime}}^{\prime}D_{4,m^{\prime}}^{4})^{*}\frac{d\Omega}{8\pi^{2}}=\frac{45}{14\pi^{2}}\beta_{2}^{4}\;.$
(24)

Lastly, I consider the case when one only keeps the zero impact parameter
requirement. Ignoring the few nucleons (or a small portion of the volume) that
may not in the overlap region, Eq. (14) is simply the average of two nuclei
with different Euler angles $\Omega_{1}$ and $\Omega_{2}$ 444I am unable to
derive an analytical formula for the contribution of the small portion of
volume not included in the overlap region. For a few special cases, not
including this volume are found to reduce slightly the coefficients of the
$\beta_{n}^{2}$ dependence.

$\displaystyle\varepsilon_{n}e^{in\Phi_{n}}=-A_{n}\frac{\beta_{n}}{2}\sum_{m}\alpha_{n,m}(D_{n,m}^{n}(\Omega_{1})+D_{n,m}^{n}(\Omega_{2}))$
(25)

From this, the mean square average similar to Eq.(15) needs to be integrated
over both $\Omega_{1}$ and $\Omega_{2}$. The crossing terms such as
$D_{n,m}^{n}(\Omega_{1})D_{n,m}^{n}(\Omega_{2})^{*}$ vanish after this
integration, and the final result is exactly half of the original value. This
argument also applies to the $\left\langle\varepsilon_{1}^{2}\right\rangle$ in
Eq.(17) and the quadrupole contribution to $\varepsilon_{4}$ in Eq. (24). So
for a more realistic selection of ultra-central collisions corresponding to
close to zero impact parameter, the coefficients of the $\beta_{n}^{2}$
dependence are a factor of two smaller. These values are listed in the bottom
row of Table 1 and they are closer to Monte Carlo Glauber model result shown
in the bottom row of Figs. 2 and 3.

##  Appendix B More detailed results

For completeness, the full set of correlations between mean square
eccentricities $\left\langle\varepsilon_{n}^{2}\right\rangle$ and deformation
parameters, $\gamma$, $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$ in Eq. (1), are
included here. The $\left\langle\varepsilon_{n}^{2}\right\rangle$ are
calculated using the nucleon Glauber model and quark Glauber model in the U+U
collisions. They are compiled in Figs. 10 and 11 as a function of two
centrality estimators, $N_{\mathrm{part}}$ and $N_{\mathrm{quark}}$,
respectively. From these plots, the ratios of
$\left\langle\varepsilon_{n}^{2}\right\rangle$ to that obtained for default
choice of each parameter are plotted in Figs. 12 and 13, respectively. For the
smaller Zr+Zr collision system, I only show the ratios in Figs. 14 and 15.

The main difference between the two centrality estimators is in the behavior
$\left\langle\varepsilon_{n}^{2}\right\rangle$ in the UCC region, more clearly
visible in the ratio plots. However, whether
$\left\langle\varepsilon_{n}^{2}\right\rangle$ themselves are calculated from
nucleons or quarks have little influences on these ratios. Another important
point is about the contribution of $\beta_{m}$ to
$\left\langle\varepsilon_{n}^{2}\right\rangle$ for $m\neq n$. Although such
mixings could in principle be used to constrain the $\beta_{3}$ using $v_{1}$
as well as $\beta_{4}$ using $v_{1}$ and $v_{3}$, this mixing also forbiddens
a straightforward disentanglement of different deformation components.
Fortunately, such mixing effects are minimal in the UCC collisions, and if one
stay in the 0–1% centrality range, each
$\left\langle\varepsilon_{n}^{2}\right\rangle$ only has one dominating
contribution:
$\left\langle\varepsilon_{n}^{2}\right\rangle_{\rm{UCC}}=a_{n}^{\prime}+b_{n}^{\prime}\beta_{n}^{2}$
for $n=2$, 3 and 4, and
$\left\langle\varepsilon_{1}^{2}\right\rangle_{\rm{UCC}}=a_{1}^{\prime}+b_{3,1}^{\prime}\beta_{3}^{2}$.
In non-central collisions, there is a modest cross correlation between
$\beta_{3}$ and $\left\langle\varepsilon_{2}^{2}\right\rangle$ and between
$\beta_{2}$ and $\left\langle\varepsilon_{3}^{2}\right\rangle$ with the former
having somewhat larger amplitudes. However, in medium-size Zr+Zr collision
system, the cross correlation between $\beta_{3}$ and
$\left\langle\varepsilon_{2}^{2}\right\rangle$ is still significant, but
almost disappears between $\beta_{2}$ and
$\left\langle\varepsilon_{3}^{2}\right\rangle$ (see Figs. 14 and 15). The
reason can be explained as follows. In large U+U collision system, one notes
that the maximum influence of $\beta_{3}$ to
$\left\langle\varepsilon_{2}^{2}\right\rangle$ appears at larger
$N_{\mathrm{part}}$ or $N_{\mathrm{quark}}$ (around 3.5% centrality) than the
location of maximum influence of $\beta_{2}$ to
$\left\langle\varepsilon_{3}^{2}\right\rangle$ (around 6% centrality). In
smaller Zr+Zr collision system, the peak location shifts towards more
peripheral region, maximum influence of $\beta_{3}$ to
$\left\langle\varepsilon_{2}^{2}\right\rangle$ shifts to around 7% centrality
and the maximum influence of $\beta_{3}$ to
$\left\langle\varepsilon_{2}^{2}\right\rangle$ shifts to around 15%
centrality. In the latter case, the fluctuation-driven
$\left\langle\varepsilon_{3}^{2}\right\rangle$ component is much more
important than enhancement from $\beta_{2}$. In the former case, the
$\left\langle\varepsilon_{2}^{2}\right\rangle$ value for undeformed case is
also larger, leading to a smaller relative increase compare to U+U for the
same $\beta_{3}$.

Focusing on the 0–1% centrality range, I then obtain the
$\left\langle\varepsilon_{n}^{2}\right\rangle$ as a function of various
$\beta_{n}^{2}$. The results are summarized in Figs. 16 and 17 for U+U and
Zr+Zr collisions, respectively. In most cases, strict linear dependencies are
observed. One noticeable exception is the relation between
$\left\langle\varepsilon_{4}^{2}\right\rangle$ and $\beta_{2}$, which is
better described by a quartic dependence $\beta_{2}^{4}$ in the UCC region of
U+U, consistent with the analytical results in Table 1. However, in the mid-
central and peripheral U+U collisions and in Zr+Zr over the full centrality
range, one finds that it is still better described by a $\beta_{2}^{2}$
dependence. Lastly, the slopes of these dependencies are nearly independent of
whether nucleons or quarks are used for $\varepsilon_{n}$ or the centrality,
with the exception of the $\gamma$ dependence of
$\left\langle\varepsilon_{2}^{2}\right\rangle$.

Section III.2 discusses briefly the influence of non-axial higher-order
deformation, analogous to the triaxiality for the quadrupole deformation. This
aspect is explored by mixing two different octupole or hexadecapole
components, while keeping the overall magnitude of the deformation to be the
same. Three cases are studied for the octupole deformation,
$1+\beta_{3}(\cos\delta Y_{3,0}+\sin\delta Y_{3,1})$, $1+\beta_{3}(\cos\delta
Y_{3,0}+\sin\delta Y_{3,2})$ and $1+\beta_{3}(\cos\delta Y_{3,0}+\sin\delta
Y_{3,3})$, the results are shown in the left three columns of Fig. 18. Only a
small, less than 3%, dependence on the mixing angle $\delta$ is observed for
$\left\langle\varepsilon_{3}^{2}\right\rangle$. The situation for hexadecapole
is a bit more involved. To simplify the discussion, I consider only the
components respecting all three reflection symmetries, $Y_{4,0}$, $Y_{4,2}$
and $Y_{4,4}$. The nuclear surface can be parametrized with two angular
variables $\gamma_{4}$ and $\delta_{4}$, in addition to $\beta_{4}$ Magierski
_et al._ (1997); Rohozinski (1997):

$\displaystyle
R(\theta,\phi)=R_{0}(1+\beta_{4}(\cos\delta_{4}Z_{0}+\sin\delta_{4}[\cos\gamma_{4}Z_{1}+\sin\gamma_{4}Z_{2}])),\;$
$\displaystyle
Z_{0}=\sqrt{\frac{7}{12}}Y_{4,0}+\sqrt{\frac{5}{12}}Y_{4,4},\;Z_{1}=\sqrt{\frac{5}{12}}Y_{4,0}-\sqrt{\frac{7}{12}}Y_{4,4},\;Z_{2}=Y_{4,2}$
(26)

The parameter $\gamma_{4}$ plays the similar role as the triaxiality parameter
$\gamma$. For example for $\delta_{4}=\mathrm{acos}(\sqrt{7/12})$,
$\gamma_{4}=0,2\pi/3,$ and $4\pi/3$ would correspond to axial-hexadecapole
shape around $z$-, $x$\- and $y$-axis, respectively. The right two columns of
Fig. 18 show results for the two mixing cases,
$1+\beta_{4}(\cos\gamma_{4}Z_{1}+\sin\gamma_{4}Z_{2})$ and
$1+\beta_{4}(\cos\delta_{4}Z_{0}+\sin\delta_{4}Z_{1})$, respectively. A clear
linear dependence on $\cos 3\gamma_{4}$ is observed in the first case. The
dependence in the second case is somewhat more complex, but I do observe it
reaches maximum when $\delta_{4}=0$ or $\pi/2$, for which hexadecapole shape
is described by $Z_{0}=\sqrt{\frac{7}{12}}Y_{4,0}+\sqrt{\frac{5}{12}}Y_{4,4}$
or $Z_{1}=\sqrt{\frac{5}{12}}Y_{4,0}-\sqrt{\frac{7}{12}}Y_{4,4}$ and positive
$\beta_{4}=|\beta_{4}|$. The minimum on the other hand corresponds to the same
shape components but with $\beta_{4}=-|\beta_{4}|$. This behavior is similar
to the influence the prolate vs oblate quadrupole deformation on the
$\left\langle\varepsilon_{2}^{2}\right\rangle$ as seen in Fig. 4.

Figure 10:  The centrality dependence, characterized by $N_{\mathrm{part}}$,
of $\left\langle\varepsilon_{n}^{2}\right\rangle$ in U+U collisions in the
presence of different quadrupole triaxiality $\gamma$ for $\beta_{2}=0.28$
(left column), different axial quadrupole deformation $\beta_{2}$ (second
column), different axial octupole deformation $\beta_{3}$ (third column) and
different axial hexadecapole deformation $\beta_{4}$ (last column) for $n=1$
(top row), $n=2$ (second row), $n=3$ (third row) and $n=4$ (bottom row). The
markers and line curves represent
$\left\langle\varepsilon_{n}^{2}\right\rangle$ calculated from nucleon Glauber
model and quark Glauber model, respectively. The functional form of the
deformation and different parameters are given in the top-row panel for each
corresponding column. Figure 11:  Same as Fig. 10 but using the
$N_{\mathrm{quark}}$ as centrality. Figure 12:  The ratio of
$\left\langle\varepsilon_{n}^{2}\right\rangle$ to the
$\left\langle\varepsilon_{n}^{2}\right\rangle$ in U+U collisions calculated
with the default choice of each parameter (indicated by marks or lines around
unity in each panel); they are obtained directly from Fig. 10. The three solid
vertical bars around unity in each one of these ratio plots indicate the
locations of 2%, 1% and 0.2% centrality. Figure 13:  Same as Fig. 12 but using
the $N_{\mathrm{quark}}$ as centrality. Figure 14:  Same as Fig. 12 but for
Zr+Zr collisions. Figure 15:  Same as Fig. 13 but for Zr+Zr collisions. Figure
16:  The parametric dependence of
$\left\langle\varepsilon_{n}^{2}\right\rangle$ in 0–1% centrality U+U
collisions on the quadrupole triaxiality $\gamma$ for $\beta_{2}=0.28$ (left
column), axial quadrupole deformation $\beta_{2}$ (second column), axial
octupole deformation $\beta_{3}$ (third column) and axial hexadecapole
deformation $\beta_{4}$ (last column) for $n=1$ (top row), $n=2$ (second row),
$n=3$ (third row) and $n=4$ (bottom row). The
$\left\langle\varepsilon_{n}^{2}\right\rangle$ and centrality can be
determined either from nucleon Glauber or quark Glauber, leading to four
different curves in each panel as indicated by the legend in the top-left
panel. They are obtained directly from plots like Figs. 10 and 11, where each
panel provides the two set of data points in the corresponding panel in this
figure. Figure 17:  Same as Fig. 16 but for 0–1% most central Zr+Zr
collisions. Figure 18:  Effects of mixing between three pairs of octupole
components in the left three columns, $1+\beta_{3}(\cos\delta
Y_{3,0}+\sin\delta Y_{3,1})$, $1+\beta_{3}(\cos\delta Y_{3,0}+\sin\delta
Y_{3,2})$ and $1+\beta_{3}(\cos\delta Y_{3,0}+\sin\delta Y_{3,3})$, and
between two pairs of hexadecapole components in the right two columns,
$1+\beta_{4}(\cos\gamma_{4}Z_{1}+\sin\gamma_{4}Z_{2})$ and
$1+\beta_{4}(\cos\delta_{4}Z_{0}+\sin\delta_{4}Z_{1})$, where
$Z_{0}=\sqrt{\frac{7}{12}}Y_{4,0}+\sqrt{\frac{5}{12}}Y_{4,4}$,
$Z_{1}=\sqrt{\frac{5}{12}}Y_{4,0}-\sqrt{\frac{7}{12}}Y_{4,4}$ and
$Z_{2}=Y_{4,2}$ (see text). The results are obtained for 0–1% central U+U
collisions and are presented separately for
$\left\langle\varepsilon_{n}^{2}\right\rangle$, $n=1$, 2, 3, and 4 from the
top to the bottom rows. The $\left\langle\varepsilon_{n}^{2}\right\rangle$ and
centrality are determined either from nucleon Glauber or quark Glauber,
leading to four different curves in each panel as indicated by the legend in
the right panel of the second row.

## References

  * Bohr and Mottelson (1998) A. Bohr and B. R. Mottelson, eds.,  _Nuclear Structure_ (World Scientific, 1998).
  * Heyde and Wood (2011) K. Heyde and J. L. Wood, Rev. Mod. Phys. 83, 1467 (2011).
  * Togashi  _et al._ (2016) T. Togashi, Y. Tsunoda, T. Otsuka, and N. Shimizu, Phys. Rev. Lett. 117, 172502 (2016), arXiv:1606.09056 [nucl-th] .
  * Heyde and Wood (2016) K. Heyde and J. L. Wood, Phys. Scripta 91, 083008 (2016).
  * Frauendorf (2018) S. Frauendorf, Phys. Scripta 93, 043003 (2018), arXiv:1710.01210 [nucl-th] .
  * Zhou (2016) S.-G. Zhou, Phys. Scripta 91, 063008 (2016), arXiv:1605.00956 [nucl-th] .
  * De Vries  _et al._ (1987) H. De Vries, C. W. De Jager, and C. De Vries, Atom. Data Nucl. Data Tabl. 36, 495 (1987).
  * Busza  _et al._ (2018) W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Heinz (2013) U. W. Heinz, J. Phys. Conf. Ser. 455, 012044 (2013), arXiv:1304.3634 [nucl-th] .
  * Nakatsukasa  _et al._ (2016) T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Rev. Mod. Phys. 88, 045004 (2016), arXiv:1606.04717 [nucl-th] .
  * Giacalone  _et al._ (2021a) G. Giacalone, J. Jia, and C. Zhang, Phys. Rev. Lett. 127, 242301 (2021a), arXiv:2105.01638 [nucl-th] .
  * C.G.Joslin and Gray (1983) C.G.Joslin and C. Gray, Molecular Physics 50, 329 (1983).
  * Teaney and Yan (2011) D. Teaney and L. Yan, Phys. Rev. C 83, 064904 (2011), arXiv:1010.1876 [nucl-th] .
  * Niemi  _et al._ (2016) H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph] .
  * Teaney and Yan (2012) D. Teaney and L. Yan, Phys. Rev. C 86, 044908 (2012), arXiv:1206.1905 [nucl-th] .
  * Jia (2014) J. Jia, J. Phys. G 41, 124003 (2014), arXiv:1407.6057 [nucl-ex] .
  * Bernhard  _et al._ (2016) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th] .
  * Bernhard  _et al._ (2019) J. E. Bernhard, J. S. Moreland, and S. A. Bass, Nature Phys. 15, 1113 (2019).
  * Nijs  _et al._ (2021) G. Nijs, W. van der Schee, U. Gürsoy, and R. Snellings, Phys. Rev. Lett. 126, 202301 (2021), arXiv:2010.15130 [nucl-th] .
  * Rosenhauer  _et al._ (1986) A. Rosenhauer, H. Stocker, J. A. Maruhn, and W. Greiner, Phys. Rev. C 34, 185 (1986).
  * Li (2000) B.-A. Li, Phys. Rev. C 61, 021903 (2000), arXiv:nucl-th/9910030 .
  * Gupta and Gale (2000) S. D. Gupta and C. Gale, Phys. Rev. C 62, 031901 (2000).
  * Heinz and Kuhlman (2005) U. W. Heinz and A. Kuhlman, Phys. Rev. Lett. 94, 132301 (2005), arXiv:nucl-th/0411054 .
  * Filip  _et al._ (2009) P. Filip, R. Lednicky, H. Masui, and N. Xu, Phys. Rev. C 80, 054903 (2009).
  * Shou  _et al._ (2015) Q. Y. Shou, Y. G. Ma, P. Sorensen, A. H. Tang, F. Videbæk, and H. Wang, Phys. Lett. B 749, 215 (2015), arXiv:1409.8375 [nucl-th] .
  * Goldschmidt  _et al._ (2015) A. Goldschmidt, Z. Qiu, C. Shen, and U. Heinz, Phys. Rev. C 92, 044903 (2015), arXiv:1507.03910 [nucl-th] .
  * Giacalone  _et al._ (2018) G. Giacalone, J. Noronha-Hostler, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 97, 034904 (2018), arXiv:1711.08499 [nucl-th] .
  * Giacalone (2019) G. Giacalone, Phys. Rev. C 99, 024910 (2019), arXiv:1811.03959 [nucl-th] .
  * Adamczyk  _et al._ (2015) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 115, 222301 (2015), arXiv:1505.07812 [nucl-ex] .
  * ALICE Collaboration (2018) ALICE Collaboration, Phys. Lett. B 784, 82 (2018), arXiv:1805.01832 [nucl-ex] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Phys. Rev. C 100, 044902 (2019), arXiv:1901.07997 [hep-ex] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), Phys. Rev. C 101, 024906 (2020), arXiv:1911.04812 [nucl-ex] .
  * Carzon  _et al._ (2020) P. Carzon, S. Rao, M. Luzum, M. Sievert, and J. Noronha-Hostler, Phys. Rev. C 102, 054905 (2020), arXiv:2007.00780 [nucl-th] .
  * Giacalone (2020a) G. Giacalone, Phys. Rev. Lett. 124, 202301 (2020a), arXiv:1910.04673 [nucl-th] .
  * Giacalone (2020b) G. Giacalone, Phys. Rev. C 102, 024901 (2020b), arXiv:2004.14463 [nucl-th] .
  * Jia  _et al._ (2021) J. Jia, S. Huang, and C. Zhang, (2021), arXiv:2105.05713 [nucl-th] .
  * (37) J. Jia, “Nuclear deformation effects via Au+Au and U+U collisions from STAR,” Contribution to the VIth International Conference on the Initial Stages of High-Energy Nuclear Collisions, January 2021, https://indico.cern.ch/event/854124/contributions/4135480/.
  * Miller  _et al._ (2007) M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007), arXiv:nucl-ex/0701025 .
  * Adler  _et al._ (2014) S. S. Adler _et al._ (PHENIX), Phys. Rev. C89, 044905 (2014), arXiv:1312.6676 [nucl-ex] .
  * Lacey  _et al._ (2018) R. A. Lacey, P. Liu, N. Magdy, M. Csanád, B. Schweid, N. N. Ajitanand, J. Alexander, and R. Pak, Universe 4, 22 (2018), arXiv:1601.06001 [nucl-ex] .
  * Loizides (2016) C. Loizides, Phys. Rev. C94, 024914 (2016), arXiv:1603.07375 [nucl-ex] .
  * Bożek  _et al._ (2016) P. Bożek, W. Broniowski, and M. Rybczyński, Phys. Rev. C94, 014902 (2016), arXiv:1604.07697 [nucl-th] .
  * Acharya  _et al._ (2019) S. Acharya _et al._ (ALICE), Phys. Lett. B790, 35 (2019), arXiv:1805.04432 [nucl-ex] .
  * Zhou and Jia (2018) M. Zhou and J. Jia, Phys. Rev. C 98, 044903 (2018), arXiv:1803.01812 [nucl-th] .
  * Mitchell  _et al._ (2016) J. T. Mitchell, D. V. Perepelitsa, M. J. Tannenbaum, and P. W. Stankus, Phys. Rev. C93, 054910 (2016), arXiv:1603.08836 [nucl-ex] .
  * De Forest and Walecka (1966) T. De Forest, Jr. and J. D. Walecka, Adv. Phys. 15, 1 (1966).
  * Myers and Schmidt (1983) W. D. Myers and K.-H. Schmidt, Nucl. Phys. A 410, 61 (1983).
  * Mazeliauskas and Teaney (2015) A. Mazeliauskas and D. Teaney, Phys. Rev. C 91, 044902 (2015), arXiv:1501.03138 [nucl-th] .
  * Mazeliauskas and Teaney (2016) A. Mazeliauskas and D. Teaney, Phys. Rev. C 93, 024913 (2016), arXiv:1509.07492 [nucl-th] .
  * Hamamoto  _et al._ (1991) I. Hamamoto, N. Xi zhen Zhang, and N. Hong-xing Xie, Phys. Lett. B 257, 1 (1991), [Erratum: Phys.Lett.B 278, 511–511 (1992)].
  * Rohozinski (1997) S. G. Rohozinski, Phys. Rev. C 56, 165 (1997).
  * Lemmon  _et al._ (1993) R. C. Lemmon, J. R. Leigh, J. X. Wei, C. R. Morton, D. J. Hinde, J. O. Newton, J. C. Mein, M. Dasgupta, and N. Rowley, Phys. Lett. B 316, 32 (1993).
  * Butler (2016) P. A. Butler, J. Phys. G 43, 073002 (2016).
  * P  _et al._ (2013) G. L. P _et al._ , Nature 497, 199 (2013).
  * Nazarewicz  _et al._ (1984) W. Nazarewicz, P. Olanders, I. Ragnarsson, J. Dudek, G. A. Leander, P. Möller, and E. Ruchowsa, Nucl. Phys. A 429, 269 (1984).
  * Raman  _et al._ (2001) S. Raman, C. W. G. Nestor, Jr, and P. Tikkanen, Atom. Data Nucl. Data Tabl. 78, 1 (2001).
  * Agbemava  _et al._ (2016) S. E. Agbemava, A. V. Afanasjev, and P. Ring, Phys. Rev. C 93, 044304 (2016), arXiv:1603.03414 [nucl-th] .
  * Libert  _et al._ (1982) J. Libert, M. Meyer, and P. Quentin, Phys. Rev. C 25, 586 (1982).
  * Möller  _et al._ (2016) P. Möller, A. J. Sierk, T. Ichikawa, and H. Sagawa, Atom. Data Nucl. Data Tabl. 109-110, 1 (2016), arXiv:1508.06294 [nucl-th] .
  * Robledo and Bertsch (2011) L. M. Robledo and G. F. Bertsch, Phys. Rev. C 84, 054302 (2011), arXiv:1107.3581 [nucl-th] .
  * Hilaire and Girod (2007) S. Hilaire and M. Girod, Eur. Phys. J. A 33, 237 (2007).
  * KIBÉDI and SPEAR (2002) T. KIBÉDI and R. H. SPEAR, Atom. Data Nucl. Data Tabl. 80, 35 (2002).
  * Giacalone  _et al._ (2021b) G. Giacalone, J. Jia, and V. Somà, Phys. Rev. C 104, L041903 (2021b), arXiv:2102.08158 [nucl-th] .
  * Hammelmann  _et al._ (2020) J. Hammelmann, A. Soto-Ontoso, M. Alvioli, H. Elfner, and M. Strikman, Phys. Rev. C 101, 061901 (2020), arXiv:1908.10231 [nucl-th] .
  * Xu  _et al._ (2021) H.-j. Xu, H. Li, X. Wang, C. Shen, and F. Wang, Phys. Lett. B 819, 136453 (2021), arXiv:2103.05595 [nucl-th] .
  * Huo  _et al._ (2014) P. Huo, J. Jia, and S. Mohapatra, Phys. Rev. C 90, 024910 (2014), arXiv:1311.7091 [nucl-ex] .
  * Alver  _et al._ (2010) B. H. Alver, C. Gombeaud, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 82, 034913 (2010), arXiv:1007.5469 [nucl-th] .
  * Bhalerao  _et al._ (2011) R. S. Bhalerao, M. Luzum, and J.-Y. Ollitrault, Phys. Rev. C 84, 054901 (2011), arXiv:1107.5485 [nucl-th] .
  * Gubser and Yarom (2011) S. S. Gubser and A. Yarom, Nucl. Phys. B 846, 469 (2011), arXiv:1012.1314 [hep-th] .
  * Staig and Shuryak (2011) P. Staig and E. Shuryak, Phys. Rev. C 84, 034908 (2011), arXiv:1008.3139 [nucl-th] .
  * Lacey  _et al._ (2016) R. A. Lacey, D. Reynolds, A. Taranenko, N. N. Ajitanand, J. M. Alexander, F.-H. Liu, Y. Gu, and A. Mwai, J. Phys. G 43, 10LT01 (2016), arXiv:1311.1728 [nucl-ex] .
  * Gardim and Ollitrault (2021) F. G. Gardim and J.-Y. Ollitrault, Phys. Rev. C 103, 044907 (2021), arXiv:2010.11919 [nucl-th] .
  * Ryssens  _et al._ (2019) W. Ryssens, M. Bender, K. Bennaceur, P. H. Heenen, and J. Meyer, Phys. Rev. C 99, 044315 (2019), arXiv:1809.04406 [nucl-th] .
  * Magierski  _et al._ (1997) P. Magierski, K. Burzynski, E. Perlinska, J. Dobaczewski, and W. Nazarewicz, Phys. Rev. C 55, 1236 (1997), arXiv:nucl-th/9611014 .