# Multi-messenger heavy-ion collision physics

Charles Gale  Department of Physics, McGill University, 3600 University
Street, Montreal, QC, Canada H3A 2T8  Jean-François Paquet  Department of
Physics, Duke University, Durham, NC 27708, USA  Björn Schenke  Physics
Department, Brookhaven National Laboratory, Upton, NY 11973, USA  Chun Shen
Department of Physics & Astronomy, Wayne State University, Detroit, MI 48201,
USA  RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, NY
11973, USA

###### Abstract

This work studies the production of direct photons in relativistic nuclear
collisions, along with the production of hadrons. Radiation from the very
first instants to the final moments of the evolution is included. The hybrid
model used here describes all stages of relativistic heavy-ion collisions.
Chronologically, those are an initial state reflecting the collision of nuclei
described within the Color Glass Condensate effective theory; a pre-
equilibrium phase based on non-equilibrium linear response; relativistic
viscous hydrodynamics, and a hadronic afterburner. The effect of the pre-
equilibrium phase on both photonic and hadronic observables is highlighted for
the first time. The potential of photon observables – spectrum, differential
elliptic and triangular flow – to reveal the chemical equilibration time is
studied. Finally, we consider “small collision systems”, including
proton+nucleus collisions and collisions of light nuclei, as probed by
hadronic and electromagnetic observables. We demonstrate how photon production
can signal the formation of quark-gluon plasma in such small systems.

###### pacs:

Valid PACS appear here

##  I Introduction

A great deal is now known about Quantum Chromodynamics (QCD), the theory of
the nuclear strong interaction; it is a gauge theory whose degrees of freedom
are quarks and gluons Yang and Mills (1954). However, its non-linear nature
makes it challenging to solve explicitly, except in some special cases. One of
these involves interactions at large momentum transfers where, because of
asymptotic freedom, the smallness of the strong coupling constant admits a
perturbative expansion. In that limit, perturbative QCD correctly predicts and
interprets a wide class of experimental results [See; forexample; ][;
andreferencestherein.]RevModPhys.67.157. Where perturbative techniques fail, a
non-perturbative realization of QCD is accomplished by discretizing the theory
on a space-time lattice and by solving numerically the path integral which
underpins the theory. Lattice QCD has had remarkable success in its analyses
of hadronic observables, both at zero Bazavov  _et al._ (2010) and finite
temperature Ratti (2018).

Studying the collisions of strongly interacting systems at relativistic
energies is currently the only practical means of investigating finite-
temperature and finite-density QCD in terrestrial settings. This is one of the
goals of the relativistic heavy-ion program, which is pursued at several large
accelerator facilities around the world. In addition to furthering knowledge
about QCD under extreme conditions, the heavy-ion program informs our
knowledge of strongly interacting systems out of equilibrium. Indeed, one of
the breakthroughs in the modeling of relativistic nuclear collisions is the
remarkable success of relativistic viscous hydrodynamics, and the realization
that the collective behavior observed in heavy-ion final states can reveal
departures from equilibrium typically characterized by transport coefficients
such as the shear viscosity to entropy density, $\eta/s$, and the bulk
viscosity to entropy density ratio, $\zeta/s$ Gale  _et al._ (2013a).
Theoretical advances such as the ones just described coupled with a vigorous
experimental program have been successful in making concrete an early
prediction of QCD: the formation of a quark-gluon plasma (QGP) Baym (2016).
The physics of the QGP has now entered a phase of characterization, and most
of the contemporary efforts in the field aim to further this goal.

Several of the observables considered in the analysis of the strongly
interacting matter in and out of equilibrium are hadrons, and their final
state behavior is often quantified using what has become known as flow
analysis Heinz and Snellings (2013). Quantitatively, in addition to familiar
transverse momentum spectra at a given rapidity $y$, the momentum anisotropy
is typically quantified in terms of Fourier coefficients as in

$\displaystyle E\frac{dN}{d^{3}p}=\left(\frac{1}{2\pi
p_{T}}\frac{dN}{dp_{T}dy}\right)\left[1+2\sum_{n=1}^{\infty}v_{n}\cos\left(n(\phi-\Psi_{n})\right)\right].$

The azimuthal angle $\Psi_{n}$ is a measure of the angle where the measured
particles are dominantly produced in the transverse plane. The coefficients
$v_{n}$ are the so-called flow coefficients.

In this work, using the same hybrid approach, we report on calculations of
hadronic and electromagnetic observables, focusing on real photons. The
ability to tackle hadronic and electromagnetic observables simultaneously
within a single theoretical framework is worthy of emphasis and it is one of
the cornerstones of our work: this is multi-messenger heavy-ion physics.
Furthermore, in addition to the viscous fluid-dynamical phase that has become
the workhorse of relativistic heavy-ion modeling, this work explores the
eventually observable consequences of a “pre-hydrodynamics” phase. The
analyses pursued here use KøMPøST Kurkela  _et al._ (2019a, b), an approach
described in more detail in the following section.

This paper is organized as follows: the next section discusses the hybrid
model used to describe the production of hadrons and of photons created during
the collision of heavy ions at RHIC and LHC energies. We then show results
obtained for charged hadrons and for the spectra and momentum anisotropies of
real photons. In the latter case, the role of chemical equilibrium is examined
in some detail. Section III ends with a discussion of “small systems”, where
predictions are made for photons produced in collisions of p+Pb and O+O at the
LHC. The paper concludes in Sec. IV.

##  II Modeling the time evolution of relativistic heavy-ion collisions

In this section, the spacetime modeling of relativistic nuclear collisions is
discussed, from the initial stage to the final configurations measured by the
detectors.

###  II.1 The pre-hydrodynamical phase of relativistic heavy-ion collisions

The theoretical modeling of the very early nuclear medium is achieved using
the IP-Glasma model. This approach is by now well-established and details are
to be found in the original publications Schenke  _et al._ (2012a);
*Schenke:2012fw. The IP-Glasma initialization is based on the collision of two
color glass condensates, which results in the generation of gluon fields,
whose time-evolution is accomplished by solving classical Yang-Mills equations
for a proper time duration which we shall label $\tau_{\rm IPG}^{\rm X}$ in
this work, where “X” is the phase immediately following the IP-Glasma time-
evolution. The phase labeled by “X” will be one of two, as detailed below.

One of the goals of this work is to study the effect of considering an out-of-
equilibrium era between IP-Glasma and viscous fluid dynamics. The pre-
hydrodynamic evolution scheme considered here is that of KøMPøST Kurkela  _et
al._ (2019a, b). In a given time interval $\Delta\tau$, causality sets a limit
on the size of the region that can affect conditions at a given spacetime
point, $(\tau,{\rm\bf x})$. Taking an average over the causally-connected
region, one can split the energy-momentum tensor into a locally homogeneous
background in the transverse plane, and perturbations

$\displaystyle T^{\mu\nu}\left(\tau,{\bf
x}\right)=\bar{T}^{\mu\nu}\left(\tau,{\rm\bf x}\right)+\delta
T^{\mu\nu}\left(\tau,{\rm\bf x}\right)\,.$ (2)

The evolution of both the background and the perturbation are performed using
response functions that scale with $\tau\epsilon^{1/4}/(\eta/s)$, where
$\epsilon$ is the energy density, $\tau$ the evolution time and $\eta/s$ the
shear viscosity to entropy density ratio. In KøMPøST, these response functions
were calculated using pure-glue QCD kinetic theory, but it is now understood
that there is a degree of universality in these functions, both for the
transversely-uniform background (see Ref. Giacalone  _et al._ (2019) and
references therein) and the linear perturbations Kamata  _et al._ (2020).

The procedure for obtaining all the elements of the energy-momentum tensor,
together with its spacetime dependence is detailed in Refs. Kurkela  _et al._
(2019a, b). It is worth noting that the scaling properties of the response
functions in KøMPøST reduce considerably its numerical cost, compared to full
simulations of the Boltzmann equation. However, there are limitations to
KøMPøST, in its current version. First, it is conformal, which inevitably
leads to discontinuities when matching to the non-conformal QCD equation of
state used in hydrodynamic simulations Kurkela  _et al._ (2019b); Nunes da
Silva  _et al._ (2021). Moreover, the lack of bulk viscosity in the conformal
case likely generates larger early radial flow. A second limitation is the
linear response approximation which limits its applicability in systems with
large fluctuations; this includes small systems such as formed in proton-
nucleus collisions.

###  II.2 The fluid-dynamical phase of relativistic heavy-ion collisions

The pre-hydrodynamical “KøMPøST phase” (or directly the IP-Glasma stage) feeds
into a fluid-dynamical evolution. This hydrodynamic phase is modelled using
music Schenke  _et al._ (2010, 2011); Paquet  _et al._ (2016), a relativistic
fluid-dynamical evolution which also takes into account transport coefficients
such as shear Schenke  _et al._ (2011) and bulk viscosity Ryu  _et al._
(2015); Paquet  _et al._ (2016); Ryu  _et al._ (2018). The specific shear
viscosity parameter111The shear viscosity per unit of entropy density. is
taken in this work to be $\eta/s=0.12$, equal to that used in the KøMPøST
phase. To make the connection with earlier calculations of hadron and photon
observables using IP-Glasma and its Yang-Mills evolution to initialize music,
and to quantitatively assess the changes introduced by the pre-hydrodynamic
phase on hadron and photon observables, we shall define and compare the
following two scenarios:

I

IP-Glasma initialization with Yang-Mills propagation for $\tau<\tau_{\rm
IPG}^{\rm EKT}$ 222EKT stands here for “Effective Kinetic Theory”., followed
by a phase where the energy-momentum tensor is evolved with non-equilibrium
linear response Kurkela  _et al._ (2019a, b). That KøMPøST phase is then
followed by a viscous fluid-dynamical evolution which begins at
$\tau=\tau_{\rm EKT}^{\rm hydro}$.

II

IP-Glasma initialization with Yang-Mills propagation for $\tau_{\rm IPG}^{\rm
hydro}$, followed directly by a viscous fluid-dynamical evolution.

In this work we adopt $\tau_{\rm IPG}^{\rm EKT}=0.1$ fm/$c$, which corresponds
to a gluon saturation scale $Q_{s}\sim 2$ GeV. For the transition from KøMPøST
to hydrodynamics, we use $\tau_{\rm EKT}^{\rm hydro}=0.8$ fm/$c$ as default,
in line with the values used in Ref. Kurkela  _et al._ (2019a). For scenario
II, $\tau_{\rm IPG}^{\rm hydro}=0.4$ fm/$c$ is used, consistent with previous
works Ryu  _et al._ (2015); Gale  _et al._ (2013b, a). To help comparisons
between scenarios I and II, a smaller number of additional calculations are
also performed in scenario I with $\tau_{\rm EKT}^{\rm hydro}=0.4$ fm/$c$.

At these switching times the complete energy-momentum tensor $T^{\mu\nu}$ is
passed on to the next stage. Schematically:

$\displaystyle T_{\rm IPG}^{\mu\nu}(\tau_{\rm IPG}^{\rm EKT})$
$\displaystyle=$ $\displaystyle T_{\rm EKT}^{\mu\nu}(\tau_{\rm IPG}^{\rm
EKT})$ $\displaystyle T_{\rm EKT}^{\mu\nu}(\tau_{\rm EKT}^{\rm hydro})$
$\displaystyle=$ $\displaystyle T_{\rm hydro}^{\mu\nu}(\tau_{\rm EKT}^{\rm
hydro})$ (3)

for scenario I, and

$\displaystyle T_{\rm IPG}^{\mu\nu}({\tau_{\rm IPG}^{\rm hydro}})$
$\displaystyle=$ $\displaystyle T_{\rm hydro}^{\mu\nu}(\tau_{\rm IPG}^{\rm
hydro})$ (4)

for scenario II. In every scenario, longitudinal boost invariance is assumed.
The fluid-dynamical phase is followed by a “hadronic afterburner”, which
ensures a dynamic kinetic freezeout, after which the hadrons free-stream to
the detectors. The Cooper-Frye prescription Cooper and Frye (1974) provides
the interface between hydrodynamics and UrQMD Shen  _et al._ (2016a). Viscous
corrections follow Ref. Schenke  _et al._ (2020): the 14-moment viscous
corrections for the shear $\delta f$ and the Chapman-Enskog form for the bulk
viscous corrections Ryu  _et al._ (2015); Schenke  _et al._ (2020). In this
work, the switching temperature to the afterburner is set to $T_{\rm
switch}=145$ MeV.

##  III Results

###  III.1  Hadrons

Figure 1: (Color online) (a) The bulk viscosity temperature profile consistent
with heavy ion data at RHIC and at the LHC, with and without the KøMPøST
phase. Panel (b) shows the corresponding bulk relaxation times and their
minimum value required by causality in the linear regime around equilibrium
Huang  _et al._ (2011).

In this subsection, the production of hadronic observables within the hybrid
model is reviewed, with focus on explicit studies of the effects of the
KøMPøST phase (by comparing scenario I and II discussed in the previous
section). The results are compared against experimental data gathered at RHIC
and at the LHC. Two parameters are adjusted to data separately for scenario I
and II: the bulk viscosity and the normalization of the energy-momentum
tensor.

It is known that the bulk viscosity can have a significant effect on hadronic
observables Ryu  _et al._ (2015). The bulk viscosity profile consistent with
the hadronic data in scenario I is different from that in scenario II, as
shown in Fig. 1. This highlights the first important conclusion of our paper:
the presence of a pre-equilibrium phase like KøMPøST can have a clear
influence on the values of the transport coefficients of QCD extracted by
interpreting the experimental data. Given the fact that KøMPøST is currently
formulated in a conformal formalism, the quantitative effects may change in a
more complete treatment Nunes da Silva  _et al._ (2021). Note that in Fig. 1,
the shaded region in yellow is where the transport coefficients are handled
dynamically by UrQMD.

The normalization of the initial energy-momentum tensor changes between
scenarios I and II because entropy production is different in our kinetic
theory and hydrodynamic phases. The KøMPøST phase in scenario I generates
$\sim 20\%$ less entropy during the preequilibrium phase compared to the
scenario II, which is corrected by a corresponding change in the normalization
factor.

Figure 2: (Color online) For Au+Au collisions at RHIC, the different panels
show (a) multiplicity values, (b) average transverse momentum, and (c)
momentum anisotropy coefficients $v_{n}\\{2\\}$ as a function of collision
centrality, for charged hadrons and identified species. See main text for
details. Results obtained with (scenario I) and without (scenario II) KøMPøST
are shown. The data are from Refs. Adler  _et al._ (2005, 2004); Adamczyk  _et
al._ (2016, 2018).

Starting with Au+Au collisions at the top RHIC energy of 200 A GeV, Fig. 2 (a)
shows multiplicity as a function of centrality class, for charged hadrons as
well as for individual identified hadrons. The multiplicities have been
integrated over a range of pseudorapidity (or rapidity, for identified
hadrons) of one unit, centered at the origin. Once the energy-momentum tensor
normalization and the bulk viscosity are tuned to data, the effect of the
preequilibrium phase on these is seen to be modest for the shape of centrality
dependence but becomes larger in peripheral collisions starting at $\sim
75\%$. The mean transverse momentum of identified hadrons generated in the
relativistic nuclear collisions, $\langle p_{T}\rangle$, is shown in Fig. 2
(b). The transverse momentum in scenario II tracks the experimentally measured
values, up to the last centrality bin, whereas that in scenario I deviates
from observed values starting at $\sim 70\%$ for mesons and somewhat earlier
for protons. Additional verifications suggest that above is predominantly a
consequence of the different hydrodynamics initialization time of KøMPøST and
IP-Glasma. Peripheral collisions generate fluid dynamical phases with a
shorter lifetime than those associated with more central collisions.
Consequently, the relative time spent in the “pre-hydrodynamics” era will
appear correspondingly longer. Since both IP-Glasma and KøMPøST are conformal
Schenke  _et al._ (2012a); Kurkela  _et al._ (2019a), this leaves a shorter
time for the bulk viscosity to generate viscous entropy and reduce the
transverse momentum, as illustrated, for example, in Refs. Ryu  _et al._
(2015); Vujanovic  _et al._ (2020); Nunes da Silva  _et al._ (2021).

This lack of viscous damping and increase of mean $p_{T}$ can also be seen in
the flow coefficients $v_{n}\\{2\\}$ of charged hadrons, obtained from two
particle correlations and shown in Fig. 2 (c), where the effect of going from
scenario I to II grows again with increasing centralities. In summary, for the
specific choice of hydrodynamic initialization times used in this work, the
difference seen between scenario I and II appears to be predominantly from the
later hydrodynamic initialization time used with KøMPøST (scenario I). This
observation is also made with photons in Section III.2.

Concentrating now on conditions prevalent at the Large Hadron Collider, the
interpretation of experimental data with our hybrid model is highlighted in
Fig. 3.

Figure 3: (Color online) For Pb+Pb collisions at the LHC, for an energy of
2760 A GeV. the different panels show (a) multiplicity values, (b) average
transverse momentum, and (c) momentum anisotropy coefficients $v_{n}\\{2\\}$
as a function of collision centrality, for charged hadrons and identified
species. See main text for details. Results obtained with (scenario I) and
without (scenario II) KøMPøST are shown. The data are from Refs. Abelev  _et
al._ (2013); Aamodt  _et al._ (2011).

The same data selection as that for RHIC is shown. Also at the LHC, the
effects of the pre-equilibrium phase are seen mostly in peripheral collisions.
They are quantitatively comparable to those observed at RHIC energies. Turning
to a slightly smaller system at a higher energy – 129Xe + 129Xe at 5440 A GeV
– yields the results shown in Fig. 4, where we also compare directly to Pb+Pb
collisions at 5020 A GeV. We do not observe a dramatic difference in the
effects of the pre-equilibrium stage between different systems. Differences
appear to set in somewhat earlier for the smaller system, as is qualitatively
expected, for the same reason that the effects are larger in peripheral
events.

Figure 4: (Color online) Comparison of results obtained for Pb+Pb collisions
at an energy of 5020 A GeV, and of 129Xe+129Xe at an energy of 5440 A GeV. The
different panels show (a) multiplicity values, (b) average transverse
momentum, and (c)+(d) momentum anisotropy coefficients $v_{n}\\{2\\}$ as
functions of collision centrality, for charged hadrons and identified
particles. See main text for details. Results obtained with (scenario I) and
without (scenario II) KøMPøST are shown. The data are from Refs. Adam  _et
al._ (2016a); Acharya  _et al._ (2019a, b); Adam  _et al._ (2016b); Acharya
_et al._ (2018).

###  III.2  Photons

Electromagnetic radiation constitutes another class of observables capable of
revealing the nature and behavior of the strongly interacting medium created
under the extreme conditions generated by heavy-ion collisions at RHIC and at
the LHC. A striking difference between electromagnetic and hadronic
observables is that the former will survive unscathed the entire spacetime
history of the collision, from the very first interaction to the stage where
final particles are recorded in the detectors. The electromagnetic signal is
penetrating, owing to $\alpha/\alpha_{\rm s}\ll 1$. A complete theoretical
understanding of the photon (real and virtual) signal will require summing the
different components of the photon spectrum, highlighting the importance of
treating electromagnetic signals with a dynamical approach which is realistic
over the entire space-time history of the collision. In this work, the
approach outlined in Ref. Paquet  _et al._ (2016) is followed, with the
exception of an additional stage to be specified shortly. The very first
nucleon-nucleon collisions will generate “prompt photons”, calculated with
next-to-leading-order perturbative QCD using INCNLO Aurenche  _et al._ (1999),
nCTEQ15-np parton distribution functions corrected for nuclear matter effects
Kovarik  _et al._ (2016), and BFG-II fragmentation functions Bourhis  _et al._
(1998). The “thermal photons” are obtained by integrating photon emission
rates Paquet  _et al._ (2016) over the entire spacetime volume occupied by the
fluid-dynamical phase. Corrections to the photon emission rate from the effect
of shear and bulk viscosities are available for many of the rates Paquet  _et
al._ (2016), and are included when available. Photon production in the
hadronic cascade — UrQMD — has been studied in the past Baeuchle and Bleicher
(2010a, b) (see also Ref. Huovinen  _et al._ (2002)), although this feature is
not included in its standard release. As in several previous works, the
photons from the late stages are computed here by using profiles from
hydrodynamics down to T = 105 MeV.333Note that profiles from hydrodynamics in
the low temperature regime 105 MeV $\leq$ T $\leq$ 145 MeV are only used for
the calculations of photons. Hadrons are treated dynamically by UrQMD. This
procedure has been found to be justified quantitatively a posteriori by recent
microscopic calculations of photon production in the hadronic transport code
SMASH Schäfer  _et al._ (2019); *Schafer:2020vvw; *Schafer:2021slz.

We also need to describe photon production from the newly introduced KøMPøST
phase. As mentioned previously, the response functions in the current
implementation of KøMPøST were evaluated with pure-glue kinetic theory. This
could appear incomplete, since photon production requires electrically charged
quarks to be present. In practice, however, there is evidence that the
response functions used in KøMPøST are similar across different microscopic
theories Giacalone  _et al._ (2019); Kamata  _et al._ (2020). In particular,
the introduction of quarks does not appear to significantly change the
response functions (Ref. Giacalone  _et al._ (2019) and references therein),
and thus is not expected to have a major effect on the _evolution_ of the
energy-momentum tensor in KøMPøST.

In principle, computing photon production in the pre-hydrodynamics phase
requires information beyond the energy-momentum tensor provided by KøMPøST:
one needs to know the photon emission rate in non-equilibrium nuclear matter.
In this work, we estimate photon emission in the KøMPøST stage in the same way
as in the hydrodynamic stage: (i) the energy-momentum tensor is decomposed
into energy density, flow velocity, bulk pressure and shear tensor; (ii) the
temperature is calculated from the energy density with the QCD equation of
state Bazavov  _et al._ (2014); and (iii) photon emission rates are folded
with the obtained temperature, flow velocity and viscous component profiles.
While this approach has limitations, it ensures that photon emission is smooth
at the transition to hydrodynamics.

Figure 5: (Color online) The suppression factor applied to the thermal photon
emission rates as a function of proper time, used to simulate the out-of-
chemical equilibrium conditions. The two curves model the effect of a chemical
equilibration time of 1 and 1.5 fm/$c$, respectively.

The establishment of chemical equilibrium over time, i.e., how fast the quarks
are produced and equilibrate in an initially purely gluonic system, can
strongly affect the production of electromagnetic probes Traxler and Thoma
(1996); Chaudhuri (2000); Gelis  _et al._ (2004); Monnai (2014); Bhattacharya
_et al._ (2016); Monnai (2015); Linnyk  _et al._ (2016); Greif  _et al._
(2017); Vovchenko  _et al._ (2016); Oliva  _et al._ (2017); Hauksson  _et al._
(2018); Monnai (2020); Churchill  _et al._ (2021); Garcia-Montero (2020). QCD
kinetic theory simulations Kurkela and Mazeliauskas (2019a, b) provide
detailed information on the approach to chemical equilibrium. For the large
coupling (low specific shear viscosity) such as those used in this work, Ref.
Kurkela and Mazeliauskas (2019a, b) found that chemical equilibration is
reached on a timescale slightly larger than the hydrodynamization time. Recall
that our hydrodynamization time is the time we switch to hydrodynamics, namely
$\tau_{\rm EKT}^{\rm hydro}$ or $\tau_{\rm IPG}^{\rm hydro}$.

We account for different possible chemical equilibration scenarios by
multiplying the photon emission rate with a suppression factor at early times,
based on the fermion energy density fraction of the equilibrium density from
Ref. Kurkela and Mazeliauskas (2019a). We investigate three scenarios of
photon production at early time: instant chemical equilibration, chemical
equilibrium444Following Ref. Kurkela and Mazeliauskas (2019a), chemical
equilibrium is defined as the time at which the suppression factor reaches
0.9. at $\tau_{\rm chem}=1$ fm/$c$, or at $\tau_{\rm chem}=1.5$ fm/$c$. The
corresponding time-dependent suppression factors are shown in Fig. 5.

Figure 6: (Color online) (a) The yield of direct photons in Au+Au collisions
at maximum RHIC energy, in the $0-20\%$ centrality class. The different
channels are enumerated in the text. Here, $\tau_{\rm chem}=1$ fm/$c$. We
compare to experimental data from the PHENIX Adare  _et al._ (2015) and STAR
Adamczyk  _et al._ (2017) Collaborations. (b) The ratio of experimental data
over the total calculated photon yield.

The direct photon yield and the individual contributions from different
sources are shown in Figure 6. We show prompt photons (dashed curve), thermal
photons (QGP + hadron phase, dashed-dotted curve), and photons from the pre-
equilibrium stage (solid yellow curve). The total signal is the dashed red
curve. The calculated photon yield is compared with measurements from the
PHENIX Adare  _et al._ (2015) and STAR Adamczyk  _et al._ (2017)
collaborations. We note the tension between the two experimental measurements
for low momenta $(p_{T}\lesssim 2$ GeV), for the same system at the same
centrality and energy.

Prompt photons that are produced earliest in the collision, dominate at high
$p_{T}\gtrsim 3\,{\rm GeV}$; thermal photons are the largest contribution at
$p_{T}\lesssim 3\,{\rm GeV}$. The pre-equilibrium contribution never
dominates, but exceeds the thermal contribution for $p_{T}\gtrsim 3\,{\rm
GeV}$. The experimental data from the PHENIX collaboration is significantly
underestimated at $p_{T}\lesssim 2$ GeV. A comparison to data from the STAR
collaboration shows better agreement.

Figure 7: (Color online) The photon yield for direct photons in Au+Au
collisions at maximum RHIC energy, in different centrality classes, calculated
for different values of the chemical equilibration time. The PHENIX data are
from Ref. Adare  _et al._ (2015) and the STAR data are from Ref. Adamczyk  _et
al._ (2017). Figure 8: (Color online) Same caption as for Figure 7, but for
the direct photon $v_{2}^{\gamma}\\{SP\\}$. Data are from Ref. Adare  _et al._
(2016). Figure 9: (Color online) Same caption as for Figure 7, but for the
direct photon $v_{3}^{\gamma}\\{SP\\}$. Data are from Ref. Adare  _et al._
(2016).

The effect of the chemical non-equilibrium on the photon spectrum is shown in
Fig. 7 for Au+Au collisions at RHIC energy. Panels (a),(c) and (e) correspond
to the $0-20\%$, $20-40\%$ and $40-60\%$ centrality class respectively. These
panels show the effect of suppressing the photon yield using the time-
dependent factor shown in Fig. 5 to reflect the dynamical buildup of the
fermion content in the model. The ratios to the instant chemical equilibration
case are shown in panels (b), (d), and (f) for clarity. The difference between
a chemical equilibration time of 1 fm/$c$ or 1.5 fm/$c$ is small, and appears
not to be resolved by the current RHIC data on photon spectra. The (chemical)
non-equilibrium effects on the yield appear in the window $1\,{\rm
GeV}\lesssim p_{T}\lesssim 7$ GeV for the physical conditions reported here.
This window of intermediate $p_{T}$ is a result of later stage emission from
the hadronic phase dominating at lower $p_{T}$, and of prompt photons
dominating at higher $p_{T}$, neither of them being affected by the chemical
equilibration.

Compared with experimental data from the PHENIX collaboration Adare  _et al._
(2015), the difference between calculation and measurement grows moving to
more peripheral collisions, where it exceeds one standard deviation for
$p_{T}\lesssim 2$ GeV in the $40-60\%$ centrality class. However, note that a
comparison between the theoretical results reported here and the data reported
by the STAR collaboration Adamczyk  _et al._ (2017), for the same system at
the same energy, show better agreement, as also shown previously in Fig. 6.
Theory and STAR measurements agree within uncertainties, for the two
centralities addressed by the experimental collaboration. The origin of the
discrepancy between the photon results reported by STAR and PHENIX at RHIC
remains unresolved to this day David (2020).

The corresponding calculations for the photon elliptic and triangular momentum
anisotropies, $v_{2}^{\gamma}$ and $v_{3}^{\gamma}$, are shown in Figures 8
and 9, respectively, and compared to experimental data from the PHENIX
collaboration.555Note that the STAR collaboration has not released
measurements of photon $v_{n}$. Importantly, photon momentum anisotropy
coefficients are not obtained through photon-photon correlations but rather
using photon-hadron correlations, because of statistics considerations. The
analysis techniques differ between PHENIX and ALICE (to be shown shortly): the
former uses the event-plane method, while the latter uses the scalar product
method. As discussed in detail in Ref. Paquet  _et al._ (2016), we calculate
the photon $v_{2}^{\gamma}$ with the scalar product method to compare with
both PHENIX and ALICE:

$\displaystyle
v_{n}^{\gamma}\left\\{SP\right\\}\left(p_{T}^{\gamma}\right)=\frac{\langle
v_{n}^{\gamma}\left(p_{T}^{\gamma}\right)v_{n}^{h}\cos\left(n\left(\Psi_{n}^{\gamma}(p_{T}^{\gamma})-\Psi_{n}^{h}\right)\right)\rangle}{\sqrt{\langle\left(v_{n}^{h}\right)^{2}\rangle}}$
(5)

where $h$ refers to hadrons. The measured event-plane flow anisotropy is close
to the scalar product results in the low resolution limit Luzum and Ollitrault
(2013).

The effect of delayed chemical equilibration on the photon spectra naturally
leads to a suppression of the number of early-stage photons as shown in Figure
7. For the $v_{2}^{\gamma}$ and $v_{3}^{\gamma}$ of photons, Figures 8 and 9
show that a delayed equilibration time increases both $v_{2/3}^{\gamma}$ in
the $1-4$ GeV $p_{T}$ range. This is a straightforward consequence of
suppressing early-stage photons that have a small $v_{2/3}^{\gamma}$ and would
dilute the larger $v_{2/3}^{\gamma}$ of photons emitted at later times. As
such, larger values of chemical equilibration time do improve slightly the
agreement with experimental measurements.

Figure 10: (Color online) The photon yield for direct photons in Au+Au
collisions at maximum RHIC energy, in three different centrality classes.
Here, $\tau_{\rm chem}=1$ fm/$c$. See main text for details.  Figure 11:
(Color online) Same caption as for Fig. 10, but for the photon elliptic flow
$v_{2}^{\gamma}(p_{T})$.

The effect of the KøMPøST phase in the photon spectrum is highlighted in Fig.
10. Starting from the top of the legend, the first curve shows photons
generated in our scenario I, with $\tau_{\rm IPG}^{\rm EKT}=0.1$ fm/$c$ and
$\tau_{\rm EKT}^{\rm hydro}=0.8$ fm/$c$. The second curve shows the direct
photon signal obtained in scenario II with $\tau_{\rm IPG}^{\rm hydro}=0.4$
fm/$c$. Finally, in order to obtain a more meaningful assessment of the role
of the KøMPøST phase with respect to scenario II, a calculation run with the
pre-equilibrium dynamics lasting only up to 0.4 fm/$c$ is done ($\tau_{\rm
IPG}^{\rm EKT}=0.1$ fm/$c$ and $\tau_{\rm EKT}^{\rm hydro}=0.4$ fm/$c$). The
ratios of both this last scenario and scenario II to scenario I are also shown
in Fig. 10. Here it is revealed that differences are maximally 10%, with the
spectra enhanced by a longer KøMPøST phase for $1.5\,{\rm GeV}<p_{T}<6\,{\rm
GeV}$ and reduced for $p_{T}\lesssim 1.5\,{\rm GeV}$. Replacing IP-Glasma with
KøMPøST between 0.1 fm/$c$ and 0.4 fm/$c$ has a small effect.

The corresponding result for the photon $v_{2}^{\gamma}$ is shown in Figure
11, for 0-20%, 20-40% and 40-60% centralities. For the most central case, all
three pre-equilibrium scenarios lead to similar $v_{2}^{\gamma}$. However, a
different result is observed for more peripheral centralities: the
$v_{2}^{\gamma}$ when transitioning from KøMPøST to hydrodynamics at
$\tau_{\rm EKT}^{\rm hydro}=0.8$ fm/$c$ is larger than the two other scenarios
that transition to hydrodynamics at $0.4$ fm/$c$. This larger dependence on
the pre-equilibrium scenarios in peripheral collisions was also observed for
hadrons. In general, peripheral collisions are expected to be more sensitive
to the pre-equilibrium phase, as a consequence of their shorter hydrodynamic
phase compared to central events. Physically, the dependence of
$v_{2}^{\gamma}$ on the transition time to hydrodynamics could be related to
the conformal nature of the pre-equilibrium models.

The overall difference between the three scenarios is small compared to the
current tension with PHENIX measurements. In this sense, current RHIC photon
data does not demand the presence of a phase like KøMPøST. However, and
importantly, the KøMPøST phase has the potential to provide a dovetail
matching Kurkela and Mazeliauskas (2019a) between IP-Glasma and relativistic
viscous fluid dynamics, from a theory point of view. Including a pre-
equilibrium phase can also avoid non-causal behavior, which can appear when
running hydrodynamics using a general off-equilibrium initial state like IP-
Glasma Chiu and Shen (2021); Plumberg  _et al._ (2021). The matching to
hydrodynamics is expected to be further improved by future non-equilibrium
models that go beyond current approximations such as conformality and linear
response.

Figure 12: (Color online) (a) Direct photon yield in Pb+Pb collisions at 2760
A GeV, in the $0-20\%$ centrality class. Also shown is the breakdown into the
different components. (b) The ratio between the experimental data and the
calculated total photon yield. Experimental data from the ALICE Collaboration
Adam  _et al._ (2016c). See main text for details.

Considering photons generated by collisions at the LHC, Fig. 12 (a) reports on
results of calculations performed for Pb+Pb at 2760 A GeV. As in the case of
RHIC, the same hydrodynamic evolution is used for photons as for hadrons. The
kinematical window where the pre-equilibrium photons contribute the most
(relative to prompt and thermal photon yields) is similar to that seen at
RHIC, around $p_{T}=4$ GeV. Experimental data from the ALICE Collaboration
Adam  _et al._ (2016c) is shown for comparison. Figure 12 (b) shows the ratio
of the experimental data to the calculated total direct photon yield.
Agreement is very good for $p_{T}\gtrsim 3\,{\rm GeV}$, with deviations,
albeit mostly within experimental errors, seen at lower $p_{T}$.

The effect of chemical non-equilibrium on the direct photon yields in two
different centrality classes is illustrated in Fig. 13. The chemical
suppression factor has the largest effects in the window $2\lesssim
p_{T}\lesssim 7$ GeV, and the ALICE data show a preference for the model with
chemical non-equilibrium.

Figure 13: (Color online) The photon yield for direct photons in Pb+Pb
collisions at $\sqrt{s_{\mathrm{NN}}}=2760$ GeV in three different centrality
classes. Here, $\tau_{\rm chem}=1$ fm/$c$. See main text for details.  Figure
14: (Color online) Direct photon elliptic flow $v_{2}^{\gamma}(p_{T})$ in
$2760\,{\rm AGeV}$ Pb+Pb collisions for different chemical equilibration times
compared to experimental data from the ALICE collaboration Acharya  _et al._
(2019c). See main text for details.

Figure 14 shows the elliptic flow of direct photons in 2760 A GeV Pb+Pb
collisions, compared to experimental data from the ALICE collaboration Acharya
_et al._ (2019c). For the measured centrality classes, the calculations
systematically undershoot the data, even if the two are mostly consistent
within uncertainties. As was the case at RHIC, a larger value of the chemical
equilibration time reduces the difference between theory and experimental
data.

Figure 15: (Color online) The direct photon yield and $v_{2}^{\gamma}$ in
0-20% centrality, for Au+Au @ 200 A GeV and Pb+Pb @ 2760 A GeV and 5020 A GeV.

Figure 15 shows the collision energy dependence of the direct photon spectra
and their elliptic flow coefficients from the top RHIC to LHC energies in
0-20% central heavy-ion collisions. We find that direct photon yield and
$v_{2}^{\gamma}$ increase with $\sqrt{s_{\mathrm{NN}}}$. Prompt photons have
flatter $p_{T}$-spectra at LHC energies. At the same time, more thermal
photons radiate from Pb+Pb collisions at the LHC than those at RHIC because of
the higher averaged fireball temperature and larger space-time volume. The
increase of direct photon $v_{2}^{\gamma}$ with increasing collision energy
suggests that the enhancement in the thermal radiation wins over the prompt
photon production at high energies.

###  III.3 Photon production in “small systems”

This section is devoted to the theoretical treatment of what has become known
in the field as “small systems”. This includes asymmetric systems involving
protons and a heavier ion, as well as symmetric collisions of two smaller
nuclei, such as O+O. The interest generated by such configurations owes much
to the fact that those systems show a similar degree of collectivity (e.g.
flow coefficients) than that observed in “large” colliding systems such as
some of the ones discussed earlier: Au+Au and Pb+Pb. This fact was met with
surprise by the community, as one of the original motivations to perform
experimental runs with proton projectiles incident on a larger target was to
provide a neutral baseline. In addition to their undeniable theoretical
interest, small systems deserve a separate section here because KøMPøST was
not designed to operate in such environments: the combination of small size
and large gradients pushes the linear treatment past its limit of validity.
Therefore, only the combination of IP-Glasma and hydrodynamics — scenario II —
is used in this section.

Figure 16: (Color online) Averaged transverse momentum and flow anisotropic
coefficients for charged hadrons in O+O and p+Pb collisions compared with
those in Pb+Pb collisions at 5020 A GeV.

In addition to investigating whether QGP is also present in small systems,
another reason to study them is that they provide another mean to explore the
mechanisms of particle production as a function of system size. For instance,
for a colliding energy of 5020 A GeV, the charged hadron multiplicity,
$dN_{\rm ch}/d\eta$,for minimum-bias O+O collisions is comparable to that for
Pb+Pb collisions in a centrality class of $70-80\%$ Huss  _et al._ (2021);
Brewer  _et al._ (2021). Similarly, O+O collisions at the same energy and in a
centrality class of $0-5\%$ can be compared to Pb+Pb collisions in a $50-70\%$
centrality class.

The calculated charged multiplicity dependence of light system (O+O and p+Pb)
observables $\left(\langle p_{T}\rangle,v_{n}\\{2\\},n=2,3\right)$ is shown in
Fig. 16. Hadrons’ averaged transverse momenta in p+Pb and O+O collisions are
larger than those in Pb+Pb collisions at the same multiplicity, which
indicates that a stronger radial flow is developed by larger pressure
gradients in the more compact systems. For $dN^{\mathrm{ch}}/d\eta>20$, the
elliptic flow coefficients have the ordering
$v_{2}(\mathrm{Pb+Pb})>v_{2}(\mathrm{O+O})>v_{2}(\mathrm{p+Pb})$, driven by
the global collision geometry. In contrast, the triangular flow $v_{3}$
coefficients do not show any system dependence at fixed charged particle
multiplicities since $v_{3}$ is mainly driven by the event-by-event initial-
state fluctuations. We note that both $v_{2}$ and $v_{3}$ increase as
$dN^{\mathrm{ch}}/d\eta$ goes below $5$. In these peripheral collisions, the
initial-state momentum anisotropy gradually begins to dominate the measured
anisotropic flow coefficients over the hydrodynamic response to the global
geometry Schenke  _et al._ (2020); Giacalone  _et al._ (2020).

Figure 17: (Color online) (a) Energy-density-weighted average temperature of
fluid dynamical simulations for central O+O and semi-peripheral Pb+Pb
collisions with the same charged hadron multiplicity at 5020 A GeV. (b) Energy
density weighted average transverse flow velocity for the two systems
described in (a).

Figure 17 shows the evolution of energy-density-weighted average temperature
(a) and flow velocity (b) in 0-5% O+O and 50-70% Pb+Pb collisions. As already
mentioned, these two collision systems produce similar charged hadron
multiplicity at midrapidity, $dN^{\mathrm{ch}}/d\eta=116$. Owing to a
relatively compact size, central O+O collisions achieve a higher average
temperature and stronger radial flow during the hydrodynamic evolution than
those in the semi-peripheral Pb+Pb collisions. The time-integrated effect of
such a stronger radial flow is reflected in the hadrons’ mean transverse
momenta in Fig. 16 (a). However, it is difficult to find experimental evidence
of the higher initial temperature in the central O+O collisions using hadronic
observables. In this case, thermal photon radiation provides a better tool to
access the early-stage dynamics.

For the same conditions – central O+O collisions and peripheral Pb+Pb
collisions – the direct photon nuclear modification factor, $R^{\gamma}_{AA}$
is shown in Fig. 18. This figure illustrates well the power of electromagnetic
radiation to report on local conditions at the time of its emission. The
difference between central O+O and peripheral Pb+Pb is now clearly visible
when comparing the photon nuclear modification ratios. For $p_{T}\sim 2$ GeV,
$R_{AA}^{\gamma}$ for O+O collisions in the $0-5\%$ centrality class is
enhanced by approximately 80% over its value in peripheral Pb+Pb events, for
$\sqrt{s_{\rm NN}}$ = 5020 A GeV.

Figure 18: (Color online) $R_{AA}^{\gamma}$ of photon in 0-5% O+O collisions
at 5020 A GeV, compared to 50-70% Pb+Pb collisions. Figure 19: (Color online)
The $R^{\gamma}_{AA}$ of direct photons in O+O collisions at 5020 A GeV. The
lower line includes only prompt photons, with isospin effect and nCTEQ15
nuclear pdf’s as described in Section III.2. The middle and top lines include
medium photons, for minimum bias collisions and 0-5% central collisions,
respectively. Figure 20: (Color online) $R^{\gamma}_{pPb}$ of direct photons
in p+Pb collisions at 5020 A GeV. The lower line includes only prompt photons,
with isospin effect and nCTEQ15 nuclear pdf’s as described in Section III.2.
The middle and top lines include medium photons, for minimum bias collisions
and 0-5% central collisions, respectively.

Finally, Figures 19 and 20 highlight the direct photon nuclear modification
factor as a function of $p_{T}$ for O+O and p+Pb collisions, respectively,
showing separately the effect of the prompt photons. The direct photon nuclear
modification factors are shown for 0-5% central collisions and for the case
where one accepts all centralities (minimum bias, 0-100%). The central
collisions show a significantly larger direct photon enhancement compared to
the minimum bias case for both O+O and p+Pb collisions. In O+O collisions
there is a 50% enhancement in minimum bias collisions (compared to the p+p
case), while in p+Pb collisions the enhancement is significantly smaller.

We find a quantitatively similar thermal photon enhancement as those reported
in Refs. Shen  _et al._ (2016b, 2017), despite the use of a different type of
initial state model (Refs. Shen  _et al._ (2016b, 2017) used a Glauber model
initial state) . The enhancement of direct photon $R^{\gamma}_{pPb}$ for
$p_{T}<5$ GeV seen here is an experimental signature of a nearly-thermalized
QGP produced in p+Pb collisions. Without thermal radiation contributions, the
direct photon $R^{\gamma}_{AA}$ is expected to be smaller than 1 because of
the nuclear shadowing effects. Note that a photon nuclear modification factor,
$R_{pA}^{\gamma}$ was measured in p+Au collisions at RHIC. That data has large
uncertainties, but is nevertheless very suggestive of an enhancement like that
described here Khachatryan  _et al._ .

##  IV Conclusion

In this work, we studied simultaneously hadron and photon production in a
state-of-the-art multistage model of heavy ion collisions, in a range of
different systems at RHIC and LHC. We studied the effect of varying the pre-
equilibrium evolution and chemical equilibration time.

The multistage model included IP-Glasma and KøMPøST to initialize the
hydrodynamic simulation music, which was then followed by the hadronic
transport UrQMD. The introduction of KøMPøST changes the entropy production of
the multistage model as well as the transport coefficients extracted from
measurements. Once these are accounted for, the influence of the pre-
equilibrium model KøMPøST on the hadronic observables was quantitatively
modest, with its largest effect being on the average transverse momentum of
charged hadrons at high centrality (peripheral events). In this work both
$\eta/s$ and $\zeta/s$ were chosen to best model the hadronic final state
data, and only $\zeta/s$ was modified when including the KøMPøST phase, while
we chose $\eta/s$ to be fixed; to explore more detailed effects on all
transport coefficients, pre-hydrodynamics can be included in future global
Bayesian analyses of relativistic nuclear collisions Pratt  _et al._ (2015);
Bernhard  _et al._ (2019); Nijs  _et al._ (2021); Everett  _et al._ (2021).

We further studied the sensitivity of photons to initial deviations from
chemical equilibrium. We introduced a suppression parameter estimated from QCD
kinetic theory studies. The effect of chemical equilibrium was found to be
strongest for the spectra around $p_{T}\approx 4\,{\rm GeV}$ and for the
$v_{n}$ for $p_{T}\gtrsim 1-1.5\,{\rm GeV}$ for heavy ion collisions at both
RHIC and LHC energies. The current photon data sets cannot resolve the
specific value of the chemical equilibration time, even though the general
trend observed for photon elliptic flow seems to favor a larger one.

Much has been written about “the photon $v_{2}$ puzzle” [See; forexample; ][;
andreferencestherein.]Gale:2012xq; David (2020). In a nutshell, the puzzle
stems from the fact that the measured direct photon elliptic flow has been
found to be as large as that of hadrons, in the region $p_{T}<4$ GeV/c, as is
also clear from the data shown in this paper. The majority of theoretical
models currently underpredict the photon spectrum and elliptic flow. No
approach with realistic dynamics can both reproduce photon spectrum and
elliptic flow, and this situation has not been modified because of the
inclusion of a pre-hydrodynamic phase like KøMPøST.

For small collision systems, such as O+O and p+Pb, we predict a strong
enhancement of photon production for transverse momenta $p_{T}\lesssim 4\,{\rm
GeV}$, compared to the expectation from prompt photon production.
Consequently, experimental measurements of photon spectra in small systems
will be essential to determine how small a QGP can be created in nuclear
collisions.

We emphasize again that photons (real and virtual) enjoy a special status in
efforts to characterize strongly interacting matter under extreme conditions.
They can be both penetrating and soft (on typical strong interaction energy
scales): hadronic observables that are penetrating are not soft, and
conversely. Measured together, strong and electroweak particles are a prime
example of multi-messenger heavy-ion physics.

###### Acknowledgements.

This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada, in part by the US Department of Energy under
Contracts DE-SC0012704, DE-SC0013460, DE-FG02-05ER41367, DE-AC02-05CH11231,
and in part by the National Science Foundation (NSF) under grant number
PHY-2012922. This research used resources of the National Energy Research
Scientific Computing Center, which is supported by the Office of Science of
the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

## References

  * Yang and Mills (1954) Chen-Ning Yang and Robert L. Mills, “Conservation of Isotopic Spin and Isotopic Gauge Invariance,” Phys. Rev. 96, 191–195 (1954).
  * Sterman  _et al._ (1995) George Sterman, John Smith, John C. Collins, James Whitmore, Raymond Brock, Joey Huston, Jon Pumplin, Wu-Ki Tung, Hendrik Weerts, Chien-Peng Yuan, Stephen Kuhlmann, Sanjib Mishra, Jorge G. Morfín, Fredrick Olness, Joseph Owens, Jianwei Qiu, and Davison E. Soper, “Handbook of perturbative QCD,” Rev. Mod. Phys. 67, 157–248 (1995).
  * Bazavov  _et al._ (2010) A. Bazavov _et al._ (MILC), “Nonperturbative QCD Simulations with 2+1 Flavors of Improved Staggered Quarks,” Rev. Mod. Phys. 82, 1349–1417 (2010), arXiv:0903.3598 [hep-lat] .
  * Ratti (2018) Claudia Ratti, “Lattice QCD and heavy ion collisions: a review of recent progress,” Rept. Prog. Phys. 81, 084301 (2018), arXiv:1804.07810 [hep-lat] .
  * Gale  _et al._ (2013a) Charles Gale, Sangyong Jeon, and Bjoern Schenke, “Hydrodynamic Modeling of Heavy-Ion Collisions,” Int. J. Mod. Phys. A28, 1340011 (2013a), arXiv:1301.5893 [nucl-th] .
  * Baym (2016) Gordon Baym, “Ultrarelativistic heavy ion collisions: the first billion seconds,”  _Proceedings, 25th International Conference on Ultra-Relativistic Nucleus-Nucleus Collisions (Quark Matter 2015): Kobe, Japan, September 27-October 3, 2015_ , Nucl. Phys. A956, 1–10 (2016), arXiv:1701.03972 [nucl-ex] .
  * Heinz and Snellings (2013) Ulrich Heinz and Raimond Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part. Sci. 63, 123–151 (2013), arXiv:1301.2826 [nucl-th] .
  * Kurkela  _et al._ (2019a) Aleksi Kurkela, Aleksas Mazeliauskas, Jean-François Paquet, Sören Schlichting, and Derek Teaney, “Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory,” Phys. Rev. Lett. 122, 122302 (2019a), arXiv:1805.01604 [hep-ph] .
  * Kurkela  _et al._ (2019b) Aleksi Kurkela, Aleksas Mazeliauskas, Jean-François Paquet, Sören Schlichting, and Derek Teaney, “Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions,” Phys. Rev. C99, 034910 (2019b), arXiv:1805.00961 [hep-ph] .
  * Schenke  _et al._ (2012a) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Fluctuating Glasma initial conditions and flow in heavy ion collisions,” Phys. Rev. Lett. 108, 252301 (2012a), arXiv:1202.6646 [nucl-th] .
  * Schenke  _et al._ (2012b) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions,” Phys. Rev. C86, 034908 (2012b), arXiv:1206.6805 [hep-ph] .
  * Giacalone  _et al._ (2019) Giuliano Giacalone, Aleksas Mazeliauskas, and Sören Schlichting, “Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions,” Phys. Rev. Lett. 123, 262301 (2019), arXiv:1908.02866 [hep-ph] .
  * Kamata  _et al._ (2020) Syo Kamata, Mauricio Martinez, Philip Plaschke, Stephan Ochsenfeld, and Sören Schlichting, “Hydrodynamization and nonequilibrium Green’s functions in kinetic theory,” Phys. Rev. D 102, 056003 (2020), arXiv:2004.06751 [hep-ph] .
  * Nunes da Silva  _et al._ (2021) Tiago Nunes da Silva, David Chinellato, Mauricio Hippert, Willian Serenone, Jun Takahashi, Gabriel S. Denicol, Matthew Luzum, and Jorge Noronha, “Pre-hydrodynamic evolution and its signatures in final-state heavy-ion observables,” Phys. Rev. C 103, 054906 (2021), arXiv:2006.02324 [nucl-th] .
  * Schenke  _et al._ (2010) Bjoern Schenke, Sangyong Jeon, and Charles Gale, “(3+1)D hydrodynamic simulation of relativistic heavy-ion collisions,” Phys. Rev. C82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Schenke  _et al._ (2011) Björn Schenke, Sangyong Jeon, and Charles Gale, “Elliptic and triangular flow in event-by-event (3+1)D viscous hydrodynamics,” Phys. Rev. Lett. 106, 042301 (2011), arXiv:1009.3244 [hep-ph] .
  * Paquet  _et al._ (2016) Jean-François Paquet, Chun Shen, Gabriel S. Denicol, Matthew Luzum, Björn Schenke, Sangyong Jeon, and Charles Gale, “Production of photons in relativistic heavy-ion collisions,” Phys. Rev. C93, 044906 (2016), arXiv:1509.06738 [hep-ph] .
  * Ryu  _et al._ (2015) S. Ryu, J. F. Paquet, C. Shen, G. S. Denicol, B. Schenke, S. Jeon, and C. Gale, “Importance of the Bulk Viscosity of QCD in Ultrarelativistic Heavy-Ion Collisions,” Phys. Rev. Lett. 115, 132301 (2015), arXiv:1502.01675 [nucl-th] .
  * Ryu  _et al._ (2018) Sangwook Ryu, Jean-François Paquet, Chun Shen, Gabriel Denicol, Björn Schenke, Sangyong Jeon, and Charles Gale, “Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider,” Phys. Rev. C97, 034910 (2018), arXiv:1704.04216 [nucl-th] .
  * Gale  _et al._ (2013b) Charles Gale, Sangyong Jeon, Björn Schenke, Prithwish Tribedy, and Raju Venugopalan, “Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics,” Phys. Rev. Lett. 110, 012302 (2013b), arXiv:1209.6330 [nucl-th] .
  * Cooper and Frye (1974) Fred Cooper and Graham Frye, “Comment on the Single Particle Distribution in the Hydrodynamic and Statistical Thermodynamic Models of Multiparticle Production,” Phys. Rev. D10, 186 (1974).
  * Shen  _et al._ (2016a) Chun Shen, Zhi Qiu, Huichao Song, Jonah Bernhard, Steffen Bass, and Ulrich Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun. 199, 61–85 (2016a), arXiv:1409.8164 [nucl-th] .
  * Schenke  _et al._ (2020) Bjoern Schenke, Chun Shen, and Prithwish Tribedy, “Running the gamut of high energy nuclear collisions,” Phys. Rev. C 102, 044905 (2020), arXiv:2005.14682 [nucl-th] .
  * Huang  _et al._ (2011) Xu-Guang Huang, Takeshi Kodama, Tomoi Koide, and Dirk H. Rischke, “Bulk Viscosity and Relaxation Time of Causal Dissipative Relativistic Fluid Dynamics,” Phys. Rev. C 83, 024906 (2011), arXiv:1010.4359 [nucl-th] .
  * Adler  _et al._ (2005) S. S. Adler _et al._ (PHENIX), “Systematic studies of the centrality and s(NN)**(1/2) dependence of the d E(T) / d eta and d (N(ch) / d eta in heavy ion collisions at mid-rapidity,” Phys. Rev. C 71, 034908 (2005), [Erratum: Phys.Rev.C 71, 049901 (2005)], arXiv:nucl-ex/0409015 .
  * Adler  _et al._ (2004) S. S. Adler _et al._ (PHENIX), “Identified charged particle spectra and yields in Au+Au collisions at S(NN)**1/2 = 200-GeV,” Phys. Rev. C 69, 034909 (2004), arXiv:nucl-ex/0307022 .
  * Adamczyk  _et al._ (2016) L. Adamczyk _et al._ (STAR), “Beam Energy Dependence of the Third Harmonic of Azimuthal Correlations in Au+Au Collisions at RHIC,” Phys. Rev. Lett. 116, 112302 (2016), arXiv:1601.01999 [nucl-ex] .
  * Adamczyk  _et al._ (2018) L. Adamczyk _et al._ (STAR), “Harmonic decomposition of three-particle azimuthal correlations at energies available at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 98, 034918 (2018), arXiv:1701.06496 [nucl-ex] .
  * Vujanovic  _et al._ (2020) Gojko Vujanovic, Jean-François Paquet, Chun Shen, Gabriel S. Denicol, Sangyong Jeon, Charles Gale, and Ulrich Heinz, “Exploring the influence of bulk viscosity of QCD on dilepton tomography,” Phys. Rev. C 101, 044904 (2020), arXiv:1903.05078 [nucl-th] .
  * Abelev  _et al._ (2013) Betty Abelev _et al._ (ALICE), “Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV,” Phys. Rev. C 88, 044910 (2013), arXiv:1303.0737 [hep-ex] .
  * Aamodt  _et al._ (2011) K. Aamodt _et al._ (ALICE), “Higher harmonic anisotropic flow measurements of charged particles in Pb-Pb collisions at $\sqrt{s_{NN}}$=2.76 TeV,” Phys. Rev. Lett. 107, 032301 (2011), arXiv:1105.3865 [nucl-ex] .
  * Adam  _et al._ (2016a) Jaroslav Adam _et al._ (ALICE), “Centrality dependence of the charged-particle multiplicity density at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV,” Phys. Rev. Lett. 116, 222302 (2016a), arXiv:1512.06104 [nucl-ex] .
  * Acharya  _et al._ (2019a) Shreyasi Acharya _et al._ (ALICE), “Centrality and pseudorapidity dependence of the charged-particle multiplicity density in Xe–Xe collisions at $\sqrt{s_{\rm NN}}$ =5.44TeV,” Phys. Lett. B 790, 35–48 (2019a), arXiv:1805.04432 [nucl-ex] .
  * Acharya  _et al._ (2019b) Shreyasi Acharya _et al._ (ALICE), “Transverse momentum spectra and nuclear modification factors of charged particles in Xe-Xe collisions at $\sqrt{s_{\rm NN}}$ = 5.44 TeV,” Phys. Lett. B 788, 166–179 (2019b), arXiv:1805.04399 [nucl-ex] .
  * Adam  _et al._ (2016b) Jaroslav Adam _et al._ (ALICE), “Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,” Phys. Rev. Lett. 116, 132302 (2016b), arXiv:1602.01119 [nucl-ex] .
  * Acharya  _et al._ (2018) S. Acharya _et al._ (ALICE), “Anisotropic flow in Xe-Xe collisions at $\mathbf{\sqrt{s_{\rm{NN}}}=5.44}$ TeV,” Phys. Lett. B 784, 82–95 (2018), arXiv:1805.01832 [nucl-ex] .
  * Aurenche  _et al._ (1999) P. Aurenche, M. Fontannaz, J. Ph. Guillet, Bernd A. Kniehl, E. Pilon, and M. Werlen, “A Critical phenomenological study of inclusive photon production in hadronic collisions,” Eur. Phys. J. C 9, 107–119 (1999), arXiv:hep-ph/9811382 .
  * Kovarik  _et al._ (2016) K. Kovarik _et al._ , “nCTEQ15 - Global analysis of nuclear parton distributions with uncertainties in the CTEQ framework,” Phys. Rev. D 93, 085037 (2016), arXiv:1509.00792 [hep-ph] .
  * Bourhis  _et al._ (1998) L. Bourhis, M. Fontannaz, and J. P. Guillet, “Quarks and gluon fragmentation functions into photons,” Eur. Phys. J. C 2, 529–537 (1998), arXiv:hep-ph/9704447 .
  * Baeuchle and Bleicher (2010a) Bjoern Baeuchle and Marcus Bleicher, “Hybrid model calculations of direct photons in high-energy nuclear collisions,” Phys. Rev. C 81, 044904 (2010a), arXiv:0905.4678 [hep-ph] .
  * Baeuchle and Bleicher (2010b) Bjoern Baeuchle and Marcus Bleicher, “Direct photon calculations in heavy-ion collisions at $\sqrt{s_{NN}}$ = 62.4 - 200 AGeV in a (3+1) dimensional hybrid approach,” Phys. Rev. C 82, 064901 (2010b), arXiv:1008.2332 [nucl-th] .
  * Huovinen  _et al._ (2002) P. Huovinen, M. Belkacem, P. J. Ellis, and Joseph I. Kapusta, “Dileptons and photons from coarse grained microscopic dynamics and hydrodynamics compared to experimental data,” Phys. Rev. C 66, 014903 (2002), arXiv:nucl-th/0203023 .
  * Schäfer  _et al._ (2019) Anna Schäfer, Juan M. Torres-Rincon, Jonas Rothermel, Niklas Ehlert, Charles Gale, and Hannah Elfner, “Benchmarking a nonequilibrium approach to photon emission in relativistic heavy-ion collisions,” Phys. Rev. D 99, 114021 (2019), arXiv:1902.07564 [nucl-th] .
  * Schäfer  _et al._ (2021a) Anna Schäfer, Juan M. Torres-Rincon, Charles Gale, and Hannah Elfner, “A Non-Equilibrium Approach to Photon Emission from the Late Stages of Relativistic Heavy-Ion Collisions,” Nucl. Phys. A 1005, 121772 (2021a), arXiv:2001.03378 [nucl-th] .
  * Schäfer  _et al._ (2021b) Anna Schäfer, Oscar Garcia-Montero, Jean-François Paquet, Hannah Elfner, and Charles Gale, “Out-of-Equilibrium Photon Production in the Late Stages of Relativistic Heavy-Ion Collisions,” (2021b), arXiv:2111.13603 [hep-ph] .
  * Bazavov  _et al._ (2014) A. Bazavov _et al._ (HotQCD), “Equation of state in ( 2+1 )-flavor QCD,” Phys. Rev. D90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Traxler and Thoma (1996) Christoph T. Traxler and Markus H. Thoma, “Photon emission from a parton gas at chemical nonequilibrium,” Phys. Rev. C53, 1348–1352 (1996), arXiv:hep-ph/9507444 [hep-ph] .
  * Chaudhuri (2000) A K Chaudhuri, “Chemical equilibration in a viscous quark-gluon plasma and electromagnetic signals,” Journal of Physics G: Nuclear and Particle Physics 26, 1433–1445 (2000).
  * Gelis  _et al._ (2004) F. Gelis, H. Niemi, P. V. Ruuskanen, and S. S. Rasanen, “Photon production from nonequilibrium QGP in heavy ion collisions,”  _Ultra-relativistic nucleus-nucleus collisions. Proceedings, 17th International Conference, Quark Matter 2004, Oakland, USA, January 11-17, 2004_ , J. Phys. G30, S1031–S1036 (2004), arXiv:nucl-th/0403040 [nucl-th] .
  * Monnai (2014) Akihiko Monnai, “Thermal photon $v_{2}$ with slow quark chemical equilibration,” Phys. Rev. C 90, 021901 (2014), arXiv:1403.4225 [nucl-th] .
  * Bhattacharya  _et al._ (2016) Lusaka Bhattacharya, Radoslaw Ryblewski, and Michael Strickland, “Photon production from a nonequilibrium quark-gluon plasma,” Phys. Rev. D 93, 065005 (2016), arXiv:1507.06605 [hep-ph] .
  * Monnai (2015) Akihiko Monnai, “Chemically non-equilibrated QGP and thermal photon elliptic flow,” in  _7th International Conference on Hard and Electromagnetic Probes of High-Energy Nuclear Collisions_ (2015) arXiv:1510.00539 [nucl-th] .
  * Linnyk  _et al._ (2016) O. Linnyk, E. L. Bratkovskaya, and W. Cassing, “Effective QCD and transport description of dilepton and photon production in heavy-ion collisions and elementary processes,” Prog. Part. Nucl. Phys. 87, 50–115 (2016), arXiv:1512.08126 [nucl-th] .
  * Greif  _et al._ (2017) Moritz Greif, Florian Senzel, Heiner Kremer, Kai Zhou, Carsten Greiner, and Zhe Xu, “Nonequilibrium photon production in partonic transport simulations,” Phys. Rev. C95, 054903 (2017), arXiv:1612.05811 [hep-ph] .
  * Vovchenko  _et al._ (2016) V. Vovchenko, Iu. A. Karpenko, M. I. Gorenstein, L. M. Satarov, I. N. Mishustin, B. Kämpfer, and H. Stoecker, “Electromagnetic probes of a pure-glue initial state in nucleus-nucleus collisions at energies available at the CERN Large Hadron Collider,” Phys. Rev. C 94, 024906 (2016), arXiv:1604.06346 [nucl-th] .
  * Oliva  _et al._ (2017) L. Oliva, M. Ruggieri, S. Plumari, F. Scardina, G. X. Peng, and V. Greco, “Photons from the Early Stages of Relativistic Heavy Ion Collisions,” Phys. Rev. C 96, 014914 (2017), arXiv:1703.00116 [nucl-th] .
  * Hauksson  _et al._ (2018) Sigtryggur Hauksson, Sangyong Jeon, and Charles Gale, “Photon emission from quark-gluon plasma out of equilibrium,” Phys. Rev. C97, 014901 (2018), arXiv:1709.03598 [nucl-th] .
  * Monnai (2020) Akihiko Monnai, “Prompt, pre-equilibrium, and thermal photons in relativistic nuclear collisions,” J. Phys. G 47, 075105 (2020), arXiv:1907.09266 [nucl-th] .
  * Churchill  _et al._ (2021) Jessica Churchill, Li Yan, Sangyong Jeon, and Charles Gale, “Emission of electromagnetic radiation from the early stages of relativistic heavy-ion collisions,” Phys. Rev. C 103, 024904 (2021), arXiv:2008.02902 [hep-ph] .
  * Garcia-Montero (2020) Oscar Garcia-Montero, “Non-thermal photons from the pre-equilibrium stages,” Journal of Physics: Conference Series 1667, 012008 (2020).
  * Kurkela and Mazeliauskas (2019a) Aleksi Kurkela and Aleksas Mazeliauskas, “Chemical Equilibration in Hadronic Collisions,” Phys. Rev. Lett. 122, 142301 (2019a), arXiv:1811.03040 [hep-ph] .
  * Kurkela and Mazeliauskas (2019b) Aleksi Kurkela and Aleksas Mazeliauskas, “Chemical equilibration in weakly coupled QCD,” Phys. Rev. D 99, 054018 (2019b), arXiv:1811.03068 [hep-ph] .
  * Adare  _et al._ (2015) A. Adare _et al._ (PHENIX), “Centrality dependence of low-momentum direct-photon production in Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV,” Phys. Rev. C 91, 064904 (2015), arXiv:1405.3940 [nucl-ex] .
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), “Direct virtual photon production in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV,” Phys. Lett. B 770, 451–458 (2017), arXiv:1607.01447 [nucl-ex] .
  * Adare  _et al._ (2016) A. Adare _et al._ (PHENIX), “Azimuthally anisotropic emission of low-momentum direct photons in Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV,” Phys. Rev. C94, 064901 (2016), arXiv:1509.07758 [nucl-ex] .
  * David (2020) Gabor David, “Direct real photons in relativistic heavy ion collisions,” Rept. Prog. Phys. 83, 046301 (2020), arXiv:1907.08893 [nucl-ex] .
  * Luzum and Ollitrault (2013) Matthew Luzum and Jean-Yves Ollitrault, “Eliminating experimental bias in anisotropic-flow measurements of high-energy nuclear collisions,” Phys. Rev. C 87, 044907 (2013), arXiv:1209.2323 [nucl-ex] .
  * Chiu and Shen (2021) Cheng Chiu and Chun Shen, “Exploring theoretical uncertainties in the hydrodynamic description of relativistic heavy-ion collisions,” Phys. Rev. C 103, 064901 (2021), arXiv:2103.09848 [nucl-th] .
  * Plumberg  _et al._ (2021) Christopher Plumberg, Dekrayat Almaalol, Travis Dore, Jorge Noronha, and Jacquelyn Noronha-Hostler, “Causality violations in realistic simulations of heavy-ion collisions,” (2021), arXiv:2103.15889 [nucl-th] .
  * Adam  _et al._ (2016c) Jaroslav Adam _et al._ (ALICE), “Direct photon production in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=$ 2.76 TeV,” Phys. Lett. B 754, 235–248 (2016c), arXiv:1509.07324 [nucl-ex] .
  * Acharya  _et al._ (2019c) Shreyasi Acharya _et al._ (ALICE), “Direct photon elliptic flow in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV,” Phys. Lett. B 789, 308–322 (2019c), arXiv:1805.04403 [nucl-ex] .
  * Huss  _et al._ (2021) Alexander Huss, Aleksi Kurkela, Aleksas Mazeliauskas, Risto Paatelainen, Wilke van der Schee, and Urs Achim Wiedemann, “Predicting parton energy loss in small collision systems,” Phys. Rev. C 103, 054903 (2021), arXiv:2007.13758 [hep-ph] .
  * Brewer  _et al._ (2021) Jasmine Brewer, Aleksas Mazeliauskas, and Wilke van der Schee, “Opportunities of OO and $p$O collisions at the LHC,” in  _Opportunities of OO and pO collisions at the LHC_ (2021) arXiv:2103.01939 [hep-ph] .
  * Giacalone  _et al._ (2020) Giuliano Giacalone, Björn Schenke, and Chun Shen, “Observable signatures of initial state momentum anisotropies in nuclear collisions,” Phys. Rev. Lett. 125, 192301 (2020), arXiv:2006.15721 [nucl-th] .
  * Shen  _et al._ (2016b) C. Shen, J. F. Paquet, G. S. Denicol, S. Jeon, and C. Gale, “Thermal photon radiation in high multiplicity p+Pb collisions at the Large Hadron Collider,” Phys. Rev. Lett. 116, 072301 (2016b), arXiv:1504.07989 [nucl-th] .
  * Shen  _et al._ (2017) Chun Shen, Jean-François Paquet, Gabriel S. Denicol, Sangyong Jeon, and Charles Gale, “Collectivity and electromagnetic radiation in small systems,” Phys. Rev. C 95, 014906 (2017), arXiv:1609.02590 [nucl-th] .
  * (77) Vladimir Khachatryan _et al._ (PHENIX Collaboration), “PHENIX measurement of low momentum direct photon radiation from p+p and p+A collisions,” Parallel talk presentation at Quark Matter 2018.
  * Pratt  _et al._ (2015) Scott Pratt, Evan Sangaline, Paul Sorensen, and Hui Wang, “Constraining the Eq. of State of Super-Hadronic Matter from Heavy-Ion Collisions,” Phys. Rev. Lett. 114, 202301 (2015), arXiv:1501.04042 [nucl-th] .
  * Bernhard  _et al._ (2019) Jonah E. Bernhard, J. Scott Moreland, and Steffen A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma,” Nature Phys. 15, 1113–1117 (2019).
  * Nijs  _et al._ (2021) Govert Nijs, Wilke van der Schee, Umut Gürsoy, and Raimond Snellings, “Bayesian analysis of heavy ion collisions with the heavy ion computational framework Trajectum,” Phys. Rev. C 103, 054909 (2021), arXiv:2010.15134 [nucl-th] .
  * Everett  _et al._ (2021) D. Everett _et al._ (JETSCAPE), “Multisystem Bayesian constraints on the transport coefficients of QCD matter,” Phys. Rev. C103, 054904 (2021), arXiv:2011.01430 [hep-ph] .
  * Gale (2013) Charles Gale, “Electromagnetic radiation in heavy ion collisions: Progress and puzzles,”  _Proceedings, 5th International Conference on Hard and Electromagnetic Probes of High-Energy Nuclear Collisions (Hard Probes 2012): Cagliari, Italy, May 27-June 1, 2012_ , Nucl. Phys. A910-911, 147–154 (2013), arXiv:1208.2289 [hep-ph] .