[a]Takeshi Nakamori

# Collaboration between high schools in Japan and Argentina for cosmic-ray
research using CosmicWatches

Mika Takanashi  Rikako Kono  Ryuta Saito  Kyoka Ozaki  Haruki Nakadai  Kazuo
S. Tanaka  Ana Beatriz Prieto  Marianela Pepe  Juan Wehinger  Lucio Martinez
Chizuru Nose  Haruhi Enomoto  Ryosuke Kita  Kyoka Maruta

###### Abstract

Cosmic rays are ubiquitous and readily available, making them a good teaching
tool for particle and astrophysics by young students. Tan-Q is an inclusive
outreach and educational project, providing students in Japanese junior-high
or high schools with research opportunities to join cosmic-ray and particle
physics. In the Tan-Q framework, the students in each school conduct their
research with help from mentors who are mainly undergraduate students.
Researchers are also extensively involved through regular Zoom meetings and
continuous communication on Slack. Some cases are inter-school, and some are
international. This paper presents one of the Tan-Q activities of joint
research between high schools in Japan and Argentina to observe cosmic-ray
muons using CosmicWatches. Our primary goal is to investigate the muon flux
differences due to the differences in circumstances like altitudes and
geomagnetic field strengths. Those involved learn not only particle physics
but also statistical data analysis methods.

##  1 Introduction

In recent years, the importance of learning through scientific research
activities has been widely recognized in school science education in Japan,
and such inquiry activities are being conducted in schools across the country.
There are multiple perspectives, such as fostering logical thinking skills and
active learning. The research themes and forms of activities vary from school
to school, but they often depend on the school, teachers, and other factors.
Moreover, there is sometimes a disparity between schools that have been
selected for the Super Science High School Program and those that have not.
There is no doubt that the field of particle and astrophysics plays a major
part in modern physics. Even though Japan has several accelerators and
observatories, it is difficult for high school students to choose their
research theme. Some of the reasons for this may include the detectors and
data acquisition systems required to measure particles and that it is not easy
for high school students to use large facilities independently. As well as the
Japanese situation, Argentina is promoting Science, Technology, Engineering,
and Mathematics (STEM) education through research-based learning, which has
recently been incorporated into the official curriculum of the country.
Moreover, the nation hosts one of the world’s greatest cosmic-ray experiments,
the Pierre Auger Observatory [1]. Although such major facilities often promote
educational programs in both countries, the opportunity is relatively limited
for nationwide students with potential interests.

On the other hand, especially secondary muons, cosmic rays can be observed
anywhere and easily studied by junior and high school students with proper
detectors and guidance. Thus, cosmic rays are a good teaching tool for
particle and astrophysics for young students. There have been some successful
examples, such as (a part of) Quarknet [2]. However, only a few cases have
existed in Japan [3], where researchers collaborate with junior and senior
high school students to observe cosmic rays. This case was categorized as
citizen science and was positioned as a collaborative project on a scientists’
theme. In Japan, citizen science projects have been conducted in recent years,
along with intelligent devices and network communication technologies. One
recent and successful example involving particle measurements is the
Thundercloud project [4] to observe gamma rays generated in thunderclouds. In
these cases, scintillation detectors are provided by scientists and installed
in schools for observation together with students.

Researches by junior and senior high school students can be divided into two
major categories: research on their themes and so-called citizen science, as
introduced above, where the researchers bring the objectives. Suppose a
framework can comprehensively support both directions; it should be possible
to provide opportunities for particle space research to a broader range of
junior and senior high school students.

##  2 Tan-Q framework

Tan-Q 111Tan-Q’s pronunciation is the same as the Japanese word corresponding
to quest, seek and research, inclusive outreach and educational project, was
triggered along with such context and situation (PI: K. S. Tanaka, 2019–).
Tan-Q consists of junior and senior high school students, mentor students from
universities, and researchers. As shown in Figure 1, schools from nationwide
Japan participate in the program and research a wide variety of topics. Some
as individual schools and some as teams of students from different schools,
including international, work on a single topic. We assign a dedicated
mentor(s) or scientist(s) for each research group to support the school
students in their research activities. We contact each other through Slack
daily and provide regular mentoring through Zoom to continuously support the
students’ research. In addition, we share files and meeting logs via Notion
and Goole Drive. By making full use of these online tools, we realize
sustainable activities without depending on the region.

Figure 1: Map showing Tan-Q’s participating institutions. Circles indicate the
locations of high schools and colleges of technology, and squares and
triangles indicate those of the institutions to which the researchers and
mentor students belong, respectively.

Many research groups in Tan-Q use CosmicWatch (CW), a small and inexpensive
desktop cosmic ray muon detector developed by Axani et al. [7] (see also
section 3.2). Through the measurement workshops organized by Tan-Q, we make
the CW known to the students and schools. And then, the CW is distributed to
schools that are interested in CW research activities. Some schools are using
CW in its original form, for example, to search for a correlation between
weather conditions and solar activity and cosmic ray rates. In other cases,
they have modified the scintillator and data acquisition programs to research
their topics.

In this paper, we present a case that is one of the Tan-Q activities with CWs,
a collaboration between Japan and Argentina.

##  3 Case: AR-JP joint research

###  3.1 Research objective

The Earth’s magnetic field has been monitored extensively on Earth’s surface
since the 1800s and also by satellites these days. There is an area known as
the South Atlantic Anomaly (SAA), with a weak geomagnetic field strength
around the South Atlantic Ocean. The SAA area is continuously growing, and its
center is moving on a 100-year scale. This phenomenon reflects the dynamics of
the Earth’s interior and is a sign of a reversal of the Earth’s magnetic
field, which is being actively studied [5]. The geomagnetic field blocks
cosmic rays (CRs) arriving at the Earth, and CRs with energies lower than the
cutoff rigidity (typically a few GeV) cannot penetrate the magnetosphere. As a
result, the magnetic field over the SAA is weaker than in other regions, and
then the CRs, having a power-law spectrum, come into the atmosphere more
frequently.

Since Argentina is close to the SAA, the arrival frequency of cosmic-ray muons
is expected to be higher than in Japan. This study aims to detect the excess
in the rate of Argentina due to the weak geomagnetic field strength. Since the
rate of muons varies with altitude, it is necessary to correct this effect and
then compare the results. It should be noted that Cecire et al. [6] performed
cosmic ray measurements associated with SAA in Chile within the framework of
QuarkNet.

###  3.2 Students’ activity : Measurements

The CW [7] consists of a plastic scintillator of $\sim 5\times 5\times 1$ cm3
and a silicon photomultiplier. Arduino Nano handles the data acquisition,
which is triggered when the pulse height exceeds a threshold value. Then pulse
height, time stamp, and environmental data such as temperature are recorded
for each event. PCB pattern and Arduino image files are available as open
source. Our project provided the CWs assembled by the mentors to the
observation teams in Japan and Argentina. They connected the CWs to a Windows
PC via USB to read the data.

Figure 2: Observation sites indicated in Yamagata prefecture, Japan (left),
and Argentina (right). Note that these two maps are on a different scale.

Figure 2 shows the locations of the observation sites. We measured the average
rate of muons at three sites in Yamagata prefecture, the northeastern district
in Japan, to derive the correlation between altitude and frequency. Two sites
of the three are high schools with relatively low altitudes, and the other is
a kind of a pension at the ski resort with a high altitude above 1300 m. By
measuring within the minimal latitude and longitude range of Yamagata
Prefecture, we can accurately track the differences in altitude.

In order to have the common solar condition which may affect the primary
cosmic ray rate, we performed simultaneous measurements in Japan and Argentina
during the same period. We took data continuously for several weeks from
January to March of 2021, every 24 hours without interruption if possible.
Then we shared the acquired data files on Google Drive both domestically and
internationally. The mentors prepared such environments and the analysis
template source codes on Google Colaboratory. Then the students analyzed the
data with Python, mounting the Google Drive for efficient data handling.

###  3.3 Students’ activity : Analyses

The flow of the analysis is as follows. The first step is to apply a cut to
the pulse-height ADC to select the qualified muon events. Figure 3(a) shows an
example spectrum. We know that the event rate near the trigger threshold
varies with temperature. If we choose the component around the peak at the
higher ADC value, the rate stabilizes independent of temperature. Another team
studied this finding in detail, working within the Tan-Q framework (Akita High
School and Yamagata Higashi High School). It is also expected from the
separate beta irradiation experiment that muons passing near the edge of the
plastic scintillator also produce a low ADC component. In other words, muons
striking near the edges are difficult to distinguish from thermal noise
triggers. Since in this study, only one CW, or one scintillator, was used at
each measurement location, we applied tighter cuts to the data aiming for
accurate measurements of muons. In parallel, we cut the periods when data
acquisition failed accidentally, mainly due to communication problems between
the CW and the PC. Figure 3(b) shows the event rates per 10 minutes in a time
series. We exclude the periods with low rates deviating from the statistical
fluctuations. We applied these data cleaning procedures to all data sets for
each location.

Next, we obtained the count rate per minute from the surviving events and
produced count rate histograms. Then we derived the average rate of muons by
fitting with a Gaussian, as shown in Figure 3(c). The rest of the project is
in progress, and the students are currently learning statistics and error
analysis. They are trying to fit a rate-altitude curve from data from three
sites in Yamagata Prefecture, Japan, and test whether the rate in Argentina is
significantly higher than that expected. Currently, we proceed to produce a
figure like Figure 3(d). The rate-alt dependence is expected to be exponential
[8], but since the altitude range is narrow, a linear fit could also be
likely.

Figure 3: Summary of the analysis flow. (a) The ADC spectrum for a 1-day run.
Two components are identified below/above the ADC value of 200. Apparent
periodic spikes are artificial. (b) Recorded event rate without the ADC cut. A
slight drop can be found before 12:00, and such time intervals are removed
from the analysis. (c) Histogram of the muon event rate derived from all-
integrated data set of Argentina. With sufficient statistics, the distribution
is well consistent with the Gaussian. (d) Sketch of the near future analysis
(see text).

Further studies are worth consideration. In order to improve the accuracy of
the rate-alt curve within Japan, it is conceivable to add data measured by
other Tan-Q teams with CWs. Data from a very close area may be favorable, and
the data collecting area is still limited to the relatively small region,
Japan. The benefit of having more data points would be even significant. CW
can connect two units and operate in coincidence mode [8]. This operation will
significantly improve the signal-to-noise ratio of the muons at the hardware
level. The more accurate the rate measurement will become, the more enhanced
the sensitivity to the excess detection will be. Comparison of the muon rates
with the models using EXPACS [9, 10] may be interesting. Long-term monitoring,
as long as a few years, could find an evolution of SAA. The CW can be a potent
tool to see the global scale phenomena through easy desktop experiments.

##  4 Concluding remark

Although Tan-Q is still a young project, it functions on a wide scale in Japan
with various activities supported by many mentors and researchers.
Furthermore, through advanced research activities, which are sometimes beyond
the standard curriculum, the project has proven highly educational for junior
and senior high school students and even mentor university students. From the
students’ point of view, Tan-Q promotes their research and, at the same time,
allows them to see other junior and senior high school students engaged in
similar research activities. They also connect with experts (sometimes outside
Tan-Q) through regular debriefing sessions held twice a year. Most students
who are conducting research activities in a single school cannot gain such
experience. Furthermore, one of the proofs of the program’s success is that
some students went on to university through the activities of the Tan-Q and
came back as mentors. In evaluating the results of the educational project, it
is desirable to collect objective indicators and data soon. Interviews are a
possible method to manage the outcomes;

  1. 1.

How do junior and senior high school students change their interests in
physics-related fields (not necessarily limited to particle and astrophysics)?

  2. 2.

How the school teachers evaluate from their point of view?

  3. 3.

What the mentors gain both as learners and educators through the Tan-Q
activities?

As for this Argentina-Japan project, it is rare both for Argentine and
Japanese high school students to research with a contact outside the country.
Therefore, this project is unique and noteworthy in the education in both
countries, as it fosters an international perspective through the global-scale
theme of studying the effects of geomagnetism.

## Acknowledgment

Tan-Q and this work is financially supported by Mitsubishi Memorial Foundation
for Educational Excellence (Category 3, 2020).

## References

  * Pierre Auger Collaboration [2015] A. Aab, et al. (Pierre Auger Collaboration), The Pierre Auger Cosmic Ray Observatory, Nuclear Instruments and Methods A, 798, 172-213, 2015
  * Bardeen et al. [2018] M. Bardeen, M. Wayne, and M. Jean Young, QuarkNet: A Unique and Transformative Physics Education Program, Educ. Sci., 8, 17(10pp), 2018
  * Teramoto et al. [2003] Y. Teramoto, Trial for constructing a cosmic-ray observation network within science museums and high schools via internet, Daigaku no butsuri kyouiku (Physics Education in University), 2003.2, 65-67, 2003 (Japanese)
  * Yuasa et al. [2020] T. Yuasa, Y. Wada, T. Enoto, Y. Furuta, H. Tsuchiya, S. Hisadomi, et al., Thundercloud Project: Exploring high-energy phenomena in thundercloud and lightning, Progress of Theoretical and Experimental Physics, 10k id.103H01, 2020 [arXiv:2007.13618]
  * Pavón-Carrasco et al. [2018] F. J. Pavón-Carrasco, and A. De Santis, The South Atlantic Anomaly: The Key for a Possible Geomagnetic Reversal, Frontiers in Earth Science, 4, 40, 2018
  * Cecire et al. [2019] K. Cecire, S. Carrilo, J. Chorny, D. Gayoso, R. González, M. Haacke, S. Olivares, J. Wegner, Outreach, Investigation, Muons, and the South Atlantic Anomaly, PoS, ICRC2019, 043, 2019
  * Axani et al. [2018] S. N. Axani, K. Frankiesivcs, and J. M. Conrad, The cosmicwatch desktop muon detector: self-contained, pocket sized particle detector, Journal of Instrumentation, 13, P03918, 2018
  * Axani et al. [2019] S. N. Axani, The Physics Behind the CosmicWatch Desktop Muon Detectors, arXiv:1908.00146, 2019
  * Sato et al. [2016] T. Sato, Analytical Model for Estimating the Zenith Angle Dependence of Terrestrial Cosmic Ray Fluxes, PLOS ONE, 11, e0160390, 2016
  * Sato et al. [2015] T. Sato, Analytical Model for Estimating Terrestrial Cosmic Ray Fluxes Nearly Anytime and Anywhere in the World: Extension of PARMA/EXPACS, PLOS ONE, 10, e0144679. 2015


# Comparisons and Predictions for Collisions of deformed 238U nuclei at
$\sqrt{s_{NN}}=\qty{193}{\giga}$

Nicolas Miro Fortier<EMAIL_ADDRESS>Department of Physics,
McGill University, 3600 University Street, Montreal, QC, H3A 2T8, Canada
Sangyong Jeon<EMAIL_ADDRESS>Department of Physics, McGill
University, 3600 University Street, Montreal, QC, H3A 2T8, Canada  Charles
Gale<EMAIL_ADDRESS>Department of Physics, McGill University, 3600
University Street, Montreal, QC, H3A 2T8, Canada

###### Abstract

We present comparisons to experimental data along with predictions of
observables for U+U collisions at 193 GeV using a multistage theoretical and
computational framework consisting of boost-invariant IP-Glasma initial state,
MUSIC hydrodynamics, and a hadronic transport cascade generated by iS3D &
SMASH. Our results show great agreement with existing anisotropic flow
measurements from RHIC [1, 2]. We provide predictions for differential flow
observables as well as multiparticle correlations and transverse-momentum-flow
correlations. When possible, we compare our predictions to results from Au+Au
collisions at 200 GeV to properly outline the effects of deformation in the
initial state on final state observables.

##  I Introduction

Heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the
Large Hadron Collider (LHC) have provided remarkable insights into the
fundamental properties of matter under extreme conditions [3]. One of the most
intriguing phenomena observed in these collisions is the creation of quark
gluon plasma (QGP), a state of matter characterized by the deconfinement of
quarks and gluons [4]. The QGP, produced in the aftermath of collisions,
exhibits collective behaviour reminiscent of a nearly perfect fluid. This
remarkable hydrodynamic behaviour allows for the conversion of initial state
anisotropies into final state observables, providing a unique window into not
only the dynamics of the QGP, but also nuclear structure and its limits.

In studies involving spherically symmetric nuclei, researchers have
extensively investigated the effects of specific subsets of initial state
anisotropies resulting from nuclear geometry fluctuations. However, a
significant gap remains in our understanding of the impact of initial state
anisotropies originating from deformed nuclei, such as 238U. These deformed
nuclei introduce a wider range of anisotropies, presenting an opportunity to
explore unique combinations of fluctuations and test certain phenomenological
hypotheses within the QGP [5].

Notably, our study predicts that central collisions involving prolate-shaped
nuclei (such as 238U) should lead to a negative correlation between $v_{2}$
and $p_{T}$. In contrast, this correlation is observed to be positive at all
centralities in collisions of spherical nuclei [6, 7]. Indeed, the prolate
geometry of 238U causes higher eccentricity ($\varepsilon_{2}$) events to
generate lower $\left\langle p_{T}\right\rangle$ and vice-versa, inducing an
anti-correlation of the two observables. These collisions also generate higher
energy densities in ultra-central configurations compared to spherical nuclei.
This extra energy density is expected to have observable effects on various
crucial observables, including elliptic flow and jet quenching [8], which
serve as key probes to characterize the properties of the QGP, as well as
collective properties and momentum correlations [9, 10, 11].

While considerable efforts have been made in studying heavy-ion collisions
involving spherical nuclei [12, 13, 14, 15], experimental data for collisions
involving deformed nuclei (such as 238U) remains limited. Collisions of
deformed nuclei, however, present an exciting opportunity for testing our
understanding of QGP dynamics. To achieve this, it is crucial to develop a
comprehensive framework that can incorporate all types of fluctuations at all
stages.

In this study, we employ an up-to-date, comprehensive and well-motivated
theoretical framework to both compare to available data from STAR [1] and
extract predictions for a wide range of observables in U+U collisions at
$\sqrt{s_{NN}}=\qty{193}{\giga}$. Simulations begin with the IP-Glasma model,
which is based on the Colour Glass Condensate (CGC) framework [16, 17] and
provides realistic event-by-event colour fluctuations as well as pre-
equilibrium flow. IP-Glasma has been successful at reproducing key observables
across both the energy and collision system spectra [18, 19, 20, 21, 22, 23,
24]. MUSIC, a relativistic viscous hydrodynamic simulation which incorporates
bulk and shear viscosities [25, 26], then takes over the evolution of the
system. Once the density and temperatures of the QGP drop sufficiently, the
system is ‘frozen-out’, meaning that fluid cells which were previously
governed by hydrodynamics are converted into hadrons before being propagated.
This final stage of the simulations is taken up by iS3D [27] & SMASH [28], two
models which combine to generate hadronic cascades. The inner workings of the
various phases of our simulations will be discussed in detail in section II.

This paper follows the following structure: section II examines the different
stages of our physical model in varying detail, ranging from an in-depth
discussion of IP-Glasma in subsection II.1 to a more cursory look at hadronic
transport model SMASH II.3. Section III details the different studied
observables and their working definitions, as well as theoretical expectations
for their sensitivity to the initial state when possible. These include
centrality selection (III.1), flow analysis methods (III.3.1) and transverse-
momentum-flow correlations (III.3.3). We then show our model’s comparative and
predictive capabilities and results in section IV. The paper will end with a
summary and conclusion, presented in section V.

##  II Theoretical Model

###  II.1 IP-Glasma

Historically, the first hydrodynamic models of heavy-ion collisions were
developed with an emphasis on the fluid dynamics of QGP [29, 30, 31, 32, 33,
34, 35, 36, 37, 38, 39, 40]. The initial states used by these first models
were mainly geometric in nature [41]. Once the relevance and success of early
heavy-ion collision simulations were established, the need for a more
physically accurate and detailed initial state model became apparent. IP-
Glasma, a QCD- and saturation-based model, was first introduced in 2012 [42]
and quickly became the standard in the field. It is based on the Colour Glass
Condensate (CGC) effective field theory [43, 44, 45] and classical gluon
production [46, 47, 48, 49, 50].

To model heavy-ion collisions, one must first generate the nuclei. In this
study, we use the deformed Woods-Saxon distribution, given by

$\displaystyle\rho(r,\theta)=\frac{\rho_{0}}{1+\exp\left(\frac{r-R(\theta)}{a}\right)}$
(1) $\displaystyle
R(\theta)=R_{0}\left(1+\beta_{2}Y_{2}^{0}(\theta)+\beta_{4}Y_{4}^{0}(\theta)\right)$
(2)

to generate nucleon configurations. Here, $\rho_{0}$ denotes the nuclear
density, $R_{0}$ is the unmodified nuclear radius and $a$ is the nuclear skin
depth. Parameters $\beta_{l}$ multiply the spherical harmonic functions
$Y_{l}^{0}(\theta)$ and generate deformation about the $x$ and $y$ axes. The
parameters used in this study of 238U are $R_{0}=\qty{6.784}{\femto}$,
$a=\qty{0.556}{\femto}$, $\rho_{0}=1$, $\beta_{2}=0.28$ and
$\beta_{4}=-0.0035$ and were taken from [51]. A comparison between an
undeformed and a deformed nucleus using the parameters used in this study is
presented in Fig. 1. Other parameters, such as the triaxial deformation
parameter $\gamma$ and the hextupole deformation parameter $\beta_{3}$ may be
factored into $R(\theta)$ to obtain a whole suite of nuclear deformations.
However, these parameters are yet to be measured for 238U.

The Woods-Saxon distribution is a simplistic yet effective model for sampling
nucleon positions, as evidenced by its widespread use in the field. However,
as with any simple model, one must understand its limitations. The $\beta_{l}$
deformation parameters are model-dependent, in so far as their extraction must
figure in the Woods-Saxon distribution. Also, the distribution itself does not
account for important nucleon-nucleon correlations, which more recent studies
have shown to have sizable effects in generating physically accurate nuclei
configurations [52, 53]. Therefore, it should be understood that all of the
results presented in this paper carry an inherent systematic uncertainty.

Figure 1:  Comparing a regular and deformed Woods-Saxon distribution with same
base radius $R_{0}$ and skin depth $a$. The deformed distribution generates
and oblong (or pill-shaped) profile. The rotational symmetry axis is the long
horizontal axis.

Following nucleon sampling, the impact parameter for the event is sampled from

$\displaystyle P(b)db=\frac{2b}{b^{2}_{\mathrm{max}}-b^{2}_{\mathrm{min}}}db$
(3)

where $b_{\mathrm{max}}=\qty{8}{\femto}$ and
$b_{\mathrm{min}}=\qty{0}{\femto}$. The boundaries were fixed with the goal of
rejecting as few events as possible due to their peripherality since
deformation effects are only tangible in central collisions of deformed
nuclei. The nuclei configurations are then shifted symmetrically by ${b}/{2}$
in the $x$-direction. The spatial distribution of nucleons is then projected
into the transverse plane.

At this point, IP-SAT [54], the impact parameter dependent dipole saturation
model, takes over. Its contribution will be in providing the saturation scale
$Q_{s}$ at all points in the transverse plane, which then allows us to sample
colour charges and initialize the system’s colour gauge fields. To do so, IP-
SAT first models the nuclear density as

$\displaystyle T(\boldsymbol{x})=\frac{e^{-\boldsymbol{x}^{2}/2B_{G}}}{2\pi
B_{G}}$ (4) $\displaystyle
T_{A}(\boldsymbol{x})=\sum_{i=1}^{A}T(\boldsymbol{x}-\boldsymbol{x_{i}})$ (5)

where $A$ represents the current nucleus’ number of nucleons and
$B_{G}=\qty{4.0}{\giga^{-2}}$ is extracted from a fit to DIS data [55]. With
the thickness function for each nucleus in hand, we solve

$\displaystyle\frac{2\pi^{2}}{N_{c}}T_{A,B}(\boldsymbol{x})r_{s}^{2}xg(x,\mu^{2}(r_{s}^{2}))\alpha_{s}(\mu^{2}(r_{s}^{2}))=1$
(6)

for $Q_{s}^{2}={2}/{r_{s}^{2}}$, where $N_{c}=3$ is the number of colours
permitted and $T_{A,B}(\boldsymbol{x})$ is the thickness function described in
Eq.(5) for the projectile (A) and target (B) nuclei. $xg(x,\mu^{2})$, the
density of gluons at a given scale $\mu$ and momentum fraction $x$, is
initialized as

$\displaystyle xg(x,\mu_{0}^{2})=A_{g}x^{\lambda_{g}}(1-x)^{5.6}$ (7)

with $A_{g}=2.308$, $\lambda_{g}=0.058$ and
$\mu_{0}^{2}=\qty{1.51}{\giga^{2}}$. Eq.(7) is then evolved to all other
values of $\mu^{2}$ using the leading-order DGLAP equation [56, 57, 58]. The
scale $\mu$ itself is related to the saturation dipole size (and scale) by

$\displaystyle\mu^{2}=\frac{4}{r_{s}^{2}}+\mu_{0}^{2}=2Q_{s}^{2}+\mu_{0}^{2}$
(8)

which sets the scale in the leading-order QCD running coupling constant given
by

$\displaystyle\alpha_{s}(\mu^{2})=\frac{12\pi}{(33-2N_{f})\ln\left(\tfrac{\mu^{2}}{\Lambda_{QCD}}\right)}.$
(9)

where $N_{f}$ is the number of quark flavours, set to $4$ in our simulation.
Eq.(6) is itself extracted from the Glauber-Mueller dipole cross-section [59].
Solving for $Q_{s}$ must be done iteratively given the intricate
interdependence of the various functions ($xg$, $\alpha_{s}$) and variables
($x$, $r_{s}$, $\mu$, $Q_{s}$).

Once $Q_{s}^{2}$ is determined, the colour charge distribution for the
projectile nucleus, for instance, can be sampled from the following colour
correlator

$\displaystyle\left\langle\rho^{a}_{A}(\boldsymbol{x})\rho^{b}_{A}(\boldsymbol{y})\right\rangle=g^{2}\mu^{2}_{A}(x,\boldsymbol{x})\delta^{ab}\delta^{2}(\boldsymbol{x}-\boldsymbol{y}),$
(10)

where $x$ is the momentum fraction currently considered and
$Cg^{2}\mu_{A}=Q_{s}$111It is important to note here that $\mu_{A}\neq\mu$.
$\mu$ is the intrinsic energy scale at hand, while $\mu_{A}$ is the scale of
colour charge fluctuations, which is related to the saturation scale.. We are
therefore sampling from a Gaussian of width proportional to the saturation
scale. The proportionality constant $C$ is determined phenomenologically and
was set to $0.505$ in this study.

This sampled colour charge distribution is of great importance, as it acts as
a source for the small-x gluon fields which comprise the CGC. The CGC action,

$\displaystyle S_{CGC}=\int
d^{4}x\left(-\frac{1}{4}F^{a}_{{\mu\nu}}F^{{\mu\nu}}_{a}+J^{\mu
a}A^{a}_{\mu}\right)$ (11)

contains a current term $J^{\mu a}$ which sources the colour gauge fields
$A_{\mu}^{a}$. The convention for the covariant derivative and, therefore, the
field strength tensor used throughout this paper is

$\displaystyle D_{\mu}=\partial_{\mu}+igA_{\mu}$ $\displaystyle
F_{{\mu\nu}}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}+ig\left[A_{\mu},A_{\nu}\right].$

The corresponding Classical Yang-Mills (CYM) equation is

$\displaystyle\left[D_{\mu},F^{{\mu\nu}}\right]=J^{\nu}$ (12)

where

$\displaystyle J^{\nu}=\rho_{A}(\boldsymbol{x})\delta^{\nu+}\delta(x^{-})$
(13)

Here, the two $\delta$-functions indicate a right-moving source on the light
cone. Note the move to light-cone coordinates $x^{\pm}=(t\pm z)/\sqrt{2}$,
meaning that our sources (large-x partons) are travelling at the speed of
light. In the pre-collision phase, the CYM equations in the $A_{-}=0$ gauge
are

$\displaystyle\nabla^{2}_{\perp}A^{a}_{+}=-\rho^{a}$ (14)

which are Poisson equations for each colour index. The more physical light-
cone gauge gluon fields can be obtained by a gauge transformation

$\displaystyle
A^{A}_{\mu}=V^{A}A_{\mu}V^{A\dagger}+\frac{-i}{g}V^{A}\partial_{\mu}V^{A\dagger}$
(15)

where

$\displaystyle
V(x^{-},\boldsymbol{x})=\mathcal{P}\exp\left(-ig\int_{-\infty}^{x^{-}}dy^{-}A_{-}(y^{-},\boldsymbol{x})\right)$
(16)

In this 2D setting, only the transverse components $A^{A}_{x}$ and $A^{A}_{y}$
are non-zero.

Once the pre-collision colour gauge fields have been determined, the collision
takes place; both fields are combined such that the forward light cone has the
following initial fields

$\displaystyle A_{i}=A^{A}_{i}+A^{B}_{i}$ (17) $\displaystyle
E^{\eta}=ig\left[A_{i}^{A},A_{i}^{B}\right],$ (18)

where $i=x,y$. The dynamics inside the forward light cone are best described
using $\tau-\eta$ coordinates, which we will use in the rest of this work. To
obtain the other components of the chromo-electric field $E^{i}$, we must
solve the covariant form of Gauss’ law,

$\displaystyle\left[D_{i},E^{i}\right]+\left[D_{\eta},E^{\eta}\right]=0.$ (19)

Since derivatives in $\eta$ vanish in our boost-invariant implementation, we
are left with $\left[D_{i},E^{i}\right]=0$, which has trivial solution
$E^{i}=0$: the initial transverse chromo-electric fields are therefore set to
0 [60, 61].

Once the initial post-collision fields are settled, we evolve the whole system
using the sourceless CYM equations,

$\displaystyle[D_{\mu},F^{{\mu\nu}}]=0$ (20)

until $\tau_{\rm hyd}=\qty{0.52}{\femto}$. At that time, the CYM stress-energy
tensor, given by

$\displaystyle
T^{{\mu\nu}}=\textrm{Tr}\left(-g^{\mu\gamma}g^{\nu\alpha}g^{\beta\delta}F_{\gamma\beta}F_{\alpha\delta}+\frac{1}{4}g^{\mu\nu}g^{\gamma\beta}g^{\alpha\delta}F_{\gamma\alpha}F_{\beta\delta}\right)$
(21)

is constructed. It is symmetric and gauge invariant, and is the bridge that
connects the pre-equilibrium dynamics of the glasma to the relativistic
hydrodynamics of the QGP. We refer to [62] for an in-depth discussion about
the properties and development of the tensor. To find the energy density
$\varepsilon$ and flow velocity $u^{\mu}$, we diagonalize $T^{{\mu\nu}}$ and
preserve the timelike eigenvalue. The flow velocity is normalized to
$u_{\mu}u^{\mu}=1$ throughout.

The shear-stress tensor $\pi^{{\mu\nu}}$, which is needed to initialize
viscous hydrodynamics, is given by

$\displaystyle\pi^{{\mu\nu}}=T^{{\mu\nu}}_{\mathrm{IPG}}-T^{{\mu\nu}}_{\mathrm{ideal}}$
(22) $\displaystyle
T^{{\mu\nu}}_{\mathrm{ideal}}=(\varepsilon+P)u^{\mu}u^{\nu}-Pg^{{\mu\nu}}$
(23)

In IP-Glasma, the pressure $P$ is given by ${\varepsilon}/{3}$ due to the
conformality of the classical gluon system. The conformal nature of the pre-
equilibrium phase, originating from the pure gluon and classical features of
the CGC, also justifies the absence of a bulk pressure $\Pi$ at this stage. On
the hydro side, the pressure is dictated by the EoS at use, which in this
study is HotQCD [63]. An issue arises when one realises that
$P_{\mathrm{IPG}}={\varepsilon}/{3}$ and $P_{\mathrm{EoS}}(\varepsilon)$ may
not match, leading to a discontinuity in our transition to hydrodynamics. This
issue is handled by initializing the QGP with a bulk pressure $\Pi$, which, at
this point, is required and given by the difference between the IPG and EoS
pressures, i.e.

$\displaystyle\Pi=P_{\mathrm{IPG}}(\varepsilon)-P_{\mathrm{EoS}}(\varepsilon)=\frac{\varepsilon}{3}-P_{\mathrm{EoS}}(\varepsilon)$
(24)

We therefore fully conserve energy and momentum in our transition to
hydrodynamics, which allows for proper tracking and accounting of quantities
such as ${dE}/{d\eta}$ from the initial state into the hydrodynamics phase.

###  II.2 MUSIC

While the CGC is an effective field theory, hydrodynamics is a long-wavelength
effective theory [64, 65]. Both, however, have the same objective: to describe
reality as efficiently and accurately as possible, while providing pragmatic
theories which allow for the extraction of usable results. In the
hydrodynamics phase, the evolution shifts from more fundamental degrees of
freedom such as the gluon fields to coarse-grained thermodynamic ensemble
averages, such as pressures and temperatures.

MUSIC [25, 66] is the numerical implementation of the following theoretical
concepts. Formally, it is a second-order relativistic viscous hydrodynamics
simulation. The fundamental quantity being evolved in hydrodynamics is the
energy-momentum tensor

$\displaystyle
T^{{\mu\nu}}_{\mathrm{Hydro}}=T^{{\mu\nu}}_{\mathrm{ideal}}+\Pi\left(u^{\mu}u^{\nu}-g^{{\mu\nu}}\right)+\pi^{{\mu\nu}}$
(25)

In second order viscous hydrodynamics, the conservation laws

$\displaystyle\partial_{\mu}T^{{\mu\nu}}_{\mathrm{Hydro}}=0$ (26)

are supplemented by equations of motion for the the viscous tensor
$\pi^{{\mu\nu}}$ and the bulk pressure $\Pi$

$\displaystyle\dot{\Pi}=\frac{1}{\tau_{\Pi}}\left(-\Pi-\zeta\Theta-\delta_{\Pi\Pi}\Pi\Theta\right)$
(27)
$\displaystyle\dot{\pi}^{\langle{\mu\nu}\rangle}=\frac{1}{\tau_{\pi}}\left(-\pi^{{\mu\nu}}+2\eta\sigma^{{\mu\nu}}-\delta_{\pi\pi}\pi^{{\mu\nu}}\Theta\right)$
(28)

Here, $\Theta=\partial_{\mu}u^{\mu}$ is the scalar expansion rate and
$\sigma^{{\mu\nu}}=\nabla^{\langle\mu}u^{\nu\rangle}$ is the velocity shear
tensor. The overdot represents the co-moving frame time derivative
$u^{\mu}\partial_{\mu}$ and the angular bracket around the indices indicate
the transverse, symmetric and traceless part of the tensor. The coefficients
$\delta_{\Pi\Pi}/\tau_{\Pi}$ and $\delta_{\pi\pi}/\tau_{\pi}$ are derived
using the 14-moment approximation, whose values are ${2}/{3}$ and ${4}/{3}$
respectively, according to [67]. The shear viscosity $\zeta$ and the bulk
viscosity $\eta$ have the following form

$\displaystyle\frac{\zeta}{s}(T)=\frac{0.282\Lambda^{2}}{\Lambda^{2}+\left(T-0.311\right)^{2}}$
(29)
$\displaystyle\Lambda=0.029\left[1-0.970\mathrm{sign}\left(T-0.311\right)\right]$
(30) $\displaystyle\frac{\eta}{s}=0.136$ (31)

where $s$ is the entropy density and the energy unit is GeV. In this study,
following the thorough Bayesian analysis performed in [21], we have used a
temperature-dependent bulk viscosity $\zeta/s$, while keeping shear viscosity
$\eta/s$ constant. With time, the QGP grows in volume while its temperature
drops. At a certain temperature, the QGP will ‘hadronize’, i.e. turn into
hadrons. The exact value of this temperature is not sharply defined [68, 69,
70]. In this study, the switching temperature is set to $T_{\rm
sw}=\qty{155}{\mega}$. Once a specific fluid cell cools down to temperature
$T_{\rm sw}$, its spatiotemporal state is saved. Once every fluid cell has
reached $T_{\rm sw}$, all of the 4-dimensional states are combined to generate
a constant-temperature hyper-surface, which terminates the hydrodynamic stage
of the simulation.

###  II.3 iS3D & SMASH

The freeze-out hyper-surface generated by MUSIC is fed into iS3D [27], a
particlization code which implements Cooper-Frye sampling [71], i.e.

$\displaystyle
E\frac{dN_{i}}{d^{3}p}=\frac{d_{i}}{\left(2\pi\right)^{3}}\int_{\Sigma}f_{i}(x,p)p_{\mu}d\sigma^{\mu}(x)$
(32)

where $\Sigma$ is an isothermal hyper-surface and $\sigma^{\mu}(x)$ is its
normal vector, $\frac{dN_{i}}{d^{3}p}$ is the momentum spectrum of particle
species $i$, $f_{i}(x,p)$ is its phase-space distribution and $d_{i}$ is the
degeneracy factor. To ensure a smooth transition between hydrodynamics and
particlization, the energy-momentum tensor $T^{{\mu\nu}}$ must be reproduced
everywhere on the hyper-surface. Therefore,

$\displaystyle
T^{{\mu\nu}}_{\mathrm{kin}}=\sum_{i}d_{i}\int\frac{d^{3}p}{\left(2\pi\right)^{3}E}p^{\mu}p^{\nu}\left(f_{{\rm
eq},i}(x,p)+\delta f_{i}(x,p)\right)$ (33)

where $f_{i}(x,p)$ has been decomposed into an equilibrium distribution (which
follows either Bose-Einstein or Fermi-Dirac statistics depending on the
species $i$) and an out-of-equilibrium correction $\delta f_{i}$. This
correction is necessary to account for the viscous nature of our hydrodynamic
evolution. Indeed, these lead to a medium which is out-of-equilibrium at the
time of sampling which in turn produces slight deviations from the equilibrium
distributions $f_{{\rm eq},i}$. To match the nature of our shear stress
$\pi^{{\mu\nu}}$ and bulk viscous pressure $\Pi$, we use the 14-moment $\delta
f_{i}$ corrections. It is an expansion of $\delta f_{i}$ which is truncated at
terms of first- and second-order in momentum ($p^{\mu}$ and $p^{\mu}p^{\nu}$)
[72, 73, 74].

An important concept, which will be revisited in more detail in subsection
III.2, is the fact that the freeze-out hyper-surface from a single
hydrodynamic event is oversampled hundreds of times [20, 22, 21]. Indeed,
because the hyper-surface stems from a hydrodynamic treatment of the QGP,
which itself is an ensemble average, the sampling of particles will converge
to the hydrodynamic value of all observables (multiplicity, momenta, flow)
only once a sufficient number of samplings are done. Therefore, a single IP-
Glasma and MUSIC event, comprised of a unique collision, impact parameter and
nuclei configuration, can generate hundreds of iS3D events, each consisting of
its own particle list containing specific species and momenta. In this work,
we test two different ways of treating these oversampled events. The
oversampled events from a single hydrodynamic event can be averaged over and
treated as a single event, relating it to its initial and hydrodynamic stages
uniquely. Alternatively, each oversampled event can be regarded as a distinct
event fitting the general prescription provided by the ensemble average hyper-
surface. Choosing one method over the other has tangible effects on computed
observables, as will be evidenced in our results.

Once the hadrons are sampled, they are evolved kinetically using SMASH [28], a
hadronic cascade code. It implements inter-particle interactions and
scatterings, as well as resonance and decays via coupled Boltzmann equations
of most of the known hadrons. We used SMASH Version 1.8 in this work.

##  III Methods and Observables

###  III.1 Centrality

Proper centrality selection is key to ensuring that the theoretical and
computational models of HICs are comparable to available experimental data. In
this study, matching multiplicities across centrality classes served as the
sole calibration of our model: the proportionality constant $C=0.505$,
described in section II.1, was calibrated to reproduce multiplicity
distributions following the approach advocated in [22] and described below.
This is the only calibration we made. All observables extracted thereafter
used this calibrated value.

When all events have gone through all stages of our framework, they are sorted
by charged particle multiplicity $dN_{\textrm{ch}}/d\eta$, then separated into
a sufficiently many number of bins, all of which contain the same number of
events. The $C$ parameter is calibrated such that the average multiplicity of
the most central bin matches that of experiment. Two experimental centrality
bins and their respective multiplicities are then selected. These are the most
central and the most peripheral centrality we would like to analyze (in this
study, $0-1\%$ and $27-28\%$ respectively). The experimental multiplicity
ratios of these two bins is then computed. The same is done with the first and
last bins from the simulations. While the ratio of the average multiplicities
of our two selected bins exceeds that of the corresponding experimental ratio,
we drop the lowest-multiplicity event from our consideration, recalculate the
ratio, and keep repeating this until the desired ratio is achieved, allowing
us to reject a negligible number of events compared to what would have to be
rejected in a fully minimally-biased study. The calibration and peripheral
event dropping steps are sufficient to ensure that we reproduce available
charged hadron multiplicity curves (see Fig. 2 below), which then allows for a
thorough comparison to all available observables.

Another important aspect of this study was the faithful emulation of Zero-
Degree Calorimeter (ZDC) binning to compare to data from STAR [1]. The ZDC is
a calorimeter that resides at $0^{\circ}$ from the detector’s point of
collision. It aims to measure neutrons which were a part of the collision
systems but didn’t participate in an event. Their lack of electric charge
means that once they are free from the confines of their nucleus, the
collider’s electromagnetic fields do not affect them, leading them to follow
straight paths directly into the calorimeters.

Experimentalists assess the centrality of a collision by counting the number
of detected neutrons: the more neutrons were found, the higher the chance that
the collision was peripheral. To emulate ZDCs within our framework, we
calculate the total number of participating nucleons from a given collision
event and subtract it from the total number of nucleons available to give us
the number of spectator nucleons $S$, i.e.

$\displaystyle S=2A-N_{\mathrm{Participants}}$ (34)

where $A=238$ for a U+U collision. To obtain the number of neutrons out of the
total number of spectator nucleons, we sample a binomial distribution

$\displaystyle
P(N)=\binom{S}{N}\left(1-\frac{Z}{A}\right)^{N}\left(\frac{Z}{A}\right)^{S-N}$
(35)

where Z is the atomic number ($92$ for U) and we aim to sample $N$, the number
of neutrons, as done in [75]. We average $20$ samplings of this distribution
per event to reduce variability, giving us the number of spectator neutrons
for each event. This method does overlook some points, such as the fact that
atomic nuclei have neutron skins (an outer shell where only neutrons are
found) [76, 77] which lead to higher probabilities of having spectator
neutrons than protons which aren’t encapsulated within our simple binomial
distribution. These considerations, however, were outside of the scope of this
study.

###  III.2 Averaging

In section II.3, we briefly discussed the concept of oversampling the
hydrodynamical hyper-surface. Given the ensemble average definition of this
hyper-surface, sampling it multiple times is an important part of ensuring
that the final state observables associated with an event converge to their
hydrodynamic values. Previously [20, 22, 21], this oversampling procedure
(which, in this work, we will call ‘oversampled average’) is followed by an
averaging of the oversampled events to create a single set of observables
related to a single hyper-surface. By design, this method smoothes out
fluctuations and accentuates the features associated with a given hyper-
surface. If, however, each oversampled event is regarded as an independent
event and averaging is only performed in a given centrality class, the effects
of short-range fluctuations and correlations will remain (we will call this
method ‘SMASH sub-event averaging’). However, events that came from the same
hyper-surface event may not be fully independent. Not all observables are
sensitive to these differing averaging methods. But, as will become evident in
our results section, many observables studied in this work are sensitive to
the averaging method.

An ideal simulation of one experimental event would be the chain of one IP-
Glasma, one hydro, and one SMASH event. This, however, requires many orders of
magnitude more computing resources than are currently available. As such, it
is important to analyze which observables are sensitive and for what reasons,
as we have done in this study. For this purpose, we also compare the above two
averaging methods to the ‘mixed events’ in Section IV, where we group all
oversampled events from a given centrality class, mix all of their particles
and create new ‘mixed’ events. Only un-correlated fluctuations should survive
in these mixed events.

###  III.3 Defining Selected Observables

The following parts will briefly introduce and define a selection of the
observables which will figure in our results in section IV.

####  III.3.1 Flow analysis

The n-th anisotropic flow coefficient $v_{n}$ is by now generally accepted as
one of the primary evidence of QGP undergoing fluid-like behaviour in
relativistic heavy-ion collisions. In this study, we will be interested in the
2- and 4-particle cumulants of various components of the flow harmonics. To
start, we define the flow vector $Q_{n}$ for each event [78],

$\displaystyle Q_{n}=\sum^{N_{\mathrm{ch}}}_{j=1}e^{in\phi_{j}}$ (36)

where $N_{\mathrm{ch}}$ is the event’s multiplicity, $j$ runs over all of the
particles of the event with transverse momentum restricted to
$\qty{0.2}{GeV}<p_{T}<\qty{2.0}{GeV}$ to conform with the STAR acceptance
window, and $\phi_{j}$ is the azimuthal angle of the $j^{\mathrm{th}}$
particle. Then, the $2^{\mathrm{nd}}$ order azimuthal correlation is given by

$\displaystyle\left\langle
2\right\rangle=\frac{\left|Q_{n}\right|^{2}-N_{\mathrm{ch}}}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)}$
(37)

while the $4^{\rm th}$ order azimuthal correlation is

$\left\langle
4\right\rangle=\frac{\left|Q_{n}\right|^{4}+\left|Q_{2n}\right|^{4}-2\operatorname{Re}\left[Q_{2n}Q_{n}^{*}Q_{n}^{*}\right]}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)(N_{\mathrm{ch}}-2)(N_{\mathrm{ch}}-3)}\\\
-2\frac{2(N_{\mathrm{ch}}-2)\cdot\left|Q_{n}\right|^{2}-N_{\mathrm{ch}}(N_{\mathrm{ch}}-3)}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)(N_{\mathrm{ch}}-2)(N_{\mathrm{ch}}-3)}.$
(38)

where $Q_{2n}$ is to be understood as the flow vector associated with the
2n-th harmonic if we are calculating the $4^{\mathrm{th}}$ order azimuthal
correlation of the n-th harmonic (i.e. if $n=2$, then $Q_{2n}=Q_{4}$). We then
take an average of these correlations over the entirety of events in their
centrality class, which finally allows us to compute the 2- and 4-particle
cumulants, i.e.

$\displaystyle v_{n}\\{2\\}=\sqrt{\left\langle\left\langle
2\right\rangle\right\rangle}$ (39) $\displaystyle
v_{n}\\{4\\}=\sqrt[4]{-\left(\left\langle\left\langle
4\right\rangle\right\rangle-2\cdot\left\langle\left\langle
2\right\rangle\right\rangle^{2}\right)}$ (40)

where $\left\langle\left\langle\cdot\right\rangle\right\rangle$ denotes
$\left\langle\cdot\right\rangle$ averaged over the given centrality.

Finally, the 2-particle scalar product $p_{T}$-differential flow is given by

$\displaystyle v_{n}\\{2\\}(p_{T})=\frac{\operatorname{Re}(\left\langle
Q_{n}^{\rm PI}(p_{T})\cdot(Q_{n}^{\rm ref})^{*}\right\rangle)}{\left\langle
N_{\rm ch}^{\rm PI}(p_{T})N_{\rm ch}^{\rm ref}\right\rangle v_{n}^{\rm
ref}\\{2\\}}$ (41)

where the superscript ‘PI’ denotes the particle species of interest, while the
superscript ‘ref’ denotes the reference flow vector. To avoid self-
correlations being represented in this observable, $Q_{n}^{\rm PI}$ is taken
from the $|\eta|<0.5$ rapidity window, while the reference flow vector is
taken from $0.5<\eta<2$.

For a more thorough treatment and discussion of these quantities, along with
their respective errors, see [79].

####  III.3.2 Multi-particle Transverse Momentum Correlators

We will be presenting results for 2- and 3-particle $p_{T}$ correlations,
which are sometimes referred to as ‘variance’ and ‘skewness’ respectively, the
difference being that the correlators do not consider self-correlations. The
event-averaged 2-particle transverse momentum correlator is defined as

$\displaystyle\left\langle\delta p\delta
p\right\rangle=\left\langle\frac{\sum_{i\neq j}(p_{i}-\left\langle
p_{T}\right\rangle)(p_{j}-\left\langle
p_{T}\right\rangle)}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)}\right\rangle$ (42)

where the sum is over particles in a given event, the averaging is over the
given centrality class and $\left\langle p_{T}\right\rangle$ denotes the
average transverse momentum in the centrality class being analyzed. The
3-particle version is

$\left\langle\delta p\delta p\delta p\right\rangle=\\\
\left\langle\frac{\sum_{i\neq j\neq k}(p_{i}-\left\langle
p_{T}\right\rangle)(p_{j}-\left\langle p_{T}\right\rangle)(p_{k}-\left\langle
p_{T}\right\rangle)}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)(N_{\mathrm{ch}}-2)}\right\rangle$
(43)

where $\left\langle p_{T}\right\rangle$ denotes the average transverse
momentum over the centrality class being analyzed. Implementing these formulae
numerically as they are presented would be unwise, as they would run at least
in $O(N_{\mathrm{events}}\cdot N_{\mathrm{ch}}^{2})$ and
$O(N_{\mathrm{events}}\cdot N_{\mathrm{ch}}^{3})$ respectively. To avoid such
computationally taxing and redundant computations, we have implemented a
modified version of a framework presented by Giacalone et al. [80]. We start
by defining $P_{n}$, the modified moments of the $p_{T}$ distributions,

$\displaystyle P_{n}=\sum_{i}^{N_{\mathrm{ch}}}(p_{i}-\left\langle
p_{T}\right\rangle)^{n}$ (44)

where $p_{i}$ is the transverse momentum of the $i^{\mathrm{th}}$ particle in
the given event. Then one can easily show

$\displaystyle\sum_{i\neq j}(p_{i}-\left\langle
p_{T}\right\rangle)(p_{j}-\left\langle p_{T}\right\rangle)=(P_{1})^{2}-P_{2}$
(45) $\sum_{i\neq j\neq k}(p_{i}-\left\langle
p_{T}\right\rangle)(p_{j}-\left\langle p_{T}\right\rangle)(p_{k}-\left\langle
p_{T}\right\rangle)=\\\ (P_{1})^{3}-3P_{2}P_{1}+2P_{3}$ (46)

allowing for the following redefinitions of the 2- and 3-particle transverse
momentum correlators

$\displaystyle\left\langle\delta p\delta
p\right\rangle=\left\langle\frac{(P_{1})^{2}-P_{2}}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)}\right\rangle$
(47) $\displaystyle\left\langle\delta p\delta p\delta
p\right\rangle=\left\langle\frac{(P_{1})^{3}-3P_{2}P_{1}+2P_{3}}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)(N_{\mathrm{ch}}-2)}\right\rangle$
(48)

All of the modified moments are computable in linear time on an event-by-event
basis, greatly reducing the computational stress required to extract such
observables from large particle lists. Also, since both quantities have been
reduced to a single term compared to the moments presented in [80], the
calculation of their respective errors has been greatly simplified.

####  III.3.3 Transverse-momentum-flow correlations

The final selected observable integrates both the flow harmonics and the
2-particle transverse momentum correlator. It is a correlator between
2-particle flow harmonics and average transverse momentum which was developed
in [81]. It is defined as

$\displaystyle\rho(v_{n}\\{2\\}^{2},\left\langle
p_{T}\right\rangle)=\frac{\operatorname{cov}(v_{n}\\{2\\}^{2},\left\langle
p_{T}\right\rangle)}{\sqrt{\operatorname{var}\left(v_{n}^{2}\right)\cdot\left\langle\delta
p\delta p\right\rangle}}$ (49)

where

$\operatorname{cov}(v_{n}\\{2\\}^{2},\left\langle p_{T}\right\rangle)=\\\
\left\langle\frac{\left|Q_{n}\right|^{2}-N_{\mathrm{ch}}}{N_{\mathrm{ch}}(N_{\mathrm{ch}}-1)}\cdot\left(\frac{\sum_{i=1}^{N_{\mathrm{ch}}}p_{i}}{N_{\mathrm{ch}}}-\left\langle
p_{T}\right\rangle\right)\right\rangle$ (50)

and

$\displaystyle\operatorname{var}\left(v_{n}^{2}\right)=v_{n}\\{2\\}^{4}-v_{n}\\{4\\}^{4}$
(51)

This correlator will be important in highlighting specific properties of
central collisions of deformed nuclei. Indeed, it should show marked
differences when compared to results from collisions of spherically symmetric
nuclei.

##  IV Results and Discussion

Our results section will be divided into two subsections; the first contains
comparisons of our model to existing data from two RHIC detectors (STAR and
PHENIX) for U+U collisions at $\sqrt{s_{NN}}=\qty{193}{\giga}$ when available
and Au+Au collisions at $\sqrt{s_{NN}}=\qty{200}{\giga}$ otherwise, while the
second will focus on predictions of our model regarding multiparticle
correlations and transverse-momentum-flow correlations in U+U collisions.

###  IV.1 Descriptions of Existing Data

####  IV.1.1 Charged Hadron Multiplicity

We begin by verifying that our model can reproduce charged hadron multiplicity
at midrapidity. As mentioned in section III.1, this observable serves as the
sole calibration tool for the proportionality constant between the saturation
scale $Q_{s}$ and colour charge fluctuations $\mu_{A}$.

Figure 2:  Charged hadron multiplicity in $|\eta|<0.5$ as a function of
centrality in our model, compared to results for $\qty{193}{\giga}$ U+U and
$\qty{200}{\giga}$ Au+Au collisions at STAR [1].

In Fig. 2, we show the number of charged particles per unit pseudo-rapidity
$dN_{\rm ch}/d\eta$ from our model compared to data from STAR at
$\qty{193}{\giga}$ (U+U) and $\qty{200}{\giga}$ (Au+Au). The experimental data
stems from a parametrization undertaken in [1]. Our model’s agreement with the
U+U data is excellent throughout. The most peripheral point, at $28\%$
centrality, dips slightly compared to the rest of our curve. This is due to
our centrality selection procedure outlined in section III.1. Indeed, our
procedure is bound to allow for events which are ‘too’ peripheral to be found
in our most peripheral bin, given that we reject events based on multiplicity
(and, therefore, peripherality) until the multiplicity ratio matches that of
the experiment. Compared to the Au curve, we find similar features, and the
fact that U collisions yield more hadrons, which is sensible given their
larger nucleonic content (and, therefore, larger total collision energy).

Our relatively narrow centrality window ($0-28\%$) is due to a conscious
choice and focus: the differences between collisions of deformed nuclei and
collisions of spherical nuclei are most prominent in central collisions. We
sought to limit our scope to more central collisions to generate sufficient
statistics in the region of interest without using computational resources to
simulate more peripheral events where the differences are not particularly
noticeable.

Figure 3:  Identified particle multiplicity in $|y|<0.5$ as a function of
centrality in our model, compared to results for $\qty{193}{\giga}$ U+U
collisions at STAR [82].

Fig. 3 shows identified particle yields as a function of centrality compared
to data from STAR [82]. For kaons and pions, the agreement is great across the
entire centrality window. Our model underestimates the proton yield by more
than $50\%$ throughout. It is important to note here that our model does not
include a baryon chemical potential. At this collision energy, $\mu_{B}$ is
small but non-zero. This may influence the proton yield.

Figure 4:  Identified particle multiplicity in $|y|<0.5$ scaled by average
number of participant nucleon pairs in the centrality class $\left\langle
N_{\rm part}\right\rangle/2$ as a function of $\left\langle N_{\rm
part}\right\rangle$ in our model, compared to results for $\qty{193}{\giga}$
U+U collisions at STAR [82].

Finally, Fig. 4 shows identified particle yields scaled by the average number
of participant nucleon pairs in a given centrality class $\left\langle N_{\rm
part}\right\rangle/2$ as a function of the number of participants. This set of
results is highly dependent on the results shown in Fig. 3, as the number of
participant nucleons and centrality are highly correlated. However, this
specific observable looks to identify where particle production comes from at
a given centrality, and how it progresses across the spectrum. Because it
increases with the average number of participant nucleons, we determine that
particle production is guided by a combination of soft and hard productions
that scale differently with $N_{\rm part}$.

####  IV.1.2 Average Transverse Momentum

Figure 5:  Identified particle mean transverse momentum $\left\langle
p_{T}\right\rangle$ in $|y|<0.5$ as a function of centrality in our model,
compared to results for $\qty{193}{\giga}$ U+U collisions at STAR [82]. Figure
6:  Two- and four-particle cumulants of elliptic flow ($v_{2}\\{2\\}$ and
$v_{2}\\{4\\}$) as functions of (left) charged particle multiplicity and
(right) centrality, compared to results for $\qty{193}{\giga}$ U+U collisions
at STAR [1]. The shaded bands represent statistical errors.

We now move to the average transverse momentum of identified particles. Fig. 5
shows identified particle mean transverse momentum $\left\langle
p_{T}\right\rangle$ as a function of centrality. The respective masses of the
identified particles are responsible for the ordering. Again, as with yields,
kaons and pions show excellent agreement, while our model overestimates
$\left\langle p_{T}\right\rangle$ for protons. This overestimation makes sense
given our underestimation of the yields: the sampled protons need to
compensate for their lack of numbers with an increase in momentum. Finally, it
is important to note that given our use of hydrodynamics, $\left\langle
p_{T}\right\rangle$ puts strong constraints on the $p_{T}$ spectrum, which
entails that a good agreement with $\left\langle p_{T}\right\rangle$ is
equivalent to reproducing the spectrum [83]. Therefore, given our excellent
agreement with experimental data, the shear and bulk viscosity
parametrizations obtained via Bayesian analysis [21] should be favoured over
other available parametrizations.

####  IV.1.3 Anisotropic Flow

We now shift our attention to integrated elliptic and triangular flow and
anisotropic flow in general. Fig. 6 shows the two- and four-particle cumulants
of elliptic flow ($v_{2}\\{2\\}$ and $v_{2}\\{4\\}$) as functions of charged
hadron multiplicity and centrality. We see that our model reproduces both
observables very well. At smaller multiplicities (or more peripheral
collisions), our model overestimates elliptic flow. This is due to the same
effect apparent in Fig. 2, namely that our peripheral centrality class may
contain more peripheral events which weren’t rejected by our centrality
selection process. It should be emphasized that no additional adjustments of
parameters were made to produce our results.

One may notice that our model does not reach as high in multiplicity as the
experimental data does. This is because the experimental points beyond
$dN_{\textrm{ch}}/d\eta\approx 850$ constitute ultra-ultra-central collisions.
That is, they are in the top $0.01\%$ of events registered at STAR; these
experimental events are rare and would require a much larger number of runs on
our part to reproduce. Nevertheless, the right side panels show that, in terms
of centrality, we are essentially reproducing the entire spectrum of central
events.

These right-side panels also serve to compare Au and U data more accurately.
We see that experimental $v_{2}\\{2\\}$’s are extremely similar throughout the
centrality spectrum, except in the very central collisions ($0-5\%$). There, U
data shows a marked increase before coming back to Au levels between $0-1\%$.
This is explained by the deformed geometry of U. In central collisions of
spherically symmetric nuclei, only one overlap shape is generated in the
transverse plane, namely a circle, which has small (or 0) eccentricity. This
in turn leads to a small elliptic flow in the final state. In central
collisions of deformed nuclei such as U, many different transverse cross-
sections are possible. Indeed, looking at Figs. 7 and 8 we see that a nucleus
travelling with its long axis parallel to the beam direction will have a
circular cross-section in the transverse plane. Similarly, if its short axis
is parallel to the beam direction, it will have an elliptical transverse
cross-section. Therefore, in collisions of deformed nuclei with near-zero
impact parameters, we can have both extremely eccentric (body-body) and
circular (tip-tip) cross-sections. The body-body collisions will have a
smaller energy density (due to their larger transverse overlap area) than
their tip-tip counterparts, which will in turn lead to slightly smaller
multiplicities. Therefore, the marked increase in $v_{2}\\{2\\}$ in $1-3\%$
centrality is due to body-body collisions, and its sharp decrease in ultra-
central ($0-1\%$) collisions is due to tip-tip collisions.

Figure 7:  Schematic representation of the asymmetry between a) the short-axis
and b) the long-axis directions. In the short axis direction, the nucleus is
not deformed (has constant R in WS parametrization). Figure 8:  Schematic
representation of the difference in nuclei overlap area between body-body and
tip-tip collisions. The former presents large eccentricity while the latter
presents near-zero eccentricity.

This figure also introduces an alternative averaging method described in
section III.2, namely the SMASH sub-event averaging. This sub-event averaging
method leads to the expression of short-range correlations which are usually
suppressed by oversampled averaging. We note that this method leads to an
overestimation of both $v_{2}\\{2\\}$ and $v_{2}\\{4\\}$ across the entire
centrality range. As it stands, neither the oversampling averaging nor the
sub-event averaging faithfully follows what takes place in a real heavy-ion
collision because hydrodynamics is an inherently coarse-grained theory of
ensemble-averaged quantities. As such, one should regard the results of the
two different averaging procedures as a part of theoretical uncertainty ($\sim
10\,\%$).

In Fig. 9 we added a mixed event curve which was left out of Fig. 6 for
clarity. By mixing the events, only the average effect of the collision
geometry should survive, whereas the effects due to the deformation will be
washed out. This is indeed what can be observed in Fig. 9: it clearly shows
that deformation effects are crucial in understanding flow in U+U collisions.

Fig. 10 shows the two-particle cumulant of elliptic flow as a function of
scaled multiplicity in two ultra-central ZDC bins: $0-0.125\%$ and $0-1\%$.
These events were selected based on their respective number of sampled
spectator neutrons, as described in section III.1. The most central bin
suffers from a small number of events, which in turn affects our statistics.
Nevertheless, the scale of the experimental data is reproduced. In the broader
ultra-central bin, our model reproduces the general shape and trend of
experimental data. It overestimates elliptic flow at smaller multiplicities
while underestimating flow at higher multiplicities. This is likely due to the
effects discussed earlier in this section regarding high-multiplicity
experimental events. Indeed, the overestimation of elliptic flow in body-body
events (scaled $N_{\rm ch}<1$) indicates that some of these events should have
generated fewer particles, while the underestimation of elliptic flow in tip-
tip events (scaled $N_{\rm ch}>1$) should have generated more particles. We
leave a study focused on improving statistics for the future.

Fig. 11 shows the two-particle cumulant of triangular flow $v_{3}\\{2\\}$ as a
function of charged particle multiplicity, which is a fluctuation-driven
observable. Looking at Fig. 11, we see that experimental data for U and Au are
practically overlapping in central collisions, confirming that initial global
geometry plays little role in this observable. Our model underestimates
triangular flow across our chosen range, which indicates that it
underestimates initial state fluctuations. This could potentially be mended by
the addition of sub-nucleonic degrees of freedom (i.e. valence quark
configurations) in our initial construction of the nuclear thickness function
$T_{A}(\boldsymbol{x})$, such as those described in [84]. We also see that,
while our SMASH sub-event average curve is slightly higher than our ensemble
average curve, the addition of short-range correlations isn’t sufficient to
cover the gap between our model and experimental data.

Figure 9:  Two-particle cumulant of elliptic flow ($v_{2}\\{2\\}$) as
functions of charged particle multiplicity, emphasizing the addition of a
‘mixed event’ curve, compared to results for $\qty{193}{\giga}$ U+U collisions
at STAR [1]. Figure 10:  Two-particle cumulant of elliptic flow
($v_{2}\\{2\\}$) as functions of scaled charged particle multiplicity for
(top) $0-0.125\%$ and (bottom) $0-1\%$ most cetral events, compared to results
for $\qty{193}{\giga}$ U+U collisions at STAR [1]. Figure 11:  Two-particle
cumulant of triangular flow ($v_{3}\\{2\\}$) as a function of charged particle
multiplicity, compared to results for $\qty{193}{\giga}$ U+U collisions at
STAR [2].

###  IV.2 Predictions

We now move on to predictions of our model for multi-particle correlations and
transverse-momentum-flow correlations for $\sqrt{s_{NN}}=\qty{193}{\giga}$ U+U
collisions.

####  IV.2.1 Differential flow

Figure 12:  Charged hadron differential anisotropic flow coefficients
$v_{2}\\{2\\}$, $v_{3}\\{2\\}$ and $v_{4}\\{2\\}$ as functions of transverse
momentum $p_{T}$ for various centrality classes for U+U at $\qty{193}{\giga}$,
compared to results for $\qty{200}{\giga}$ Au+Au collisions at PHENIX [14].

Fig. 12 displays our predictions for the differential $v_{n}$ for U+U at
$\sqrt{s_{NN}}=193\,\hbox{GeV}$ compared to the experimental results from
Au+Au at $\sqrt{s_{NN}}=200\,\hbox{GeV}$. For $v_{3}$, We find that our
model’s calculations of $v_{3}\\{2\\}(p_{T})$ for U+U is similar to the event-
plane $v_{n}$ data from PHENIX [14] for Au+Au222The PHENIX data here does not
use the scalar particle method.. This is expected because $v_{3}$ depend
mainly on local fluctuations which are similar in the two systems.
Differential elliptic flow ($v_{2}$) for U+U is larger across all centrality
classes than $v_{2}$ for Au+Au. For the $0-10\%$ centrality class, this is
consistent with the deformation effects discussed in section IV.1.3: we expect
the elliptic flow to be enhanced in this region of the centrality spectrum
because of the elliptic shape of 238U nuclei. In the $10-20\%$ and $20-30\%$
classes, the differences are smaller and are themselves consistent with
elliptic flow being noticeably larger for U+U collisions compared to Au+Au
throughout the collision spectrum, as evidenced in Fig. 6. Our model’s $v_{4}$
progressively deviates from the Au+Au results going from the central to
peripheral collisions, but the deviation is less noticeable than $v_{2}$.

Figure 13:  Identified particle differential elliptic flow coefficients
$v_{2}\\{2\\}$ as a function of transverse momentum $p_{T}$ for various
centrality classes for U+U at $\qty{193}{\giga}$, compared to results for
$\qty{200}{\giga}$ Au+Au collisions at STAR [15].

Finally, Fig. 13 shows our model predicts a larger differential elliptic flow
for identified particles than what was measured in Au+Au collisions at STAR
[15]. In the ultra-central regions ($0-5\%$ & $5-10\%$), the effect is clear
and crosses hadronic lines. However, in the more peripheral regions ($10-20\%$
& $20-30\%$), this difference becomes much smaller, and varies considerably
from one hadron to another; our model’s predictions for anti-protons is
similar to the Au experimental data. Higher-$p_{T}$ differential flow seems to
converge consistently across all identified particles.

####  IV.2.2 Multi-particle Momentum Correlations

Multi-particle correlations, when compared to future experimental data and
combined with primary observables such as elliptic flow, will help constrain
238U’s deformation parameters and will allow for further analysis using
various sets of deformation parameters [85].

In the top panel of Fig. 14, we show the 2-particle $p_{T}$ correlator. We
have included a mixed event curve to ensure that no underlying $p_{T}$
correlations exist - its position on the plot confirms this. A sizeable
difference exists between our oversampled event average and SMASH sub-event
average curves, which are diametrically opposed with respect to our mixed
event curve. This indicates that the inclusion (or exclusion) or short-range
correlations are a key determinant of the behaviour of this observable.
Experimental data will be key in determining whether or not these short-range
correlations play an important role in observed 2-particle correlations.

The bottom panel of Fig. 14 shows the 3-particle $p_{T}$ correlator. In
contrast to the upper panel, both curves are extremely similar, with the
oversampled average only slightly larger than the SMASH sub-event averaging
values. The oversampled averaging techniques give out larger correlations,
which serve as a further contrast to the behaviour of the curves in the upper
panel. Once again, the mixed event curve is included to ensure that no
implicit correlations exist for this observable.

####  IV.2.3 Transverse-Momentum-Flow Correlation

We end our results section with plots showing predictions for the correlations
between integrated anisotropic flow and mean transverse momentum. This
observable has garnered interest as a tell-tale sign of deformation [80].

In Fig. 15, we see that both averaging methods present very similar behaviour
across both plots. The mixed event curve could not be included in this plot as
it is undefined across our centrality interval. Indeed, referring to Eq.(49),
the correlator $\rho(v^{2}_{n},\left\langle p_{T}\right\rangle)$ requires a
division by $\left\langle\delta p_{T}\delta p_{T}\right\rangle$ which is
evenly $0$ across our centrality interval, as evidenced by Fig. 14.
Triangular-flow-momentum (bottom panel) correlations seem to be dominated by
fluctuations, with the only real clear trend being that the correlator remains
positive throughout the centrality range. Elliptic-flow-momentum correlations,
on the other hand, have a clear and relatively stable trend, with a crossover
from positive to negative correlation happening at around $7\%$ centrality.

In central collisions of spherical nuclei [6, 7], we observe a dip in the
correlator due to the correlation between the inverse of the transverse
overlap area (larger mean $p_{T}$) and initial state eccentricity
$\epsilon_{2}$ (which gets smaller in more central collisions). However, no
anti-correlation is observed. Our model predicts that for U+U an anti-
correlation should be observed in central collisions. We can make sense of
this prediction by using the correlation between inverse overlap area and
eccentricity. Indeed, tip-tip collisions will generate relatively small
overlap areas that are high in energy density (leading to high $\left\langle
p_{T}\right\rangle$) while generating almost no eccentricity and therefore,
elliptic flow. Therefore, when contrasted with other events (such as body-body
events) in a given central centrality class, the correlator finds that
$\left\langle p_{T}\right\rangle$ and $v_{2}\\{2\\}$ are anti-correlated, as
events in these classes having lower $\left\langle p_{T}\right\rangle$ and
higher eccentricities contrast other events in the same centrality class
having larger $\left\langle p_{T}\right\rangle$ and smaller eccentricities.

##  V Summary and Conclusion

We have presented a detailed synthesis of our multi-phase model consisting of
IP-Glasma initial state, MUSIC viscous relativistic hydrodynamics, and iS3D +
SMASH particle sampling and transport, and have shown its results for a wide
variety of observables for U+U collisions at $\sqrt{s_{NN}}=\qty{193}{\giga}$.
The purpose of this work is to describe the published data and to present
timely predictions for observables specifically relevant to deformed nuclei.

Our model shows good agreement across all available experimental data, which
was composed of charged hadron multiplicity, identified particle yields and
integrated anisotropic flow. It underestimated the fluctuation-driven two-
particle cumulant of triangular flow but was still well within $15\%$ of the
provided experimental data across our entire centrality range. It also
underestimated proton yield and overestimated proton $\left\langle
p_{T}\right\rangle$ for protons as a result. This provides an opportunity to
revisit the parameters of iS3D in order to determine whether a different
$\delta f$-correction could help fix the yields. All-in-all, however, our
model performed extremely well against a limited set of experimental data.

Our model’s physics-based approach allowed us to make predictions regarding
flow and bulk observables and interpret them phenomenologically. Furthermore,
we were able to use it to test different averaging techniques and determine
their effects on observables. In the case of 2-particle momentum correlations,
these two methods led to widely different values. This clear difference
presents itself as a golden opportunity to determine which types of
fluctuations and correlations (either local or global) dominate the $p_{T}$
spectrum in U+U collisions at $\sqrt{s_{NN}}=\qty{193}{\giga}$. This
prediction could also be of use in constraining 238U deformation parameters,
as momentum correlators are sensitive to nuclear deformation [85]. Another
compelling prediction of our model was that of a definite cross-over towards
anti-correlation of $v_{2}\\{2\\}$ and $\left\langle p_{T}\right\rangle$ in
central collisions, at about the $7\%$ mark. This is in stark contrast to
reported correlations of spherically symmetric nuclei [7], as well as results
obtained in other works using our model [21]. If such an anti-correlation is
detected in future experimental work, it would validate our phenomenological
reasoning as well as our model. Finally, our predictions regarding the
differential flow for charged hadrons and identified particles were consistent
across centrality classes when compared to data for Au+Au collisions at RHIC.
Indeed, the observables acted as expected given available integrated flow
observables.

Figure 14:  (top) 2- and (bottom) 3-particle momentum correlators as functions
of centrality, with $\qty{0.2}{\giga}\leq p_{T}\leq\qty{2.0}{\giga}$. Figure
15:  (top) Elliptic and (bottom) triangular flow and $\left\langle
p_{T}\right\rangle$ correlations as functions of centrality.

Our work here sets the stage for further calculations regarding deformed
nuclei, be it 238U specifically or others, such as 129Xe [86]. Indeed, while
our current work seems to suggest a good understanding of the structure of
238U, many details remain to be explored. As our models become better at
describing heavy-ion collisions, we must prove increasingly diligent regarding
the nuclear parametrizations we use. While our 238U parametrization was
satisfactory, more details are certainly required to represent physical
reality adequately; we think here of different and more recent nuclear
parametrization paradigms, such as Nuclear Density Functional Theory (NDFT)
[87], which provide deeper and more physical constraints on nucleon positions
and correlations. In the end, we would want to use sensitive observables
available experimentally to constrain nuclear deformation, allowing us to
further perfect our state-of-the-art simulations and provide a deeper
understanding of strongly interacting matter at various stages.

## Acknowledgments

We would like to acknowledge the support of the entirety of our research group
at McGill University. We also acknowledge insightful conversations with R.
Modarresi Yazdi, M. Heffernan, S. McDonald, S. Shi, B. Schenke C. Shen, J. Jia
and C. Zhang. This work was funded by the Natural Sciences and Engineering
Research Council of Canada (NSERC) [SAPIN-2018-00024 ; SAPIN-2020-00048].
Computations were made on the Béluga supercomputer system from McGill
University, managed by Calcul Québec (calculquebec.ca) and Digital Research
Alliance of Canada (alliancecan.ca). The operation of this supercomputer is
funded by the Canada Foundation for Innovation (CFI), Ministère de l’Économie,
des Sciences et de l’Innovation du Québec (MESI) and le Fonds de recherche du
Québec - Nature et technologies (FRQ-NT).

## References

  * Adamczyk  _et al._ [2015] L. Adamczyk, J. Adkins, G. Agakishiev, M. Aggarwal, Z. Ahammed, I. Alekseev, J. Alford, A. Aparin, D. Arkhipkin, E. Aschenauer, _et al._ , Azimuthal Anisotropy in U+U and Au+Au Collisions at RHIC, Physical Review Letters 115, 10.1103/physrevlett.115.222301 (2015).
  * Adam  _et al._ [2019] J. Adam, L. Adamczyk, _et al._ , Azimuthal harmonics in small and large collision systems at RHIC top energies, Physical Review Letters 122, 10.1103/physrevlett.122.172301 (2019).
  * Busza  _et al._ [2018] W. Busza, K. Rajagopal, and W. van der Schee, Heavy ion collisions: The big picture and the big questions, Annual Review of Nuclear and Particle Science 68, 339 (2018).
  * Pasechnik and Šumbera [2017] R. Pasechnik and M. Šumbera, Phenomenological review on quark–gluon plasma: Concepts vs. observations, Universe 3, 7 (2017).
  * Shuryak [2000] E. V. Shuryak, High energy collisions of strongly deformed nuclei: An old idea with a new twist, Phys. Rev. C 61, 034905 (2000).
  * Schenke  _et al._ [2020a] B. Schenke, C. Shen, and D. Teaney, Transverse momentum fluctuations and their correlation with elliptic flow in nuclear collisions, Physical Review C 102, 10.1103/physrevc.102.034905 (2020a).
  * Aad  _et al._ [2019] G. Aad, B. Abbott, D. Abbott, A. Abed Abud, K. Abeling, _et al._ , Measurement of flow harmonics correlations with mean transverse momentum in Pb-Pb and p-Pb collisions at 5.02 TeV with the ATLAS detector, The European Physical Journal C 79, 10.1140/epjc/s10052-019-7489-6 (2019).
  * Heinz and Kuhlman [2005] U. Heinz and A. Kuhlman, Anisotropic flow and jet quenching in ultrarelativistic U+U collisions, Physical Review Letters 94, 10.1103/physrevlett.94.132301 (2005).
  * Giacalone  _et al._ [2021a] G. Giacalone, J. Jia, and C. Zhang, Impact of nuclear deformation on relativistic heavy-ion collisions: Assessing consistency in nuclear physics across energy scales, Physical Review Letters 127, 242301 (2021a).
  * Jia [2022a] J. Jia, Shape of atomic nuclei in heavy ion collisions, Physical Review C 105, 014905 (2022a).
  * Bally  _et al._ [2022] B. Bally, M. Bender, G. Giacalone, and V. Somà, Evidence of the triaxial structure of xe 129 at the large hadron collider, Physical Review Letters 128, 082301 (2022).
  * Adam  _et al._ [2016] J. Adam _et al._ , Centrality dependence of the charged-particle multiplicity density at midrapidity in Pb-Pb collisions at 5.02 TeV, Physical Review Letters 116, 10.1103/physrevlett.116.222302 (2016).
  * Acharya  _et al._ [2020] S. Acharya _et al._ , Production of charged pions, kaons, and (anti-)protons in Pb-Pb and inelastic pp collisions at 5.02 TeV, Physical Review C 101, 10.1103/physrevc.101.044907 (2020).
  * Adare  _et al._ [2011] A. Adare _et al._ , Measurements of higher order flow harmonics in Au+Au collisions at 200 GeV, Physical Review Letters 107, 10.1103/physrevlett.107.252301 (2011).
  * Adams  _et al._ [2005] J. Adams _et al._ , Azimuthal anisotropy in Au+Au collisions at 200 GeV, Physical Review C 72, 10.1103/physrevc.72.014904 (2005).
  * Gelis  _et al._ [2010] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, The color glass condensate, Annual Review of Nuclear and Particle Science 60, 463 (2010).
  * Iancu and Venugopalan [2004] E. Iancu and R. Venugopalan, The Color Glass Condensate and high energy scattering in QCD, Quark-Gluon Plasma 3 , 249 (2004).
  * Gale  _et al._ [2013a] C. Gale, S. Jeon, B. Schenke, P. Tribedy, and R. Venugopalan, Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics, Phys. Rev. Lett. 110, 012302 (2013a), arXiv:1209.6330 [nucl-th] .
  * Gale  _et al._ [2013b] C. Gale, S. Jeon, and B. Schenke, Hydrodynamic Modeling of Heavy-Ion Collisions, Int. J. Mod. Phys. A28, 1340011 (2013b).
  * Schenke  _et al._ [2020b] B. Schenke, C. Shen, and P. Tribedy, Running the gamut of high energy nuclear collisions, Physical Review C 102, 10.1103/physrevc.102.044905 (2020b).
  * Heffernan  _et al._ [2023a] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Bayesian quantification of strongly-interacting matter with color glass condensate initial conditions (2023a), arXiv:2302.09478 [nucl-th] .
  * McDonald  _et al._ [2017] S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Hydrodynamic predictions for Pb+Pb collisions at 5.02 TeV, Physical Review C 95, 10.1103/physrevc.95.064913 (2017).
  * Heffernan  _et al._ [2023b] M. R. Heffernan, C. Gale, S. Jeon, and J.-F. Paquet, Early-times yang-mills dynamics and the characterization of strongly interacting matter with statistical learning (2023b), arXiv:2306.09619 [nucl-th] .
  * McDonald  _et al._ [2023] S. McDonald, S. Jeon, and C. Gale, The 3+1D initialization and evolution of the Glasma (2023), arXiv:2306.04896 [hep-ph] .
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, (3+1)D hydrodynamic simulation of relativistic heavy-ion collisions, Physical Review C 82, 10.1103/physrevc.82.014903 (2010).
  * Ryu  _et al._ [2015] S. Ryu, J.-F. Paquet, C. Shen, G.-S. Denicol, B. Schenke, S. Jeon, and C. Gale, Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions, Physical Review Letters 115, 10.1103/physrevlett.115.132301 (2015).
  * McNelis  _et al._ [2020] M. McNelis, D. Everett, and U. Heinz, Particlization in fluid dynamical simulations of heavy-ion collisions: The iS3D module (2020), arXiv:1912.08271 [nucl-th] .
  * Weil  _et al._ [2016] J. Weil, V. Steinberg, J. Staudenmaier, L. G. Pang, D. Oliinychenko, J. Mohs, M. Kretz, T. Kehrenberg, A. Goldschmidt, B. Bäuchle, J. Auvinen, M. Attems, and H. Petersen, Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions, Physical Review C 94, 10.1103/physrevc.94.054905 (2016).
  * Venugopalan and Prakash [1990] R. Venugopalan and M. Prakash, Flow effects on transverse momentum spectra in ultrarelativistic nuclear collisions, Phys. Rev. C 41, 221 (1990).
  * Sollfrank  _et al._ [1997] J. Sollfrank, P. Huovinen, M. Kataja, P. V. Ruuskanen, M. Prakash, and R. Venugopalan, Hydrodynamical description of 200-A/GeV/c S + Au collisions: Hadron and electromagnetic spectra, Phys. Rev. C 55, 392 (1997), arXiv:nucl-th/9607029 .
  * Kolb  _et al._ [1999] P. F. Kolb, J. Sollfrank, and U. W. Heinz, Anisotropic flow from AGS to LHC energies, Phys. Lett. B 459, 667 (1999), arXiv:nucl-th/9906003 .
  * Morita  _et al._ [2000] K. Morita, S. Muroya, H. Nakamura, and C. Nonaka, Numerical analysis of two pion correlation based on a hydrodynamical model, Phys. Rev. C 61, 034904 (2000), arXiv:nucl-th/9906037 .
  * Nonaka  _et al._ [2000] C. Nonaka, E. Honda, and S. Muroya, (3+1)-dimensional relativistic hydrodynamical expansion of hot and dense matter in ultrarelativistic nuclear collision, Eur. Phys. J. C 17, 663 (2000), arXiv:hep-ph/0007187 .
  * Kolb  _et al._ [2001] P. F. Kolb, P. Huovinen, U. W. Heinz, and H. Heiselberg, Elliptic flow at SPS and RHIC: From kinetic transport to hydrodynamics, Phys. Lett. B 500, 232 (2001), arXiv:hep-ph/0012137 .
  * Kolb  _et al._ [2000] P. F. Kolb, J. Sollfrank, and U. W. Heinz, Anisotropic transverse flow and the quark hadron phase transition, Phys. Rev. C 62, 054909 (2000), arXiv:hep-ph/0006129 .
  * Huovinen  _et al._ [2001] P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen, and S. A. Voloshin, Radial and elliptic flow at RHIC: Further predictions, Phys. Lett. B 503, 58 (2001), arXiv:hep-ph/0101136 .
  * Hirano  _et al._ [2002] T. Hirano, K. Morita, S. Muroya, and C. Nonaka, Hydrodynamical analysis of hadronic spectra in the 130 GeV/nucleon Au+Au collisions, Phys. Rev. C 65, 061902 (2002), arXiv:nucl-th/0110009 .
  * Morita  _et al._ [2002] K. Morita, S. Muroya, C. Nonaka, and T. Hirano, Comparison of space-time evolutions of hot / dense matter in $\sqrt{s_{NN}}=17$ GeV and $130$-GeV relativistic heavy ion collisions based on a hydrodynamical model, Phys. Rev. C 66, 054904 (2002), arXiv:nucl-th/0205040 .
  * Hirano and Tsuda [2002] T. Hirano and K. Tsuda, Collective flow and two pion correlations from a relativistic hydrodynamic model with early chemical freezeout, Phys. Rev. C 66, 054905 (2002), arXiv:nucl-th/0205043 .
  * Huovinen and Ruuskanen [2006] P. Huovinen and P. V. Ruuskanen, Hydrodynamic Models for Heavy Ion Collisions, Ann. Rev. Nucl. Part. Sci. 56, 163 (2006), arXiv:nucl-th/0605008 .
  * Miller  _et al._ [2007] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Glauber modeling in high-energy nuclear collisions, Annual Review of Nuclear and Particle Science 57, 205 (2007).
  * Schenke  _et al._ [2012] B. Schenke, P. Tribedy, and R. Venugopalan, Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions, Physical Review C 86, 10.1103/physrevc.86.034908 (2012).
  * McLerran and Venugopalan [1994a] L. McLerran and R. Venugopalan, Computing quark and gluon distribution functions for very large nuclei, Physical Review D 49, 2233 (1994a).
  * McLerran and Venugopalan [1994b] L. McLerran and R. Venugopalan, Gluon distribution functions for very large nuclei at small transverse momentum, Phys. Rev. D 49, 3352 (1994b).
  * Iancu  _et al._ [2002] E. Iancu, A. Leonidov, and L. McLerran, The colour glass condensate: An introduction, QCD Perspectives on Hot and Dense Matter , 73 (2002).
  * Gyulassy and McLerran [1997] M. Gyulassy and L. D. McLerran, Yang-Mills radiation in ultrarelativistic nuclear collisions, Phys. Rev. C 56, 2219 (1997), arXiv:nucl-th/9704034 .
  * Kovchegov and Rischke [1997] Y. V. Kovchegov and D. H. Rischke, Classical gluon radiation in ultrarelativistic nucleus-nucleus collisions, Phys. Rev. C 56, 1084 (1997), arXiv:hep-ph/9704201 .
  * Krasnitz and Venugopalan [1999] A. Krasnitz and R. Venugopalan, Non-perturbative computation of gluon mini-jet production in nuclear collisions at very high energies, Nuclear Physics B 557, 237 (1999).
  * Krasnitz  _et al._ [2003] A. Krasnitz, Y. Nara, and R. Venugopalan, Gluon production in the color glass condensate model of collisions of ultrarelativistic finite nuclei, Nuclear Physics A 717, 268 (2003).
  * Dumitru and McLerran [2002] A. Dumitru and L. D. McLerran, How protons shatter colored glass, Nucl. Phys. A 700, 492 (2002), arXiv:hep-ph/0105268 .
  * De Vries  _et al._ [1987] H. De Vries, C. W. De Jager, and C. De Vries, Nuclear charge-density-distribution parameters from elastic electron scattering, At. Data Nucl. Data Tables; (United States) 36, 10.1016/0092-640X(87)90013-1 (1987).
  * Alvioli  _et al._ [2009] M. Alvioli, H. Drescher, and M. Strikman, A monte carlo generator of nucleon configurations in complex nuclei including nucleon-nucleon correlations, Physics Letters B 680, 225 (2009).
  * Dalal and MacGregor [2022] R. Dalal and I. J. D. MacGregor, Nucleon-nucleon correlations inside atomic nuclei: synergies, observations and theoretical models (2022), arXiv:2210.06114 [nucl-th] .
  * Kowalski and Teaney [2003] H. Kowalski and D. Teaney, Impact parameter dipole saturation model, Physical Review D 68, 10.1103/physrevd.68.114005 (2003).
  * Rezaeian  _et al._ [2013] A. H. Rezaeian, M. Siddikov, M. Van de Klundert, and R. Venugopalan, Analysis of combined HERA data in the impact-parameter dependent saturation model, Phys. Rev. D 87, 034002 (2013).
  * Altarelli and Parisi [1977] G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126, 298 (1977).
  * Dokshitzer [1977] Y. L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics., Sov. Phys. JETP 46, 641 (1977).
  * [58] F. Gelis, Dglap equations solver.
  * Mueller [1990] A. H. Mueller, Small x Behavior and Parton Saturation: A QCD Model, Nucl. Phys. B 335, 115 (1990).
  * Kovner  _et al._ [1995a] A. Kovner, L. McLerran, and H. Weigert, Gluon production at high transverse momentum in the mclerran-venugopalan model of nuclear structure functions, Physical Review D 52, 3809 (1995a).
  * Kovner  _et al._ [1995b] A. Kovner, L. McLerran, and H. Weigert, Gluon production from non-abelian weizsäcker-williams fields in nucleus-nucleus collisions, Phys. Rev. D 52, 6231 (1995b).
  * Blaschke  _et al._ [2016] D. N. Blaschke, F. Gieres, M. Reboud, and M. Schweda, The energy–momentum tensor(s) in classical gauge theories, Nuclear Physics B 912, 192 (2016), mathematical Foundations of Quantum Field Theory: A volume dedicated to the Memory of Raymond Stora.
  * Bazavov  _et al._ [2018] A. Bazavov, P. Petreczky, and J. Weber, Equation of state in 2+1 flavor QCD at high temperatures, Physical Review D 97, 10.1103/physrevd.97.014510 (2018).
  * Kolb and Heinz [2003] P. F. Kolb and U. W. Heinz, Hydrodynamic description of ultrarelativistic heavy ion collisions, Quark Gluon Plasma 3 , 634 (2003), arXiv:nucl-th/0305084 .
  * Jeon and Heinz [2015] S. Jeon and U. Heinz, Introduction to hydrodynamics, International Journal of Modern Physics E 24, 1530010 (2015).
  * Schenke  _et al._ [2011] B. Schenke, S. Jeon, and C. Gale, Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics, Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Denicol  _et al._ [2014] G. S. Denicol, S. Jeon, and C. Gale, Transport coefficients of bulk viscous pressure in the 14-moment approximation, Physical Review C 90, 10.1103/physrevc.90.024912 (2014).
  * Pandav  _et al._ [2022] A. Pandav, D. Mallick, and B. Mohanty, Search for the QCD critical point in high energy nuclear collisions, Progress in Particle and Nuclear Physics 125, 103960 (2022).
  * Yin [2018] Y. Yin, The QCD critical point hunt: emergent new ideas and new dynamics (2018), arXiv:1811.06519 [nucl-th] .
  * Luo and Xu [2017] X. Luo and N. Xu, Search for the QCD critical point with fluctuations of conserved quantities in relativistic heavy-ion collisions at RHIC: an overview, Nuclear Science and Techniques 28, 10.1007/s41365-017-0257-0 (2017).
  * Cooper and Frye [1974] F. Cooper and G. Frye, Single-particle distribution in the hydrodynamic and statistical thermodynamic models of multiparticle production, Phys. Rev. D 10, 186 (1974).
  * Grad [1949] H. Grad, On the kinetic theory of rarefied gases, Communications on Pure and Applied Mathematics 2, 331 (1949).
  * Monnai and Hirano [2009] A. Monnai and T. Hirano, Effects of bulk viscosity at freezeout, Phys. Rev. C 80, 054906 (2009).
  * Israel and Stewart [1979] W. Israel and J. Stewart, Transient relativistic thermodynamics and kinetic theory, Annals of Physics 118, 341 (1979).
  * Schenke  _et al._ [2014] B. Schenke, P. Tribedy, and R. Venugopalan, Initial-state geometry and fluctuations in Au + Au, Cu + Au, and U + U collisions at energies available at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 89, 064908 (2014).
  * Liu  _et al._ [2022] L.-M. Liu, C.-J. Zhang, J. Xu, J. Jia, and G.-X. Peng, Free spectator nucleons in ultracentral relativistic heavy-ion collisions as a probe of neutron skin, Physical Review C 106, 10.1103/physrevc.106.034913 (2022).
  * Giacalone  _et al._ [2023] G. Giacalone, G. Nijs, and W. van der Schee, Determination of the neutron skin of 208Pb from ultrarelativistic nuclear collisions (2023), arXiv:2305.00015 [nucl-th] .
  * Borghini  _et al._ [2001] N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, New method for measuring azimuthal distributions in nucleus-nucleus collisions, Physical Review C 63, 10.1103/physrevc.63.054906 (2001).
  * Bilandzic  _et al._ [2011] A. Bilandzic, R. Snellings, and S. Voloshin, Flow analysis with cumulants: Direct calculations, Physical Review C 83, 10.1103/physrevc.83.044913 (2011).
  * Giacalone  _et al._ [2021b] G. Giacalone, F. G. Gardim, J. Noronha-Hostler, and J.-Y. Ollitrault, Skewness of mean transverse momentum fluctuations in heavy-ion collisions, Physical Review C 103, 10.1103/physrevc.103.024910 (2021b).
  * Bożek [2016] P. Bożek, Transverse-momentum–flow correlations in relativistic heavy-ion collisions, Physical Review C 93, 10.1103/physrevc.93.044908 (2016).
  * Abdallah  _et al._ [2023] M. Abdallah, B. Aboona, J. Adam, _et al._ , Pion, kaon, and (anti)proton production in U+U collisions at 193 GeV measured with the STAR detector, Physical Review C 107, 10.1103/physrevc.107.024901 (2023).
  * Pratt  _et al._ [2015] S. Pratt, E. Sangaline, P. Sorensen, and H. Wang, Constraining the equation of state of superhadronic matter from heavy-ion collisions, Physical Review Letters 114, 10.1103/physrevlett.114.202301 (2015).
  * Kumar and Toll [2022] A. Kumar and T. Toll, Investigating the structure of gluon fluctuations in the proton with incoherent diffraction at HERA, The European Physical Journal C 82, 10.1140/epjc/s10052-022-10774-3 (2022).
  * Jia [2022b] J. Jia, Probing triaxial deformation of atomic nuclei in high-energy heavy ion collisions, Physical Review C 105, 10.1103/physrevc.105.044905 (2022b).
  * Acharya  _et al._ [2019] S. Acharya, F. Tolares-Acosta, D. Adamová, J. Adolfsson, _et al._ , Centrality and pseudorapidity dependence of the charged-particle multiplicity density in xe-xe collisions at 5.44 tev, Physics Letters B 790, 35 (2019).
  * Colò [2020] G. Colò, Nuclear density functional theory, Advances in Physics: X 5, 1740061 (2020), https://doi.org/10.1080/23746149.2020.1740061 .