UNIVERSIDADE ESTADUAL DE CAMPINAS
Instituto de Física “Gleb Wataghin”

João Paulo Picchetti

Estudando o impacto do tamanho do nucleon em colisões de íons pesados
relativísticos

Studying the impact of the nucleon size in relativistic heavy-ion collisions

Campinas

2022

João Paulo Picchetti

Studying the impact of the nucleon size in relativistic heavy-ion collisions

Estudando o impacto do tamanho do nucleon em colisões de íons pesados
relativísticos

Dissertação apresentada ao Instituto de Física

“Gleb Wataghin” da Universidade Estadual de Campinas

como parte dos requisitos exigidos para obtenção do

título de Mestre em física, na área de física.

Dissertation presented to the “Gleb Wataghin”

Institute of Physics of the University of Campinas

in partial fulfillment of the requeriments for the

degree of Master in Physics, in the area of Physics.

Supervisor: Jun Takahashi

ESTE TRABALHO CORRESPONDE À VERSÃO FINAL
DA DISSERTAÇÃO DEFENDIDA PELO ALUNO
JOÃO PAULO PICCHETTI, E ORIENTADO
PELO PROF. DR. JUN TAKAHASHI.

Campinas

2022

### Abstract

Under the extreme conditions of temperature generated in relativistic heavy-
ion collisions, a fascinating fluid-like state of matter where quarks and
gluons are no longer confined is formed, the Quark Gluon Plasma (QGP). The
most modern computational approaches are multi-stage (hybrid) simulations, in
which different models are used in a chain structure, each one dedicated to
the description of a specific stage of the collision.

The hydrodynamic stage of the simulation requires an energy density profile of
the system as an initial condition. In the process of converting the two
colliding nuclei in such an energy distribution, some specification about the
nucleon size inevitably has to be made. Nucleons are usually modeled as
bidimensional Gaussians, and the Gaussian width (the nucleon-width) is a free
parameter of the simulation. A best-fit value of the nucleon-width can be
inferred by Bayesian Analyses, where the model is confronted with experimental
data.

Some of the most recent analyses have obtained surprisingly large values for
the nucleon width parameter, exceeding in over 50 % the current value for the
proton charged radius. This motivates the development of a better
understanding of the role played by this parameter inside the simulation.

In this work, we perform simulations of relativistic heavy-ion collisions
using a state-of-the-art hybrid simulation chain, using three different values
of the nucleon width inside the initial condition generator
T${}_{\text{R}}$ENTo, and systematically investigate its effects on the
initial condition characteristics and observables. The nucleon-width strongly
affects the eccentricity harmonics and the gradients in the initial condition.
The mean $p_{T}$ of particles in the simulation using $w$ = 0.5 fm is much
larger than experimental data. We associate this to the combination of
stronger gradients in the initial condition and the coupling of a conformal
pre-equilibrium dynamics to the hydrodynamic simulation.

Keywords: Heavy-ion phenomenology, Heavy-ion collisions, High-energy nuclear
physics.

###### Contents

  1. 1 Introduction
    1. 1.1 The Standard Model
    2. 1.2 Quantum Chromodynamics (QCD)
    3. 1.3 The QCD Phase Diagram
  2. 2 Relativistic heavy-ion collisions
    1. 2.1 Historical overview
    2. 2.2 Space-time evolution of relativistic heavy-ion collisions
    3. 2.3 Basic kinematic variables, coordinate system and centrality
    4. 2.4 Experimental observables and evidences of QGP formation
      1. 2.4.1 Multiplicity
      2. 2.4.2 Transverse momentum distributions
      3. 2.4.3 Anisotropic flow
    5. 2.5 The Initial Condition
      1. 2.5.1 The Glauber Model
      2. 2.5.2 Eccentricity harmonics
  3. 3 Simulating heavy-ion collisions
    1. 3.1 The initial condition: T${}_{\text{R}}$ENTo
    2. 3.2 Pre-equilibrium dynamics: K$\o{}$MP$\o{}$ST
    3. 3.3 Hydrodynamic evolution: MUSIC
    4. 3.4 Particlization: iSS
    5. 3.5 Hadronic phase: UrQMD
  4. 4 Motivation
    1. 4.1 The nucleon size
    2. 4.2 Constraining the hybrid simulations: Bayesian Analysis
    3. 4.3 The nucleon-width parameter
  5. 5 Results
    1. 5.1 Characterizing the initial condition
      1. 5.1.1 Participants and binary collisions
      2. 5.1.2 Eccentricity harmonics
    2. 5.2 Final state observables
      1. 5.2.1 Charged particle multiplicity density at mid-rapidity
      2. 5.2.2 Mean transverse momentum
      3. 5.2.3 Anisotropic flow
      4. 5.2.4 Blast Wave analysis
  6. 6 Conclusions

##  Introduction

While at Earth-like conditions quarks are confined inside the protons and
neutrons, which are in their turn bound together in the atomic nuclei, at
sufficiently high energies a phase transition occurs to a fluid-like state of
matter where quarks and gluons are no longer confined, the Quark Gluon Plasma
(QGP).

This chapter is structured as an introduction to high-energy physics, focusing
on the topics which are relevant for the physics of relativistic heavy-ion
collisions and for the development of this work. Section 1.1 is an
introduction to the Standard Model of particle physics, while Section 1.2
provides an overview of the most important aspects of Quantum Chromodynamics
(QCD), the theory of the strong interaction. Section 1.3 is dedicated to a
brief discussion about the QCD phase diagram and phase transitions in QCD.

###  1.1 The Standard Model

Particle physics studies the fundamental building blocks of the Universe and
the interactions between them. While the gravitational interaction is the
subject of Einstein’s Theory of General Relativity, the other three
fundamental forces can be studied in a single unified framework, which is the
Standard Model of Particle Physics (SM) [1]. These are the electromagnetic
force, which is the only long range SM interaction (and therefore, the only
one that we are able to actually experience daily), and also the strong and
weak forces, which act on a short range and are relevant only in the scale of
quantum mechanics. The strong force is responsible for holding the protons and
neutrons together in atomic nuclei, while the weak force is related to
processes such as the $\beta$-decay.

Figure 1.1: Standard Model particles. Figure taken from [2].

The formal description of particles and their interactions is the object of
Quantum Field Theory (QFT) [3], which combines quantum mechanics, special
relativity and classical field theory. In QFT, particles are viewed as excited
states of quantum fields, and the interactions are mediated by the exchange of
virtual bosons. Historically, QFT was first developed as a quantum theory of
the electromagnetic interaction [4], and the term Quantum Electrodynamics
(QED) was first used by Paul Dirac in 1927.

In the SM, matter is made of twelve elementary spin-1/2 fermions which
interact through the fundamental forces. To interact via one of the forces, a
particle must posses the respective interaction charge. All the elementary
fermions interact through the weak force, while particles with electric charge
participate in the electromagnetic interaction. The charge of the strong
interaction is called color charge. The twelve fundamental fermions are
separated in two groups of six based on whether they have color charge or not.
The six fermions which are color charged are called quarks, while the ones
which are not are known as leptons.

The elementary fermions interact by exchanging spin-1 bosons (the gauge
bosons). For this reason, the gauge bosons are usually referred to as the
force carriers (or mediators) of the SM. The photon ($\gamma$) is the force
carrier of the electromagnetic interaction, the gluon ($g$) is the mediator of
the strong interaction, and the weak interaction is mediated by the neutral
$Z$ boson and also by the charged $W^{\pm}$ bosons. In addition, all massive
SM particles are coupled to the Higgs field. The Higgs boson was theoretically
proposed by Peter Higgs in 1964 [5] and experimentally discovered in 2012 both
by the ATLAS [6] and CMS [7, 8] Collaborations. Table 1.1 shows a schematic
view of which of the elementary fermions participate in each interaction.

| Strong | Electromagnetic | Weak
---|---|---|---
Quarks (u, d, c, s, t, b) | ✓ | ✓ | ✓
Charged leptons ($e^{-}$, $\mu^{-}$, $\tau^{-}$) |  | ✓ | ✓
Neutrinos ($\nu_{e}$, $\nu_{\mu}$, $\nu_{\tau}$) |  |  | ✓
Table 1.1: Standard Model particles and interactions in which they participate

There is a mass hierarchy between the elementary fermions in the SM, that
separates them in three generations. The first generation consists on the u
and d quarks, which form the protons (uud) and neutrons (udd), and also the
electron ($e^{-}$) and the electron neutrino ($\nu_{e}$). These particles
dominate the low-energy world in which we live in and, with the exception of
the electron neutrino, are the building blocks of atoms. Being the
constituents of atomic nuclei, protons and neutrons are also known as
nucleons. For each particle in the first generation, there is a heavier one,
which differs from its first generation partner only by the mass. The second
generation is formed by the c and s quarks, the muon ($\mu^{-}$) and the muon
neutrino ($\nu_{\mu}$), while the third generation is composed by the t and b
quarks, the tau ($\tau^{-}$) and the tau neutrino ($\nu_{\tau}$). Each
elementary fermion has also a corresponding anti-particle, which is identical
except for the electric charge, which has the opposite sign.

###  1.2 Quantum Chromodynamics (QCD)

Quantum Chromodynamics (QCD) is the QFT that studies the strong interaction
[9, 10]. In opposition to the electric charge, which comes in two types
(positive and negative), color charge comes in six types: red, green, blue,
and the corresponding anti-colors. In QCD, quarks interact by exchanging
gluons, which are themselves color charged and also interact with each other
through the same mechanism. This is an expression of the fact that QCD is a
non-Abelian field theory, which means that the generators of the theory’s
symmetry group, in the case SU(3), don’t commute. These gluon-gluon
interactions appear in Feynman diagrams as three-gluon and four-gluon
vertices, for which there are no analogues in QED. The fact that the force
carriers of QCD are themselves color charged is one the things that make it so
distinct from the electroweak part of the SM, and also so complicated. There
are $N_{c}^{2}-1=8$ “types” of gluons (where $N_{c}=3$ is the number of
colors), meaning eight possible color states a gluon may have, known as the
“color octet” [1].

Figure 1.2: The basic QCD interaction vertices: the quark-gluon vertex, and
the three-gluon and four-gluon vertices, for which there are no analogues in
QED. Figure taken from from [1].

A remarkable feature of QCD is that the coupling constant of the theory
($\alpha_{s}$) is strongly dependent on the energy scale of the interaction
or, in other words, the four-momentum of the exchanged gluon ($Q$). The energy
dependence of the coupling constant is dictated by the $\beta$-function [11]:

$Q^{2}\frac{\partial}{\partial
Q^{2}}\frac{\alpha_{s}}{4\pi}=\beta(\alpha_{s}),$ (1.1)

which can be expressed as a perturbative series. The first term of the series
is:

$\beta_{0}=11-\frac{2}{3}n_{f},$ (1.2)

where $n_{f}$ is the number of quark flavors considered at the energy scale in
question. Then, at $\beta_{0}$ order, there is an analytic solution to
Equation (1.1):

$\alpha_{s}(Q^{2})=\frac{4\pi}{\beta_{0}\ln{(Q^{2}/\Lambda_{QCD}^{2})}},$
(1.3)

where $\Lambda_{\text{QCD}}$ = 220 MeV. Figure 1.3 shows different
measurements of the strong coupling constant $\alpha_{s}(Q^{2})$ as a function
of $Q$. The coupling constant is larger for low $Q^{2}$ and decreases
significantly with increasing energy.

Figure 1.3: The QCD coupling constant $\alpha_{s}(Q^{2})$ as a function of $Q$
from different experiments. The order of perturbation theory used in the
extraction of $\alpha_{s}$ is also indicated. From [12].

Due to this behaviour of the coupling constant, the characteristic energy
scale $\Lambda_{\text{QCD}}$ roughly separates QCD in two regimes. At higher
energies, the coupling constant is small, and perturbation theory can be
applied (perturbative QCD), while at lower energies perturbation theory is not
applicable (non-perturbative QCD). In general, interactions with low momentum
transfer are called soft, while interactions whit high momentum transfer are
referred to as hard interactions.

This strong energy dependence of the coupling constant (known as running of
$\alpha_{s}$) leads to two main properties of QCD:

  * •

Color confinement: color charged particles are never observed as free
particles. Quarks and gluons are always observed as confined bound states of
three quarks (baryons) or a quark and an anti-quark (mesons), which are the
hadrons. In nature, hadrons are always found in a color singlet state, and it
is not possible to separate them into their colored constituents. Although
color confinement is qualitatively well established, an analytical
demonstration of it is still an open problem.

  * •

Asymptotic freedom: at very high energies, quarks and gluons interact weakly.
This means that inside the hadrons, when they are very close to each other,
quarks coexist as almost free particles. Asymptotic freedom was discovered by
David Gross and Frank Wilczek [13], and independently by David Politzer [14]
in 1973, which awarded them the Nobel Prize in 2004.

It seems reasonable to say that QCD began with the invention of the quark
model in 1964 by Murray Gell-Mann [15], with the proposition that baryons and
mesons were not elementary, but rather bound states of elementary particles,
precisely the quarks. In 1969, Richard Feynman introduced the ideia that at
very high energies, hadrons are not only formed by its valence quarks, but
contain a collection of gluons and virtual pairs of quarks, which came to be
known as partons [16]. Later on in 1973 Gell-Mann, Harald Fritzsch and
Heinrich Leutwyler introduced color as the source of the “strong fields” [17].

###  1.3 The QCD Phase Diagram

While at low temperatures and densities quarks are confined inside hadrons, at
very high energies or densities, lattice calculations QCD predict the
existence of a phase transition to a state where quarks and gluons are
deconfined [18, 19]. This state of matter, where quarks and gluons themselves
are the degrees of freedom, is called Quark Gluon Plasma (QGP).

The idea of a QCD phase diagram was first proposed by N. Cabbibo and G. Parisi
in 1975 [20]. It consists of a diagram that shows in one axis the baryonic
density (or alternatively, the baryonic chemical potential $\mu_{B}$) and in
the other, the temperature. Figure 1.4 shows a more modern view of the QCD
phase diagram.

Figure 1.4: The QCD Phase Diagram: baryonic density is shown in the horizontal
axis, and temperature in the vertical axis. The orange band represents the
phase transition region from hadronic matter to QGP. Figure taken from [21].

It is known since 2006 from lattice calculations that the phase transition
between hadronic matter and QGP in the region of high temperature and low
baryonic density is a smooth crossover [18]. On the other hand, in the region
of low temperature and non-zero baryonic potential it is believed to exist a
first order phase transition from hadronic matter to QGP, where the first
derivatives of the thermodynamic fields are discontinuous [19]. An active
topic of research is the search for a critical end point (CEP) above the phase
transition line, where the phase transition changes its order, which is the
objective of the Beam Energy Scam (BES) program at RHIC, for example [22].

Figure 1.5 shows pressure, energy density and entropy density as a function of
temperature (results from lattice QCD at zero baryonic chemical potential).
There is an abrupt (yet smooth) increase in the number of degrees of freedom
of the system around the temperature of 150 MeV, which corresponds to the
crossover phase transition between nuclear matter and quark matter.

Figure 1.5: Pressure, energy density and entropy density as a function of
temperature, obtained from lattice calculations. Results obtained from a
hadron resonance gas (HRG) model calculation are shown as solid lines. Figure
obtained from [23].

Lattice QCD calculations are more accurately performed at $\mu_{B}$ = 0. At
non-vanishing baryonic density, lattice calculations are not so accurate,
which accounts for the broadening of the orange band in Figure 1.4 near the
horizontal axis of the phase diagram. The region of the diagram of (almost)
vanishing baryonic chemical potential and high temperature (above the phase
transition temperature) can be accessed in heavy ion collisions experiments.
These were the conditions believed to exist in the whole Universe in the first
moments after the Big Bang. On the other hand, the conditions of low
temperature and high baryonic potential are believed to exist in the interior
of very dense neutron stars.

##  Relativistic heavy-ion collisions

The extreme temperatures necessary for the QCD phase transition to occur can
be achieved in relativistic heavy-ion collisions. Chapter 2 aims to provide an
introduction to the field of relativistic heavy-ion collisions. Section 2.1
briefly discusses the history of relativistic heavy-ion collisions
experiments. Section 2.2 gives an overview of the space-time evolution of a
heavy-ion collision, and its main stages. Section 2.3 presents the main
kinematic variables and the coordinate system used in heavy-ion experiments,
and introduces most of the basic vocabulary which will be used in the rest of
this work. Section 2.4 is devoted the the main experimental observables used
to characterize the QGP, and evidences of the QGP formation. Section 2.5 talks
about the modelling of the initial state of the collision, and introduces
concepts that will be necessary for what follows.

###  2.1 Historical overview

The field of relativistic heavy-ion collisions is strongly interdisciplinary.
It makes connections between nuclear physics, which relies mostly on effective
models, and particle physics, which can be studied from first principles under
the formalism of QFT. As a system of thousands of particles is formed,
thermodynamic quantities such as temperature, pressure and entropy are used in
the theoretical description of heavy-ion collisions, providing also a
connection to the world of statistical physics. Moreover, the conditions
created in heavy-ion accelerators were believed to exist in the early stages
of the Universe, making the field of heavy-ion collisions also relevant for
cosmology [24].

Nuclear physics went through a revolution in the 1970s and early 1980s, as the
first particle accelerators were adapted to accelerate heavy nuclei. The most
famous example is probably the coupling of the SuperHilac linear accelerator
to the Bevatron in Berkeley, forming the Bevalac [25]. At the same time, the
Dubna Syncrophasotron was also modified to accelerate heavy ions [26].

Collisions with energy exceeding 10 GeV per nucleon were first made at the
Alternating Gradient Synchrotron (AGS) at the Brookhaven National Laboratory
(BNL) in the year of 1986, and in 1987 at the Super Proton Synchrotron (SPS)
at the European Organization for Nuclear Research (CERN). Later on, an
experimental breakthrough took place with the construction of the Relativistic
Heavy Ion (RHIC) [27], designed to accelerate gold nuclei at the center-of-
mass energy of $\sqrt{s_{\text{NN}}}$ = 200 GeV. After a first run in the year
of 2000, where the maximum energy of $\sqrt{s_{\text{NN}}}$ = 130 GeV was
achieved, three runs were performed at the full capacity of
$\sqrt{s_{\text{NN}}}$ = 200 GeV between 2001 and 2004. Contradicting the
expectation of the community that the deconfined state formed in such
collisions would be a weakly interacting system of quarks and gluons (similar
to a gas), RHIC data showed that it was indeed much more similar to a strongly
interacting fluid. The findings of the first runs at RHIC were summarized in
the famous “White Papers” [28, 29, 30, 31], published by each of RHIC’s
collaborations (PHOBOS, BRAHMS, STAR and PHENIX). Together, these four
articles have over 11.500 citations.

Between 1998 and 2008 the Large Hadron Collider (LHC) was built by CERN near
the border of France and Switzerland. With a 27 kilometers circumference, it
is the most energetic particle collider built to this day, achieving the
center-of-mass energy of $\sqrt{s_{\text{NN}}}$ = 5.02 TeV in Pb-Pb collisions
[32], and reaching the maximum energy of $\sqrt{s}$ = 13 TeV in p-p collisions
[33].

###  2.2 Space-time evolution of relativistic heavy-ion collisions

After being accelerated almost at the speed of light in circular trajectories
and in opposite directions, two highly relativistic nuclei collide inside the
accelerator beam pipe. Much less than a second later, an enormity of particles
is detected, containing many hadron species, leptons and photons. What
happened between the collision and the observation of this final-state
particles? A heavy-ion collision is a complicated process, which can be
divided in a few main stages:

  * •

Pre-equilibrium: in the first moments of the collision, just after the two
Lorentz contracted nuclei overlap and quarks and gluons become deconfined, the
system is far from equilibrium. At this point, the energy density across all
of the system is much larger than 500 $\text{MeV}/\text{fm}^{3}$, the typical
energy density inside hadrons [34]. Most of the initial interactions between
partons are soft, but rare scatterings with very large momentum transfer also
occur, leading to the formation of jets [35, 36]. In this initial moments
after the collision, the system expands violently, close to the speed of
light.

  * •

Hydrodynamic evolution: approximately 1 fm/c after the two nuclei collided,
the system has been driven to local equilibrium. At this point, the QGP
behaves as a strongly interacting, almost ideal fluid, with an extremely low
shear viscosity to entropy ratio of $\eta/s\approx 1/4\pi$ [37]. This is the
hydrodynamic stage of the collision. In general, hydrodynamics is an effective
theory based on conservation laws which describes the dynamics of a system
over long times and long distances, and is applicable when the microscopic
distance between particles ($\ell$) is much smaller than the macroscopic
length scale of the system ($L$).

  * •

Particlization: about 10 fm/c after the beginning of the hydrodynamic
evolution, the QGP has expanded and consequently, cooled down. At some point,
the energy density of the system is small enough so that quarks and gluons
recombine into hadrons.

  * •

Hadronic phase: in the final stage of the collision, the system is in a hadron
resonance gas (HRG) phase. The final state hadrons, which are mainly pions
(the lightest hadron), are boosted in the radial direction by the expansion of
the system, and propagate until they reach the detectors. In this process,
resonances might decay into more stable particles and the particles interact
with each other, both elastically and inelastically. At some point, the decays
and inelastic scatterings cease, and each particle species abundance stays
approximately stable. This stage of the hadronic phase is known as chemical
freeze-out. Later, as the system continues to expand and becomes more dilute,
the elastic scatterings cease, and the final momentum distribution of the
particles is fixed, which defines the kinetic freeze-out. After kinetic
freeze-out, the final state hadrons propagate as free particles until they
reach the detectors of the experiment.

Figure 2.1: Space-time evolution of a heavy ion collision and its stages:
initial scatterings and thermalization are followed by the hydrodynamic stage,
where the QGP expands as a fluid. At some point, quarks and gluons recombine
into hadrons (hadronization) which, after chemical and kinetic freeze-out, are
observed.

###  2.3 Basic kinematic variables, coordinate system and centrality

The coordinate system usually employed in studies of relativistic heavy-ion
collisions is shown in Figure 2.2, on top of a sketch of a typical collider
type experiment, in this case, the ALICE experiment at the LHC. The $z$-axis
lies is the direction of the beam. In symmetric collisions, the two colliding
nuclei travel at the same speed in the positive and negative $z$ directions,
so that it is specially convenient to perform calculations in the center-of-
mass frame (which, in this case, coincides with the laboratory frame). The
$x$-axis is oriented parallel to the ground, an the $y$-axis is placed along
the vertical direction. In analogy with spherical coordinates, a polar angle
$\theta$ and an azimuthal angle $\phi$ are defined.

Figure 2.2: Coordinate system used in this work, shown on top of a sketch of
the ALICE experiment. Figure obtained from [38].

The QGP itself is extremely short-lived and never observed, so it is necessary
to extract information about its characteristics solely based on the particles
that reach the experiment’s detector. A few kinematic variables are usually
employed to characterize the final state particles. Transverse momentum
($p_{T}$) is the projection of the particle’s momentum in the transverse plane
of the collision:

$p_{T}=\sqrt{p_{x}^{2}+p_{y}^{2}},$ (2.1)

High-$p_{T}$ particles in the final state emerge from the QGP travelling fast
and almost perpendicular to the beam direction. The rapidity variable ($y$)
can be used to characterize the particle’s momentum in the longitudinal
direction:

$y=\frac{1}{2}\ln{\left(\frac{E+p_{z}}{E-p_{z}}\right)}$ (2.2)

Rapidity is close to zero when a particle has little longitudinal momentum,
that is, when it is moving approximately perpendicular to the direction of the
beam. Alternatively, when a particle has large $p_{z}$, then
$|y|\rightarrow\infty$. Rapidity has the convenient property of being additive
under longitudinal boosts. Alternatively, the pseudo-rapidity variable
($\eta$) can be used:

$\eta=\frac{1}{2}\ln{\left(\frac{|\vec{p}|+p_{z}}{|\vec{p}|-p_{z}}\right)}=-\ln{\left(\tan{\frac{\theta}{2}}\right)}$
(2.3)

Pseudo-rapidity is zero perpendicular to the beam axis ($\theta=90^{\circ}$)
and $|\eta|\rightarrow\infty$ close to the beam axis ($\theta=0^{\circ}$ or
$\theta=180^{\circ}$). Pseudo-rapidity is specially useful because it is only
necessary to know the particle’s angle of emission $\theta$, while more
details about the particle, such as mass and momentum, are necessary to know
its rapidity. For this reason, rapidity is used when dealing with identified
particle species, and pseudo-rapidity is used when dealing with charged
particles. In the relativistic limit, where $E\approx|\vec{p}|$, the rapidity
and pseudo-rapidity coincide: $y\approx\eta$. Experimental measurements are
usually performed in a particular region of pseudo-rapidity or rapidity, as
$|y|<0.5$ or $|\eta|<0.8$, for example.

Being extended objects, the two nuclei might collide “head-on”, with their
centers aligned, or they might collide in a way such that their centers are
dislocated and the overlap between them is only partial. The variable that
quantifies the amount of overlap between the projectiles is the impact
parameter ($b$), defined as the distance between the centers of the two
colliding nuclei. In practice, the impact parameter is not known for each
collision, and it is necessary to use another quantity, analogous to the
impact parameter, to classify collisions. For this purpose, the concept of
centrality is used.

Figure 2.3: Illustration of the concept of centrality. The impact parameter,
number of participants, and number of charged particles are shown in the
horizontal axis. The interaction cross-section is shown in the vertical axis,
in arbitrary units. From [39].

The concept of centrality is based on the intuitive idea that collisions with
greater overlap area between the colliding nuclei (small impact parameter)
should produce more particles, while in collisions with a smaller overlap
area, less particles should be observed in the final state. In this way, it is
possible to classify collisions (events) based on the number of detected
particles, and arrange them in centrality classes. Out of all events, the 10 %
which produced more particles constitute the 0 - 10 % centrality class. The 20
% ones that produced more particles, removing those from the 0 - 10 % class,
constitute the 10 - 20 % class, and so on.

Events where the nuclei collide “head-on” are called central collisions, and
as the impact parameter grows the collisions are said to be more peripheral.
Collisions where the impact parameter is greater than twice the nuclear
radius, so that there is no overlap of hadronic matter, are called ultra-
peripheral [40]. Figure 2.3 illustrates the concept of centrality: the blue
line shows the nucleus-nucleus cross-section distribution as a function of the
number of charged particles (the impact parameter is also indicated) and the
dashed lines show the centrality percentiles. Nucleons that collide with at
least another (the darker circles in Figure 2.11) are called participants,
while the remaining are known as spectators. Central collisions are
characterized by a large number of participants, while in more peripheral
collisions, there are less participant nucleons.

###  2.4 Experimental observables and evidences of QGP formation

####  2.4.1 Multiplicity

Multiplicity is a very straightforward observable, defined simply as the
number of charged particles detected in the final state of the collision,
usually measured in a given pseudo-rapidity region. Nevertheless, multiplicity
in the central pseudo-rapidity region is an important observable to study the
bulk properties of the QGP formed in heavy-ion collisions. Table 2.1 shows,
for each centrality class, the charged-particle multiplicity at mid-rapidity
and the mean number of participants for Pb-Pb collisions at
$\sqrt{s_{\text{NN}}}$ = 2.76 TeV, as measured by the ALICE Collaboration
[41]. The number of participant nucleons in each centrality class was obtained
from Glauber Model estimates (see Section 2.5.1).

Centrality | $dN_{\text{ch}}/d\eta$ | $N_{\text{part}}$
---|---|---
0 - 5 % | 1601 $\pm$ 60 | 382.8 $\pm$ 3.1
5 - 10 % | 1294 $\pm$ 49 | 329.7 $\pm$ 4.6
10 - 20 % | 966 $\pm$ 37 | 260.5 $\pm$ 4.4
20 - 30 % | 649 $\pm$ 23 | 186.4 $\pm$ 3.9
30 - 40 % | 426 $\pm$ 15 | 128.9 $\pm$ 3.3
40 - 50 % | 261 $\pm$ 9 | 85.0 $\pm$ 2.6
50 - 60 % | 149 $\pm$ 6 | 52.8 $\pm$ 2.0
60 - 70 % | 76 $\pm$ 4 | 30.0 $\pm$ 1.3
70 - 80 % | 35 $\pm$ 2 | 15.8 $\pm$ 0.6
Table 2.1: Charged-particle multiplicity density at mid-rapidity and mean
number of participant nucleons (obtained from Glauber Model estimates) for
each centrality class in Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76
TeV, as measured by the ALICE experiment. Data from [41].

The number of charged particles observed in the final state is much larger
than the original number of participant nucleons, and decreases as the
collisions become more peripheral. This should be expected, since for more
peripheral collisions there is a decrease in the overlap area between the two
nuclei. To study bulk particle production, it is more convenient to calculate
the multiplicity per participant pair, as in Figure 2.4. This also allows
comparisons with particle production in other collision systems, such as p-p
and p-Pb. Figure 2.4 shows that multiplicity is not merely proportional to the
number of participants: in the most central collisions, approximately 10
charged particles are produced for each pair of participants, while for the
most peripheral ones, this number drops to about 4 particles.

Figure 2.4: Charged particle multiplicity density in $|\eta|<0.5$ divided by
the number of participant pairs as a function of the number of participant
nucleons for different collision systems. From [32].

Figure 2.4 also shows, as should be expected, that multiplicity increases with
increasing collision energy, which can be seen by noticing that the data of
Pb-Pb at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV is multiplied by 1.2. From the
theoretical point of view, particle production is usually described by two
categories of models: two-component models which combine soft interactions and
perturbative QCD processes [42, 43] and saturation models [44, 45, 46].

####  2.4.2 Transverse momentum distributions

Transverse momentum distributions ($p_{T}$-spectra) are also among the most
commonly measured observables. A transverse momentum distribution is defined
simply as a histogram counting the number of particles detected in each
$p_{T}$ bin, per unit rapidity. $p_{T}$-spectra contain information about the
kinetic properties of the final state particles, and therefore are a powerful
tool to study the bulk properties of the QGP. Figure 2.5 shows, for each
centrality class, the $p_{T}$-spectra of pions, kaons and protons produced in
Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV, measured by the ALICE
Collaboration.

Figure 2.5: Transverse momentum distributions of pions, kaons and protons in
Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV measured by ALICE. From
[47].

Most particles are produced at low-$p_{T}$ ($p_{T}\lesssim$ 1 GeV) and the
soft region of the spectra, with $p_{T}\lesssim$ 3 GeV, is well described by a
thermal distribution. The high-$p_{T}$ (hard) part of the spectra, on the
other hand, exhibits a power-law behavior. Many statistical models have been
employed to extract physical parameters by fitting $p_{T}$ spectra [48, 49].
These include the non-extensive Tsallis statistics [50], the QCD-inspired
Hagedorn inverse power law [51, 52] and the Pearson distribution [53]. A
phenomenological model widely used to characterize $p_{T}$-spectra and obtain
information about the kinetic freeze-out is the Blast Wave (BW) model [54],
which is used in this work (see Section 5.2.4).

Figure 2.6: Mean $p_{T}$ of charged pions (left) and protons (right) in Pb-Pb
collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV and Au-Au collisions at
$\sqrt{s_{\text{NN}}}$ = 200 GeV, as a function of multiplicity. From [47].

A straightforward quantity to extract from $p_{T}$-spectra is the mean
transverse momentum of the detected particles (mean $p_{T}$):

$\langle p_{T}\rangle=\frac{\int d^{2}p_{T}p_{T}\frac{dN}{d^{2}p_{T}dy}}{\int
d^{2}p_{T}\frac{dN}{d^{2}p_{T}dy}}$ (2.4)

Figure 2.6 shows that the mean $p_{T}$ increases with multiplicity, suggesting
that the transverse expansion of the system is somewhat more “violent” in
central collisions.

####  2.4.3 Anisotropic flow

Unlike gases, where particles are far apart from each other and rarely meet,
in fluids particles are constantly interacting with their neighbors, so that
fluids present collective behavior. The experimental observation of collective
behavior in relativistic heavy-ion collisions is probably the most compelling
evidence that indeed a QGP is formed in such experiments. The QGP itself is
never observed: collective behavior manifests itself as anisotropy in the
momentum distribution of the final state particles. Due to the fact that in a
non-central collision the overlap region of the two nuclei has an
approximately elliptical (“almond”) shape, greater pressure gradients develop
in the $x$-direction as compared to the $y$-direction. In the hydrodynamic
evolution, the QGP flows preferentially in the $x$-direction, as the fluid
expands and the elliptical shape of the system becomes more circular. In this
process, it is said that the spatial anisotropy is transferred to momentum
space. Figure 2.7 illustrates, for a non-central collision, the time evolution
of the shape of the system in the transverse plane:

Figure 2.7: Time evolution of the shape of the system in the transverse plane
of the collision in a non-central event. The initially elliptical shape turns
more circular as the system expands, and particles flow preferably in the
$x$-direction. Figure from [55].

As a result, more particles are detected close to $\phi$ = 0 and $\phi$ =
$\pi$ and less particles are detected near $\phi$ = $\pi/2$ and $\phi$ =
$3\pi/2$. This momentum anisotropy in the transverse plane can be quantified
by expanding the azimuthal distribution of particles in a Fourier series [56]:

$\frac{dN}{d\phi}\propto
1+2\sum_{n=1}^{\infty}v_{n}\cos[n(\phi-\Psi_{\text{RP}})]$ (2.5)

The coefficients of the series (flow harmonics) are given by:

$v_{n}=\langle\cos[n(\phi-\Psi_{\text{RP}})]\rangle,$ (2.6)

where the expectation value symbol denotes an average over all particles and
$\Psi_{\text{RP}}$ is the reaction plane angle [57]. The series has only the
cosine terms, due to the reflection symmetry with respect to the reaction
plane. The first coefficient of the series is known as directed flow, while
the harmonics $v_{2}$ and $v_{3}$ are called elliptic and triangular flow,
respectively.

Figure 2.8: Illustration of particle distributions in the transverse plane of
the collision with (a) $v_{2}>0$ and (b) $v_{2}$ = 0. Figure from [58].

Figure 2.9 shows, on the left, the flow harmonics integrated over transverse
momentum as a function of centrality, and on the right, the flow harmonics as
a function of $p_{T}$ (usually called differential flow) for three centrality
classes.

Figure 2.9: Flow measurements in Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ =
2.76 TeV. Left: integrated flow harmonics as a function of centrality. Right:
$p_{T}$-differential flow for three centrality classes. From [59].

There is little flow in central collisions, where the initial spatial
anisotropy is small. Flow grows for non-central collisions, as the overlap
region between the two colliding nuclei acquires an anisotropic shape. For
very peripheral collisions, flow decreases once again. It can also be seen
that, inside a given centrality class, elliptic flow increases with transverse
momentum. Figure 2.9 also shows, on the bottem left panel, the ratio between
the flow harmonics and the corresponding eccentricity harmonics (see Section
2.5.2).

###  2.5 The Initial Condition

####  2.5.1 The Glauber Model

In general, the distinction between participants and spectators is not as
straightforward as Figure 2.3 might suggest. Moreover, some participants
collide only once, while some might suffer multiple collisions, and there is
also the possibility of two nucleons “passing by” each other without
interacting at all. A quantitative estimate of the number of participant
nucleons ($N_{\text{part}}$) and nucleon-nucleon (or binary) collisions
($N_{\text{coll}}$) can be provided by the Glauber Model [60, 61]. In the
1950’s Roy Glauber pioneered the use of quantum scattering theory for
composite systems such as heavy nuclei, allowing the calculation of geometric
quantities which can be related to the experimental observables. The Glauber
formalism treats collisions of composite systems, such as p-A and A-A as a
superposition of collisions between their constituent nucleons. The model uses
the eikonal approximation for quantum scattering, which assumes that the
exchanged momentum between the colliding projectiles is very small compared to
their total longitudinal momentum, so that the nucleons travel in straight
trajectories.

There are two most relevant inputs for any Glauber based calculation. The
first of them is the nuclear density distribution, for which the most natural
choice is a Woods-Saxon distribution:

$\rho(r)=\rho_{0}\left(1+\exp{\left(\frac{r-R}{a}\right)}\right)^{-1}$ (2.7)

In Equation (2.7), $\rho_{0}$ is the nucleon density in the center of the
nucleus, $R$ corresponds to the nuclear radius, $a$ is the skin depth and $r$
is the distance from the center of the nucleus. The second one is the
inelastic nucleon-nucleon cross section, which is assumed to be dependent only
on energy, and not on the collision system, the position of the nucleon inside
the nucleus, or any other characteristic of the nuclear environment. Table 2.2
shows some measured values of the nucleon-nucleon inelastic scattering cross-
section ($\sigma_{\text{NN}}$) for different collision energies.

$\sqrt{s}$ (TeV) |  $\sigma_{\text{NN}}$ ($\text{fm}^{2}$) | Reference
---|---|---
0.2 | 4.23 | [62]
2.76 | 6.28 | [63]
5.02 | 7.0 | [64]
7 | 7.32 | [63]
Table 2.2: Nucleon-nucleon inelastic scattering cross-section for different
collision energies

In a collision of two nuclei $A$ and $B$ at impact parameter $b$, the
transverse density of each nucleus is represented by the thickness function,
obtained by integrating the nuclear density in the beam direction:

$T_{A}(x,y)=\int_{-\infty}^{+\infty}\rho_{A}(x,y,z)dz,$ (2.8)

and analogously for B. The density of binary collisions at point $(x,y)$ is
given by the product of the thickness functions of the two nuclei and the
inelastic nucleon-nucleon cross-section:

$n_{\text{BC}}(x,y;b)=\sigma_{\text{NN}}T_{A}(x+b/2,y,z)T_{B}(x-b/2,y,z),$
(2.9)

from which the total number of binary collisions is obtained by integration in
the transverse plane:

$N_{\text{BC}}(b)=\int\sigma_{\text{NN}}T_{A}(x+b/2,y,z)T_{B}(x-b/2,y,z)dxdy.$
(2.10)

The number of participants (or wounded nucleons), on its turn, is given by the
integral [65]:

$\begin{split}N_{\text{WN}}(b)=\int
T_{A}(x+b/2,y,z)\left(1-\left(1-\frac{\sigma_{\text{NN}}T_{B}(x-b/2,y,z)}{B}\right)^{B}\right)\\\
+T_{B}(x-b/2,y,z)\left(1-\left(1-\frac{\sigma_{\text{NN}}T_{A}(x+b/2,y,z)}{A}\right)^{A}\right)dxdy\end{split}$
(2.11)

In modern computational simulations, the Glauber Model is usually implemented
through Monte Carlo calculations. In the Monte Carlo Glauber approach, the two
colliding nuclei are assembled according to a nuclear density distribution on
an event-by-event basis. Then, the binary collisions are performed.

Figure 2.10: Event generated using a Monte Carlo Glauber approach.
Participants are represented by darker circles. Left: transverse plane view.
Right: view along the beam axis. Figure from [60].

The simplest criterion that can be employed to count the number of binary
collisions is simply geometric, and justified by the black disk approximation
[66]. It consists in assuming that a nucleon-nucleon collision occurs every
time that:

$d\leq\sqrt{\frac{\sigma_{\text{inel}}^{NN}}{\pi}},$ (2.12)

where $d$ is the distance between two nucleons in the transverse plane. It is
well known that the nucleons themselves are composite particles, and it should
be noted that there has also been some progress in performing Glauber
calculations at the sub-nucleonic level [67].

####  2.5.2 Eccentricity harmonics

In a collision with $b=0$, it is natural to imagine that the transverse
geometry of the system should be approximately circular. In a non-central
collision, on the other hand, the overlap region of the two colliding nuclei
should have a rather elliptic shape, being elongated in the $y$ direction.
Besides this evident “almond” shape, more complex geometries of the system are
generated by event-by-event fluctuations on the spatial distribution of the
nucleons inside the nuclei.

Figure 2.11: Transverse view of two colliding nuclei. In central events the
geometry is approximately circular, and in non-central events it is
approximately elliptical, but more complex irregularities are also always
presents. Figure from [68].

In practice, the shape of the overlap area of the two colliding nuclei is not
perfectly elliptical, and Figure 2.11 shows that even for $b=0$ the overlap
region is not perfectly circular: there is always some spatial anisotropy.
This spatial anisotropy can be quantified by the eccentricity harmonics:

$\varepsilon_{n}=\frac{\int r^{n}e^{in\varphi}s(r,\varphi)rdrd\varphi}{\int
r^{n}s(r,\varphi)rdrd\varphi},$ (2.13)

which are calculated using the transverse entropy density distribution of the
system $s(x,y)$ as a weight function. The eccentricity harmonics have a
geometric interpretation: each harmonic quantifies the spatial anisotropy
associated with a particular geometric shape in the initial state (Figure
2.12). In this way, $\varepsilon_{2}$ is appropriately called ellipticity,
while $\varepsilon_{3}$ and $\varepsilon_{4}$ are referred to as triangularity
and quadrangularity, respectively.

Figure 2.12: Eccentricity harmonics $\varepsilon_{n}$ for $n$ = 2, 3, 4.

On an event-by-event basis, there is a strong relation between the initial
state geometry and final state flow observables. The initial geometry of the
system dictates the direction of the pressure gradients in the beginning of
the hydrodynamic evolution, so that is said that the final state momentum
anisotropy is a hydrodynamic response to the initial state anisotropy. For $n$
= 2, 3 there is a mapping between the eccentricity harmonics and its
corresponding flow harmonics [69, 70, 71], which can be expressed as:

$v_{n}=f(\varepsilon_{n})+\delta_{n},$ (2.14)

where $\delta_{n}$ is an error. As a first approximation, which is excellent
for central events, a linear mapping can be used:

$f(\varepsilon_{n})=\kappa_{n}\varepsilon_{n},$ (2.15)

but for more peripheral events a linear + cubic mapping is a better predictor:

$f(\varepsilon_{n})=\kappa_{n}\varepsilon_{n}+\kappa_{n}^{{}^{\prime}}|\varepsilon_{n}|^{2}\varepsilon_{n}$
(2.16)

Figure 2.13 shows the mapping between ellipticity and elliptic flow, and
between triangularity and triangular flow for Pb-Pb collisions at
$\sqrt{s_{\text{NN}}}$ = 2.76 TeV in the 45 - 50 % centrality class.

Figure 2.13: Mapping between eccentricity harmonics and flow harmonics for
$n=2,3$. Each blue dot corresponds to an event. The dotted black line
represents the linear estimator, and the full red line represents the cubic
estimator. From [69].

##  Simulating heavy-ion collisions

Contradicting the expectation that the deconfined state formed in heavy-ion
collisions would be similar to a “gas” of quarks and gluons, experimental data
collected at RHIC showed that instead a strongly interacting state is formed,
similar to an almost ideal fluid. This discovery motivated the employment of
relativistic hydrodynamics modelling to simulate heavy-ion collisions, which
had great success in describing a variety of soft hadronic observables.
Hydrodynamics has passed a series of tests and showed to be a good effective
description of the bulk evolution of the QGP. More recently, hybrid
simulations, where a different model is used to simulate each specific stage
of the collision, are considered the most modern way to simulate heavy-ion
collisions.

Chapter 3 provides an introduction to the models which constitute the hybrid
simulation chain used to generate the data presented in this work. Section 3.1
introduces T${}_{\text{R}}$ENTo, the initial condition generator. Section 3.2
and Section 3.3 present K$\o{}$MP$\o{}$ST and MUSIC, which are responsible for
simulating the pre-hydrodynamic phase and hydrodynamic phase of the collision,
respectively. Section 3.4 and Section 3.5 are devoted to the late stages of
the collision: particlization, and final state resonance decays and hadronic
interactions.

###  3.1 The initial condition: T${}_{\text{R}}$ENTo

T${}_{\text{R}}$ENTo (Reduced Thickness Event-by-event Nuclear Topology) is a
Glauber inspired model used to generate an initial entropy density profile of
two colliding projectiles (proton-proton, proton-nucleus or nucleus-nucleus)
in the transverse plane of the collision [72, 73]. It is an effective model,
meaning that it is not based in first principles, and no assumptions are made
about the specific mechanisms of entropy production.

Figure 3.1: Transverse entropy distributions from T${}_{\text{R}}$ENTo for
three values of the entropy deposition parameter $p$. Figure from [73].

In a collision between two protons $A$ and $B$, with impact parameter $b$ in
the $x$ direction, their nuclear densities are given by:

$\rho_{A,B}=\rho_{\text{proton}}(x\pm b/2,y,z).$ (3.1)

In T${}_{\text{R}}$ENTo, the thickness function of each nucleon:

$T_{A,B}(x,y)=\int dz\rho_{A,B}(x,y,z)$ (3.2)

is given by a boosted two-dimensional Gaussian:

$T_{A,B}(x,y)=\frac{1}{2\pi
w^{2}}\exp\left(-\frac{x^{2}+y^{2}}{2w^{2}}\right)$ (3.3)

A fluctuated thickness is assigned to each one of the protons $A$ and $B$:

$\tilde{T}_{A,B}(x,y)=w_{A,B}.T_{A,B}(x,y),$ (3.4)

where $w_{A,B}$ are independent random weights sampled from a gamma
distribution with unit mean:

$P_{k}(w)=\frac{k^{k}}{\Gamma(k)}w^{k-1}e^{-kw}.$ (3.5)

These weights introduce additional multiplicity fluctuations, which are
necessary to reproduce the large multiplicity fluctuations observed
experimentally in proton-proton collisions [74]. T${}_{\text{R}}$ENTo proposes
a function $f(\tilde{T}_{A},\tilde{T}_{B})$ to convert projectile thickness
into entropy deposition:

$f(\tilde{T}_{A},\tilde{T}_{B})\propto\frac{dS}{d^{2}x_{\perp}d\eta}=s(x,y)$
(3.6)

This function is the reduced thickness, which is a generalized mean controlled
by a real parameter $p$:

$f=\tilde{T}_{R}(p;\tilde{T}_{A},\tilde{T}_{B})=\left(\frac{\tilde{T}_{A}^{p}+\tilde{T}_{B}^{p}}{2}\right)^{1/p}.$
(3.7)

Entropy is determined by the reduced thickness up to an overall normalization
factor. For certain values of $p$, the reduced thickness simplifies to the
arithmetic, geometric and harmonic means:

$\tilde{T}_{R}=\small\begin{dcases}\text{min}(\tilde{T}_{A},\tilde{T}_{B})&\quad
p\rightarrow-\infty\\\
2\tilde{T}_{A}\tilde{T}_{B}/(\tilde{T}_{A}+\tilde{T}_{B})&\quad
p=-1\quad\text{(harmonic)}\\\ \sqrt{\tilde{T}_{A}\tilde{T}_{B}}&\quad
p=0\quad\text{(geometric)}\\\ (\tilde{T}_{A}+\tilde{T}_{B})/2&\quad
p=1\quad\text{(arithmetic)}\\\ \text{max}(\tilde{T}_{A},\tilde{T}_{B})&\quad
p\rightarrow+\infty\\\ \end{dcases}$ (3.8)

Figure 3.2: Reduced thickness function for three values of the entropy
deposition parameter $p$ in a non-central collision. Figure taken from [72].

Composite systems are treated as superpositions of proton-proton collisions.
In a nucleus-nucleus collision, a set of nucleon positions is sampled for each
projectile using an uncorrelated Woods-Saxon distribution. In this process,
nucleons are forbidden to be placed closer than a distance $d_{\text{min}}$
from another one that was previously positioned. For each pair of nucleons,
the collision probability:

$P_{\text{coll}}(b)=1-\exp\left[-\sigma_{gg}T_{AB}(b)\right]$ (3.9)

is sampled once to decide if the nucleons collide. In Equation (3.9),
$T_{AB}(b)$ is the overlap integral of the two nucleons $A$ and $B$:

$T_{AB}(b)=\int dxdyT_{A}(x-b/2,y)T_{B}(x+b/2,y)$ (3.10)

The nucleons which interact at least once are labelled as participants and
assigned a fluctuated thickness, and the others are discarded. Then, in a
collision between two nuclei $A$ and $B$, the fluctuated thickness for nucleus
$A$ is the sum of the individual fluctuated thickness of each of its
participants:

$\tilde{T}_{A}=\sum_{i=1}^{N_{\text{part}}}w_{i}\int
dz\rho_{\text{proton}},(x-x_{i},y-y_{i},z-z_{i})$ (3.11)

and analogously for nucleus $B$. $\sigma_{gg}$ is an effective parton-parton
cross-section tuned so that the total proton-proton cross section matches the
experimental inelastic nucleon-nucleon cross-section:

$\int 2\pi bdbP_{\text{coll}}(b)=\sigma_{\text{NN}}^{\text{inel.}}$ (3.12)

Following the Glauber model assumptions, the inelastic nucleon-nucleon cross-
section remains constant, and the probability of two nucleons interacting
depends solely on the impact parameter. It is independent of the number of
previous collisions suffered by each nucleon, and the location of the nucleon
inside the nucleus.

###  3.2 Pre-equilibrium dynamics: K$\o{}$MP$\o{}$ST

The initial state formed just after the collision is far from equilibrium, and
complicated. In the time scale of approximately 1 fm/c, this out of
equilibrium initial state evolves to a state of local equilibrium, which
defines the beginning of the hydrodynamic evolution. K$\o{}$MP$\o{}$ST [75]
proposes an effective macroscopic description of the dynamics of the system in
the early moments of the collision, providing a map of the form:

$T^{\mu\nu}(\tau_{0},\mathbf{x})|_{\text{out-of-equilibrium}}\longrightarrow
T^{\mu\nu}(\tau_{\text{hydro}},\mathbf{x}),$ (3.13)

relating the energy-momentum tensor in an early time $\tau_{0}\sim$ 0.1 fm/c
to a later time where hydrodynamics should become applicable,
$\tau_{\text{hydro}}\sim$ 1 fm/c.

Figure 3.3: Ilustration of the time evolution of the simulation: an initial
condition passes through K$\o{}$MP$\o{}$ST pre-equilibrium before entering
relativistic hydrodynamics. Figure from [76].

K$\o{}$MP$\o{}$ST initially divides the energy-momentum tensor into a local
homogeneous background and small perturbations around it:

$T^{\mu\nu}(\tau_{0},\mathbf{x^{\prime}})=\overline{T}^{\mu\nu}(\tau_{0})+\delta
T_{\mu\nu}(\tau_{0},\mathbf{x^{\prime}})$ (3.14)

As a first approximation, these perturbations can be studied in the framework
of linear response theory. Then the energy-momentum tensor at
$\tau_{\text{hydro}}$ can be calculated by adding to the evolved background
the response to the initial perturbations from equilibrium:

$T^{\mu\nu}(\tau_{\text{hydro}},\mathbf{x^{\prime}})=\overline{T}_{\mathbf{x}}^{\mu\nu}(\tau_{\text{hydro}})+\frac{\overline{T}_{\mathbf{x}}^{\tau\tau}(\tau_{\text{hydro}})}{\overline{T}_{\mathbf{x}}^{\tau\tau}(\tau_{0})}\int
d^{2}\mathbf{x^{\prime}}G_{\alpha\beta}^{\mu\nu}(\mathbf{x},\mathbf{x^{\prime}},\tau_{\text{hydro}},\tau_{0})\delta
T_{\mathbf{x}}^{\alpha\beta}(\tau_{0},\mathbf{x^{\prime}})$ (3.15)

The first term is the (non-linear) equilibration of the background and the
second term is a convolution of the initial perturbations and the response
functions, which is responsible for populating the non-diagonal terms of the
energy-momentum tensor. The second term is the convolution of the initial
deviations from equilibrium and the response functions, normalized by the
background energy density. It is important to note that, due to causality, the
energy-momentum tensor at ($\tau_{\text{hydro}},\mathbf{x}$) is only affected
by cells that are within the causal neighborhood of $\mathbf{x}$:

$|\mathbf{x^{\prime}}-\mathbf{x}|<c(\tau_{\text{hydro}}-\tau_{0}),$ (3.16)

which is represented by the white circle in Figure 3.4.

The non-equilibrium evolution of the background and the response functions
have to be calculated according to an underlying microscopic description. A
microscopic description of the out-of-equilibrium system formed in the
collision can be provided by QCD effective kinetic-theory (EKT) [77] :

$p^{\mu}\partial_{\mu}f(x,p)=\mathcal{C}_{2\leftrightarrow
2}[f]+\mathcal{C}_{1\leftrightarrow 2}[f]$ (3.17)

In this relativistic Boltzmann equation, $\mathcal{C}_{2\leftrightarrow 2}[f]$
is the collision integral for elastic scatterings at leading order and
$\mathcal{C}_{1\leftrightarrow 2}[f]$ is the collision integral for inelastic,
particle number changing processes. K$\o{}$MP$\o{}$ST can also be used in a
free streaming limit. In the free streaming scenario the partons, which are
taken to be massless, don’t interact, so that the evolution of the system is
dictated by a homogeneous Boltzmann equation:

$p^{\mu}\partial_{\mu}f(x,p)=0$ (3.18)

In this case, the perturbations around equilibrium can be obtained
analytically and expressed in terms of the Bessel functions [75].

Figure 3.4: Transverse energy density distribution evolved from
$\tau_{\text{EKT}}$ = 0.2 fm/c to $\tau_{\text{hydro}}$ = 1.2 fm/c using
K$\o{}$MP$\o{}$ST EKT. The white circle represents the causal circle of the
white dot. From [75].

The energy-momentum tensor coming from T${}_{\text{R}}$ENTo (or from any
Glauber-based calculation) initially has the simple, diagonal form:

$T^{\mu\nu}(\tau_{0},\mathbf{x})=\begin{pmatrix}e(\mathbf{x})&0&0&0\\\
0&\frac{1}{2}e(\mathbf{x})&0&0\\\ 0&0&\frac{1}{2}e(\mathbf{x})&0\\\
0&0&0&0\end{pmatrix}$

At the end of the simulation, the energy-momentum tensor can be decomposed in
the Landau frame and written in the usual form used in hydrodynamic
simulations:

$T^{\mu\nu}=eu^{\mu}u^{\nu}+P(e)\Delta^{\mu\nu}+\pi^{\mu\nu},$ (3.19)

where $\Delta^{\mu\nu}=g^{\mu\nu}+u^{\mu}u^{\nu}$, and $P=P(e)$ is the
equation of state (see Section 3.3). At this point, the energy density
distribution of the system should be smoother, closer to the assumptions of
local equilibrium required by hydrodynamics (Figure 3.4).

###  3.3 Hydrodynamic evolution: MUSIC

After pre-equilibrium and thermalization, the dynamics of the QGP is well
described by relativistic hydrodynamics [78, 79], which is based on local
conservation laws:

$\partial_{\mu}T^{\mu\nu}=0$ (3.20)

For a relativistic ideal fluid, the energy-momentum tensor takes the form:

$T^{\mu\nu}=eu^{\mu}u^{\nu}-P(g^{\mu\nu}-u^{\mu}u^{\nu}),$ (3.21)

where $e$ is the energy density and $P$ is the pressure. $u^{\mu}$ is the
fluid velocity field:

$u^{\mu}=\frac{dx^{\mu}}{d\tau},$ (3.22)

which is normalized so that $u^{\mu}u_{\mu}=1$, and the derivative is taken
with respect to proper time. Including the viscous corrections, the energy-
momentum tensor can be written as:

$T^{\mu\nu}=eu^{\mu}u^{\nu}-(P+\Pi)\Delta^{\mu\nu}+\pi^{\mu\nu}$ (3.23)

In Equation (3.23) $\pi^{\mu\nu}$ is the shear stress tensor, $\Pi$ is the
bulk viscous pressure and $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ is
called the projection operator. Shear viscosity is a measure of the resistance
between layers of a fluid, while the bulk viscosity is associated with the
fluid expansion rate (Figure 3.5). Pressure is not independent, it is related
to the energy density by the equation of state $P=P(e)$.

Figure 3.5: Illustration of shear viscosity (left) and bulk viscosity (right).

MUSIC [80] is a publicly available C++ code which performs 2+1D relativistic
hydrodynamics simulations. Given an initial condition (which is the energy
density distribution of the system at $\tau_{\text{hydro}}$) and an equation
of state, MUSIC evolves the energy-momentum tensor by solving the hydrodynamic
equations of motion [81, 82] numerically, while respecting the conservation
laws in Equation (3.20).

MUSIC uses temperature-dependent parametrizations for the shear viscosity to
entropy density ratio $\eta/s$ and also for the bulk viscosity to entropy
density ratio $\zeta/s$. The entropy density can be used as a substitute for
number density, so that these quantities are an attempt of capturing the
“viscosity per unit”. For the shear viscosity, which should reach a minimum
near the critical temperature $T_{c}$, a modified linear parametrization is
used:

$(\eta/s)(T)=(\eta/s)_{\text{min}}+(\eta/s)_{\text{slope}}(T-T_{c})\left(\frac{T}{T_{c}}\right)^{(\eta/s)_{\text{curve}}},$
(3.24)

which contains three free parameters: a minimum value, a slope above $T_{c}$
and a curvature parameter. When the curvature parameter
$(\eta/s)_{\text{curve}}$ is set to zero, Equation (3.24) reduces to a simple
linear parametrization. For the bulk viscosity, an unnormalized Cauchy
distribution is used:

$(\zeta/s)(T)=\frac{(\zeta/s)_{\text{max}}}{1+\left(\frac{T-(\zeta/s)T_{0}}{(\zeta/s)_{\text{width}}}\right)^{2}},$
(3.25)

which also has three tunable parameters: a maximum value, the width of the
peak and its location $T_{0}$. Figure 3.6 shows the functional form of
Equations (3.8) and (3.9) for some values of the free parameters.

Figure 3.6: Shear and bulk viscosity parametrizations given in Equations (3.8)
and (3.9) for different choices of the free parameters. Figure taken from
[68].

In MUSIC, reaching a switching temperature $T_{\text{switch}}$ is the
criterion used for stopping hydrodynamics. Connecting all the cells at which
the temperature reached the switching temperature forms a 4-dimensional
spacetime hypersurface (Figure 3.8). This hypersurface is complete when the
temperature of all cells drops to this switching temperature, and hydrodynamic
evolution reaches an end.

###  3.4 Particlization: iSS

Once all points of the grid have cooled down to the switching temperature, it
is necessary to change the description of the system from the hydrodynamic
picture to a hadron gas description, so that each cell must be converted into
discrete particles. This change from hydrodynamics to a microscopic
description is known as particlization [83].

Figure 3.7: Illustration of the particlization procedure: a fluid cell is
converted into hadrons.

The momentum distribution of particles of species $i$ (with degeneracy factor
$g_{i}$) emerging from particlization can be calculated using the Cooper-Frye
prescription [84]:

$E\frac{dN_{i}}{dp^{3}}=\frac{g_{i}}{(2\pi)^{3}}\int_{\Sigma}f(p)p^{\mu}d\sigma_{\mu}$
(3.26)

The integral in Equation (3.26) is performed on the isothermal hypersurface
$\Sigma$ of temperature $T_{\text{switch}}$ and $d\sigma_{\mu}$ is a volume
element of the four-dimensional surface normal to the surface. The integral is
calculated using the one-particle distribution function.

Figure 3.8: Two-dimensional visualization of a freeze-out hypersurface at
$T_{\text{switch}}$ = 150 MeV (red line). Figure from [80].

If the system was at perfect thermal equilibrium, the distribution function
would simply be the Bose-Einstein or Fermi-Dirac distribution:

$f_{0}(p)=\frac{1}{e^{p\cdot u/T}\mp 1}$ (3.27)

This is, however, not the case and out-of-equilibrium corrections to the
distribution function must be made. These corrections can be divided in two
main categories. One possibility is to use small linear corrections to the
equilibrium distribution function $f_{0}$ [85, 86]:

$f=f_{0}+\delta f$ (3.28)

Another strategy is to perform a transformation in the momentum vector inside
the distribution function [87] using a linear transformation matrix
$\lambda_{ij}$:

$p_{i}\rightarrow p_{i}^{\prime}+\sum_{j}\lambda_{ij}p_{j}$ (3.29)

In this work, we use a linear correction derived from the Boltzmann equation
using relaxation time approximation (RTA) [88, 89, 90]:

$\delta f=f_{0}(1\pm
f_{0})\frac{\tau}{ET}\left[\frac{1}{2\eta}p^{i}p^{j}\pi_{ij}+\frac{1}{\zeta}\left(\frac{p^{2}}{3}-c_{s}^{2}E^{2}\right)\Pi\right],$
(3.30)

where $\pi_{ij}$ is the shear stress tensor and $\tau$ is a shear and bulk
relaxation time, which is taken to be constant and the same for all particle
species.

iSS (iSpectraSampler) [91, 92] is a Monte Carlo particle sampler that uses the
Cooper-Frye formula to calculate the momentum distribution of particles coming
from a freeze-out hypersurface. iSS generates sets of momenta and positions
for particles emitted at each cell at the end of the hydrodynamic simulation.
The procedure of doing this many times on the same hypersurface is known as
oversampling. Oversampling is often used due to the fact that the hydrodynamic
part of the simulation is the most expensive one computationally, and one
would like to “reuse” the same hydrodynamic evolution to generate more events.

###  3.5 Hadronic phase: UrQMD

Once the hydrodynamic cells have been converted into particles, it is
necessary to simulate the last stage of the collision, in which these hadrons
travel to the detector. In the way, they interact and resonances might decay
into more stable species. A microscopic description of this hadronic phase is
provided by the Ultra-Relativistic Quantum Molecular Dynamics Model (UrQMD)
[93, 94], which was developed in the late 90’s.

UrQMD is a microscopic transport model, in which particles travel in classical
trajectories and interact stochastically according to experimental hadron-
hadron cross-sections. The UrQMD package [95] solves a Boltzmann equation for
each hadron species:

$\frac{df_{i}(x,p)}{dt}=\mathcal{C}_{i}(x,p),$ (3.31)

where the collision term accounts for $2\leftrightarrow 2$ and
$1\leftrightarrow 2$ processes, such as binary collisions and decays. UrQMD
contains 53 baryon species and 24 different meson species. Two particles
collide if the distance between them satisfies:

$d\leq d_{0}=\sqrt{\frac{\sigma_{\text{tot}}}{\pi}},$ (3.32)

where the cross-section is taken to be the free cross-section of the
interaction in question. The interactions and decays performed by UrQMD
constitute the last stage of the hybrid simulation.

##  Motivation

Chapter 4 presents the motivation of this work: evidences that the size of
nucleons, which is a free parameter in the initial condition, has been
overestimated in most recent Bayesian Analyses, in which the model is
confronted with experimental data to obtain a best-fit value of each
parameter. Section 4.1 brings a short review about the current understanding
about the nucleon size. Section 4.2 provides an overview of the general
structure of Bayesian Analyses used to constrain the parameters of hybrid
simulations of heavy-ion collisions. In Section 4.3, the motivation of the
work is exposed in detail, which sets the ground for Chapter 5.

###  4.1 The nucleon size

Experimentally, probing inside hadrons is a huge challenge. From the
theoretical point of view, studying the hadronic structure is also a hard
task, as it is associated with the non-perturbative part of QCD. Not by
accident, details about the structure of hadrons are in general still poorly
known, including information about the nucleons.

The charge radius of the proton can be defined as the slope of its
electromagnetic form factor at zero momentum transfer:

$r_{p}^{2}=6\left.\frac{dG_{p}^{E}}{dQ^{2}}\right|_{Q^{2}=0}$ (4.1)

The first indirect measurement of the proton charge radius was made by
Hofstader et al. in the famous electron scattering experiment [96, 97] which
later on won the Nobel Prize in 1961. In this experiment, the proton radius
was estimated by fitting the electric form factor with a dipole form. The
charge radius of the proton can also be determined by Lamb shift measurements,
performed for both electronic and muonic hydrogen atoms [98]. Most of the
first Lamb shift experiments performed with electronic atoms gave as a result
for the proton charge radius roughly the same value as electron scattering
experiments, $r_{p}$ = 0.88 fm, which came to be known as the large radius.
Lamb shift results with muonic atoms, on the other hand, gave the so-called
small radius, $r_{p}$ = 0.84184(67) fm, which differed by 5$\sigma$ from the
large radius. This discrepancy came to be known as the proton radius puzzle
[99, 100]. As both theory and experiment advanced, electronic and muonic Lamb
shift experiments came to an agreement, and the present CODATA value for the
proton charge radius is $r_{p}$ = 0.8414(19) fm [101], closer to the small
radius. Table 4.1 shows the results of some of the most modern measurements of
the proton charge radius performed with different experimental techniques. The
Paris electronic Lamb shift measurement from 2018 [102] remains an exception,
with an extracted charge radius that is actually closer to the large radius.

$r_{p}$ (fm) | Year | Method | Reference
---|---|---|---
0.8335(95) | 2017 | Electronic Lamb shift | [103]
0.877(13) | 2018 | Electronic Lamb shift | [102]
0.833(10) | 2019 | Electronic Lamb shift | [104]
0.831(7)(12) | 2019 |  $e^{-}p$ scattering | [105]
Table 4.1: Recent extractions of the proton charge radius from electronic Lamb
shift experiments and $e^{-}p$ scattering.

The strong radius of the proton should be smaller than its charge radius, once
it is associated with the short-range strong interaction. The electron is
electrically charged, and penetrates in the nucleus, which makes it a good
probe for the electronic structure of the proton. On the other hand, it is not
sensitive to the color fields in the nuclear interior. Low-energy scattering
of protons on nuclei provides information about both the electronic and strong
interaction. However, unlike electrons, the protons are extended, composite
objects, and interact in a much more complicated way as compared to point-like
particles. The gluonic structure of nuclei can be studied via $J/\Psi$
scattering. The $J/\Psi$ meson is formed by a charm quark-antiquark pair and
has a naturally small dipole, which should scatter on the nucleons
individually. The $J/\Psi$ meson interacts mainly by two-gluon exchange, so
that scattering off a nuclei should really probe the spatial distribution of
the color fields of the nucleons. The proton strong radius has been determined
from analysis of diffractive $J/\Psi$ photoproduction in $e^{-}p$ collisions
at the HERA collider [106]. Results suggest that the transverse strong size of
the proton is approximately 0.5 fm.

The mass radius of the proton is a quantity that has not yet been determined
experimentally. The gravitational field created by a single proton is
extremely weak, so that a direct measurement is currently very difficult.
Nevertheless, the proton mass radius was recently estimated by Kharzeev, who
defines the mass radius of the proton through the form factor of the trace of
the energy-momentum tensor of QCD in the weak gravitational field
approximation [107]. This form factor is extracted from data on
photoproduction of $J/\Psi$ and $\Upsilon$ quarkonia from the GlueX
Collaboration [108]. The extracted mass radius of the proton was
$R_{\text{m}}$ = 0.55 $\pm$ 0.03 fm.

###  4.2 Constraining the hybrid simulations: Bayesian Analysis

Hybrid simulations such as the one presented in Chapter 3 are mostly based on
effective models, which contain numerous free parameters. Each of these
parameters affects every one of the observables and, reciprocally, each
observable is in general affected by every one of the parameters. Due to the
large amount of free parameters and their complex (usually non-linear)
relation to the observables, a robust quantitative tool is necessary to
constrain the values of the model parameters and guarantee that they are
connectable to actual physical properties of the QGP. Loosely speaking, this
procedure should consist in some sort of “global fit”, in which the model is
confronted with experimental data, and all parameters and observables are
simultaneously taken into account.

The framework of Bayesian statistics provides a tool for tackling this kind of
problem, and has been employed to constrain the parameters of the models which
make up hybrid simulations of heavy-ion collisions. Such analyses are based in
a powerful identity, Bayes’ theorem, which seeks to quantitatively measure a
”degree of belief” about something, given some previously available
information. Given a proposition A and some known information B (the
evidence), Bayes’ theorem states that the probability $P(\text{A}|\text{B})$
of A being true given that evidence B is known (the posterior probability) is
given by:

$P(\text{A}|\text{B})=\frac{P(\text{B}|\text{A})P(\text{A})}{P(\text{B})}$
(4.2)

$P(\text{A})$ is the prior, and represents the initial degree of belief in A.
$P(\text{B}|\text{A})$ is the likelihood for B to be true if the proposition A
holds. $P(\text{B})$, usually referred to as Bayes evidence for information B,
acts as a normalization factor:

$P(\text{B})=\int P(\text{B}|\text{A})P(\text{A})d\text{A},$ (4.3)

so that the statement of Bayes’ theorem can be expressed as:

$P(\text{A}|\text{B})\propto P(\text{B}|\text{A})P(\text{A})$ (4.4)

It states that the posterior probability distribution of A given that B is
true is proportional to the product of the previous knowledge about A and the
likelihood of B given A.

How can this formalism be applied to constrain the parameters associated with
properties of the QGP? In this case, the hybrid model is the proposition to be
tested, and the experimental data collected in colliders constitute the
available evidence. In this context, Equation (4.4) can be written as:

$\mathcal{P}(\mathbf{x}|\mathbf{y}_{\text{exp}})\propto\mathcal{P}(\mathbf{y}_{\text{exp}}|\mathbf{x})\mathcal{P}(\mathbf{x}),$
(4.5)

where $\mathbf{x}$ is a vector that stores in its entries the model parameters
and $\mathbf{y}$ is a vector containing $m$ observation points (sometimes
called calibration data).

For the prior distribution, which expresses the knowledge about the model
parameters without any comparison with data, a uniform distribution is chosen
for each parameter $x_{i}$:

$\mathcal{P}(x_{i})\propto\small\begin{dcases}1,&\text{if}\quad
x_{i,\text{min}}\leq x_{i}\leq x_{i,\text{max}}\\\ 0,&\text{else}\\\
\end{dcases}$ (4.6)

This is also sometimes referred to as the agnostic probability distribution:
each parameter is previously believed to be equally likely to assume any value
inside an interval, which is chosen based on basic physical constraints. It is
assumed that the priors for each parameter are independent, so that the joint
prior is the product of the individual priors:

$\mathcal{P}(\mathbf{x})\propto\prod_{i}\Theta(x_{i}-x_{i,\text{min}})\Theta(x_{i}-x_{i,\text{max}})$
(4.7)

The exact form of the likelihood function is rarely known. In the case where
uncertainties are normally distributed, the most natural choice is a
multivariate Gaussian:

$\mathcal{P}(\mathbf{y}_{\text{exp}}|\mathbf{x})=\frac{1}{\sqrt{(2\pi)^{m}\text{det}\Sigma}}\exp\left(-\frac{1}{2}\Delta\mathbf{y}^{T}\Sigma^{-1}\Delta\mathbf{y}\right),$
(4.8)

where:

$\Delta\mathbf{y}=\mathbf{y}_{\text{model}}-\mathbf{y}_{\text{exp}}$ (4.9)

and $\Sigma$ is the total covariance matrix, which contains the model
uncertainties, experimental uncertainties and also correlations between
uncertainties. In the case of a multi-system analysis, the joint likelihood
function is given by the product of the systems’ individual likelihood
functions. For example, in an analysis which uses simultaneously data from Pb-
Pb and p-Pb collisions:

$\mathcal{P}(\mathbf{y}_{\text{exp}}|\mathbf{x})=\mathcal{P}(\mathbf{y}_{\text{exp}}^{\text{Pb-
Pb}}|\mathbf{x})\mathcal{P}(\mathbf{y}_{\text{exp}}^{\text{p-Pb}}|\mathbf{x})$
(4.10)

With all the ingredients on the right side of Equation (4.5) in hand,
estimation of the posterior is accomplished via Markov Chain Monte Carlo
(MCMC) algorithms [109], which allow estimation of the shape of the posterior
distribution. Such methods produce a representative sample of the posterior
distribution by performing a random walk in parameter space weighted by the
posterior. To speed up the output calculation for each parameter space point,
a Gaussian process emulator is used to substitute the model. The set of
parameters which maximizes the posterior (the mode of the posterior
distribution) is called the maximum a posteriori (MAP) set of parameters:

$\mathbf{x}_{\text{MAP}}=\underset{\mathbf{x}}{\text{arg
max}}\>\mathcal{P}(\mathbf{x}|\mathbf{y}_{\text{exp}}).$ (4.11)

Since all the priors are uniform distributions, the set of parameters that
maximizes the posterior is the one which maximizes the likelihood function. In
other words, the MAP parameters are the ones which best fit the experimental
data. The MAP value of each parameter can be visualized as the maximum of its
marginalized distribution:

$\mathcal{P}(x_{i}|\mathbf{y}_{\text{exp}})=\int
dx_{1}...dx_{i-1}dx_{i+1}...dx_{n}\mathcal{P}(\mathbf{x}|\mathbf{y}_{\text{exp}}),$
(4.12)

which is obtained by integrating out all the other parameters on the posterior
distribution (shown in the diagonal panels of Figure 4.1).

Figure 4.1: Posterior distributions for the model parameters. In blue, results
obtained by calibrating to identified particles data. In red, results obtained
by calibrating the model to charged particle yields. In the diagonal, the
marginalized posterior for each parameter is shown, and correlations between
parameters are shown in the off-diagonal. From [68].

In the numerical setup used for this work (which was presented in Chapter 3),
the simulation has a total of 13 free parameters. These are:

  * •

The entropy deposition parameter ($p$) in Equation (3.7).

  * •

Fluctuation parameter ($k$) in Equation (3.5).

  * •

Gaussian nucleon-width ($w$) in Equation (3.3).

  * •

Minimum distance between nucleons ($d_{\text{min}}$).

  * •

Normalization constant of the entropy density profile, which converts reduced
thickness into deposited entropy in Equation (3.6). This normalization factor
depends on the energy of the beam.

  * •

Pre-equilibrium time. In this work, we use a free-streaming scenario.

  * •

Shear viscosity to entropy ratio minimum value $(\eta/s)_{\text{min}}$, slope
parameter $(\eta/s)_{\text{slope}}$ and curvature parameter
$(\eta/s)_{\text{crv}}$ in Equation (3.24).

  * •

Bulk viscosity to entropy ratio maximum value $(\zeta/s)_{\text{max}}$, width
$(\zeta/s)_{\text{width}}$ and location $T_{0}$ in Equation (3.25).

  * •

Switching temperature between hydrodynamics and the hadronic phase
($T_{\text{switch}}$).

Since the main goal of this work was to study the importance of the nucleon-
width, we used a set of parameters taken from the Bayesian analysis of the
DUKE group in [110] (see Section 4.3), and changed only the Gaussian nucleon
width.

###  4.3 The nucleon-width parameter

Event-by-event fluctuations on the spatial distribution of nucleons inside the
nuclei in heavy-ion collisions have a large impact on observables [111]. Monte
Carlo Glauber based approaches incorporate this idea by assembling the nuclei
on an event-by-event basis using some nuclear distribution. In the process of
assembling two nuclei in order to generate an initial condition of a
collision, it is inevitable to somehow make some specification about the size
of their constituent nucleons.

The exact form of the transverse distribution of charge inside a nucleon is
not well known, and a Gaussian ansatz has been successfully implemented as an
effective description of the transverse profile of a nucleon moments before
the collision:

$T_{\text{nucleon}}(x,y)=\frac{1}{2\pi
w^{2}}\exp\left(-\frac{x^{2}+y^{2}}{2w^{2}}\right)$ (4.13)

The Gaussian shape seems to capture the general idea that the charges of the
nucleon (both electric and color charge) are not uniformly distributed: they
are more concentrated near its center, although the exact shape of the
distribution is not known.

Since the beginning of the 2000’s, hydrodynamic simulations of heavy-ion
collisions have used the Gaussian ansatz to model the nucleons in the initial
condition, with a width of approximately $w$ = 0.4 fm. With this Gaussian
width, one can estimate the corresponding “nucleon radius” as the root mean
square (RMS) transverse radius:

$\sqrt{\langle r^{2}\rangle}=w\sqrt{2}\approx 0.56\ \text{fm},$

which falls in the interval between the proton color radius and its electric
charge radius (See Section 4.1). Consistency between the nucleon-width
parameter $w$ and experimental estimates of the nucleon size is certainly
desirable. At this point, hybrid simulations as sophisticated as the one
presented in Chapter 3 did not exist, and the Bayesian formalism discussed in
Section 4.2 had never been applied to constrain any of the available models
parameters.

In 2016, the first Bayesian analysis was performed by the Duke group [112],
using T${}_{\text{R}}$ENTo as the initial condition generator, coupled to a
hydrodynamic simulation and a hadronic afterburner. The model was compared
with Pb-Pb data at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV, and found an optimal
value for the nucleon-width parameter of $w$ = 0.43 fm for an analysis using
only identified particle yields and $w$ = 0.49 fm when considering only
charged particle yields. This results was in agreement with the latest
experimental estimates of the proton size. Moreover, hydrodynamic simulations
using IP-Glasma initial conditions with a nucleon-width of 0.4 fm had a great
success in quantitatively describing experimental data [113].

In 2019, the Duke group performed a second global analysis [110], with
considerable differences from the first one. This time, T${}_{\text{R}}$ENTo
was used to initialize the energy density profile, and a free-streaming pre-
equilibrium dynamics was introduced between the initial condition and
hydrodynamics. There were new parametrizations for the temperature dependence
of specific shear and bulk viscosity (Equations (3.8) and (3.9)), and the
model was compared with Pb-Pb data at both $\sqrt{s_{\text{NN}}}$ = 2.76 TeV
and 5.02 TeV. This second Bayesian analysis by the Duke group returned a MAP
value of 0.956 fm for the nucleon-width parameter, which came as a striking
surprise. The corresponding RMS transverse radius is approximately
$\sqrt{\langle r^{2}\rangle}$ = 1.352 fm, which exceeds in over 60 % the
present CODATA value for the proton charge radius [101].

Figure 4.2: Marginalized posterior distributions for the nucleon-width
parameter from the 2019 DUKE analysis [68, 110] (left) and the JETSCAPE
analysis [114] (right). On the right, the result obtained with Grad viscous
corrections is shown in blue. In red, the result using Chapman-Enskog
corrections.

The fact that this result was not expected is visually evident from the
marginalized posterior distribution of the nucleon-width parameter (Left side
of Figure 4.2) in [68, 110]: the MAP value almost hits the upper bound imposed
by the prior (1 fm)!

More recently, an analogous analysis was performed by the JETSCAPE
Collaboration in 2020 [114]. The simulation chain used was in general similar
to the one employed in the second DUKE analysis. However, different
parametrizations were used for the shear and bulk viscosity, and three
different out-of-equilibrium corrections to the distribution function in the
particlization procedure were considered. The model was simultaneously
confronted with Pb-Pb data from the LHC at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV
and Au-Au at $\sqrt{s_{\text{NN}}}$ = 200 GeV data from RHIC. Surprisingly,
this analysis returned even larger values for the nucleon-width parameter, in
the range of 0.9 - 1.2 fm.

In the latest global analysis by the DUKE group, nucleon sub-structure was
enabled in T${}_{\text{R}}$ENTo, and the model was compared to data from Pb-Pb
and p-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 5.02 TeV. The analysis
returned an optimal value of $w$ = 0.88 fm for a nucleon with 6 constituents
of width $w_{c}$ = 0.43 fm each. Moreover, other similar analyses seem to
always favour such “diffus” nucleons, with $w\approx$ 0.8 fm [115, 116, 117].

Such findings seem to demand a better understanding of the results of recent
Bayesian analyses and of the role played by the nucleon-width parameter inside
the simulation as a whole. In this work, we aim not to constrain the value of
the nucleon-width, but instead seek a systematic study of its effects in the
initial condition characteristics and in final state observables. By using a
hybrid simulation chain calibrated with the parameters from [110] and changing
only the value of the nucleon-width parameter, we perform simulations for
three values of $w$. How intensely and in which way $w$ affects different
observables? Is it possible to obtain a good description of experimental data
with the calibrated chain using values of $w$ consistent with estimates of the
proton charge radius? If not, which observables(s) are poorly described in
this scenario, and why? What led the Bayesian Analyses to favour such large
values of the nucleon-width parameter? These are some of the questions that we
seek to answer.

##  Results

In this work we use the numerical setup presented in Chapter 3 (equivalent to
the one used in [110]) and perform simulations of Pb-Pb collisions at
$\sqrt{s_{\text{NN}}}$ = 2.76 TeV using the best-fit (MAP) values obtained in
the same reference, with exception of the nucleon-width parameter. For the
nucleon-width parameter, we consider three values:

  * •

$w$ = 0.5 fm - The “small nucleon”. This is closer to the values used in
simulations before the “Bayesian Era”, for which the corresponding RMS is
consistent with estimates of the proton electric charge radius and strong
force radius.

  * •

$w$ = 1.0 fm - The “medium nucleon”. Very close to the result of the Bayesian
Analysis considered in this work [110]. This should be considered to most
“well calibrated” simulation.

  * •

$w$ = 1.5 fm - The “large nucleon”. An exaggeratedly large value for the
nucleon-width parameter.

Explicitly, the set of parameters used in this work, which is the same as
[110] except for the nucleon-width, is:

Normalization (TeV) | 286.23
---|---
$p$ | 0.007
$k$ | 0.918
$w$ (fm) | 0.5, 1.0, 1.5
$d_{\text{min}}$ (fm) | 1.27
$\tau_{\text{fs}}$ (fm/c) | 1.2
$(\eta/s)_{\text{min}}$ | 0.081
$(\eta/s)_{\text{slope}}$ ($\text{GeV}^{-1}$) | 1.11
$(\eta/s)_{\text{crv}}$ | -0.48
$(\zeta/s)_{\text{max}}$ | 0.052
$(\zeta/s)_{\text{width}}$ (GeV) | 0.022
$(\zeta/s)T_{0}$ (MeV) | 183
$T_{\text{switch}}$ (MeV) | 151
Table 5.1: Set of parameters used in this work, obtained from the Bayesian
analysis in [110].

In Section 5.1, we characterize the effects of changing $w$ in the initial
condition. In this analysis, 1.000 T${}_{\text{R}}$ENTo initial conditions
were generated for 10 values of the impact parameter, with the exception of
the ellipticity fluctuations study, for which $10^{6}$ events were generated
in 8 impact parameter intervals (centrality classes). In Section 5.2, we
present the results of the full simulation and analyze observables. In this
analysis, for each value of the nucleon-width parameter, 1.000 minimum-bias
(all values of impact parameter mixed together) events were generated. The
centrality selection was made after that, based on the total entropy of the
initial conditions.

###  5.1 Characterizing the initial condition

The main objective of this work is to study the impact of the nucleon-width
parameter on final observables. We begin this task first by analyzing how the
nucleon size affects the general characteristics of the initial condition
generated by T${}_{\text{R}}$ENTo. There is a strong relation between the
initial condition characteristics and observables, so that many of the effects
of the nucleon-width parameter on final state observables can be anticipated
(at least in approximation) by analyzing its impact on the initial state.

We begin our investigation first by visualizing the initial condition. Figures
5.1 and 5.2 show examples of transverse entropy density profiles $s(x,y)$
generated by T${}_{\text{R}}$ENTo in the 0 - 5 % and 30 - 40 % centrality
classes, respectively:

Figure 5.1: Entropy density distribution in the transverse plane of Pb-Pb
collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV in the 0 - 5 % centrality
class for: $w=0.5$ fm (left), $w=1.0$ fm (center) and $w=1.5$ fm (right).

Figure 5.2: Entropy density distribution in the transverse plane of Pb-Pb
collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV in the 30 - 40 % centrality
class for: $w=0.5$ fm (left), $w=1.0$ fm (center) and $w=1.5$ fm (right).

A straightforward visual analysis of the initial conditions suggests two most
prominent effects of changing the nucleon size on the initial state:

  * •

The system formed in the collision increases with increasing nucleon size.
Furthermore, there is a decrease in system size with increasing centrality
(due to the decrease of the overlap area of the two nuclei), which is more
pronounced when using smaller nucleons.

  * •

The entropy distribution’s granularity (“lumpiness”) decreases with increasing
nucleon size. In collisions with smaller nucleons, more local high energy
density regions (hotspots) can be seen in the initial condition, which gets
smoother as the nucleon grows. When using the large nucleon, there is a single
peak in the entropy density profile, approximately in the center of the
system, which falls smoothly towards the edges of the grid.

Figure 5.3: Mean gradient (per cell) of the entropy density profile calculated
numerically as a function of impact parameter. Each point is the mean of 1.000
events, and the shaded region represents the event-by-event dispersion around
the mean.

A visibly higher granularity of the entropy density distribution should mean
stronger gradients in the initial condition. This can be verified by sweeping
over all cells of the grid an calculating the gradient numerically:

$|\vec{\nabla}s|=\sqrt{\left(\frac{\partial s}{\partial
x}\right)^{2}+\left(\frac{\partial s}{\partial y}\right)^{2}},$ (5.1)

where for each cell $(i,j)$ we have:

$\frac{\partial s}{\partial
x}\approx\frac{s(x_{i+1},y_{j})-s(x_{i},y_{j})}{x_{i+1}-x_{i}}$ (5.2)
$\frac{\partial s}{\partial
y}\approx\frac{s(x_{i},y_{j+1})-s(x_{i,}y_{j})}{y_{j+1}-y_{j}}$ (5.3)

Figure 5.3 shows the mean initial condition gradient divided by the number of
grid cells as a function of impact parameter, calculated considering 1.000
events generated using T${}_{\text{R}}$ENTo for each value of the impact
parameter. There is a strong decrease in the gradients from initial conditions
generated with $w$ = 0.5 fm to the profiles with $w$ = 1.0 fm. The difference
between the medium and the large nucleons, on the other hand, is small. There
is a decrease also as a function of impact parameter, simply due to the
decrease of the overlap area between the two nuclei, which leads to a decrease
of the formed system size. Stronger gradients in the initial condition result
in a larger mean of the transverse momentum distribution of the particles
detected at the end of the collision [118].

####  5.1.1 Participants and binary collisions

The number of participant nucleons ($N_{\text{part}}$) and the number of
binary collisions ($N_{\text{coll}}$) are the central objects of any Glauber
based calculation. How are these affected by changing the nucleon size? Figure
5.4 shows results for $\langle N_{\text{part}}\rangle$ and $\langle
N_{\text{coll}}\rangle$ as a function of impact parameter, calculated
considering 1.000 events generated using T${}_{\text{R}}$ENTo for each value
of the impact parameter.

Figure 5.4: Mean number of participant nucleons (left) and binary collisions
(right) as a function of impact parameter. Each point is the mean of 1.000
events, and the shaded region represents the event-by-event dispersion around
the mean.

The left side of Figure 5.4 shows that, at the same impact parameter, there is
a slight increase in the number of participants for collisions with more
diffuse nucleons. In general, the number of binary collisions is simply
proportional to the number of participants. On this particular case, however,
this is not always true: in central events, the small nucleons interact more
often, while in more peripheral events, there are more binary collisions when
using the the large nucleons. This is due to the fact that that the collision
probability used in T${}_{\text{R}}$ENTo:

$P_{\text{coll}}(b)=1-\exp[-\sigma_{gg}T_{AB}(b)]$ (5.4)

depends on the nucleon size. This can be seen by substituting the Gaussian
form of the thickness functions in the overlap integral:

$T_{AB}(b)=\int dxdy\,T_{A}(x-b/2)T_{B}(x+b/2)=\frac{1}{4\pi
w^{2}}\exp\left(-\frac{b^{2}}{4w^{2}}\right).$ (5.5)

Figure 5.5 shows the functional form of Equation (5.4) as a function of the
binary collision impact parameter, using a constant value for the effective
parton-parton cross-section $\sigma_{gg}=1$.

Figure 5.5: Binary collision probability as a function of impact parameter.
The collision probability is more concentrated near $b$ = 0 for the small
nucleons.

The form of the collision probability reflects the form of the nucleons
themselves: when the nucleon-width is small, two nucleons are very likely to
interact if they meet at small impact parameter, an this probability falls
quickly as the impact parameter grows. As the nucleons get larger, there is a
weaker dependence on the impact parameter. This explains why at small impact
parameter the small nucleons collide more often, while the opposite happens in
more peripheral events: the collision probability is more concentrated near
$b$ = 0 for the smaller nucleons. This also explains the slight increase in
the number of participants as the nucleon-width increases: for larger
nucleons, the collision probability is more spread, so that they are more
likely to interact (at least once) when they meet at relatively large
distances ($b$ $>$ 2 fm) from each other.

####  5.1.2 Eccentricity harmonics

A simple visual analysis of Figures 5.1 and 5.2 suggests that the geometry of
the system is sensitive to the nucleon size. More quantitatively, this should
manifest as a sensitivity of the eccentricity harmonics, as they are
calculated using precisely the transverse entropy density distribution of the
system as a weight function:

$\varepsilon_{n}=\frac{\int r^{n}e^{in\varphi}s(r,\varphi)rdrd\varphi}{\int
r^{n}s(r,\varphi)rdrd\varphi}.$ (5.6)

Figure 5.6 shows mean values of ellipticity $\varepsilon_{2}$ (left) and
triangularity $\varepsilon_{3}$ (right) of the initials conditions calculated
from Equation (5.6) as a function of impact parameter, considering 1.000
events generated using T${}_{\text{R}}$ENTo for each value of the impact
parameter.

Figure 5.6: Mean value of the eccentricity harmonics for $n=2$ (left) and
$n=3$ (right) as a function of impact parameter. Each point is the average of
1.000 events, and the shaded region represents the event-by-event dispersion
around the mean.

In general, $\varepsilon_{2}$ is small in central collisions, where the
overlap region between the two colliding nuclei is approximately circular. As
the impact parameter grows, the overlap area acquires a more elliptical shape,
and the value of $\varepsilon_{2}$ also grows. At some point, for very large
values of $b$, the overlap area between the nuclei becomes very small
(essentially created by the collision between two nucleons), and ellipticity
decreases once again. This behavior as a function of impact parameter is, in
general, true for higher order harmonics as well, although there is a weaker
dependence with the impact parameter.

Ellipticity and triangularity are strongly affected by the nucleon size: the
mean value of both eccentricity harmonics decreases as the nucleons grow
larger. The smoother entropy distributions generated when using larger
nucleons are more spatially isotropic. In particular, the triangular pattern
is almost not present in the initial conditions generated with $w$ = 1.5 fm,
and there is almost no dependence on the impact parameter.

Borrowing the concept of cumulants from the flow analysis (see Section 5.2.3),
we can calculate the analogous estimates for the eccentricity harmonics from
cumulants:

$\varepsilon_{n}\\{2\\}=\sqrt{\langle|\varepsilon_{n}|^{2}\rangle}$ (5.7)
$\varepsilon_{n}\\{4\\}=\sqrt[4]{\langle|\varepsilon_{n}|^{4}\rangle-2\langle|\varepsilon_{n}|^{2}\rangle^{2}}$
(5.8)

The ratio $\varepsilon_{n}\\{4\\}/\varepsilon_{n}\\{2\\}$ is a standard
measure of event-by-event eccentricity fluctuations. Such calculations demand
great statistics: in Figure 5.7, which shows ellipticity fluctuations as a
function of centrality, each point was calculated using $10^{6}$ events.

Figure 5.7: Event-by-event ellipticity fluctuations as a function of
centrality. Each dot was calculated using $10^{6}$ T${}_{\text{R}}$ENTo
initial conditions.

Figure 5.7 shows that although the nucleon-width strongly affects the average
value of ellipticity $\langle\varepsilon_{2}\rangle$, there is not a strong
impact on event-by-event fluctuations around the mean for more central events.
For more peripheral events (from 40 - 50 % on), there is a decrease in
fluctuations as the nucleon size increases. It seems intuitive that the degree
of fluctuation decreases for larger values of $w$ (as the initial conditions
are smoother), and that this effect is more pronounced in peripheral events.
In central events, where there are lots of binary collisions, initial
conditions are more likely to “look alike” regardless of the nucleon size,
while in peripheral collisions, where the system size is smaller, the effect
of the nucleon size is visible.

Figure 5.8 shows the results for higher order harmonics ($n=4$ and $n=5$),
calculated considering 1.000 events generated using T${}_{\text{R}}$ENTo for
each value of impact parameter. It is possible to see that the effects of
changing the nucleon-width are even more pronounced for these higher order
harmonics, for which the corresponding geometric patterns are “sharper” when
compared to ellipticity and triangularity.

Figure 5.8: Mean value of the eccentricity harmonics for $n=4$ (left) and
$n=5$ (right) as a function of impact parameter. Each point is the mean of
1.000 events, and the shaded region represents the event-by-event dispersion
around the mean.

The results for these higher order harmonics show that indeed initial
conditions generated using smaller nucleons have a more detailed and complex
geometric structure, as Figures 2.1 and 2.2 suggest.

###  5.2 Final state observables

While Section 5.1 presented results obtained by analyzing T${}_{\text{R}}$ENTo
initial conditions, this section presents results obtained after the complete
simulation, as discussed in Chapter 3. For each value of the nucleon-width
parameter, 1.000 minimum-bias (all values of impact parameter mixed together)
events were generated. The centrality selection was made after that, based on
the total entropy of the initial conditions (in the usual way, by ordering
them from the lowest to the highest entropy values, and separating them in
percentiles), which is a good predictor for the final multiplicity.

####  5.2.1 Charged particle multiplicity density at mid-rapidity

We begin the analysis of final state observables by calculating the charged
particle multiplicity density in the mid pseudo-rapidity region
($|\eta|<0.5$). Figure 5.9 shows charged particle multiplicity density in the
region $|\eta|<0.5$ as a function of centrality.

Figure 5.9: Charged particle multiplicity density in the central pseudo-
rapidity region as a function of centrality. Experimental data from the ALICE
Collaboration[41]. The error bars represent the event-by-event dispersion
around the mean.

As seen in Figure 5.9, our results show that, in more central events,
collisions with larger nucleons produce more particles, while the opposite
happens for the more peripheral classes. Although $w=1.0$ fm is essentially
the MAP value of the nucleon-width in [110], a better description of the
centrality dependence of charged particle multiplicity is actually obtained
using $w=0.5$ fm. In fact, until the 40 - 50 % class, the simulation using the
small nucleon provides a good description (within the 10 % range) of
experimental data. The simulated points were separately normalized by a
constant factor so that the simulations agree exactly with experiment (and
with each other) in the 20 - 30 % centrality interval.

####  5.2.2 Mean transverse momentum

The stronger gradients in the initial energy density profile of the system
when using smaller nucleons should have a great effect on raising the mean
transverse momentum of particles in the final state of the collision. Figure
5.10 shows the mean transverse momentum of charged pions in the mid-rapidity
region for Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV, compared
with data from ALICE. The mean transverse momentum was calculated considering
the charged pions spectra from our simulation, in the $|y|<$0.5 rapidity
interval, and with 0.1 $<p_{T}<$ 3 GeV/c (to be in agreement with the
experimental acceptance).

Figure 5.10: Mean transverse momentum of charged pions as a function of
centrality. Data from the ALICE Collaboration [47], measured in Pb-Pb
collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV. The error bars represent the
event-by-event dispersion around the mean.

For the 0 - 5 % centrality class, the mean transverse momentum is roughly the
same for the three simulations, and slightly above the experimental data
(around 10 %). The simulations using the two larger nucleons provide a good
description of data (within the 10 % range) across all centralities, but as
events get more peripheral, the simulation with $w=0.5$ fm generates particles
with too large transverse momentum. In fact, there is an increase with
centrality, which is known not to be the case.

Raising the $\langle p_{T}\rangle$ as a consequence of increasing the
granularity of the initial condition has been well established for almost 15
years [118]. Certainly this effect is present when the nucleon size is
changed, but there also seems to be more to it. To investigate this, we note
there were two main simplifications in the pre-equilibrium phase considered in
[110]:

  * •

In the pre-equilibrium phase, partons are assumed to be massless. Therefore,
free-streaming takes place with $v=c$. This is known to result in an
exaggerated large out-of-equilibrium bulk pressure when switching to
hydrodynamics, which, on its turn, is responsible for an artificially large
value of the mean transverse momentum in the final state [119].

  * •

The pre-equilibrium phase lasts the same time for all centrality classes. As
the system size decreases fore more peripheral collisions, it is reasonable to
expect that the pre-equilibrium phase should last less when compared to
central events (as is the case for the hydrodynamic evolution, for example).
The use of a constant pre-equilibrium time results in an artificial increase
of the mean transverse momentum with centrality [120], as the violent pre-
hydrodynamic expansion lasts longer than it should.

The simulation using the small nucleon is the most sensitive to this second
simplification in the free-streaming phase, due to the fact that the system
formed is considerably smaller when compared to using the two larger nucleon-
width values (Figures 5.1 and 5.2). The first simplification, on the other
hand, affects all centrality classes. Figure 5.10 shows that, in this
scenario, where the effects of the stronger gradients due to using smaller
nucleons are combined with the effects of these two simplifications made in
the pre-hydrodyamic stage of the simulation, the simulation with $w$ = 0.5 fm
produces wrong results (exceeding up to 60 % the experimental data).

####  5.2.3 Anisotropic flow

Calculating the flow harmonics explicitly as an expectation value like in
Equation (2.6) requires event-by-event knowledge of the reaction plane angle
$\Psi_{\text{RP}}$, which is not accessible experimentally. As an alternative,
flow harmonics are usually calculated using multiparticle azimuthal
correlations [121, 122], sometimes also referred to as the cumulants method.
The idea behind the method is that, as collectivity is related to correlations
between particles in momentum space, it should be possible to extract flow
from correlation functions.

Let $\langle k\rangle$ denote the azimuthal correlation function of $k$
particles in a single event. In this notation, the 2-particle and 4-particle
correlation functions are given by [123]:

$\langle 2\rangle=e^{in(\phi_{1}-\phi_{2})}=\frac{1}{P_{M,2}}\sum_{i\neq
j}^{M}e^{in(\phi_{i}-\phi_{j})}$ (5.9) $\langle
4\rangle=e^{in(\phi_{1}+\phi_{2}-\phi_{3}-\phi_{4})}=\frac{1}{P_{M,4}}\sum_{i\neq
j\neq k\neq l}^{M}e^{in(\phi_{i}+\phi_{j}-\phi{k}-\phi_{l})},$ (5.10)

where $M$ is the number of charged particles in the event, $n$ is the order of
the harmonic, $\phi$ is the azimuthal angle, and:

$P_{M,k}=\frac{M!}{(M-k)!}$ (5.11)

is the number of permutations of $k$ particles in an event of multiplicity
$M$, so that:

$P_{M,2}=M(M-1)$ (5.12) $P_{M,4}=M(M-1)(M-2)(M-3).$ (5.13)

Inside a given centrality bin, the $k$-particle correlation function can be
written as $\langle\langle k\rangle\rangle$, where the outer brackets denote
an average over all events inside the centrality class.

How can the flow coefficients be estimated from these correlation functions?
We begin by rewriting $\langle\langle 2\rangle\rangle$, adding an subtracting
$\Psi_{\text{RP}}$ on the exponential:

$\langle\langle 2\rangle\rangle=\langle\langle
e^{in[(\phi_{1}-\Psi_{\text{RP}})-(\phi_{2}-\Psi_{\text{RP}})]}\rangle\rangle.$
(5.14)

Now, assuming that the only correlation between $\phi_{1}$ and $\phi_{2}$ is
due to collective flow, the inner average can be factorized in a product:

$\langle\langle 2\rangle\rangle\approx\langle\langle
e^{in(\phi_{1}-\Psi_{\text{RP}})}\rangle\langle
e^{-in(\phi_{2}-\Psi_{\text{RP}})}\rangle\rangle=\langle v_{n}^{2}\rangle,$
(5.15)

and analogously, $\langle\langle 4\rangle\rangle\approx\langle
v_{n}^{4}\rangle$. This factorization is not exact, so that the equality is
only approximate. Denoting by $c_{n}\\{k\\}$ the $n$th-order cumulant from
correlations of $k$ particles:

$c_{n}\\{2\\}=\langle\langle 2\rangle\rangle$ (5.16)
$c_{n}\\{2\\}=\langle\langle 4\rangle\rangle-2\langle\langle 2\rangle\rangle,$
(5.17)

so that the flow coefficients $v_{n}$ estimated from the 2 and 4-particle
cumulants are given by:

$v_{n}\\{2\\}=\sqrt{c_{n}\\{2\\}}$ (5.18)
$v_{n}\\{4\\}=\sqrt[4]{-c_{n}\\{4\\}}$ (5.19)

As the calculation of correlation functions is computationally expensive, in
practice the flow harmonics are usually extracted using the $Q$-vectors, which
are defined as [123]:

$Q_{n}=\sum_{i=1}^{M}e^{in\phi_{i}}.$ (5.20)

Using the $Q$-vector, the single-event correlation functions can be expressed
analytically. For example:

$|Q_{n}|^{2}=\sum_{i,j=1}^{M}e^{in(\phi_{i}-\phi_{j})}=M+\sum_{i\neq
j}^{M}e^{in(\phi_{i}-\phi_{j})},$ (5.21)

so that:

$\langle 2\rangle=\frac{|Q_{n}|^{2}-M}{M(M-1)},$ (5.22)

and more complicated, yet analogous expressions can be derived for higher
order correlations.

Figure 5.11 shows the $p_{T}$-integrated elliptic flow calculated using two-
particle correlations as a function of centrality, compared with data from the
ALICE Collaboration [59]. Each dot is the mean value of the events in that
centrality class, and the error bars represent the event-by-event dispersion
around the mean.

Figure 5.11: Integrated elliptic flow calculated using two-particle
correlations as a function of centrality. Data from [59]. The error bars
represent the event-by-event dispersion around the mean.

There is a pronounced effect of the nucleon-width parameter on integrated
flow, which decreases as the nucleon size increases. This follows the same
trend as ellipticity in the initial condition, and is simply a consequence of
the mapping between initial state geometry and final state momentum anisotropy
mentioned in Section 2.5.2:

$v_{n}=f(\varepsilon_{n})+\delta_{n}$ (5.23)

Using $w=0.5$ fm provides a better description of data description of $v_{2}$
data (within the 10 % range), while the simulations with the larger nucleons
do not produce enough elliptic flow.

Figure 5.12: Integrated triangular flow calculated using two-particle
correlations as a function of centrality. Data from [59]. The error bars
represent the event-by-event dispersion around the mean.

In Figure 5.12, which shows the integrated triangular flow, this effect is
even more pronounced: the simulations with the larger nucleons do not produce
enough triangular flow, while the simulation with the small nucleons is able
to reproduce experimental data within the 10 % range. Once again, the mapping
between initial state eccentricity and anisotropic flow is present.

####  5.2.4 Blast Wave analysis

So far the only information we have extracted from the $p_{T}$-spectra was the
mean value. What else can we learn about the effects of the nucleon-width
parameter in the transverse momentum distributions? A hydrodynamic-inspired
phenomenological model widely used to characterize $p_{T}$-spectra is the
Blast Wave (BW) model [54]. In the BW model, particles are emitted from the
“source” with an approximately thermal distribution, at the same time that
they are boosted in the radial direction by the transverse expansion of the
system. Starting from this picture, the following analytical expression is
obtained:

$\frac{1}{p_{T}}\frac{dN}{dp_{T}}\propto\int_{0}^{R}rdrm_{T}I_{0}\left(\frac{p_{T}\sinh\rho}{T_{\text{kin}}}\right)K_{1}\left(\frac{m_{T}\cosh\rho}{T_{\text{kin}}}\right),$
(5.24)

where:

$\rho=\tanh^{-1}\beta_{T}=\tanh^{-1}\left(\left(\frac{r}{R}\right)^{n}\beta_{s}\right),$
(5.25)

and $m_{T}=\sqrt{m^{2}+p_{T}^{2}}$ is the transverse mass. $R$ is the system
radius, and $I_{0}$ and $K_{1}$ are modified Bessel functions. This expression
contains three free parameters:

  * •

$T_{\text{kin}}$ \- Kinetic freeze-out temperature

  * •

$\beta_{T}$ \- Transverse expansion velocity ($\beta_{s}$ is the surface
velocity)

  * •

$n$ \- Velocity profile exponent

which can be extracted by fitting $p_{T}$-spectra with the expression in
Equation (5.24). This is usually done in the form of a combined fit, in which
the $p_{T}$-spectra of three particle species are fitted simultaneously,
providing a more stable fit and more robust insights. We follow the
prescription used by the ALICE Collaboration [47], in which the
$p_{T}$-spectra of charged pions, kaons and protons were fitted in the
following $p_{T}$ intervals:

  * •

$\pi^{+}$ \+ $\pi^{-}$ : 0.5 - 1.0 GeV/c

  * •

$K^{+}$ \+ $K^{-}$ : 0.2 - 1.5 GeV/c

  * •

$p$ \+ $\bar{p}$ : 0.3 - 3 GeV/c

Figure 5.13: Blast Wave combined fit parameters. Each point corresponds to a
centrality class. Central events are on the right, and peripheral events on
the left. Data from [47].

The standard visualization of the results of a BW analysis is such as in
Figure 5.13: the transverse expansion velocity is displayed in the horizontal
axis and the kinetic freeze-out temperature in the vertical axis. It is clear
that the simulation with $w=0.5$ fm provides an insufficient description of
experimental data, occupying a small region on the right side of Figure 5.13,
while the experimental data is well spread across the plane. This suggests
that, besides the mean transverse momentum, the simulation with the smaller
nucleon size seems to be problematic in what concerns the description of the
shape of $p_{T}$-spectra as a whole. Also, we note that for very peripheral
collisions, the simulation using the small nucleons changes the trend of the
transverse expansion velocity, following what was already observed in the
centrality dependence of the mean transverse momentum of charged pions.

Based on our results, we conclude that the only experimental observable which
constitutes the core of the Bayesian Analyses that is actually poorly
described when using the small nucleons is the mean transverse momentum, which
is too large. We associate this to two simplifications made in the pre-
hydrodynamic stage of the simulation. In this sense, the large values of the
nucleon-width parameter might have been an artefact of the Bayesian Analyses,
necessary to lower the mean transverse momentum, which was being artificially
raised. In addition, the combined BW fit suggests that other issues are
present in the transverse momentum distributions.

##  Conclusions

This work was motivated by recent results in Bayesian Analyses used to
constrain the parameters of simulations of heavy-ion collisions, in which the
nucleon-width parameter assumed surprisingly large values, inconsistent with
current measurements of the proton size. In this work, we have performed full
hybrid simulations of Pb-Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV
using a state-of-the-art simulation chain for three different values of the
nucleon-width parameter, and systematically analyzed its effects on the
initial state characteristics and observables.

The nucleon size strongly affects the initial condition characteristics:
collisions generated using smaller nucleons have stronger entropy gradients
and are more spatially anisotropic.

Now we are ready to answer the questions that were asked in Section 4.3:

  * •

How intensely and in which way $w$ affects different observables? : The
nucleon-width visibly affects all the considered observables. As the nucleon
size increases, there is a decrease in the mean transverse momentum of final
particles (due to the weaker gradients in the initial conditions) and a
decrease on elliptic and triangular flow (due to the decrease in the initial
condition ellipticity and triangularity).

  * •

Is it possible to obtain a good description of experimental data with the
calibrated chain using values of $w$ consistent with estimates of the proton
charge radius? If not, which observables(s) are poorly described in this
scenario, and why? : Overall, it is possible to obtain a good description
(within the 10 % range) of experimental data using $w=0.5$ fm, with the
exception of the mean transverse momentum. In this case, the simulation with
the smaller nucleons result in values of $\langle p_{T}\rangle$ above data,
and a wrong centrality dependence. This was linked to simplifications made in
the pre-equilibrium phase which artificially increase the mean transverse
momentum.

  * •

What led the Bayesian Analyses to favour such large values of the nucleon-
width parameter? : Our results suggest that the large value of the nucleon-
width parameter returned by the Bayesian Analysis in [110] might have been an
“attempt” to lower the $\langle p_{T}\rangle$, which was being artificially
raised. However, it is not possible to state that this is the only effect
present, and totally accounts for the enlargement of the nucleon-width
parameter MAP value.

A combined Blast Wave analysis suggests that the shape of the $p_{T}$-spectra
as a whole is poorly described by the simulation with $w=$ 0.5 fm.

This work highlights the importance of performing computational simulations
based on physical models in order to extract information about the collision
and allow for comparisons with experimental data collected in heavy-ion
colliders. These simulations are computationally demanding and present a
complex relation between its parameters and observables, so that systematic
studies such as the ones presented in this work are necessary.

Finally, this work complements the recent findings in [124, 125], which show
that indeed the nucleon-width parameter should assume values smaller than 0.7
fm.

## References

  * [1] Mark Thomson  “Modern Particle Physics”  Cambridge University Press, 2013
  * [2] MissMJ  “Standard Model of Elementary Particles” https://commons.wikimedia.org/wiki/File:Standard_Model_of_Elementary_Particles.svg
  * [3] Matthew D. Schwartz  “Quantum Field Theory and the Standard Model”  Cambridge University Press, 2014
  * [4] Julian Schwinger  “Quantum Electrodynamics. I. A Covariant Formulation”  In _Phys. Rev._ 74 American Physical Society, 1948, pp. 1439–1461  DOI: 10.1103/PhysRev.74.1439
  * [5] Peter W. Higgs  “Broken Symmetries and the Masses of Gauge Bosons”  In _Phys. Rev. Lett._ 13 American Physical Society, 1964, pp. 508–509  DOI: 10.1103/PhysRevLett.13.508
  * [6] “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC”  In _Physics Letters B_ 716.1, 2012, pp. 1–29  DOI: https://doi.org/10.1016/j.physletb.2012.08.020
  * [7] Serguei Chatrchyan  “Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC”  In _Phys. Lett. B_ 716, 2012, pp. 30–61  DOI: 10.1016/j.physletb.2012.08.021
  * [8] Serguei Chatrchyan  “Observation of a New Boson with Mass Near 125 GeV in $pp$ Collisions at $\sqrt{s}$ = 7 and 8 TeV”  In _JHEP_ 06, 2013, pp. 081  DOI: 10.1007/JHEP06(2013)081
  * [9] F Wilczek  “Quantum Chromodynamics: The Modern Theory of the Strong Interaction”  In _Annual Review of Nuclear and Particle Science_ 32.1, 1982, pp. 177–209  DOI: 10.1146/annurev.ns.32.120182.001141
  * [10] Gerhard Ecker  “Quantum chromodynamics”  In _2005 European School of High-Energy Physics_ , 2006  arXiv:hep-ph/0604165
  * [11] Alexandre Deur, Stanley J. Brodsky and Guy F. de Téramond  “The QCD running coupling”  In _Progress in Particle and Nuclear Physics_ 90, 2016, pp. 1–74  DOI: https://doi.org/10.1016/j.ppnp.2016.04.003
  * [12] Particle Data Group  “Review of Particle Physics” 083C01  In _Progress of Theoretical and Experimental Physics_ 2020.8, 2020  DOI: 10.1093/ptep/ptaa104
  * [13] David J. Gross and Frank Wilczek  “Ultraviolet Behavior of Non-Abelian Gauge Theories”  In _Phys. Rev. Lett._ 30 American Physical Society, 1973, pp. 1343–1346  DOI: 10.1103/PhysRevLett.30.1343
  * [14] H David Politzer  “Asymptotic freedom: An approach to strong interactions”  In _Physics Reports_ 14.4, 1974, pp. 129–180  DOI: https://doi.org/10.1016/0370-1573(74)90014-3
  * [15] M. Gell-Mann  “A schematic model of baryons and mesons”  In _Physics Letters_ 8.3, 1964, pp. 214–215  DOI: https://doi.org/10.1016/S0031-9163(64)92001-3
  * [16] R.. Feynman  “The behavior of hadron collisions at extreme energies”  In _Conf. Proc. C_ 690905, 1969, pp. 237–258
  * [17] H. Fritzsch, Murray Gell-Mann and H. Leutwyler  “Advantages of the Color Octet Gluon Picture”  In _Phys. Lett. B_ 47, 1973, pp. 365–368  DOI: 10.1016/0370-2693(73)90625-4
  * [18] Y. Aoki et al.  “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics”  In _Nature_ 443, 2006, pp. 675–678  DOI: 10.1038/nature05120
  * [19] Mark G. Alford, Krishna Rajagopal and Frank Wilczek  “QCD at finite baryon density: Nucleon droplets and color superconductivity”  In _Phys. Lett. B_ 422, 1998, pp. 247–256  DOI: 10.1016/S0370-2693(98)00051-3
  * [20] N. Cabibbo and G. Parisi  “Exponential hadronic spectrum and quark liberation”  In _Physics Letters B_ 59.1, 1975, pp. 67–69  DOI: https://doi.org/10.1016/0370-2693(75)90158-6
  * [21] Antonin Maire  “Introduction to the Quark-Gluon Plasma session in RJC 2014”, 2015
  * [22] David Tlusty  “The RHIC Beam Energy Scan Phase II: Physics and Upgrades”  In _13th Conference on the Intersections of Particle and Nuclear Physics_ , 2018  arXiv:1810.04767 [nucl-ex]
  * [23] Heng-Tong Ding, Frithjof Karsch and Swagato Mukherjee  “Thermodynamics of strong-interaction matter from Lattice QCD”  In _International Journal of Modern Physics E_ 24, 2015  DOI: 10.1142/S0218301315300076
  * [24] Chin-Hao Chen  “Studying the Early Universe via Quark-Gluon Plasma” Proceedings of the 9th International Symposium on Cosmology and Particle Astrophysics  In _Nuclear Physics B - Proceedings Supplements_ 246-247, 2014, pp. 38–41  DOI: https://doi.org/10.1016/j.nuclphysbps.2013.10.063
  * [25] Reinhard Stock  “Relativistic nucleus nucleus collisions: From the BEVALAC to RHIC.”  In _J. Phys. G_ 30, 2004, pp. S633–S648  DOI: 10.1088/0954-3899/30/8/001
  * [26] A M Baldin et al.  “Beams of Highly Charged Ions at the Dubna Syncrophasotron”  In _Physica Scripta_ T3 IOP Publishing, 1983, pp. 43–44  DOI: 10.1088/0031-8949/1983/t3/009
  * [27] Gordon Baym  “RHIC: From dreams to beams in two decades”  In _Nucl. Phys. A_ 698, 2002, pp. XXIII–XXXII  DOI: 10.1016/S0375-9474(01)01342-2
  * [28] B.. Back  “The PHOBOS perspective on discoveries at RHIC”  In _Nucl. Phys. A_ 757, 2005, pp. 28–101  DOI: 10.1016/j.nuclphysa.2005.03.084
  * [29] I. Arsene  “Quark gluon plasma and color glass condensate at RHIC? The Perspective from the BRAHMS experiment”  In _Nucl. Phys. A_ 757, 2005, pp. 1–27  DOI: 10.1016/j.nuclphysa.2005.02.130
  * [30] John Adams  “Experimental and theoretical challenges in the search for the quark gluon plasma: The STAR Collaboration’s critical assessment of the evidence from RHIC collisions”  In _Nucl. Phys. A_ 757, 2005, pp. 102–183  DOI: 10.1016/j.nuclphysa.2005.03.085
  * [31] K. Adcox  “Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration”  In _Nucl. Phys. A_ 757, 2005, pp. 184–283  DOI: 10.1016/j.nuclphysa.2005.03.086
  * [32] Jaroslav Adam  “Centrality dependence of the charged-particle multiplicity density at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV”  In _Phys. Rev. Lett._ 116.22, 2016, pp. 222302  DOI: 10.1103/PhysRevLett.116.222302
  * [33] Shreyasi Acharya  “Underlying Event properties in pp collisions at $\sqrt{s}$ = 13 TeV”  In _JHEP_ 04, 2020, pp. 192  DOI: 10.1007/JHEP04(2020)192
  * [34] Wit Busza, Krishna Rajagopal and Wilke Schee  “Heavy Ion Collisions: The Big Picture and the Big Questions”  In _Annual Review of Nuclear and Particle Science_ 68.1, 2018, pp. 339–376  DOI: 10.1146/annurev-nucl-101917-020852
  * [35] Megan Connors, Christine Nattrass, Rosi Reed and Sevil Salur  “Jet measurements in heavy ion physics”  In _Rev. Mod. Phys._ 90, 2018, pp. 025005  DOI: 10.1103/RevModPhys.90.025005
  * [36] Jasmine Brewer  “Jets as a probe of the quark-gluon plasma”  In _PoS_ HardProbes2020, 2021, pp. 012  DOI: 10.22323/1.387.0012
  * [37] Ulrich Heinz and Raimond Snellings  “Collective flow and viscosity in relativistic heavy-ion collisions”  In _Ann. Rev. Nucl. Part. Sci._ 63, 2013, pp. 123–151  DOI: 10.1146/annurev-nucl-102212-170540
  * [38] Costanza Cavicchioli  “Development and Commissioning of the Pixel Trigger System for the ALICE Experiment at the CERN Large Hadron Collider”, 2010
  * [39] Barbara Betz  “Jet Propagation and Mach-Cone Formation in (3+1)-dimensional Ideal Hydrodynamics”, 2009
  * [40] “Evidence for light-by-light scattering in heavy-ion collisions with the ATLAS detector at the LHC”  In _Nature Physics_ 13, 2017, pp. 852–858  DOI: https://doi.org/10.1038/nphys4208
  * [41] Kenneth Aamodt  “Centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV”  In _Phys. Rev. Lett._ 106, 2011, pp. 032301  DOI: 10.1103/PhysRevLett.106.032301
  * [42] Fritz W. Bopp, R. Engel, J. Ranft and S. Roesler  “Inclusive distributions at the LHC as predicted from the DPMJET-III model with chain fusion”  In _J. Phys. G_ 35.5, 2008, pp. 054001.20  arXiv:0706.3875 [hep-ph]
  * [43] Wei-Tian Deng, Xin-Nian Wang and Rong Xu  “Hadron production in p+p, p+Pb, and Pb+Pb collisions with the HIJING 2.0 model at energies available at the CERN Large Hadron Collider”  In _Phys. Rev. C_ 83, 2011, pp. 014915  DOI: 10.1103/PhysRevC.83.014915
  * [44] Javier L. ALbacete and Adrian Dumitru  “A model for gluon production in heavy-ion collisions at the LHC with rcBK unintegrated gluon densities”, 2010  arXiv:1011.5161 [hep-ph]
  * [45] Dmitri Kharzeev, Eugene Levin and Marzia Nardi  “Color glass condensate at the LHC: Hadron multiplicities in pp, pA and AA collisions”  In _Nucl. Phys. A_ 747, 2005, pp. 609–629  DOI: 10.1016/j.nuclphysa.2004.10.018
  * [46] Nestor Armesto, Carlos A. Salgado and Urs Achim Wiedemann  “Relating high-energy lepton-hadron, proton-nucleus and nucleus-nucleus collisions through geometric scaling”  In _Phys. Rev. Lett._ 94, 2005, pp. 022002  DOI: 10.1103/PhysRevLett.94.022002
  * [47] Betty Abelev  “Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV”  In _Phys. Rev. C_ 88, 2013, pp. 044910  DOI: 10.1103/PhysRevC.88.044910
  * [48] Qi Wang, Pei-Pin Yang and Fu-Hu Liu  “Comparing a few distributions of transverse momenta in high energy collisions”  In _Results Phys._ 12, 2019, pp. 259–267  DOI: 10.1016/j.rinp.2018.11.067
  * [49] Rohit Gupta and Satyajit Jena  “Model comparison of the transverse momentum spectra of charged hadrons produced in $PbPb$ collision at $\sqrt{s_{NN}}=5.02$ TeV”, 2021  arXiv:2103.13104 [hep-ph]
  * [50] Constantino Tsallis  “Possible generalization of Boltzmann-Gibbs statistics”  In _Journal of Statistical Physics_ 52, 1988, pp. 479–487  DOI: 10.1007/BF01016429
  * [51] R. Hagedorn  “Statistical thermodynamics of strong interactions at high-energies”  In _Nuovo Cim. Suppl._ 3, 1965, pp. 147–186
  * [52] R. Hagedorn and J. Ranft  “Statistical thermodynamics of strong interactions at high-energies. 2. Momentum spectra of particles produced in pp-collisions”  In _Nuovo Cim. Suppl._ 6, 1968, pp. 169–354
  * [53] Rohit Gupta and Satyajit Jena  “A unified formalism to study $soft$ as well as $hard$ part of the transverse momentum spectra”  In _24th DAE-BRNS High Energy Physics Symposium_ , 2021  arXiv:2103.13896 [hep-ph]
  * [54] Philip J. Siemens and John O. Rasmussen  “Evidence for a Blast Wave from Compressed Nuclear Matter”  In _Phys. Rev. Lett._ 42 American Physical Society, 1979, pp. 880–883  DOI: 10.1103/PhysRevLett.42.880
  * [55] Peter F. Kolb and Ulrich W. Heinz  “Hydrodynamic description of ultrarelativistic heavy ion collisions”, 2003, pp. 634–714  arXiv:nucl-th/0305084
  * [56] S. Voloshin and Y. Zhang  “Flow study in relativistic nuclear collisions by Fourier expansion of Azimuthal particle distributions”  In _Z. Phys. C_ 70, 1996, pp. 665–672  DOI: 10.1007/s002880050141
  * [57] Sergei A. Voloshin, Arthur M. Poskanzer and Raimond Snellings  “Collective phenomena in non-central nuclear collisions”  In _Landolt-Bornstein_ 23, 2010, pp. 293–333  DOI: 10.1007/978-3-642-01539-7˙10
  * [58] Raimond Snellings  “Elliptic Flow: A Brief Review”  In _New J. Phys._ 13, 2011, pp. 055008  DOI: 10.1088/1367-2630/13/5/055008
  * [59] K. Aamodt  “Higher harmonic anisotropic flow measurements of charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV”  In _Phys. Rev. Lett._ 107, 2011, pp. 032301  DOI: 10.1103/PhysRevLett.107.032301
  * [60] Michael L. Miller, Klaus Reygers, Stephen J. Sanders and Peter Steinberg  “Glauber modeling in high energy nuclear collisions”  In _Ann. Rev. Nucl. Part. Sci._ 57, 2007, pp. 205–243  DOI: 10.1146/annurev.nucl.57.090506.123020
  * [61] David d’Enterria and Constantin Loizides  “Progress in the Glauber Model at Collider Energies”  In _Ann. Rev. Nucl. Part. Sci._ 71, 2021, pp. 315–344  DOI: 10.1146/annurev-nucl-102419-060007
  * [62] A. Adare  “Transverse energy production and charged-particle multiplicity at midrapidity in various systems from $\sqrt{s_{NN}}=7.7$ to 200 GeV”  In _Phys. Rev. C_ 93.2, 2016, pp. 024901  DOI: 10.1103/PhysRevC.93.024901
  * [63] Betty Abelev  “Measurement of inelastic, single- and double-diffraction cross sections in proton–proton collisions at the LHC with ALICE”  In _Eur. Phys. J. C_ 73.6, 2013, pp. 2456  DOI: 10.1140/epjc/s10052-013-2456-0
  * [64] Betty Abelev  “Pseudorapidity density of charged particles in $p$ \+ Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV”  In _Phys. Rev. Lett._ 110.3, 2013, pp. 032301  DOI: 10.1103/PhysRevLett.110.032301
  * [65] A. Bialas, M. Bleszynski and W. Czyz  “Multiplicity Distributions in Nucleus-Nucleus Collisions at High-Energies”  In _Nucl. Phys. B_ 111, 1976, pp. 461–476  DOI: 10.1016/0550-3213(76)90329-1
  * [66] Martin M. Block and Francis Halzen  “New experimental evidence that the proton develops asymptotically into a black disk”  In _Phys. Rev. D_ 86, 2012, pp. 051504  DOI: 10.1103/PhysRevD.86.051504
  * [67] Constantin Loizides  “Glauber modeling of high-energy nuclear collisions at the subnucleon level”  In _Phys. Rev. C_ 94.2, 2016, pp. 024914  DOI: 10.1103/PhysRevC.94.024914
  * [68] Jonah E. Bernhard  “Bayesian parameter estimation for relativistic heavy-ion collisions”, 2018  arXiv:1804.06469 [nucl-th]
  * [69] Jacquelyn Noronha-Hostler, Li Yan, Fernando G. Gardim and Jean-Yves Ollitrault  “Linear and cubic response to the initial eccentricity in heavy-ion collisions”  In _Phys. Rev. C_ 93.1, 2016, pp. 014909  DOI: 10.1103/PhysRevC.93.014909
  * [70] Fernando G. Gardim, Frederique Grassi, Matthew Luzum and Jean-Yves Ollitrault  “Mapping the hydrodynamic response to the initial geometry in heavy-ion collisions”  In _Phys. Rev. C_ 85, 2012, pp. 024908  DOI: 10.1103/PhysRevC.85.024908
  * [71] Maurı́cio Hippert et al.  “Momentum-dependent flow fluctuations as a hydrodynamic response to initial geometry”  In _Nucl. Phys. A_ 1005, 2021, pp. 121982  DOI: 10.1016/j.nuclphysa.2020.121982