# Thermalization in QCD:
theoretical approaches, phenomenological applications,
and interdisciplinary connections

Jürgen Berges<EMAIL_ADDRESS>Institut für Theoretische
Physik, Universität Heidelberg, Philosophenweg 16, D-69120 Heidelberg, Germany
Michal P. Heller<EMAIL_ADDRESS>_On leave from:_ National Centre
for Nuclear Research, Pasteura 7, Warsaw, PL-02093, Poland  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, D-14476 Potsdam, Germany  Aleksas Mazeliauskas
<EMAIL_ADDRESS>Theoretical Physics Department, CERN, CH-1211
Genève 23, Switzerland  Raju Venugopalan<EMAIL_ADDRESS>Physics Department,
Brookhaven National Laboratory, Upton, New York 11973-5000, USA

###### Abstract

Heavy-ion collisions at BNL’s Relativistic Heavy-Ion Collider (RHIC) and
CERN’s Large Hadron Collider (LHC) provide strong evidence for the formation
of a quark-gluon plasma, with temperatures extracted from relativistic viscous
hydrodynamic simulations shown to be well above the transition temperature
from hadron matter. How the strongly correlated quark-gluon matter forms in a
heavy-ion collision, its properties off-equilibrium, and the thermalization
process in the plasma, are outstanding problems in QCD. We review here the
theoretical progress in this field in weak coupling QCD effective field
theories and in strong coupling holographic approaches based on gauge-gravity
duality. We outline the interdisciplinary connections of different stages of
the thermalization process to non-equilibrium dynamics in other systems across
energy scales ranging from inflationary cosmology, to strong field QED, to
ultracold atomic gases, with emphasis on the universal dynamics of non-thermal
and of hydrodynamic attractors. We survey measurements in heavy-ion collisions
that are sensitive to the early non-equilibrium stages of the collision and
discuss the potential for future measurements. We summarize the current state-
of-the art in thermalization studies and identify promising avenues for
further progress.

††preprint: CERN-TH-2020-080

###### Contents

  1. I Introduction
    1. I.1 Big picture questions
    2. I.2 Outline and scope of the review
  2. II First principles approaches to the spacetime evolution of a heavy-ion collision
  3. III Hadron structure at high energies
    1. III.1 QCD at small $x$ and high parton densities
    2. III.2 Effective Field Theory for high parton densities: the Color Glass Condensate
    3. III.3 Renormalization group evolution in the CGC EFT
    4. III.4 DIS and the dipole model
    5. III.5 RG evolution and geometric scaling
    6. III.6 The state of the art in the CGC EFT
  4. IV Non-equilibrium QCD matter at high occupancy
    1. IV.1 Multi-particle production in strong fields
    2. IV.2 The LO Glasma: classical gluon fields from shockwave collisions
    3. IV.3 Non-perturbative evolution of high occupancy fields
      1. IV.3.1 Real time evolution of boost invariant fields on the lattice
      2. IV.3.2 Glasma flux tubes
      3. IV.3.3 The IP-Glasma model
    4. IV.4 The Glasma at NLO
      1. IV.4.1 Dynamics of $p^{\eta}=0$ modes: QCD factorization and energy evolution
      2. IV.4.2 Dynamics of $p^{\eta}\neq 0$ modes: plasma instabilities and the classical-statistical approximation
  5. V Far-from-equilibrium gluon and quark production: From plasma instabilities to non-thermal attractors
    1. V.1 Non-equilibrium time evolution equations from the quantum effective action
      1. V.1.1 Macroscopic field, spectral and statistical functions
      2. V.1.2 Resummed evolution equations to leading order
    2. V.2 Nonlinear evolution of plasma instabilities
      1. V.2.1 Dynamical power counting
      2. V.2.2 Classical-statistical field theory limit
    3. V.3 Non-thermal attractor
      1. V.3.1 Universal scaling far from equilibrium
      2. V.3.2 Identifying the weak-coupling thermalization scenario
      3. V.3.3 Non-thermal attractors in scalar field theories
    4. V.4 Separation of scales far from equilibrium and ultrasoft scale dynamics
      1. V.4.1 Non-equilibrium evolution of the spatial Wilson loop
      2. V.4.2 Effective condensate dynamics
    5. V.5 Early-time fermion production and quantum anomalies
      1. V.5.1 Real-time simulations for fermions and gauge fields beyond the classical-statistical approximation
      2. V.5.2 Real-time off-equilibrium dynamics of quantum anomalies
  6. VI Equilibration in QCD kinetic theory
    1. VI.1 The quasi-particle description of QCD plasmas
      1. VI.1.1 Chiral kinetic theory
    2. VI.2 Leading order kinetic theory
      1. VI.2.1 Elastic two-body scattering
      2. VI.2.2 Fokker-Planck limit of elastic scatterings
      3. VI.2.3 Effective collinear one-to-two processes
      4. VI.2.4 Bethe-Heitler and LPM limits of collinear radiation
    3. VI.3 Bottom-up thermalization
      1. VI.3.1 Stage one: collisional broadening
      2. VI.3.2 Stage two: collinear cascade
      3. VI.3.3 Stage three: mini-jet quenching
      4. VI.3.4 Numerical realization of bottom-up thermalization
    4. VI.4 Self-similar evolution in the high-occupancy regime
      1. VI.4.1 Self-similar scaling
      2. VI.4.2 Pre-scaling phenomenon
    5. VI.5 Extrapolation to stronger couplings
      1. VI.5.1 Hydrodynamic attractors in QCD kinetic theory
      2. VI.5.2 Entropy production during equilibration
      3. VI.5.3 Chemical equilibration of QGP
      4. VI.5.4 Equilibration of spatially inhomogeneous systems
  7. VII Ab initio holographic description of strong coupling phenomena
    1. VII.1 Holography and heavy-ion collisions
    2. VII.2 Controlled strong coupling regime
    3. VII.3 Early time physics of Bjorken flow in ${\cal N}=4$ SYM and other strongly-coupled CFTs
    4. VII.4 Hydrodynamic attractors in holography
    5. VII.5 Shockwave collisions in ${\cal N}=4$ SYM and other strongly-coupled CFTs
      1. VII.5.1 Planar shocks
      2. VII.5.2 Transverse dynamics in holography
    6. VII.6 Other topics in holography at its intersection with thermalization in heavy-ion collisions
      1. VII.6.1 Far from equilibrium physics in non-conformal strongly-coupled QFTs
      2. VII.6.2 Away from the strong coupling regime
      3. VII.6.3 Non-local correlators
  8. VIII Signatures of non-equilibrium QCD
    1. VIII.1 Hard Probes
    2. VIII.2 Long-range rapidity correlations
    3. VIII.3 Bulk observables
    4. VIII.4 Future prospects
  9. IX Interdisciplinary connections
    1. IX.1 Strong interactions: Unitary Fermi gas
    2. IX.2 Highly occupied systems I: Pre-heating in the early universe
    3. IX.3 Highly occupied systems II: Bose gases far from equilibrium
    4. IX.4 Highly occupied systems III: Classicalization and unitarization of gravitational amplitudes
    5. IX.5 Anomalous currents in non-equilibrium QED: Condensed matter systems and strong laser fields
    6. IX.6 Thermalization and entanglement
  10. X Summary and Outlook

##  I Introduction

###  I.1 Big picture questions

Ultrarelativistic collisions of heavy nuclei at the Relativistic Heavy-Ion
Collider (RHIC) and the Large Hadron Collider (LHC) produce several thousand
particles in each event generating the hottest and densest matter on Earth
Arsene  _et al._ (2005); Adcox  _et al._ (2005); Alver  _et al._ (2007); Adams
_et al._ (2005); Muller  _et al._ (2012); Roland  _et al._ (2014); Foka and
Janik (2016a, b). At the highest LHC energies, temperatures of the order of
five trillion Kelvin are attained Adam  _et al._ (2016). Temperatures on this
scale only previously existed at the earliest instants of our universe, a 10th
of a microsecond after the Big Bang. Lattice gauge theory studies Bazavov  _et
al._ (2019a) show strongly interacting matter at these temperatures to be well
over a crossover temperature from hadron matter to a regime where the degrees
of freedom describing bulk thermodynamic quantities are the fundamental quark
and gluon fields of Quantum Chromodynamics (QCD). The results of experimental
and theoretical studies indicate that shortly after the heavy-ion collision
the produced quark-gluon fields form a strongly correlated state of matter,
widely known as the quark-gluon plasma (QGP) Shuryak (1980).

The heavy-ion experiments at RHIC and LHC therefore provide us with a unique
opportunity to study terrestrially the spacetime evolution of this non-Abelian
QGP. A striking finding from the RHIC and LHC experiments is that the
experimental data are consistent with a description of the QGP as a nearly
perfect fluid with a very low value of shear viscosity to entropy density
ratio of $\upeta/s\leq 0.2$ (in natural units) Romatschke and Romatschke
(2019). These values are very close to $\upeta/s=1/(4\pi)$, a universal
property of a class of gauge theories with a large number of degrees of
freedom at infinite coupling Policastro  _et al._ (2001); Buchel and Liu
(2004); Kovtun  _et al._ (2005); Iqbal and Liu (2009).

While our understanding of the thermal properties of QGP matter and the flow
of the nearly perfect fluid has developed significantly, progress in
theoretical descriptions of the early stages of heavy-ion collisions is
relatively recent. In particular, there is a growing realization that the far-
from-equilibrium dynamics that characterizes early time physics is extremely
important in understanding collective phenomena in the heavy-ion experiments.

Since the thermalization process represents an initial value problem in
quantum field theory, one needs to understand the problem all the way from the
many-body correlations in the colliding hadrons, to how multi-particle
production occurs as they overlap in the collision, to the subsequent loss of
information in the thermalization of the matter produced. We will here outline
key questions prompted by the dynamics of each stage of the spacetime
evolution of quark-gluon matter as it thermalizes111For a complementary
perspective on open questions in heavy-ion collisions, we refer the reader to
Ref. Busza  _et al._ (2018).:

  * •

What are the many-body correlations of the strongly interacting matter in the
colliding nuclei?

In QCD, a proton (or any other nucleus) is not a state with a fixed number of
constituents. Instead, the proton must be view as a collection of short or
long lived configurations of partons (quarks, antiquarks, and gluons) each of
which carries the quantum numbers of the proton. At low energies, a probe of
the proton’s internal structure will see three quasi-particle “valence” quarks
carrying the quantum numbers of the proton. Configurations carrying larger
numbers of parton constituents are too short lived on the timescales over
which the proton interacts with the probe.

The picture changes significantly when the proton or nucleus is boosted to
high energies. Because of time dilation, short lived configurations containing
large numbers of partons live much longer; the probe interacting with the
hadron in this state will see a many-body configuration of partons. Indeed the
likelihood the probe will scatter off such a state is greater than if it were
to scatter off longer lived states with fewer constituents.

For this very reason, one expects such many-body configurations of partons to
dominate the physics of multi-particle production in the collisions of
ultrarelativistic nuclei. To understand how precisely this occurs requires
deep knowledge of the spatial and momentum distributions of partons in the
boosted nuclei, the nature of their correlations, and how these correlations
change with system size and with the collision energy.

An important complication is that each of the configurations that the proton
fluctuates into knows about confinement and other non-trivial features of the
QCD vacuum. The many-body parton description of the proton is robust for very
fine resolutions of its structure. As the resolution of the probe decreases,
this quasi-particle picture will break down completely and confining/vacuum
effects will be manifest. At what scale this occurs remains an open problem.

An interesting possibility we will discuss is the emergence of dynamical
scales, at shorter distances than the confining scale, that may dictate how
the proton scatters off the external probe.

  * •

What is the physics of the first Yocto-second of the collision?

The parton configurations that describe the substructure of a boosted nucleus
have their momenta distributed between fewer fast modes that carry the quantum
numbers of the nucleus and more plentiful soft modes. In this picture of
ultrarelativistic heavy-ion collisions, the fast modes in each of the two
nuclei interact relatively weakly with the other nucleus and populate the
“fragmentation regions” of the two nuclei corresponding to polar angles very
close to the beam axes Van Hove and Pokorski (1975). The slower degrees of
freedom interact more strongly with each other and produce strongly
interacting gluon matter.

This picture of initial conditions was proposed in a seminal paper by Bjorken
to describe the subsequent hydrodynamic flow of the quark-gluon plasma Bjorken
(1983) albeit he did not address how thermalization occurs in this spacetime
scenario. An interesting question in this regard is whether the strong
interactions of the soft modes with each other is due to strong coupling or
whether it can be due to the large occupancy of these soft modes. The answer
to this question may also influence the degree of transparency of the fast
modes, in particular a “limiting fragmentation” scaling phenomenon seen in
data.

A spacetime scenario in which both both soft and hard modes in the nuclei
interact very strongly and generate hydrodynamic flow was suggested by Landau.
Is there a transition between these two spacetime pictures with energy Gelis
_et al._ (2006a); Casalderrey-Solana  _et al._ (2013); if so, can they be
distinguished by phenomena such as limiting fragmentation Gonçalves  _et al._
(2019)?

  * •

Is there a unifying theoretical description of quark-gluon matter off-
equilibrium ?

The quark-gluon matter formed in the first few yoctoseconds of the heavy-ion
collision is very far from equilibrium. A key question in its description is
whether the relevant QCD coupling is large or small. There are two systematic
approaches in this regard. One employs QCD at very weak coupling
$\alpha_{S}\rightarrow 0$ and very high gluon occupancy $f_{g}$ satisfying
$\alpha_{S}f_{g}\sim 1$. The quark-gluon matter in this limit is called the
Glasma Lappi and McLerran (2006); Gelis and Venugopalan (2006a), and as we
will describe, a rich picture of strongly correlated dynamics in the far-from-
equilibrium non-Abelian matter emerges in this limit.

The other limit is that of very strong ’t Hooft coupling of
$\alpha_{S}N_{c}\rightarrow\infty$, as the number of colors
$N_{c}\rightarrow\infty$. In this limit, holographic approaches based on
gauge-gravity duality Maldacena (1999); Gubser  _et al._ (1998a); Witten
(1998a) are robust and can be used to obtain exact results in “QCD-like” non-
Abelian gauge theories. The best understood example of a holographic quantum
field theory is the ${\cal N}=4$ superconformal Yang-Mills theory. This and
other holographic setups, including theories with less supersymmetry and/or
with running couplings, are described by solutions of classical gravity with
matter in 5-dimensional anti-de Sitter geometries.

Ab initio calculations at strong coupling using holography predict the
applicability of hydrodynamics over a time scale set by the local energy
density, when the expanding matter in heavy-ion collisions settings is
characterized by a large spatial anisotropy in its energy-momentum tensor
Chesler and Yaffe (2010, 2011); Heller  _et al._ (2012a). This is at variance
with the presumption of local thermal equilibrium for the applicability of
hydrodynamics and has led to a paradigm shift in which the transition to
hydrodynamic flow is referred to as hydrodynamization rather than
thermalization Casalderrey-Solana  _et al._ (2014a).

A potentially rich line of inquiry is to isolate what features of the
nonequilibrium evolution of strongly correlated/coupled quark-gluon matter are
universal. One example is the approach to local thermal equilibrium governed
by viscous hydrodynamics, with the differences between frameworks manifesting
themselves in the values of transport coefficients. Another example is
universality in time dependence across a class of non-equilibrium states for
certain observables. In the weak coupling scenario, at high occupancies, these
are far-from-equilibrium attractors associated with non-thermal fixed points
Berges  _et al._ (2008a). A further example, for both weakly and strongly
coupled models, is the emergence of hydrodynamic attractors Heller and
Spalinski (2015); Romatschke (2018).

A related important set of questions concerns the use of effective theories
like hydrodynamics for systems far away from equilibrium. It’s been shown
recently that the solution of kinetic equations with quarks and gluons
exhibits the emergence of early hydrodynamic-like behavior around a far-from-
equilibrium state Mazeliauskas and Berges (2019), which is qualitatively
different from the more conventional formulation of hydrodynamics Romatschke
and Strickland (2003); Martinez  _et al._ (2012); Bazow  _et al._ (2014). A
quantitative answer to these questions can better clarify the sense in which
we understand information extracted from experiment to provide evidence for
the QGP as a thermal fluid.

Yet another line of inquiry is to determine how features of the dynamics
evolve between the weakly coupled and strongly coupled regimes. In particular,
it has been shown that some features of hydrodynamization appear to evolve
smoothly with $\upeta/s$ Keegan  _et al._ (2016a); Heller  _et al._ (2018);
Grozdanov  _et al._ (2016).

The role of topology in far-from-equilibrium dynamics is a fascinating
question with wide appeal. The idea of a Chiral Magnetic Effect (CME) Kharzeev
_et al._ (2008) was first introduced in the context of heavy-ion collisions
and is now much investigated in condensed matter physics Li and Kharzeev
(2016). Outstanding questions in the heavy-ion context have to do with the
role of topological transitions Mace  _et al._ (2016), their relation to the
infrared structure of the Glasma Berges  _et al._ (2017a), the persistence of
induced magnetic fields, and the development of chiral kinetic theory and
anomalous hydrodynamics describing the evolution of the CME Kharzeev  _et al._
(2016).

  * •

Signatures of quark-gluon matter off-equilibrium

If matter in bulk equilibrates fully in heavy-ion collisions, the only
information of the nonequilibrium evolution that survives is what is imprinted
as initial conditions for its subsequent hydrodynamic evolution. The
exceptions are electroweak and so-called “hard probes”; both of these are
sensitive to the full history of the spacetime evolution of QCD matter.

A significant development in recent years is the vastly improved ability of
the RHIC and LHC experiments to perform “event engineering” whereby final
states can be studied by varying the “control parameters” corresponding to
nuclear size, centrality of collision impact and final state multiplicities
(triggering thereby on typical versus rare event configurations) across a wide
range in energy and system size Schukraft  _et al._ (2013). Can we constrain
the current state of the art computational techniques to accurately reflect
the systematics of this event engineering, and further, to use these to
isolate empirically the out-of-equilibrium dynamics?

  * •

Interdisciplinary connections

The study of the out-of-equilibrium dynamics of strongly correlated systems is
an important topic of significant contemporary interest in a number of sub-
fields of physics. The ideas and methods outlined in this review have
significant overlap with these fields. Can one exploit these interdisciplinary
connections to make progress?

A concrete example is the concept of turbulent thermalization in the reheating
of the early universe following inflation Micha and Tkachev (2004). This idea
in turn is a relativistic generalization of the concept of weak wave
turbulence in fluids Zakharov  _et al._ (2012).

Another more recent development is that of universal behavior in the evolution
of the aforementioned Glasma and overoccupied cold atomic gases Berges  _et
al._ (2015a).

The search for effective theories far from equilibrium is a major current
research direction in quantum many-body physics and is of interest to a wide
range of communities working on complex systems. Understanding the
thermalization of closed quantum many-body systems is also a vigorous research
topic in condensed matter physics, integrable systems, quantum information and
related disciplines Eisert  _et al._ (2015).

Not least, are the connections to black holes and string theory illustrated by
holography but also with respect to general questions regarding the scrambling
of information Lashkari  _et al._ (2013); Maldacena  _et al._ (2016) and the
unitary dynamics underlying black hole formation and evaporation Hawking
(1974, 1976); Penington (2019); Almheiri  _et al._ (2019, 2020a).

###  I.2 Outline and scope of the review

Our review will attempt to cover up-to-date research on these big picture
questions to provide a comprehensive perspective of our understanding of far-
from-equilibrium strongly correlated quark-gluon matter and the approach to
thermalization in heavy-ion collisions. Some of the questions outlined have
seen significant theoretical developments in the last decade, or more,
allowing for a sophisticated formulation of the extant open questions. Others
remain open to conceptual and technical refinement.

In Section II, we will outline two ab initio approaches to the problem of
thermalization in QCD. As noted, one employs the weak coupling methods of QCD
in the limit of very high gluon occupancies. The other exploits the AdS/CFT
correspondence between strongly coupled QCD-like theories with large number of
colors and classical gravity in a higher dimensional anti-de Sitter geeometry.

The weak coupling approach to thermalization is by now highly developed. In
the asymptotic limit of very high energies (the Regge limit of QCD), the
thermalization process can be followed from a description of the matter in the
hadron wavefunction, through multi-particle production in the collision, the
non-equilibrium evolution of the matter produced, and its subsequent
equilibration.

A key element, discussed in Section III.1, is the emergence an energy-
dependent close packing “saturation” scale $Q_{s}$ Gribov  _et al._ (1983)
when the phase space of quarks and gluons (partons) in the wavefunctions of
the colliding hadrons becomes large.

In Section III.2, we motivate the construction of a QCD effective field theory
(EFT) called the Color Glass Condensate (CGC), which describes many-body
parton correlations in the hadron wavefunctions Gelis  _et al._ (2010);
Kovchegov and Levin (2012). A renormalization group (RG) framework Jalilian-
Marian  _et al._ (1998a); Iancu  _et al._ (2001) allows one to study the
energy evolution of parton many-body correlations as the hadron is boosted to
higher energies. This RG framework is outlined in Section III.3.

Deeply inelastic scattering (DIS) of electrons off protons and nuclei is a
powerful tool to probe the high energy structure of the latter. A “dipole”
model is discussed in Section III.4, wherein the interplay between the size of
the DIS quark-antiquark dipole and $Q_{S}$ is manifest in a phenomenon called
“geometrical scaling”.

The RG evolution of geometrical scaling is outlined in Section III.5, which
allows one to establish a strong correspondence between high energy QCD and
universal features of reaction-diffusion processes. This correspondence may be
particularly valuable in understanding hadron structure in rare events, such
as those probed in high multiplicity hadron-hadron collisions. Finally Section
III.6 summarizes the state-of-the-art on hadron structure in the CGC EFT.

The structure of the high energy hadron in this effective field theory is a
coherent state of static color sources and dynamical gluon fields. The
formalism for multi-particle production and the very early time evolution of
high occupancy parton matter called the Glasma Lappi and McLerran (2006) is
spelt out in Section IV.

The general quantum field theory formalism for multi-particle production in
the presence of strong time-dependent sources (with strength of order of the
inverse coupling) is discussed in Section IV.1. Inclusive quantities such as
multiplicities or energy densities, and their spacetime correlations, can be
computed systematically in the Glasma in powers of the coupling $\alpha_{S}\ll
1$.

At leading order in this power counting, the Glasma fields themselves, being
nonperturbative, are of magnitude $1/\alpha_{S}$. This is discussed in Section
IV.2, while Section IV.3 outline the shockwave initial conditions and
numerical realization Krasnitz and Venugopalan (1999) of the subsequent
spacetime evolution of Yang-Mills fields. In the latter sub-section, we also
discuss emergent stringy “Glasma flux tube” structure and an IP-Glasma model
motivated by these simulations. This model is used extensively as initial
conditions in hydrodynamic simulations of heavy-ion collisions.

The role of quantum fluctuations about the LO Glasma fields is discussed in
Section IV.4. These can be classified into two types. One of these,
corresponding to zero modes of the rapidity variable,“dress” the nuclear
wavefunctions through the RG procedure we describe in Section III.3. Such zero
modes factorize–to high accuracy; there is no cross-talk between the incoming
hadron wavefunctions Gelis  _et al._ (2008a). Non-zero quantum fluctuation
modes in the Glasma are unstable and display exponential growth Romatschke and
Venugopalan (2006a). We show how the physics of these unstable modes at very
early proper times $\tau\leq\frac{1}{Q_{S}}\ln^{2}(1/\alpha_{S})$ in the
Glasma is captured in a classical-statistical approximation.

Section V.1 addresses the non-linear evolution of Glasma fields in the terms
of the quantum two-particle-irreducible (2PI) effective action, which situates
the classical-statistical results discussed in Section IV.3 in the big picture
context of the real time evolution of quantum fields. In particular, in
Section V.2, we outline the dynamical power counting for the description of
plasma instabilities and the emergence of secondary instabilities due to non-
linear interactions of unstable modes. We discuss the range of validity of
classical-statistical field theory for the subsequent evolution, which
motivates fully 3+1-dimensional numerical simulations of the expanding Glasma
fields.

Section V.3 discusses the emergence of a non-thermal attractor described by a
self-similar gluon distribution, whose dependence on momentum, and an overall
cooling rate, are characterized by universal numbers independent of the
initial conditions. Because the numerical simulations correctly describe
dynamics in the infrared, the attractor solution helps identify the right
effective kinetic theory amongst several competing options.

Kinetic theory increasingly captures the relevant dynamics of the
thermalization process as the system expands and cools. This is prefigured in
the evolution of the Glasma, as Section V.4 discusses, by the emergence of
hard, semi-hard and soft scales characteristic of a weakly coupled plasma. In
Section V.5 we discuss other striking properties of the Glasma such as an
effective condensate, quantum anomalies, and fermion pair production.

As the expanding system cools and dilutes, quantum corrections that were sub-
leading in the classical-statistical description become important. The 2PI
framework allows one in principle to quantify these contributions; for smooth
spatio-temporal gradients, one expects emergent quasi-particle dynamics to
satisfy the Boltzmann kinetic theory, the subject of Section VI.1. In Section
VI.2, we discuss leading order kinetic theory framework, progressively from
elastic $2\leftrightarrow 2$ scatterings, to effective collinear
$1\leftrightarrow 2$ processes, with special note of interference and plasma
instability effects. The classical-statistical simulations point to a kinetic
theory equilibration scenario termed “bottom-up” thermalization. This multi-
stage process and the numerical implementation thereof are discussed in
Section VI.3. The self-similar evolution in the overlap region of kinetic and
classical-statistical descriptions are further discussed in Section VI.4.

The thermalization process discussed thus far is strictly valid for values of
the QCD coupling much smaller than realized in experiment. Section VI.5
addresses the problem in the language of hydrodynamic attractors. The
dependence on the coupling is replaced with a dimensionless combination of the
kinematic viscosity, the temperature, and the entropy density; reasonable
values are extracted for entropy production, as well as thermal and chemical
equilibration times, for values of the kinematic viscosity extracted from
hydrodynamic simulations.

In Sections VII.1 and VII.2, we provide an overview of holography based strong
coupling approaches to the early stages of heavy-ion collisions. Our focus is
on the conceptual features, universal mechanisms, and predictions from these
studies222Discussions with a focus on computational techniques are reviewed
for example in Refs. Chesler and Yaffe (2014); Ammon and Erdmenger (2015);
Baggioli (2019).

Thermalization in strongly-coupled gauge theories as an initial value problem
for higher-dimensional geometries is addressed in Section VII.2. We discuss in
particular attempts to mitigate the freedom in specifying initial states at
strong coupling to better model heavy-ion collisions.

In Section VII.3, we discuss 1+1-dimensional boost invariant flow where
hydrodynamization and hydrodynamic attractors were first discovered. Section
VII.4 continues the discussion of hydrodynamic attractors addressing
complementary aspects of this construction.

In analogy to the Glasma shockwave collisions discussed in Sections IV.3 and
IV.4, the holographic shockwave collisions discussed in Section VII.5 allow
one to address the breaking of boost invariance, emergence of insensitivity to
shockwave width, and the corresponding far-from-equilibrium signatures of
excited matter. Adding a non-trivial layer of detail, the transverse structure
of shockwaves allows one to study the development of the transverse flow
essential to the description of heavy-ion collisions.

Section VII.6 discusses efforts on more realistic holographic descriptions of
heavy-ion collisions that model confinement, the breaking of conformal
invariance, the running of the coupling, and large $N_{c}$ suppressed non-
local correlations.

Section VIII is devoted to a discussion of signatures of non-equilibrium
dynamics in heavy-ion data. While electromagnetic and high transverse momentum
strongly interacting final states are sensitive to early time dynamics,
significant contributions to their rates accrue from all stages of the
spacetime evolution of the system. The chiral magnetic effect too is
especially sensitive to early-time dynamics (since the electromagnetic fields
generated by charged spectator protons diminish very rapidly); the measurement
is similarly afflicted with large backgrounds that can obscure the signal.

Correlations amongst high momentum final states, especially those long range
in rapidity, are promising, as is “event engineering” the response of these
final states to variations in energy and system size. We also discuss how bulk
observables, in combination with these final states, can constrain
thermalization scenarios.

The thermalization process in QCD has strong (and in several instances very
concrete) interdisciplinary connections to other areas of physics. Some of
these, such as the formal analogy of color memory in QCD to gravitational
memory (Section III.2), of gluon saturation in hadron wavefunctions to
reaction-diffusion processes in statistical physics (Section III.5), or the
physics of early time instabilities to Weibel instabilities in plasma physics
(Section IV.4 and Section V.2), are discussed explicitly in the material
outlined. Our discussion of holography can of course, in its entirety, be
viewed in this light.

In Section IX, we elaborate further on some of these interdisciplinary
connections. As we observe in Section IX.1, amongst the first such connections
to be noted is the strikingly similar behavior of the QGP and unitary Fermi
gases. In Section IX.2, IX.3 and IX.4, we outline the qualitatively similar
early-time Glasma dynamics to other overoccupied systems, from inflationary
dynamics in the early universe to those of overoccupied ultracold Bose gases
to a quantum portrait of black holes as highly occupied graviton states. In
the case of the correspondence of the Glasma to overoccupied Bose gases, we
discuss in Section V.3 and VI.4, quantatitatively identical non-thermal
attractors, suggestive of a classification of far-from-equilibrium systems
into universality classes analogous to those for critical phenomena Hohenberg
and Halperin (1977). An exciting development is the advent of state-of-the-art
cold atom experiments that can provide deep insight into such universal
dynamics.

Powerful interdisciplinary connections are also observed between the chiral
magnetic effect and non-equilibrium dynamics in strongly correlated systems.
The continuing development of powerful lasers also offer great promise in the
precision study of anomalous currents off-equilibrium. Finally, but not least,
are the deep emerging connections of the thermalization process to the physics
of entanglement both in holography and in table top cold atom and condensed
matter systems. These interdisciplinary connections are briefly outlined in
Section IX.5 and IX.6.

We end the review in Section X with a brief summary and outlook to future
developments in our understanding of thermalization in QCD.

This outline suggests that thermalization in QCD is a rich field with many
research directions and we have had to make choices in our presentation due to
space limitations. Since our work is focused on theory developments, we have
not been able to do justice to experimental results and phenomenological
analyses of data. An important topic we do not address is the off-equilibrium
dynamics of QCD matter in the vicinity of a critical point Akamatsu  _et al._
(2019); Bzdak  _et al._ (2020); Bluhm  _et al._ (2020). Another is the related
topic of hydrodynamic fluctuations Akamatsu  _et al._ (2017); An  _et al._
(2019). Other noteworthy omissions in our presentation include the discussion
of holographic DIS Polchinski and Strassler (2003); Hatta  _et al._ (2008);
Shuryak and Zahed (2018), holographic hard probes Herzog  _et al._ (2006);
Gubser (2006); Chesler  _et al._ (2009a, b, 2013) and features of linear
response theory Kovtun and Starinets (2005). Some aspects of holographic
approaches that we omit or treat only partially are discussed in Casalderrey-
Solana  _et al._ (2014a); DeWolfe  _et al._ (2014); Chesler and van der Schee
(2015); Heller (2016); Florkowski  _et al._ (2018a).

##  II First principles approaches to the spacetime evolution of a heavy-ion
collision

Quantum Chromodynamics (QCD), the modern theory of the strong force in nature,
is a nearly perfect theory, the only free parameters being the quark masses
Wilczek (2000). The Lagrangian of the theory can be written compactly as

${\cal L}_{\rm
QCD}=-\frac{1}{4}F_{\mu\nu}^{a}F^{\mu\nu,a}+\sum_{f}{\bar{\Psi}}_{i}^{f}\left(i\gamma^{\mu}D_{\mu,ij}-m_{f}\delta_{ij}\right)\Psi_{j}^{f}\,.$
(1)

Here
$F_{\mu\nu}^{a}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}-gf^{abc}A_{\mu}^{b}A_{\nu}^{c}$
is the QCD field strength tensor for the color gauge fields $A_{\mu}^{a}$ that
live in the adjoint representation of SU(3), with $a=1,\cdots,8$ and $f^{abc}$
are the structure constants of the gauge group. The quark fields live in the
fundamental representation of SU(3) and are labeled by their color and flavor
indices $\Psi_{i}^{f}$ where the color index $i=1,\cdots,3$ and $f$ denotes
the flavors of quarks with masses $m_{f}$. Finally, contracted with the Dirac
matrix $\gamma^{\mu}$ is the covariant derivative
$D_{\mu,ij}=\partial_{\mu}\delta_{ij}+igt_{ij}^{a}A_{\mu}^{a}$, with
$t_{ij}^{a}$ being the generators of SU(3) in the fundamental representation.

The theory is rich in symmetry. The structure of the Lagrangian is dictated by
the invariance of the quark and gluon fields under local SU(3) color gauge
transformations. In addition, for massless quarks, the theory has a global
chiral $SU(3)_{L}\times SU(3)_{R}$ symmetry, global baryon number $U(1)_{V}$
and axial charge $U(1)_{A}$ symmetries, and the quark and gluon fields are
invariant under scale transformations. Not least, the Lagrangian is invariant
under discrete parity, charge and time reversal symmetries.

All of these symmetries, but color, are broken by vacuum/quantum effects that
give rise to all the emergent phenomena in the theory such as confinement,
asymptotic freedom, quantum anomalies and the spontaneous breaking of chiral
symmetry.

Because QCD is a confining theory, it is not analytically tractable in general
and numerical methods are essential to uncover its properties. Euclidean
lattice Monte Carlo methods can be applied to compute, with good accuracy,
“static” properties of the theory such as the mass spectrum of hadrons,
magnetic moments, and thermodynamic properties of QCD at finite temperature
Lin  _et al._ (2018); Detmold  _et al._ (2019).

These methods are however not very useful in determining dynamical “real time”
features of theory because of the contributions of a large number of paths to
the QCD path integral in Minkowski spacetime. There are promising approaches
to surmount this difficulty such as steepest descent Lefshetz thimble methods
but they are currently only applicable to problems in 1+1-dimensions Alexandru
_et al._ (2017). Likewise, quantum computing offers an alternative paradigm to
compute real time dynamics, but its applicability to QCD likely remains far in
the future Preskill (2018). From this perspective, a first principles
understanding of the rich dynamics revealed in heavy-ion collision experiments
at RHIC and the LHC appears insurmountable.

As we shall discuss at length, this perspective may be overly pessimistic.
However before we explain why, one should note that the production of high
transverse momentum and massive particles (jets, heavy quarkonia being two
notable examples) can be computed with high precision in QCD Collins  _et al._
(1989). This is because these processes correspond to very short transverse
distances and asymptotic freedom tells us that the QCD coupling $\alpha_{S}$
is weak at these scales.

Remarkably, the same may be true for multi-particle production in the very
high energy (“Regge”) limit of QCD. At these very high energies, or
equivalently small values of Bjorken $x_{\rm Bj}$, the number of gluons in the
hadron proliferate rapidly333In deeply inelastic scattering of electrons off
hadrons (DIS), the Regge limit corresponds to $x_{\rm Bj}\sim Q^{2}/s$ where
$Q^{2}$ is the squared four-momentum transfer and $s$ the squared center-of-
mass energy. In the parton model, $x_{\rm Bj}\approx x$, where $x$ is the
light cone fraction of the momentum of the hadron carried by the struck
parton. In hadron-hadron collisions, it is more appropriate to speak in terms
of momentum fractions, so we shall henceforth use $x$ instead of $x_{\rm
Bj}$.. This rapid growth of gluon distributions with decreasing $x$ was
observed at the HERA DIS collider Abt  _et al._ (1993); Ahmed  _et al._
(1995); Derrick  _et al._ (1993, 1995); Martin  _et al._ (1994); Lai  _et al._
(1995). It should saturate when the phase space occupancy $f\sim
1/\alpha_{S}$, corresponding the maximal close packing of gluons permissible
in QCD Nielsen and Olesen (1978); Savvidy (2020). The search for gluon
saturation is a major motivation for an Electron-Ion Collider (EIC) Aschenauer
_et al._ (2019).

Gluon saturation predicts that for partons within the hadron resolved at the
scale $Q^{2}\gg\Lambda_{\rm QCD}^{2}$ (where $\Lambda_{\rm QCD}$ is the
intrinsic QCD scale), there is always a value of $x$ below which close packing
occurs. This corresponds to an emergent saturation scale $Q_{S}(x)$; in the
Regge limit of $x\rightarrow 0$, $Q_{S}\rightarrow\infty$, giving
$\alpha_{S}(Q_{S})\rightarrow 0$ Gribov  _et al._ (1983); Mueller and Qiu
(1986).

In this picture, which is strongly motivated by the QCD parton model, multi-
particle production at high energies occurs when two hadrons collide and
release the highly occupied gluon clouds within Bjorken (1983). If
$\alpha_{S}(Q_{S})\ll 1$, one can compute key features of this nontrivial
dynamics ab initio, following the decoherence of the strong correlated gluons
(with $\alpha_{S}f\sim 1$) in the hadron wave functions through a
nonequilibrium temporal evolution, to thermalization, hydrodynamic evolution
and the formation of a quark-gluon plasma (QGP).

The properties of strongly correlated saturated gluons in nuclear
wavefunctions at high energies are described in a Color Glass Condensate (CGC)
effective field theory (EFT) McLerran and Venugopalan (1994a, b, c); Iancu and
Venugopalan (2003); Gelis  _et al._ (2010); Kovchegov and Levin (2012);
Blaizot (2017). In this framework, which we will elaborate on in the next
section, a systematic treatment of coherent multiple scattering, and the
energy evolution of these scatterings, is feasible to high orders in
$\alpha_{S}$.

Hadron-hadron collisions in the CGC framework can be described as the
collision of two gluon shockwaves. Further, factorization theorems allow one
to separate quantum fluctuations that describe the energy evolution of the
shockwaves from those that appear after the collision. Multi-particle
production in the collision generates, on time scales $\tau\sim 1/Q_{S}$,
overoccupied non-equilibrium matter called the “Glasma” Lappi and McLerran
(2006); Gelis and Venugopalan (2006a).

The decoherence and subsequent self-similar spacetime evolution of the Glasma
can be described using classical-statistical methods Romatschke and
Venugopalan (2006a, b). This approach breaks down when the occupancies in the
expanding fluid fall to $f\sim 1$. However, there is a regime of overlap
between classical-statistical evolution and kinetic theory Mueller and Son
(2004); Jeon (2005) which allows one to identify the right kinetic framework
for the subsequent evolution of the Glasma. As we will demonstrate in Section
VI, the QCD effective kinetic theory Arnold  _et al._ (2003a) successfully
describes the equilibration from the over-occupied initial state described by
classical fields to the near equilibrium final state of a fluid described by
viscous hydrodynamics.

Such a “bottom-up” thermalization scenario first outlined in Baier  _et al._
(2001) has been now confirmed by detailed numerical simulations of kinetic
theory Kurkela and Zhu (2015). While a clean separation of scales, and regimes
of applicability of the different stages of evolution, are only transparent at
very weak coupling, one can extrapolate the results to the realistic couplings
corresponding to heavy-ion collisions at RHIC and the LHC, obtaining sensible
results consistent with heavy-ion phenomenology.

In particular, one finds that viscous hydrodynamics is applicable at times
$\tau\sim 1\,\text{fm}/c$, approximately 3 Yocto-seconds after the collision.
We note that, thanks to a hydrodynamic attractor Heller  _et al._ (2013a),
relativistic viscous fluid dynamics can be applicable even far off-equilibrium
and it remains valid even if $\alpha_{S}$ becomes large enough to invalidate
the assumptions underlying kinetic theory.

Since the hydrodynamic description of heavy-ion collisions is quite advanced,
our treatment will not be comprehensive and we will refer the reader to
existing reviews on hydrodynamics for details Gale  _et al._ (2013a); Jeon and
Heinz (2015); Florkowski  _et al._ (2018a); Romatschke and Romatschke (2019).
Further, our discussion of signatures in Section VIII will focus on extant and
possible future signatures of non-equilibrium early time dynamics in
experiment.

There are a number of caveats to the weak coupling picture we have outlined.
Firstly, though gluon saturation is plausible, there is no smoking gun to date
from high energy experiments. Secondly, for the center-of-mass energies of
RHIC and the LHC, it is not clear that the saturation scales in the nuclei are
large enough for weak coupling to apply. A corollary is that the role of
strong coupling dynamics at distance scales $\sim 1/\Lambda_{\rm QCD}$ remains
to be quantified as energy scales are changed. This may be process dependent;
some quantities such as total cross-sections may be more sensitive to the
strongly coupled dynamics of confinement even at the very high energies where
saturation scales are large.

Further, the factorization theorems developed for A+A collisions have only
been proven444As we will argue, this may not impact our study of
thermalization; one cannot however claim the same of more differential many-
body correlations that are long range in rapidity Dusling  _et al._ (2010);
Iancu and Triantafyllopoulos (2013). to leading logarithmic accuracy in $x$
Gelis  _et al._ (2008a, b). Not least, it is unclear whether the extrapolation
from very weak couplings to realistic ones can be performed without
encountering additional dynamics that can modify the conclusions drawn from
naive extrapolations555It should be noted that weak coupling approaches are
non-perturbative at high occupancy and contain non-trivial dynamics in the
infrared. How these scale with coupling is an interesting problem..

This motivates our looking to strong coupling holography for lessons in the
context of ultra-relativistic heavy-ion collisions where it is not obvious
whether the relevant coupling is large or small. The holographic paradigm is
under full theoretical control when strong coupling governs the dynamics of
certain QCD-like gauge theories at all scales. Note that this does not
necessarily imply scale invariance. While real time dynamics in QCD is at
present completely intractable in the regime of strong coupling, for
holographic gauge theories it reduces to a problem of finding time dependent
solutions of higher dimensional Einstein’s equations with negative
cosmological constant and, possibly, matter fields. The latter is a technical
problem that existing numerical relativity techniques allow one to solve.

From the numerical relativity point of view, the key novelty of holography was
the fact that one then solves Einstein’s equations for spacetimes that
effectively look like a box with a Lorentzian boundary, an example being
Minkowski space. This boundary is the locus where the corresponding quantum
field theory lives and asymptotic boundary conditions arise as values of
sources for associated local gauge invariant operators such as the energy-
momentum tensor.

Therefore the asymptotic boundary conditions are motivated in holography by
the quantum field theory physics one is interested in. Such sources can be
completely trivial, in which case the strongly coupled quantum field theory is
conformally invariant and lives in Minkowski space. Alternately, one can
trigger a renormalization group flow by providing a constant asymptotic value
to a scalar field corresponding to a relevant operator; further, one can turn
on nontrivial sources only for a finite duration of time to create an
interesting far-from-equilibrium state from the vacuum, thermal, or any other
desired state.

From the unifying perspective of the present review, the higher dimensional
metric is the analog of the distribution function in the kinetic theory
description of weakly coupled quantum field theory and the Einstein equations
are the analog of the Boltzmann equation. While these are two very different
mathematical notions, explicit comparisons Keegan  _et al._ (2016a); Heller
_et al._ (2018) can be made at the level of the spatio-temporal dependence of
the expectation value of the energy-momentum tensor. For kinetic theory, it
appears as the second moment of the distribution function. In holography,
standard techniques allow one to reconstruct 1-point functions of local gauge
invariant operators (such as the energy-momentum tensor) in quantum field
theory states using higher dimensional geometries.

The applications of holography to non-equilibrium processes in quantum field
theories is a vast subject. We will focus in this review on its applications
to the early stages of heavy-ion collisions and to the pre-equilibrium
dynamics of the expanding plasma formed. The key quantity of interest will be
the expectation value of the energy-momentum tensor because it can be directly
compared to the predictions of weak coupling approaches to QCD. It can also be
employed as an initial condition for hydrodynamic codes which generate the
bulk spectra of particles, the systematics of which can be compared to results
from heavy-ion experiments.

A powerful guiding principle potentially unifying both weak and strong
coupling approaches is universality, wherein systems across energy scales with
vastly different microscopic reactions can demonstrate dynamical features on
larger scales that are independent of initial conditions and can be
characterized by the same scaling functions and scaling exponents. A deeper
understanding is provided by the renormalization group classification of this
dynamics into universality classes. While this classification is a powerful
and highly developed paradigm in critical phenomena, its application as a
governing principle in off-equilibrium real time dynamics is relatively
recent.

In particular, in addition to the AdS/CFT correspondence, we will discuss a
remarkable universality between weakly coupled (albeit highly occupied,
$\alpha_{S}f\sim 1$) Yang-Mills fields in QCD and highly occupied ultracold
Bose gases Berges  _et al._ (2015a). Since the latter are described by the
infrared dynamics of self-interacting scalar field theories, such a
correspondence, if robust, might allow one to employ the sophisticated many-
body techniques developed for highly occupied scalar field theories to obtain
insight into self-similar dynamics in gauge theories Berges  _et al._ (2015b).

##  III Hadron structure at high energies

The thermalization process in hadron-hadron collisions, viewed as an initial
value problem, requires a deep understanding of the structure of QCD matter in
the wavefunctions of the colliding hadrons. This was already appreciated in
the early days of QCD, where it was argued that the highly Lorentz contracted
large $x$ valence partons in the ultrarelativistic hadron wavefunctions go
through unscathed in the collision, while their accompanying small $x$ “fur
coat of wee-parton vacuum fluctuations” Bjorken (1976) interacts strongly to
form hot and dense matter Bjorken (1983). In this section, we will discuss
significant developments in this picture of the hadron wavefunction and in the
next, we shall discuss how this picture translates ab initio into multi-
particle production after a heavy-ion collision.

###  III.1 QCD at small $x$ and high parton densities

A great success of perturbative QCD (pQCD) is the QCD parton model Bjorken and
Paschos (1969), wherein the complex dynamics of quark and gluon fields in
hadrons can, at high energies and large momentum resolutions, be viewed as
that of a weakly interacting gas of partons (single-particle quark, antiquark
and gluon states). More specifically, it is argued that if one picks a
lightcone666Lightcone coordinates are $k^{\pm}=(k^{0}\pm k^{3})/\sqrt{2}$ and
lightcone fields are defined as $A^{\pm}=(A^{0}\pm A^{z})/\sqrt{2}$; we will
work in the metric $g_{\pm,\mp}=1;g_{i,j}=-1$, where $i,j$ represent the two
transverse coordinates. gauge $A^{+}=0$, and quantizes the quark and gluon
fields of QCD along a light front surface (say, $x^{+}=0$), the Hamiltonian of
free quark and gluon fields shares the same vacuum777This argument is flawed
upon closer examination because Haag’s theorem requires the vacua of free and
interacting relativistic theories to have inequivalent commutation relations
of field operators. In lightcone quantization, this requires a careful
treatment of $k^{+}=0$ vacuum modes Nakanishi and Yamawaki (1977). However for
the purposes of a perturbative treatment of lightcone wavefunctions, it may be
sufficient to assume one can project out such modes Collins (2018);
Fitzpatrick  _et al._ (2018). as the full interacting theory Brodsky  _et al._
(1998). This allows one to construct the hadron wavefunction as a linear
combination of a complete set of multi-parton eigenstates, each of which is an
eigenstate of the free QCD Hamiltonian.

In this lightcone framework, the parton distribution functions measured in DIS
experiments can be interpreted as one-body states of quarks and gluons that
carry a lightcone momentum fraction $x=k^{+}/P^{+}$, where $k^{+}$ is the
parton’s lightcone momentum and $P^{+}$ the lightcone momentum of the hadron.
As noted previously in Section II, the HERA DIS experiments revealed that
gluon distributions grow very rapidly with decreasing $x$ at fixed large
resolution scales $Q^{2}$. This is consistent with the predictions of the
DGLAP Gribov and Lipatov (1972); Lipatov (1975); Altarelli and Parisi (1977);
Dokshitzer (1977) evolution equations of pQCD; in the so-called double
logarithmic limit of small $x$ and large $Q^{2}$, the change in the nuclear
gluon distribution $xG_{A}(x,Q^{2})$ with $\log(Q^{2})$ is proportional to the
gluon distribution itself, leading to its rapid growth888Here and henceforth,
we will identify the Bjorken variable of DIS $x_{\rm Bj}$ with $x$. The
identification of the two is robust only in the simple parton model. There are
corrections to this identity in pQCD that can be computed systematically..

Limiting the discussion to parton distributions is to focus only on one-body
lightcone distributions. However as first argued in Gribov  _et al._ (1983)
and Mueller and Qiu (1986), two-body “higher twist” distributions in a
lightcone operator product expansion999In OPE language, these higher twist
contributions are suppressed by powers of $1/Q^{2}$. (OPE) corrections grow as
$(xG_{A}(x,Q^{2}))^{2}$. For fixed $Q^{2}$, these become as large as the
leading twist one-body parton distribution as $x\rightarrow 0$.

More importantly, the net contribution of such many-body
contributions101010These include the screening of bremsstrahlung gluons by
real and virtual gluons, and the recombination of softer gluons into harder
gluons. is opposite to that of the leading term, softening the growth in the
gluon distribution. When the gluon phase space density is maximal, of order
$1/\alpha_{S}$, all n-body lightcone distributions contribute equally. This
saturation of gluon distributions in a nucleus of radius $R_{A}$, corresponds
to the generation of the saturation scale $Q_{S}$, where parametrically, for
$Q^{2}=Q_{S}^{2}$ the maximal occupancy on the l.h.s is equated to the gluon
phase space density:

$\frac{1}{\alpha_{S}(Q_{S})}=\frac{xG_{A}(x,Q_{S}^{2})}{2(N_{c}^{2}-1)\pi
R_{A}^{2}Q_{S}^{2}}\,.$ (2)

Fig. 1 illustrates the gluon saturation phenomenon and the interpretation of
$Q_{S}$ as the emergent “close packing” scale.

Figure 1: Transverse hadron profile resolved in scattering with fixed squared
momentum transfer $Q^{2}$ and increasing center-of-mass energy $\sqrt{s}$. The
requirement that proliferating soft gluons have maximal occupancy
$1/\alpha_{S}$ generates the close packing saturation scale $Q_{S}$. Figure
adapted from Iancu and Venugopalan (2003).

###  III.2 Effective Field Theory for high parton densities: the Color Glass
Condensate

Since the usual formalism of pQCD relies on two-body and higher twist
distributions being small, an alternative framework is necessary to understand
the physics of gluon saturation. Fortuitously, the problem of high parton
densities can be formulated as a classical effective field theory (EFT) on the
light front. We will outline here an explicit construction performed for
nuclei with large atomic number $A\gg 1$ McLerran and Venugopalan (1994a, b,
c).

An important ingredient in this construction (in the infinite momentum frame
(IMF) $P^{+}\rightarrow\infty$ of the nucleus) is a Born-Oppenheimer
separation in time scales between the Lorentz contracted large $x$ ($k^{+}\sim
P^{+}$) “valence” modes and the noted “wee fur” of small $x$ ($k^{+}\ll
P^{+}$) gluons and “sea” quark-antiquark pairs. For partons of transverse
momentum $k_{\perp}$, their lightcone lifetimes are given by

$\displaystyle\tau_{\rm wee}$ $\displaystyle=$ $\displaystyle{1\over
k^{-}}={2k^{+}\over k_{\perp}^{2}}\equiv{2xP^{+}\over k_{\perp}^{2}}\,\,\,{\rm
and}$ $\displaystyle\tau_{\rm valence}$ $\displaystyle\approx$
$\displaystyle{2P^{+}\over k_{\perp}^{2}}\,\longrightarrow\,\tau_{\rm
wee}\ll\tau_{\rm valence}\,,$ (3)

suggesting that the valence parton modes are static over the times scales over
which wee modes are probed. However one cannot integrate out the valence
sources completely out of the theory because they are sources of color charge
for wee partons and must couple to these in a gauge invariant manner.

Note further that since wee partons have large lightcone wavelengths
($\lambda_{\rm wee}\sim 1/k^{+}=1/xP^{+}$), they can resolve a lot of color
charge provided their transverse wavelength is not too large. The inequality

$\lambda_{\rm wee}\sim{1\over k^{+}}\equiv{1\over xP^{+}}\gg\lambda_{\rm
valence}\equiv{R_{A}\,m_{N}\over P^{+}}\,,$ (4)

where on the r.h.s the Lorentz contraction factor is $P^{+}/m_{N}$ (with
$m_{N}$ the nucleon mass), suggests that wee partons with $x\ll A^{-1/3}$
resolve partons111111Wee partons with wavelength $k_{\perp}\leq\Lambda_{\rm
QCD}\sim 1$ fm-1, see no color charge at all since color is confined (in
nucleons!) on this scale. It is only wee partons with
$k_{\perp}\gg\Lambda_{\rm QCD}$ that see color charges from different nucleons
along the longitudinal direction. all along the longitudinal extent
$2\,R_{A}\sim A^{1/3}$ in units of the inverse nucleon mass.

These charges will be random since they are confined to different nucleons and
do not know about each other. A wee parton with momentum $k_{\perp}$ resolves
an area in the transverse plane $(\Delta x_{\perp})^{2}\sim 1/k_{\perp}^{2}$.
The number of valence partons it interacts simultaneously with is

$k\equiv k_{(\Delta x_{\perp})^{2}}={N_{\rm valence}\over\pi
R_{A}^{2}}\,(\Delta x_{\perp})^{2}\,,$ (5)

which is proportional to $A^{1/3}$ since $N_{\rm valence}=3\cdot A$ in QCD.
For a large nucleus with $k\gg 1$, one can show for $N_{c}\geq 2$ that the
most likely color charge representation that the wee gluons couple to is a
higher dimensional classical representation of order $\sqrt{k}$ Jeon and
Venugopalan (2004).

Thus wee partons couple to $\rho$, the classical color charge charge per unit
transverse area of large $x$ sources. On average, since the the charge
distributions are random, the wee partons will couple to zero charge; however,
fluctuations locally can be large. These conditions can be represented as

$\langle\rho^{a}(x_{\perp})\rangle=0\,;\,\langle{\rho^{a}(x_{\perp})\rho^{b}(y_{\perp})}\rangle=\mu_{A}^{2}\,\delta^{ab}\,\delta^{(2)}(x_{\perp}-y_{\perp})\,,$
(6)

where $a=1,\cdots,N_{c}^{2}-1$ and $\mu_{A}^{2}={g^{2}A\over 2\pi R_{A}^{2}}$
is the color charge squared per unit area. For a large nucleus ($A\gg 1$),
$\mu_{A}^{2}\propto A^{1/3}\gg\Lambda_{\rm QCD}^{2}$ is a large scale. Since
it is the largest scale in the problem, $\alpha_{S}(\mu_{A}^{2})\ll 1$. This
result is remarkable because it provides a concrete example suggesting that
QCD at small $x$ is a weakly coupled EFT wherein systematic computations of
its many-body properties are feasible.

We can now put together the kinematic and dynamical arguments above and write
down the generating functional for the small $x$ effective action as

${\cal
Z}[j]=\int[d\rho]\,W_{\Lambda^{+}}[\rho]\,\left\\{{\int^{\Lambda^{+}}[dA]\delta(A^{+})e^{iS_{\Lambda^{+}}[A,\rho]-\int
j\cdot
A}}\over{\int^{\Lambda^{+}}[dA]\delta(A^{+})e^{iS_{\Lambda^{+}}[A,\rho]}}\right\\}\,.$
(7)

Here $\Lambda^{+}$ denotes the longitudinal momentum scale that separates the
static color sources from the dynamical gauge fields and the gauge invariant
weight functional $W_{\Lambda^{+}}[\rho]$ describes the distribution of these
sources at the scale $\Lambda^{+}$, with its path integral over $\rho$
normalized to unity.

The CGC effective action can be written in terms of the sources $\rho$ and the
fields $A$ as

$\displaystyle S_{\Lambda^{+}}[A,\rho]={1\over 4}\int
d^{4}x\,F_{\mu\nu}^{a}\,F^{\mu\nu,a}$ (8) $\displaystyle+$
$\displaystyle{i\over N_{c}}\int d^{2}x_{\perp}dx^{-}\delta(x^{-}){\rm
Tr}\left(\rho\,U_{-\infty,\infty}[A^{-}]\right)\,.$

The first term here is the Yang-Mills action in the QCD Lagrangian given in
Eq. (1). Since wee gluons in the CGC are treated in fully generality, their
dynamics is specified by the Yang-Mills action. The second term denotes the
coupling of the wee gluon fields to the large $x$ color charge densities
$\rho$, which we have argued are static lightcone sources. Because the sources
are eikonal sources along the lightcone, their gauge invariant coupling to the
wee fields is described by the path ordered exponential along the lightcone
time direction121212The second term in the effective action can alternatively
Jalilian-Marian  _et al._ (2001) be written as ${\rm
Tr}\left(\rho\log(U_{-\infty,\infty})\right)$; for the discussion of interest,
they provide identical results. $U_{-\infty,\infty}={\cal P}\exp\left(ig\int
dx^{+}A^{-,a}T^{a}\right)$. Physically, $U$ corresponds to the color rotation
of the color sources in the background of wee gluon fields.

The weight functional in the effective action (for the Gaussian random color
charges in Eq. (6)), in what is now called the McLerran-Venugopalan (MV) model
McLerran and Venugopalan (1994a, b); Kovchegov (1999), can equivalently be
written as131313The next sub-leading term (parametrically suppressed as
$A^{-1/6}$) is proportional to the SU(3) cubic Casimir
$d_{abc}\rho^{a}\rho^{b}\rho^{c}$ Jeon and Venugopalan (2005) and further
corrections are discussed in Dumitru and Petreska (2012).

$W_{\Lambda^{+}}[\rho]=\exp\left(-\int
d^{2}x_{\perp}{\rho^{a}\rho^{a}(x_{\perp})\over 2\mu_{A}^{2}}\right)\,.$ (9)

For the static configuration of $\rho$’s in Eq. (7), the saddle point of the
effective action is given by the Yang-Mills (YM) equations:

$D_{\mu}F^{\mu\nu,a}=\delta^{\nu+}\,\delta(x^{-})\,\rho^{a}(x_{\perp})\,,$
(10)

whose solution is the non-Abelian analog of the Weizäcker–Williams fields in
classical electrodynamics. The chromo-electromagnetic gluon field strengths
are singular on the nuclear sheet of width $\Delta x^{-}\sim 2Rm_{N}/P^{+}$
and pure gauge outside. The gauge fields in lightcone gauge are given by
$A^{-}=0$ and

$A_{\rm cl}^{k}={1\over
ig}V(x^{-},x_{\perp})\nabla^{k}V^{\dagger}(x^{-},x_{\perp})\,,$ (11)

where $k=1,2$ are the transverse coordinates and $V={\cal
P}\exp\left(\int_{-\infty}^{x^{-}}dz^{-}{1\over\nabla_{\perp}^{2}}{\tilde{\rho}}(z^{-},x_{\perp})\right)$.
The path ordering here is in $x^{-}$. This solution of the equations of motion
requires smearing of the sources in $x^{-}$ Jalilian-Marian  _et al._ (1997);
Kovchegov (1999). Further, ${\tilde{\rho}}$ that appears in the solution is
the color charge density in Lorenz gauge $\partial_{\mu}{A^{\prime}}^{\mu}=0$,
where one has the solution ${A_{\rm
cl}^{\prime}}^{+}={1\over\nabla_{\perp}^{2}}{\tilde{\rho}}(x^{-},x_{\perp})$,
${A_{\rm cl}^{\prime}}^{-}={A_{\rm cl}^{\prime}}_{\perp}=0$.

The explicit solution of the gauge field in lightcone gauge is therefore in
terms of ${\tilde{\rho}}$. However since we are interested in color averages
over products of gauge fields, one can replace
$[d\rho]\rightarrow[d{\tilde{\rho}}]$ in the path integral because the
Jacobian in the transformation is simple Jalilian-Marian  _et al._ (1997) and
does not contribute to correlation functions. Therefore many-body
distributions in lightcone gauge can be computed straightforwardly by
expressing them in terms of color charges in covariant gauge.

As a simple example, the number distribution of gluons in a large nucleus can
be computed by averaging the solution in Eq. (11) with the weight functional
$W$,

$\langle AA\rangle_{\rho}=\int[d{\tilde{\rho}}]A_{\rm cl.}[\rho]A_{\rm
cl.}[{\tilde{\rho}}]W_{\Lambda^{+}}[{\tilde{\rho}}]\,.$ (12)

For the Gaussian weight in the MV model, one obtains an analytical
solution141414In the MV model, $Q_{s}^{2}=c_{A}\mu_{A}^{2}$, where the
coefficient $c_{A}$ is be determined numerically Lappi (2008).. Defining the
occupation number to be $\phi=\frac{(2\pi)^{3}}{2(N_{c}^{2}-1)}\frac{dN}{\pi
R^{2}d^{2}k_{\perp}dy}$, one obtains at large transverse momenta
($k_{\perp}\gg Q_{s}$), $\phi\propto\frac{Q_{s}^{2}}{k_{\perp}^{2}}$, where
the $k_{\perp}$ dependence is that of Weizsäcker-Williams fields. However at
smaller transverse momenta ($k_{\perp}\ll Q_{s}$), the distribution is
modified substantially: $\phi\sim\frac{1}{\alpha_{S}}\log(Q_{s}/k_{\perp})$.

Thus in lightcone gauge, the strongly non-linear behavior of the fields is
responsible for the softening of the infrared behaviour of the classical
fields. This non-linearity is responsible for the phenomenon of saturation and
is seen already at the classical level in in the EFT construction of the MV
model.

We are now in a position to understand the term Color Glass Condensate (CGC)
Iancu and Venugopalan (2003); Gelis  _et al._ (2010) used to describe the
ground state properties of a hadron/nucleus at very high energies. Color is
obvious since the state is comprised primarily of a large number of gluons and
“sea” quark-antiquark pairs. It is a condensate because the gluons have
occupation numbers $\phi\sim 1/\alpha_{S}$ and have momenta peaked at
$k_{\perp}\sim Q_{s}$. Finally, it is a glass because small $x$ gluons and sea
quarks are generated by random sources with lifetimes much longer than the
characteristic time scales of the scattering. This explains the structure of
the path integral in Eq. (7), where the path integral over the curly brackets
is performed first for fixed color charge distributions and then averaged over
an ensemble of such diistributions.

To take a specific example, consider the inclusive cross-section in the DIS
scattering of a virtual photon on the nucleus. This is illustrated in Fig. 2,
and in the CGC EFT is expressed as the cross-section for a fixed distribution
of sources convoluted with an ensemble of such sources:

$\langle\mathrm{d}\sigma\rangle=\int[\mathcal{D}{\tilde{\rho}}_{A}]\,W_{\Lambda^{+}}[{\tilde{\rho}}_{A}]\,\mathrm{d}\hat{\sigma}[{\tilde{\rho}}_{A}]\,.$
(13)

The interpretation of this expression is that on the time scale ($t\sim 1/Q$,
where $Q$ is the virtuality of the scattering), the probe resolves a colored
condensate of gluons with a well-defined number density of longitudinal modes
down to $x\sim x_{\rm Bj}\ll 1$. Due to time dilation, the averaging over
$\rho_{A}$ with $W$ takes place on much larger time scale, as indicated by Eq.
(3). This two-stage averaging process clarifies the putative conundrum of how
to reconcile gauge invariance with the presence of a colored condensate.

Figure 2: DIS in the dipole picture. The virtual photon emitted by the
electron splits into a $q{\bar{q}}$ dipole which scatters off dynamical small
$x$ gauge fields coupled to the static large $x$ lightcone sources. Figure
from Iancu and Venugopalan (2003).

The CGC classical equations possess a “color memory” effect Pate  _et al._
(2017) corresponding to the large gauge transformation $V$ of a quark after
interacting with the gluon shockwave. This color memory effect generates a
transverse momentum kick $p_{\perp}\sim Q_{S}$ to the quark, and could be
measured in DIS experiments Ball  _et al._ (2019). Remarkably, this is exactly
analogous to the inertial displacement of detectors after the passage of a
gravitational shockwave Strominger and Zhiboedov (2016). In the latter case,
this gravitational memory is deeply related to asymptotic spacetime symmetries
of gravity at null infinity (in a Penrose spacetime diagram) and likewise
related to soft graviton theorems151515An “infrared triangle” between
asymptotic symmetries, memory and soft theorems Strominger (2017) allows for
an elegant interpretation of the infrared structure of QED Bieri and Garfinkle
(2013); Kapec  _et al._ (2017).. One might imagine such infrared relations are
not applicable to QCD because of color confinement. However because
$Q_{S}\gg\Lambda_{\rm QCD}$, QCD is weakly coupled in the Regge limit. The
identification of color memory in the CGC may therefore help isolate universal
features of the infrared structure of gauge theories and gravity.

###  III.3 Renormalization group evolution in the CGC EFT

We discussed thus far a classical EFT for large nuclei and Gaussian sources
where the separation between fields (wee partons) and sources (valence
sources) was picked randomly to be at the momentum scale $\Lambda^{+}$.
Physical observables such as the inclusive cross-section in Eq. (13) should
not depend on $\Lambda^{+}$. This invariance is the essence of the
renormalization group and is realized as follows in the EFT.

Consider the NLO contributions to Eq. (13) illustrated in Fig. 3. The lower
horizontal dashed line represents the LO separation of color sources in the
target from the fields at an $x=\Lambda^{+}/P^{+}$. Quantum fluctuations in
the classical background field of the target, illustrated by the one loop real
(top) and virtual (bottom) diagrams in Fig. 3, while apparently of
${\mathcal{O}}(\alpha_{S})$ are actually
$\sim\alpha_{S}\log(\Lambda^{+}/{\Lambda^{\prime}}^{+})$ from the phase space
integration of these modes between the two horizontal lines, and
${\mathcal{O}}(1)$ when ${\Lambda^{\prime}}^{+}=\Lambda^{+}e^{-1/\alpha_{S}}$
(or equivalently, when $x_{\rm wee}=x_{\rm val.}\,e^{-1/\alpha_{S}}$).

Figure 3: Renormalization group evolution of the hadron wavefunction in DIS at
high energies. Top figure: Real quantum fluctuations of gluons with
longitudinal momenta ${\Lambda^{\prime}}^{+}\ll k^{+}\ll\Lambda^{+}$ in a high
energy hadron with lightcone momentum $P^{+}$. The dark blobs denote dressed
gluon propagators including multiple coherent scattering off the target;
likewise, the $\otimes$ symbol denotes multiple scattering of target gluons
off the quark-antiquark pair. Bottom figure: Virtual quantum correcctions.
${\tilde{\rho}}$ denotes the color charge density of color sources at the
scale $\Lambda^{+}$ and ${\tilde{\rho}}^{\prime}$ is the charge density of
sources after evolution in $x$ to the scale ${\Lambda^{\prime}}^{+}$. The
vertical dashed line represents the cut separating the eamplitude from the
complex conjugate amplitude. Figure adapted from Roy and Venugopalan (2020).

However these large NLO contributions can be absorbed in the LO form of the LO
cross-section in Eq. (13) at the scale ${\Lambda^{\prime}}^{+}$ in Fig. 3 by
redefining the weight functional $W_{\Lambda^{+}}[\rho]\rightarrow
W_{{\Lambda^{\prime}}^{+}}[\rho^{\prime}]$. Here
${\tilde{\rho}}^{\prime}={\tilde{\rho}}+\delta{\tilde{\rho}}$, is new color
source density at ${\Lambda^{\prime}}^{+}$ that incorporates the color charge
density $\delta{\tilde{\rho}}$ induced by quantum fluctuations between
$\Lambda^{+}$ and ${\Lambda^{\prime}}^{+}$. One can thus write

$\langle\mathrm{d}\sigma_{\text{LO+NLO}}\rangle=\int[\mathcal{D}{\tilde{\rho}}_{A}]\,W_{{\Lambda^{\prime}}^{+}}[{\tilde{\rho}}_{A}]\,\mathrm{d}\hat{\sigma}_{\text{LO}}[{\tilde{\rho}}_{A}]\,,$
(14)

where

$W_{{\Lambda^{\prime}}^{+}}[{\tilde{\rho}}_{A}]=\Big{(}1+\log({\Lambda^{+}}/{\Lambda^{\prime}}^{+})\mathcal{H}_{\text{LLx}}\Big{)}W_{\Lambda+}[{\tilde{\rho}}_{A}]\,,$
(15)

with the quantum fluctuations absorbed into the “JIMWLK Hamiltonian”
$\mathcal{H}_{\text{LLx}}$ we shall discuss further shortly.

Since the l.h.s of Eq. (14) should not depend on the arbitrary “factorization
scale” $\Lambda^{+}$, the derivative of both l.h.s and r.h.s with respect to
it should be zero. From Eq. (15), one can therefore deduce the
JIMWLK161616JIMWLK stands for the last names of the principal authors. RG
equation Jalilian-Marian  _et al._ (1998b, a); Iancu  _et al._ (2001)

$\frac{\partial}{\partial
Y}W_{Y}[{\tilde{\rho}}_{A}]=\mathcal{H}_{\text{LLx}}\,W_{Y}[{\tilde{\rho}}_{A}]\,,$
(16)

where the JIMWLK Hamiltonian Weigert (2002)

${\mathcal{H}}_{\rm LLx}={1\over
2}\,\int_{x_{\perp},y_{\perp}}{\delta\over\delta{\tilde{\rho}}^{a}(x_{\perp})}\chi^{ab}(x_{\perp},y_{\perp})[{\tilde{\rho}}]{\delta\over{\delta{\tilde{\rho}}^{b}(y_{\perp})}}\,,$
(17)

describes the evolution of the gauge invariant weight functional $W$ with
rapidity $Y=\log(\Lambda_{0}^{+}/\Lambda^{+})\equiv\log(x_{0}/x)$, once the
non-perturbative initial conditions for $W$ are specified at an initial
$x_{0}$.

The Hamiltonian is computed in the CGC EFT, with
$\chi^{ab}(x_{\perp},y_{\perp})[{\tilde{\rho}}]=\langle\delta{\tilde{\rho}}^{a}(x_{\perp})\delta{\tilde{\rho}}^{b}(y_{\perp})\rangle_{{\tilde{\rho}}}$
the two-point function of induced charge densities171717Note that here and
henceforth in this section, $\int_{x_{\perp}}=\int d^{2}x_{\perp}$ and
$\int_{x_{\perp},y_{\perp}}=\int d^{2}x_{\perp}d^{2}y_{\perp}$. in the
classical background field of the hadron. Note that with this computation of
${\mathcal{H}}_{\rm LLx}$, the solution of Eq. (16) resums leading logarithms
$\alpha_{S}\log(x_{0}/x)$ (LLx) to all orders in perturbative theory. Thus
this powerful RG procedure extends the accuracy of computations of the cross-
section from
$\langle\mathrm{d}\sigma_{\text{LO+NLO}}\rangle\rightarrow\langle\mathrm{d}\sigma_{\text{LO+LLx}}\rangle$.

The JIMWLK RG equation can equivalently be expressed as a hierarchy of
equations (the Balitsky-JIMWLK hierarchy independently derived in Balitsky
(1996)) for the expectation value of an operator $O$:

${\partial\langle O\rangle_{Y}\over dY}=\Bigg{<}{1\over
2}\int_{x_{\perp},y_{\perp}}\\!{\delta\over\delta\alpha^{a}(x_{\perp})}\chi^{ab}(x_{\perp},y_{\perp}){\delta\over\delta\alpha^{b}(y_{\perp})}O[\alpha]\Bigg{>}_{Y},$
(18)

where $\alpha^{a}=\frac{1}{\nabla_{\perp}^{2}}{\tilde{\rho}}^{a}$. Remarkably,
Eq. (18) has the form of a generalized Fokker-Planck equation in functional
space, where $Y$ is “time” and $\chi$ is the diffusion coefficient Weigert
(2002).

There is no known analytical solution to the JIMWLK equation; as we shall
discuss, it can be solved numerically. However good approximations exist in
different limits. In a “weak field” (and leading twist) limit $g\alpha\ll 1$,
one recovers for the number distribution (and the corresponding occupation
number $\phi$) extracted from Eq. (12), the celebrated LLx BFKL equation
Kuraev  _et al._ (1977); Balitsky and Lipatov (1978) of pQCD. Another mean
field “random phase” approximation Iancu and McLerran (2001); Weigert (2002)
allows one to evaluate the occupation number $\phi$ in the “strong field”
limit of $g\alpha\sim 1$.

The analytical approximations in the different limits can, for large
rapidities $Y$, be summarized as Iancu  _et al._ (2003)

$\displaystyle\phi={1\over\pi\gamma_{s}c{\bar{\alpha}}_{S}}\,\log\left(1+\left({Q_{s}^{2}\over
k_{\perp}^{2}}\right)^{\gamma_{s}}\right)\,,$ (19)

where $c=4.88$, ${\bar{\alpha}}_{S}=\alpha_{S}N_{c}/\pi$, and the anomalous
dimension $\gamma_{s}$ is varied from $0.63$ in the all twist saturation
regime to $1$ corresponding to the leading twist pQCD “DGLAP” double log
regime. Intermediate between the two, as we shall discuss shortly, is
“shadowed” leading twist dynamics in a so-called geometric scaling window.

The longitudinal extent of the wee gluon cloud generated by the RG evolution
has a width $x^{-}=\frac{1}{k^{+}}\sim\frac{1}{Q_{s}}$. This is much more
diffuse relative to the width $e^{-1/\alpha_{S}}\frac{1}{Q_{s}}$ of valence
modes. The RG evolution also predicts that the width of the wee gluon cloud
shrinks with increasing boost (or rapidity) relative to an “observer” quark-
antiquark pair, albeit at a slower rate than their larger $x$ counterparts.

###  III.4 DIS and the dipole model

In this sub-section, and the next, we will concretely relate the CGC EFT to
the structure functions that are measured in DIS. These comparisons are
essential for precision tests of the CGC EFT picture of high energy nuclear
wavefunctions. They also play an important role in constraining the saturation
scale and the shadowing of nuclear distributions that are key to determining
the initial conditions for early time dynamics in heavy-ion collisions. These
connections will become more evident in Section IV C.

Specifically, we will now show how one starts with the inclusive DIS inclusive
cross-section we discussed previously and systematically derive the QCD
Glauber model that describes the multiple scattering of the quark-antiquark
DIS probe off gluons in the target nucleus.

The inclusive cross-section can be expressed in full generality as
$\langle\mathrm{d}\sigma\rangle=L_{\mu\nu}W^{\mu\nu}$ where $L_{\mu\nu}$ is
the well-known lepton tensor Peskin and Schroeder (1995) representing the
squared amplitude for the emission of a virtual photon with four-momentum
$q^{\mu}$ and $W^{\mu\nu}$ is the spin-averaged DIS hadron tensor,

$W^{\mu\nu}={\rm Im}\frac{i}{2\pi}\int d^{4}x\,\,e^{iq\cdot x}\langle
P|T\left(j^{\mu}(x)j^{\nu}(0)\right)|P\rangle\,,$ (20)

where the r.h.s contains the expectation value of the time-ordered product (at
two spacetime points) of electromagnetic currents
$j^{\mu}={\bar{\psi}}\gamma^{\mu}\psi$ in the ground state of the proton.

For a nucleus in the IMF, this can be reexpressed as McLerran and Venugopalan
(1999); Venugopalan (1999)

$\displaystyle W^{\mu\nu}=\frac{1}{2\pi}\frac{P^{+}}{m_{N}}{\rm Im}\int
d^{2}X_{\perp}dX^{-}\int d^{4}x\,e^{iq\cdot x}\times$
$\displaystyle\langle{\rm
Tr}\big{(}\gamma^{\mu}S_{A}(X+\frac{x}{2},X-\frac{x}{2})\gamma^{\nu}S_{A}(X^{-}\frac{x}{2},X^{+}\frac{x}{2})\big{)}\rangle\,,$
(21)

where $S_{A}(x,y)=-i\langle\psi(x){\bar{\psi}}(y)\rangle_{A}$ is the quark
propagator in the gauge fields $A^{\mu}$ of the nucleus181818We emphasize that
the second average in Eq. (21) corresponds to the averaging over the static
color sources ${\tilde{\rho}}$. Note further that in obtaining this result we
have employed a relativistic normalization of the nuclear wavefunction
$\langle
P|P\rangle=\frac{P^{+}}{m_{N}}(2\pi)^{3}\delta^{3}(0)\equiv\frac{P^{+}}{m_{N}}\int
d^{2}X_{\perp}dX^{-}$..

In the CGC, the leading contribution is obtained by replacing the full QCD
background field by the saturated classical background field:
$A^{\mu}\rightarrow A_{\rm cl}^{\mu}$, where $A_{\rm cl}^{\mu}$ are the non-
Abelian Weizsäcker-Williams (WW) fields in Eq. (11). In $A^{-}=0$
gauge191919The solution of the YM equations is identical in this case to the
solution in Lorenz gauge., the momentum space quark propagator in the
classical background field is remarkably simple, given by McLerran and
Venugopalan (1999)

$S_{\rm A_{\rm cl}}(p,q)=S_{0}(p){\cal T}_{q}(p,q)S_{0}(q)\,,$ (22)

where the free Dirac propagator is $S_{0}=\frac{i\not{p}}{p^{2}+i\varepsilon}$
and
${\mathcal{T}}_{q}(q,p)=\pm(2\pi)\delta(p^{-}-q^{-})\gamma^{-}\int_{z_{\perp}}e^{-i(\mathbf{q}_{\perp}-\mathbf{p}_{\perp})\cdot\mathbf{z}_{\perp}}V^{\pm
1}(\mathbf{z}_{\perp})$ is the effective vertex corresponding to the multiple
scattering of the quark (or antiquark) off the shock wave background field.
The latter is represented by the eikonal path ordered phase $V=P\exp\left(\int
dz^{-}{1\over\nabla_{\perp}^{2}}{\tilde{\rho}}(z^{-},x_{\perp})\right)$ we
introduced previously after Eq. (11).

If we plug Eq. (22) into Eq. (21), with a little work detailed in McLerran and
Venugopalan (1999), we find, to this order of accuracy202020The DIS structure
function $F_{2}(x_{\rm Bj},Q^{2})={\tilde{\Pi}}_{\mu\nu}W^{\mu\nu}$, where the
projector ${\tilde{\Pi}}_{\mu\nu}=\frac{3P\cdot
q}{2a}\left[\frac{P^{\mu}P^{\nu}}{a}-\frac{g^{\mu\nu}}{3}\right]$, with
$a=\frac{P\cdot q}{2x_{\rm Bj}}+m_{N}^{2}$. The Bjorken variable $x_{\rm
Bj}=Q^{2}/2P\cdot q$ and $Q^{2}=-q^{2}>0$. As previously, we will replace
$x_{\rm Bj}$ by $x$ henceforth., that McLerran and Venugopalan (1999);
Venugopalan (1999),

$F_{2}(x,Q^{2})=\frac{Q^{2}}{4\pi^{2}\alpha_{\rm
em}}\int_{0}^{1}\\!\\!\\!dz\\!\int_{r_{\perp}}\\!\\!\\!|\Psi_{\gamma^{*}\rightarrow
q{\bar{q}}}|^{2}\sigma_{q{\bar{q}}A}(x,Q^{2})\,.$ (23)

This expression for the structure function has the nice interpretation of the
convolution of the probability of the virtual photon to split into a quark-
antiquark pair (which can be computed in QED) with the “dipole” scattering
cross-section of the quark-antiquark pair to scatter off the nucleus (that
must be computed in QCD).

Specifically, in this equation $|\Psi_{\gamma^{*}\rightarrow
q{\bar{q}}}|^{2}=|\Psi_{\gamma^{*}\rightarrow
q{\bar{q}}}^{T}|^{2}+|\Psi_{\gamma^{*}\rightarrow q{\bar{q}}}^{L}|^{2}$, is
the sum of the probabilities for transversely and longitudinally polarized
virtual photons to split into quark-antiquark pairs (with longitudinal
momentum fractions $z$ and $1-z$ respectively of the virtual photon) and are
given by Bjorken  _et al._ (1971),

$\displaystyle|\Psi_{\gamma^{*}\rightarrow q{\bar{q}}}^{T}|^{2}$
$\displaystyle=$ $\displaystyle\frac{3\alpha_{\rm
em}}{2\pi^{2}}\sum_{f}e_{f}^{2}\Big{[}(z^{2}+(1-z)^{2})\epsilon_{f}^{2}K_{1}^{2}(\epsilon_{f}r_{\perp})$
$\displaystyle+$ $\displaystyle
m_{f}^{2}K_{0}^{2}(\epsilon_{f}r_{\perp})\Big{]}\,,$
$\displaystyle|\Psi_{\gamma^{*}\rightarrow q{\bar{q}}}^{L}|^{2}$
$\displaystyle=$ $\displaystyle\frac{3\alpha_{\rm
em}}{2\pi^{2}}\sum_{f}e_{f}^{2}\,4\,Q^{2}z^{2}(1-z)^{2}K_{0}^{2}(\epsilon_{f}r_{\perp})\,,$
(24)

where the sum runs over quark flavors $f$ with mass $m_{f}$, $K_{0,1}$ are the
modified Bessel functions, and $\epsilon_{f}^{2}=z(1-z)Q^{2}+m_{f}^{2}$.

The dipole cross-section $\sigma_{q{\bar{q}}A}(x,Q^{2})$, for impact parameter
$b_{\perp}=(x_{\perp}+y_{\perp})/2$, is given by

$\sigma_{q{\bar{q}}A}=2\int d^{2}b_{\perp}{\cal
N}_{Y}(b_{\perp},r_{\perp})\,,$ (25)

where the forward scattering amplitude ${\cal
N}_{Y}(b_{\perp},r_{\perp})=1-{\cal S}_{Y}(b_{\perp},r_{\perp})$, with the
S-matrix

${\cal S}_{Y}(r_{\perp})=\frac{1}{N_{c}}\langle{\rm
Tr}\left(V(x_{\perp})V^{\dagger}(y_{\perp})\right)\rangle_{Y}\,.$ (26)

Since the weight functional in $\langle\cdots\rangle$ is Gaussian (see Eq.
(9)) in the MV model, one can compute the S-matrix in Eq. (26) explicitly. One
obtains McLerran and Venugopalan (1999); Kovchegov (1999); Venugopalan (1999);
Iancu and Venugopalan (2003),

${\cal
S}_{Y}(r_{\perp})=\exp\left[-\alpha_{S}\frac{\pi^{2}}{2N_{c}}\frac{r_{\perp}^{2}AxG_{N}(x,1/r_{\perp}^{2})}{\pi
R_{A}^{2}}\right]\,,$ (27)

where $G_{N}$ denotes the gluon distribution212121In the MV model, the
expression in the exponential is
$\frac{\alpha_{S}C_{F}}{4}r_{\perp}^{2}\mu_{A}^{2}\log(\frac{1}{r_{\perp}^{2}\Lambda_{\rm
QCD}^{2}})$. This is reexpressed in terms of the gluon distribution which is
self-consistently computed in this model from the Fourier transform of Eq.
(12). One further promotes this leading log definition of the gluon
distribution by employing the DGLAP RG equation to obtain the $x$ dependence
in the gluon distribution. in the proton at the scale
$\frac{1}{r_{\perp}^{2}}$.

For very small values of $r_{\perp}$, one can expand out the exponential and
one obtains the leading order result for the dipole cross-section in pQCD. As
$r_{\perp}$ grows, the S-matrix decreases; the saturation scale is defined as
the value of $r_{\perp}$ at which the S-matrix has a value that is
significantly smaller than what one anticipates in pQCD. One choice in the
literature is $S=e^{-\frac{1}{4}}$ Kowalski  _et al._ (2008); while this
implies that there is some freedom in setting the value of the saturation
scale, its growth with decreasing $x$ is determined by the growth in the gluon
distribution.

The MV result in Eq. (27) is the QCD Glauber model Mueller (1990) which gives
the survival probability of a dipole after multiple independent scatterings
off the nucleus. It can be refined by introducing an impact parameter
distribution inside the proton Bartels  _et al._ (2002), the so-called IP-Sat
model, which can be further extended to model the S-matrix for nuclei Kowalski
and Teaney (2003); Kowalski  _et al._ (2008).

The IP-Sat model provides very good agreement with a wide range of small $x$
DIS data on e+p scattering at HERA Rezaeian  _et al._ (2013). The latter
constrains the parameters of this model, which in turn is an essential
ingredient of the IP-Glasma model of the initial conditions for heavy-ion
collisions. We will discuss the IP-Glasma model in Sec. IV.3.3.

An advantage of the MV model formulation is that one can compute with relative
ease Fujii (2002); Blaizot  _et al._ (2004a); Dominguez  _et al._ (2011);
Dusling  _et al._ (2018a); Fukushima and Hidaka (2017) not just the dipole
Wilson line correlator but quadrupole and higher point correlators that appear
in semi-inclusive final states in e+A and p+A collisions.

###  III.5 RG evolution and geometric scaling

The MV model of Gaussian random distributions is valid for a large nucleus at
rapidities when bremsstrahlung of soft gluons is not significant, namely, for
$\alpha_{s}Y\leq 1$. The classical expressions we derived have no $x$
dependence. For moderate $x$, one can introduce $x$ dependence in framework
along the lines of the IP-Sat model we discussed. However when $\alpha_{s}Y\gg
1$, the weight functional is qualitatively modified on account of significant
gluon radiation.

In this regime, RG evolution of the S-matrix in Eq. (26) is described by the
Balitsky-JIMWLK hierarchy in Eq. (18). In addition to the coherent multiple
scattering effects in the MV model, this framework captures the the real
bremsstrahlung and virtual quantum corrections we discussed previously and
sketched in Fig. 3.

Substituting the expectation value of the correlator of Wilson lines in Eq.
(26) into the Balitsky-JIMWLK hierarchy in Eq. (18), leads to the closed
form222222In general, the evolution of the two-point correlator in Eq. (26)
also depends on the correlator $\langle{\rm Tr}(V_{x}^{\dagger}V_{z}){\rm
Tr}(V_{z}^{\dagger}V_{y})\rangle$. This factorizes into $\langle{\rm
Tr}(V_{x}^{\dagger}V_{z})\rangle\langle{\rm Tr}(V_{z}^{\dagger}V_{y})\rangle$
for a large nucleus and for $N_{c}\gg 1$. The numerical simulations of these
hierarchy that we will soon discuss show these corrections to be of the order
of $10\%$. Balitsky-Kovchegov (BK) Balitsky (1996); Kovchegov and Mueller
(1998) equation for the RG evolution in rapidity of the dipole scattering
amplitude:

$\displaystyle{\partial{\cal N}_{Y}({x_{\perp},y_{\perp}})\over\partial
Y}={\bar{\alpha}_{s}}\int_{z_{\perp}}{(x_{\perp}-y_{\perp})^{2}\over(x_{\perp}-z_{\perp})^{2}(z_{\perp}-y_{\perp})^{2}}$
(28) $\displaystyle\times$ $\displaystyle\Big{[}{\cal
N}_{Y}(x_{\perp},z_{\perp})+{\cal N}_{Y}(y_{\perp},z_{\perp})-{\cal
N}_{Y}(x_{\perp},y_{\perp})$ $\displaystyle-{\cal
N}_{Y}(x_{\perp},z_{\perp}){\cal N}_{Y}(z_{\perp},y_{\perp})\Big{]}\,.$

The BK equation is the paradigmatic equation for gluon saturation in high
energy QCD. For ${\cal N}_{Y}\ll 1$, the non-linear term in the last line
above can be ignored and the equation reduces to the linear BFKL equation as
anticipated previously. In this limit, the amplitude has the solution,

$\displaystyle{\cal
N}_{Y}(r_{\perp})\approx\exp\left(\omega{\bar{\alpha}}_{s}Y-{\rho\over
2}-{\rho^{2}\over 2\beta{\bar{\alpha}}_{S}Y}\right)\,,$ (29)

where $\omega=4\log 2\approx 2.77$, $\beta=28\,\zeta(3)\approx 33.67$ and
$\rho=\log(1/r_{\perp}^{2}\Lambda_{\rm QCD}^{2})$. This solution gives the
rapid “Markovian” growth of the dipole cross-section in rapidity due to the
copious production of softer and softer gluons.

However when ${\cal N}_{Y}\sim 1$, the non-linear term arising from the fusion
and screening of soft gluons completely saturates the growth of the dipole
cross-section. If we impose a saturation condition ${\cal N}_{Y}=1/2$, for
$r_{\perp}=2/Q_{S}$, on Eq. (29), the argument of the exponential vanishes for
$\rho_{s}=\log(Q_{S}^{2}/\Lambda_{\rm QCD}^{2})$, with

$Q_{s}^{2}=\Lambda_{\rm QCD}^{2}\,\,e^{c{\bar{\alpha}}_{S}Y}\,\,\,\,{\rm
where}\,\,\,\,c=4.88\,.$ (30)

If we now write $\rho=\rho_{S}+\delta\rho$, where
$\delta\rho=\log(1/r_{\perp}^{2}Q_{S}^{2})$, one finds that Iancu  _et al._
(2002)

${\cal N}_{Y}\approx\left(r_{\perp}^{2}Q_{S}^{2}\right)^{\gamma_{s}}\,,$ (31)

for $Q^{2}<Q_{S}^{4}/\Lambda_{\rm QCD}^{2}$. Here $\gamma_{s}$ is the
anomalous dimension we introduced in Eq. (19).

This “geometrical scaling” of the forward scattering amplitude means that Eq.
(25) scales with $Q^{2}/Q_{S}^{2}(x)$ alone instead of $x$ and $Q^{2}$
separately. Remarkably, this phenomenon was observed at HERA, providing a
strong hint for the saturation picture Stasto  _et al._ (2001). Further, the
wider scaling window $Q^{2}<Q_{S}^{4}/\Lambda_{\rm QCD}^{2}$ stretching beyond
$Q_{S}$ provides a first principles explanation for a so-called “leading twist
shadowing” of nuclear parton distributions relative to those in the proton
Frankfurt  _et al._ (2012).

Figure 4: Evolution of the unintegrated gluon distribution as a function of
the squared transverse momentum from solution of the Balitsky-Kovchegov
equation. The different curves represent increasing rapidities (left to right)
for fixed and running coupling. Figure from Dusling  _et al._ (2010).

The BK equation, in a reaction-diffusion approximation, can be formally mapped
into a well-known equation in statistical physics, the Fischer-Kolmogorov-
Petrovsky-Piscounov (FKPP) equation describing such processes Munier and
Peschanski (2003). In this context, geometrical scaling appears as a late-time
solution of a non-linear equation describing a traveling wavefront of constant
velocity. In Fig. 4, we show numerical results for the unintegrated gluon
distribution $\phi(k_{\perp}^{2})=\frac{\pi
N_{c}k_{\perp}^{2}}{2\alpha_{s}}\int_{0}^{+\infty}d^{2}r_{\perp}e^{ik_{\perp}\cdot
r_{\perp}}\left[1-{\cal N}_{Y}(r_{\perp})\right]^{2}$, which displays this
traveling wave front structure.

The correspondence of high energy QCD to reaction-diffusion processes is very
rich; advances in the dynamics of the latter can provide deeper insight into
the stochastic dynamics of the former. Specific applications to DIS have been
discussed recently Mueller and Munier (2018a, b).

It is important to note that $Q_{S}^{2}$ in Eq. (30) (and the amplitude in Eq.
(31)) grows very rapidly with rapidity, much faster than seen in the HERA
data. However this is significantly modified by running coupling corrections,
which are part of the next-to-leading-logs in $x$ (NLLx) contributions to QCD
evolution. The significant effect of these running coupling corrections is
clearly seen in Fig. 4.

These give Mueller and Triantafyllopoulos (2002),

$Q_{s,{\rm running}\,\alpha_{s}}^{2}=\Lambda_{\rm
QCD}^{2}\exp\left(\sqrt{2b_{0}c(Y+Y_{0})}\right)\,,$ (32)

where $b_{0}$ is the coefficient of the logarithm in the one loop QCD
$\beta$-function 232323Sub-leading corrections in $Y$ to $Q_{S}$ have been
computed to high order and are discussed in Beuf (2010).. The running coupling
results give a power law increase of the amplitude consistent with the HERA
data. Further, the qualitative features of geometric scaling persist even
after including running coupling effects and next-to-leading-logarithms (NLLx)
corrections to the BFKL kernel are accounted for, albeit the window for
geometrical scaling is significantly smaller Triantafyllopoulos (2008).

For a large nucleus at the saturation boundary
$Y_{0}\propto\log^{2}(A^{1/3})$, one recovers the $A^{1/3}$ scaling of the
saturation scale in the MV model from Eq. (32) for $Y_{0}\gg Y$. A striking
result, for $Y\gg Y_{0}$, is that the saturation scale for fixed impact
parameter becomes independent of $A$. Strongly correlated gluons in the
asymptotic Regge limit, lose memory of the initial conditions whereby they
were generated.

###  III.6 The state of the art in the CGC EFT

In previous sub-sections, we outlined a description of the wavefunction of a
high energy nucleus in the CGC EFT, emphasizing a qualitative understanding of
gluon saturation and key related analytical results. There have been
significant developments since in the CGC EFT.

On the formal side, the Balitsky-JIMWLK framework for the LLx evolution of
n-point Wilson line correlators, has been extended to NLLx Balitsky and
Chirilli (2013a); Kovner  _et al._ (2014a, b); Balitsky and Grabovsky (2015);
Caron-Huot (2018). For the 2-point dipole correlator, which satisfies the LLx
BK equation, the formalism has been extended to NLLx Balitsky and Chirilli
(2008) and even (for ${\cal N}=4$ supersymmetric Yang-Mills) to NNLLx in a
recent tour de force computation Caron-Huot and Herranen (2018). The BFKL/BK
kernel however receives large collinear contributions that need to be resummed
in so-called small $x$ resummation schemes for quantitative predictions Salam
(1998); Ciafaloni  _et al._ (1999); Iancu  _et al._ (2015); Ducloué  _et al._
(2019).

While as we have discussed, there are good analytical approximations, a full
analytical solution of the BK equation does not exist. Numerical simulations
have however been known for some time for the LLx BK equation Albacete  _et
al._ (2004), the LLx+running coupling BK equation Albacete  _et al._ (2005);
Albacete and Kovchegov (2007), and even more recently the full NLLx equation
implementing collinear resummation Lappi and Mäntysaari (2016); Ducloué  _et
al._ (2020). In particular, it is shown in Ducloué  _et al._ (2020) that this
NLLx framework provides very good agreement with the HERA data.

Figure 5: Solution of JIMWLK equation for the correlator of Wilson lines
$V(x_{\perp})V^{\dagger}(y_{\perp})$ probed by the DIS dipole Dumitru  _et
al._ (2011a). As the nucleus is boosted from low energy (or rapidity) to high
energy, the regions with large values of these correlator shrink spatially,
corresponding to larger values of $Q_{S}$.

Numerical simulations have also been performed of higher point correlators in
the Balitsky-JIMWLK hierarchy. As noted, Eq. (18) has the form of a functional
Fokker-Planck equation. This can therefore be reexpressed as a Langevin
equation in the space of Wilson lines Weigert (2002); Blaizot  _et al._
(2003a), allowing one to simulate the rapidity evolution of two-point Wilson
line correlators Rummukainen and Weigert (2004) as well as 4-point quadrupole
and sextupole242424These are probed in semi-inclusive DIS Dominguez  _et al._
(2011) and in proton-nucleus collisions Kovner and Lublinsky (2011); Dusling
_et al._ (2018a, b). correlators Dumitru  _et al._ (2011a); Lappi and Ramnath
(2019); Lappi and Mäntysaari (2013). Fig. 5 shows a result for the dipole
correlator from these simulations. Unfortunately, a similar Langevin
representation is not known at present for the NLLx JIMWLK Hamiltonian.

Precision computations require not just higher order computations of the
JIMWLK kernel but higher order computations of so-called process dependent
“impact factors” as well. This is analogous to pQCD computations of
coefficient functions that are convoluted, order-by-order, with the DGLAP
splitting functions Vermaseren  _et al._ (2005). For inclusive DIS, analytical
expressions exist for the virtual photon impact factor
$|\Psi_{\gamma^{*}\rightarrow q{\bar{q}}}|^{2}$ in Eq. (23) Balitsky and
Chirilli (2013b). More recently, NLO impact factors have been computed for DIS
exclusive diffractive light vector meson production Boussarie  _et al._ (2017)
and DIS inclusive photon+dijet production Roy and Venugopalan (2019, 2020).
Numerical implementation of these results to compute cross-sections remains a
formidable task and an essential component of precision studies of gluon
saturation at a future Electron-Ion Collider (EIC). A summary of extant EIC
studies can be found in Accardi  _et al._ (2016); Aschenauer  _et al._ (2019).

An outstanding problem at small $x$ is the impact parameter dependence of
distributions. The BFKL kernel at large impact parameters contributes a
Coulomb tail $\sim 1/b_{\perp}^{2}$; the conformal symmetry of the kernel and
geometric scaling suggest a particular dependence of the saturation scale on
the impact parameter Gubser (2011). The Coulomb tail is however not regulated
by saturation and violates the Froissart bound on the asymptotic behavior of
total cross-sections Kovner and Wiedemann (2003). This is only cured non-
perturbatively by the generation of a mass gap in QCD. How to address this
problem in the context of scattering, and quantifying the importance of the
Coulomb tail for different processes, remain as open problems despite much
work. This may be less of a problem in large nuclei with $\Lambda_{\rm
QCD}R_{A}\gg 1$ (until asymptotically high energies) because the contribution
of the Coulomb tail may be suppressed relative to protons, for which
$\Lambda_{\rm QCD}R_{A}\sim 1$.

##  IV Non-equilibrium QCD matter at high occupancy

The CGC EFT provides us with a powerful tool to address multi-particle
production in heavy-ion collisions from first principles. As noted, a key
organizing principle in developing an understanding of hadron structure at
high energies, is the kinematic separation in the hadron wavefunction between
static color sources at large $x$ and small $x$ gauge fields. In the CGC EFT,
the latter are pure gauge coherent classical fields (of ${\mathcal{O}}(1/g)$)
and are static in $x^{+}$ for nuclei with $P^{+}\rightarrow\infty$; after the
collision, they acquire a time dependence that leads to multi-particle
production.

The corresponding ab initio problem in quantum field theory then is the
computation of multi-particle production in the presence of strong fields. One
well-known example is $e^{+}e^{-}$ pair production in strong electromagnetic
fields Gelis and Tanji (2016); another is Hawking radiation from the Black
Hole horizon Parikh and Wilczek (2000).

In the following, we will sketch the elements of a formalism to compute
inclusive multi-particle production in the collision of two CGCs. The non-
equilibrium matter formed at very early times in the collision is the Glasma
Lappi and McLerran (2006); Gelis and Venugopalan (2006a), a state with high
occupancy $f\sim{\mathcal{O}}(1/\alpha_{S})$. Unlike the CGC, this state
decays and eventually thermalizes. Besides its relevance in the thermalization
process, the Glasma is a strongly correlated state of matter with distinct and
universal properties. These features of the Glasma will be discussed at length
in Section V.

###  IV.1 Multi-particle production in strong fields

To compute multi-particle production systematically in the collision of the
CGC gluon “shock waves”, we will begin with the first principles Lehmann-
Symanzik-Zimmerman (LSZ) formalism in quantum field theory (QFT). In the LSZ
formalism252525For simplicity, we consider here a self-interacting $\phi^{3}$
scalar theory; our discussion extends straightforwardly to the Yang-Mills
case., the amplitude for $n$-particles in the “out” state generated from the
“in-vacuum” can be expressed as

$\displaystyle\langle$ $\displaystyle\\!p_{1,{\rm out}}\cdots p_{n,{\rm
out}}|0_{\rm
in}\rangle=\frac{1}{Z^{n/2}}\int\Bigg{[}\prod_{i=1}^{n}d^{4}x_{i}e^{ip_{i}\cdot
x_{i}}$ (33) $\displaystyle\times$
$\displaystyle\left(\partial_{x_{i}}^{2}+m^{2}\right){\delta\over\delta
J(x_{i})}\Bigg{]}\exp\left(i{\cal V}\right)\,.$

Here $p_{1},\cdots p_{n}$ denote the momenta of the produced particles and the
“in-out” vacuum-amplitude $\langle 0_{\rm out}|0_{\rm in}\rangle=\exp(i{\cal
V})$, where ${\cal V}$ is the sum of all connected vacuum-vacuum diagrams
coupled to external sources. An illustration of multi-particle production for
the problem at hand is shown in Fig. 6.

In QFT computations, one usually sets $J=0$ after the functional
differentiation and $\langle 0_{\rm out}|0_{\rm in}\rangle$ is a pure phase.
When $J$ is physical, $|\langle 0_{\rm out}|0_{\rm
in}\rangle|^{2}=\exp(-2\,{\rm Im\,}{\cal V})\neq 1$. In computing multi-
particle production in this context, it is useful to employ262626For other
discussions of the SK formalism in the context of the CGC and the Glasma, see
Jeon (2014); Wu and Kovchegov (2018); Leonidov and Radovskaya (2019). For a
recent discussion in the context of thermal field theory, see Ghiglieri  _et
al._ (2020). the Schwinger-Keldysh (SK) QFT formalism Schwinger (1961);
Keldysh (1964). One introduces $+$ and $-$ vertices with opposite signs of the
coupling in Feynman diagrams, and likewise for the sources $J_{\pm}$. The
corresponding “$+$” and “$-$” fields live on the upper and lower segments of a
closed time contour ranging forward in time from $t=-\infty$ on the upper
contour and back to $-\infty$ on the lower contour, as shown in Fig. 6. Time
ordered “$++$” (anti-time ordered “$--$”) Green’s functions “live” on the
upper (lower) contour, and the mixed $+-$ “Wightman” functions connect the
upper and lower contours.

Figure 6: a) Multi-particle production from cut “vacuum-vacuum” graphs
connecting time dependent sources of the two nuclei after the collision.
Figure from Gelis  _et al._ (2010). b) The Schwinger-Keldysh closed time
contour on which the sources and fields are defined.

Following LSZ, the probability to produce $n$-identical particles is

$P_{n}=\frac{1}{n!}\prod_{i=1}^{n}{d^{3}p_{i}\over(2\pi)^{3}2E_{p_{i}}}|\langle
p_{1,{\rm out}}\cdots p_{n,{\rm out}}|0_{\rm in}\rangle|^{2}\,,$ (34)

where $E_{p_{i}}^{2}=p_{i}^{2}+m^{2}$. Plugging the expression for the
amplitude in Eq. (33) into the r.h.s, one can express the result as Gelis and
Venugopalan (2006b)

$\displaystyle P_{n}=\frac{1}{n!}{\cal D}^{n}\exp\left(i{\cal V}[J_{+}]-i{\cal
V}[J_{-}]\right)|_{J_{+}=J_{-}=J}\,,$ (35)

with

${\cal
D}=\int_{x,y}\\!\\!\\!\\!\\!\\!Z\,G_{+-}^{0}(x,y)\frac{\left(\partial_{x_{i}}^{2}+m^{2}\right)}{Z}\frac{\left(\partial_{y_{i}}^{2}+m^{2}\right)}{Z}{\delta\over\delta
J_{+}(x)}{\delta\over\delta J_{-}(y)}\,.$ (36)

Here $\int_{x}=d^{4}x$,
$G_{+-}^{0}(x,y)=\int{d^{3}p_{i}\over(2\pi)^{3}2E_{p_{i}}}e^{ip\cdot(x-y)}\equiv\theta(p^{0})\delta^{(3)}(x-y)$
and $Z$ is the residue of the pole of the renormalized propagator.

The action of the operator ${\cal D}$ can be understood as follows. The “$+$”
piece with
$\frac{\left(\partial_{x_{i}}^{2}+m^{2}\right)}{Z}{\delta\over\delta
J_{+}(x)}$ acts on a particular diagram in the connected sum of vacuum-vacuum
connected diagrams ${\cal V}[J_{+}]$ by removing a source $J_{+}$ and then
amputating the renormalized propagator to which it is attached. The same
procedure is followed for the “$-$” piece; the two amputated propagators are
then sewn together by the renormalized “cut” propagator $ZG_{+-}^{0}$.

Computing $P_{n}$ in a theory with physical sources is hard272727One might
imagine it sufficient to compute ${\rm Im\,}{\cal V}$ to $n={\rm few}$ since
particle multiplicity is naively suppressed by $\alpha_{S}$. However since the
effective expansion parameters in the Glasma are $\alpha_{S}\log(1/x)$ and
$\alpha_{S}^{2}A^{1/3}$, both of ${\mathcal{O}}(1)$ in heavy-ion collisions,
$n\gg 1$ is typical. because one also has to compute the disconnected vacuum-
vacuum graphs for each $n$. However if we define a generating functional
$F(z)=\sum_{n}z^{n}P_{n}$, Eq. (35) gives

$F(z)=\exp\left(z{\cal D}\right)\exp\left(i{\cal V}[J_{+}]-i{\cal
V}[J_{-}]\right)|_{J_{+}=J_{-}=J}\,,$ (37)

and successive differentiation of this equation with respect to $z$ (and
setting $z=1$), generates the $n$-particle correlators $\langle
n(n-1)(n-2)\cdots\rangle$. These moments do not require one compute the
disconnected vacuum-vacuum graphs, since they also appear in the normalization
of $P_{n}$ and therefore cancel out282828Such cancellations are seen in the
Abramovsky-Gribov-Kancheli (AGK) rules Abramovsky  _et al._ (1973) that
implement the combinatorics of cut/uncut vacuum-to-vacuum graphs in Reggeon
field theory Gelis and Venugopalan (2007). in the moments.

This is illustrated by expressing the r.h.s of Eq. (37) for $z=1$ as

$\exp\left(i{\cal V_{\rm SK}}[J_{+},J_{-}]\right)=\exp\left({\cal
D}\right)\exp\left(i{\cal V}[J_{+}]-i{\cal V}[J_{-}]\right)\,,$ (38)

where now $i{\cal V_{\rm SK}}[J_{+},J_{-}]$ represents the sum over all
vacuum–to–vacuum connected graphs that live on the SK closed time contour. One
can then express the inclusive multiplicity as Gelis and Venugopalan (2006b)

$\langle
N\rangle=\int_{x,y}ZG_{+-}^{0}(x,y)\left[\Gamma_{+}(x)\Gamma_{-}(y)+\Gamma_{+-}(x,y)\right]_{J_{\pm}=J}\,,$
(39)

with the amputated one-point and two-point Green’s functions in the Schwinger-
Keldysh formalism defined respectively as

$\Gamma_{\pm}(x)=\Delta_{x}^{\rm R}\frac{\delta i{\cal V}_{\rm SK}}{\delta
J_{\pm}(x)}\,;\,\Gamma_{+-}(x,y)=\Delta_{x}^{\rm R}\Delta_{y}^{\rm
R}\frac{\delta^{2}i{\cal V}_{\rm SK}}{\delta^{2}J_{+}(x)J_{-}(y)}\,,$ (40)

with $\Delta_{x}^{\rm R}=\frac{\partial_{x}^{2}+m^{2}}{Z}$.

In summing over all the nodes of all the trees connecting $\Gamma_{+}(x)$ to
the sources, the time (anti-time) ordered Feynman propagators in each tree on
the upper (lower) SK contour are recursively converted to retarded
propagators: $G_{R}=G_{++}-G_{+-}\equiv G_{-+}-G{--}$. This is equivalent to
solving the classical equations of motion with retarded boundary conditions
when $J_{\pm}=J$ ! A further important result is that the renormalized cut
propagator $\Gamma_{+-}$ is obtained by solving the small fluctuation
equations of motion in the classical background, also as an initial value
problem with retarded boundary conditions.

As we discussed previously, the classical fields, and sources thereof, of the
colliding CGC’s are static shock waves; as such, they do not spontaneously
decay and are thus part of the nuclear wavefunction. After the collision, the
colored sources become time dependent. Thus $\Gamma_{\pm}$ in Eq. (39)
corresponds to $\partial_{x}^{2}{\cal A}_{\pm,{\rm cl.}}^{\mu}$ where ${\cal
A}_{\pm,{\rm cl.}}^{\mu}$ is the time dependent ${\mathcal{O}}(1/g)$ Glasma
field in the forward lightcone. The two-point function $\Gamma_{+-}(x,y)$ in
Eq. (39) is ${\mathcal{O}}(1)$ and therefore NLO in the power counting for the
inclusive multiplicity in the Glasma. The formalism can be extended to higher
orders in $\alpha_{S}$. Its generalization to higher multiplicity moments was
developed in Gelis and Venugopalan (2006c).

###  IV.2 The LO Glasma: classical gluon fields from shockwave collisions

Since at LO in our power counting only the product
$\Gamma_{+}(x)\Gamma_{-}(y)\equiv\partial_{x}^{2}{\cal
A}_{+}^{\mu}\partial_{x}^{2}{\cal A}_{-}^{\nu}$ in Eq. (39) contributes, one
obtains for a fixed distribution of lightcone sources
$\rho_{\pm,1,2}=\rho_{1,2}$ (where $1,2$ denote the two nuclei) Gelis  _et
al._ (2007)

${d\langle N\rangle_{\rm LO}\over
dYd^{2}p_{\perp}}[\rho_{1},\rho_{2}]=\frac{1}{16\pi^{3}}\int_{x,y}\\!\\!\\!\Delta_{x}^{\rm
R}\Delta_{y}^{\rm
R}\varepsilon_{\lambda}^{\mu}\varepsilon_{\lambda}^{\nu}{\cal A}_{\mu}(x){\cal
A}_{\nu}(y)\,,$ (41)

where repeated indices are summed over. Note too that ${\cal A}(x)\equiv{\cal
A}_{\mu}[\rho_{1},\rho_{2}](x)$ and $m=0$ in $\Delta_{x,y}^{\rm R}$. An
integration by parts,

$\int d^{4}x\,e^{ip\cdot x}\partial_{x}^{2}{\cal
A}_{\mu}(x)=\int_{x^{0}\rightarrow+\infty}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!d^{3}x\,e^{ip\cdot
x}\left(\partial_{0}-iE_{p}\right){\cal A}_{\mu}(x)\,,$ (42)

shows that Eq. (41) can be computed by solving the classical YM equations in
Eq. (10) (with
$J^{\mu}=\delta^{\mu+}\delta(x^{-})\rho_{1}(x_{\perp})+\delta^{\mu-}\delta(x^{+})\rho^{2}(x_{\perp})$
and ${\cal A}_{\mu}(x)|_{x^{0}=-\infty}=0$) to determine ${\cal A}_{\mu}(x)$.

In the discussion to follow, it will be convenient to introduce the
$(\tau,\eta,x_{\perp})$ coordinate system, where the proper time
$\tau=\sqrt{(x^{0})^{2}-(x^{3})^{2}}$ and the spacetime rapidity
$\eta=\frac{1}{2}\log(\frac{x^{0}+x^{3}}{x^{0}-x^{3}})$, and $g_{\mu\nu}={\rm
diag}(1,-\tau^{2},-1,-1)$. A convenient gauge to solve the YM equations in the
forward lightcone is the Fock-Schwinger gauge ${\cal A}^{\tau}\equiv
x^{+}{\cal A}^{-}+x^{-}{\cal A}^{+}=0$. In this gauge292929A perturbative
solution was also found in Lorenz gauge $\partial_{\mu}{\cal A}^{\mu}=0$, but
is not easily extended to discuss the full non-perturbative solution to the YM
equations Kovchegov and Rischke (1997)., the solution to the YM equations are
manifestly boost invariant: ${\cal A}^{\mu}(\tau,\eta,x_{\perp})\equiv{\cal
A}^{\mu}(\tau,x_{\perp})$ and one obtains Kovner  _et al._ (1995a, b);
Gyulassy and McLerran (1997),

${\cal A}^{i}=A_{1,{\rm cl.}}^{i}+A_{2,{\rm cl.}}^{i}\,;\,{\cal
A}^{\eta}=\frac{ig}{2}[A_{1,{\rm cl.}}^{i},A_{2,{\rm cl.}}^{i}]\,,$ (43)

with $\partial_{\tau}{\cal A}^{i}=0$ and $\partial_{\tau}{\cal A}^{\eta}=0$ at
$\tau=0^{+}$. This solution is obtained by matching the delta-functions on the
lightcone wedges shown in Fig. 7.

Since the gauge fields are functionals of $\rho_{1,2}$, the full average
inclusive multiplicity in the Glasma is obtained by averaging over many
nuclear collisions, each with its distribution of color sources in the two
nuclei303030Due to color confinement at distances scales $1/\Lambda_{\rm
QCD}$, one requires $\int_{0}^{1/\Lambda_{\rm
QCD}}d^{2}x_{\perp}\rho_{1,2}^{a}=0$ for each such configuration.. This can be
expressed as

$\displaystyle{d\langle\langle N\rangle\rangle_{\rm LO}\over
dYd^{2}p_{\perp}}$ $\displaystyle=$
$\displaystyle\int[D\rho_{1}][D\rho_{2}]W_{Y_{{\rm beam}-Y}}^{\rm
MV}[\rho_{1}]W_{Y_{{\rm beam}+Y}}^{\rm MV}[\rho_{2}]$ (44)
$\displaystyle\times$ $\displaystyle{d\langle N\rangle_{\rm LO}\over
dYd^{2}p_{\perp}}[\rho_{1},\rho_{2}]\,,$

where the $W^{\rm MV}$’s for each of the nuclei at LO are the weight
functionals in the MV model in Eq. (9) and are independent, to this order, of
$Y_{{\rm beam}\mp Y}$, where $Y_{\rm beam}=\log(\sqrt{s}/m_{N})$ is the beam
rapidity.

Figure 7: Gauge field configurations in the spacetime diagram of the
collision. Before the collision, the gauge fields are pure gauge,
corresponding to zero field strength outside the gluon shockwaves. After the
collision, the solution of the Yang-Mills generate finite field strengths in
the Glasma leading to multiparticle production. Figure from Lappi and McLerran
(2006).

With the initial conditions in Eq. (43), the YM equations for $\tau=0^{+}$ can
be solved perturbatively to lowest non-trivial order in
${\mathcal{O}}(\frac{\rho_{1}}{\nabla_{\perp}^{2}}\frac{\rho_{2}}{\nabla_{\perp}^{2}})$;
in this “dilute-dilute” approximation, one obtains for identical nuclei,

${d\langle\langle N\rangle\rangle_{\rm LO}\over dYd^{2}p_{\perp}}={\pi
R_{A}^{2}}\frac{g^{6}\mu_{A}^{4}}{(2\pi)^{4}}\frac{2N_{c}(N_{c}^{2}-1)}{p_{\perp}^{4}}{\cal
L}(p_{\perp},\Lambda)\,.$ (45)

This result, which agrees with the pQCD bremsstrahlung formula first derived
by Gunion and Bertsch Gunion and Bertsch (1982) is valid for $p_{\perp}\gg
Q_{S}$ (recall $Q_{S}\propto g^{2}\mu_{A}$) and ${\cal L}(p_{\perp},\Lambda)$
is a logarithmically divergent function, screened at
$\Lambda\approx\Lambda_{\rm QCD}$.

From our dipole model discussion (see Eq. (27) and related discussion),
$Q_{S}^{2}\propto G_{A}(x,p_{\perp}^{2})$, where $p_{\perp}$ is the momentum
conjugate to the dipole size. This suggests that Eq. (45) can be generalized
to a “$k_{\perp}$ factorization” form ${d\langle\langle N\rangle\rangle_{\rm
LO}\over dYd^{2}p_{\perp}}\propto\alpha_{S}\int
dk_{\perp}^{2}\phi_{A}(x_{1},k_{\perp}^{2})\phi_{B}(x_{2},(k_{\perp}-p_{\perp})^{2})$.
Here $\frac{\phi_{A,B}(x,k_{\perp}^{2})}{k_{\perp}^{2}}$ is the Fourier
transform of the dipole scattering amplitude313131This distribution is
distinct from the WW-distribution and coincides with it only for large
$k_{\perp}$ Kharzeev  _et al._ (2003); Blaizot  _et al._ (2004b). in the each
of the hadrons we discussed previously in Sec. III.5. This $k_{\perp}$
factorization formula Gribov  _et al._ (1983) is widely used in
phenomenological studies of hadron-hadron collisions.

The dilute-dilute analytical approximation for shockwave collisions can be
generalized to compute the inclusive multiplicity to lowest order
${\mathcal{O}}(\frac{\rho_{1}}{\nabla_{\perp}^{2}})$ in one of the sources but
to all orders ${\mathcal{O}}((\frac{\rho_{2}}{\nabla_{\perp}^{2}})^{n})$ in
the other. In this “dilute-dense” case as well, the inclusive gluon
multiplicity can be expressed as a $k_{\perp}$-factorized convolution of the
unintegrated gluon distributions in the projectile and target. It is valid for
$Q_{S,1}^{2}(x_{1})\ll Q_{S,2}^{2}(x_{2})$, corresponding to the forward (or
backward) kinematic regions of the shockwave collision where the parton
momentum fractions are $x_{1}\gg x_{2}$ . Alternately, it can be a good
approximation in proton-nucleus collisions, where $Q_{S,A}^{2}\sim
A^{1/3}Q_{S,p}^{2}$ Kovchegov and Mueller (1998); Dumitru and McLerran (2002).

###  IV.3 Non-perturbative evolution of high occupancy fields

####  IV.3.1 Real time evolution of boost invariant fields on the lattice

While analytical results for the inclusive multiplicity are available only in
limited kinematic regions, the YM equations for shockwave collisions can be
solved numerically to all orders
${\mathcal{O}}((\frac{\rho_{1,2}}{\nabla_{\perp}^{2}})^{n})$ Krasnitz and
Venugopalan (1998, 1999) to obtain the full non-perturbative result to Eq.
(44) Krasnitz and Venugopalan (2000, 2001); Krasnitz  _et al._ (2001); Lappi
(2003); Krasnitz  _et al._ (2003a). Hamilton’s equations are solved in Fock-
Schwinger gauge ${\cal A}^{\tau}=0$ with the initial conditions at $\tau=0$
specified by Eq. (43). To preserve gauge invariance, lattice gauge theory
techniques can be adapted to this problem. The boost invariance of the LO
shockwave gauge fields provides a significant simplification whereby the 3+1-D
Kogut-Susskind QCD lattice Hamiltonian Kogut and Susskind (1975) can be
“dimensionally reduced” to the 2+1-D form Krasnitz and Venugopalan (1999)

$\displaystyle aH$
$\displaystyle=\sum_{\mathbf{x}}\left[\frac{g^{2}a}{\tau}{\rm
tr}\,E^{i}E^{i}+\frac{2\tau}{g^{2}a}\left(N_{c}-{\rm Re}\,{\rm
tr}\,U_{1,2}\right)\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.+\frac{\tau}{a}{\rm
tr}\,\pi^{2}+\frac{a}{\tau}\sum_{i}{\rm
tr}\left(\Phi-\tilde{\Phi}_{i}\right)^{2}\right]\,.$ (46)

Here the sum is over all discretized cells in the transverse plane and for
clarity, we have omitted the cell index $j$ for all quantities in this
expression. Further, $E^{i}$ with $i\in\\{1,2\\}$ are the components of the
transverse electric field living on each site; discretizing the initial
conditions gives $E^{i}=0$ at $\tau=0$. The spatial plaquette of link
variables $U^{i}_{j}$,

$U_{1,2}^{j}=U^{1}_{j}\,U^{2}_{j+\hat{e}_{1}}\,U^{1{\dagger}}_{j+\hat{e}_{2}}\,U^{2{\dagger}}_{j}\,,$
(47)

(where $+\hat{e_{i}}$ indicates a shift from $j$ by one lattice site in the
$i=1,2$ transverse direction) represents the squared longitudinal magnetic
fields in the Glasma. In Eq. (IV.3.1), we have represented
$A_{\eta}(\tau,x_{\perp})$ as an adjoint scalar field $\Phi$ because, as a
result of boost invariance, it transforms covariantly under $\eta$-dependent
gauge transformations:

$\tilde{\Phi}_{i}^{j}=U_{j}^{i}\Phi_{j+\hat{e}_{i}}U_{j}^{i{\dagger}}\,.$ (48)

Finally, $\pi=E_{\eta}=\dot{\Phi}/\tau$ in Eq. (IV.3.1) represents the
longitudinal electric field.

The details of the numerical simulations of the real time evolution of gauge
fields can be found in Krasnitz and Venugopalan (1999); Lappi (2003). In the
early work, only uniform sheets of nuclei were considered with constant ($x$
independent) values of $Q_{S}$. These were subsequently relaxed to consider
finite nuclei Krasnitz  _et al._ (2003b, c); more realistic simulations with
event-by-event simulations of RHIC and LHC collisions were developed later in
the IP-Glasma model we shall discuss shortly Schenke  _et al._ (2012a).

As anticipated, the numerical results reproduce the perturbative result in Eq.
(45) at large $k_{\perp}\gg Q_{S}$. However, unlike that expression, there is
no logarithmic factor ${\cal L}(k_{\perp},\Lambda_{\rm QCD})$. At momenta
$k_{\perp}<Q_{S}$, the $1/k_{\perp}^{4}$ distribution is modified to a form
that is well fit by a Bose-Einstein exponential distribution Krasnitz  _et
al._ (2003a). Even more remarkably, the non-linear dynamics generates a
plasmon mass323232This plasmon mass is parametrically larger than the
confining scale; its properties have been investigated recently in a number of
approaches Dumitru  _et al._ (2014); Boguslavski  _et al._ (2019). that
screens the momentum distribution in the infrared Krasnitz and Venugopalan
(2001); Lappi and Peuron (2018). The energy density is therefore well-defined
at all proper times without infrared or ultraviolet divergences Lappi (2006).

####  IV.3.2 Glasma flux tubes

An interesting consequence of the LO Glasma solution is that the Weizäcker-
Williams plane polarized $E$ and $B$ fields in the colliding CGCs become
purely longitudinal immediately after the collision at $\tau=0^{+}$;
$E_{\eta},B_{\eta}\neq 0$ and $E_{i},B_{i}=0$. It was pointed out in Kharzeev
_et al._ (2002) that this configuration satisfies the identity

$Q_{\rm CS}=\frac{\alpha_{S}}{2\pi}\int d^{4}x\,{\rm Tr}\,E_{\eta}\cdot
B_{\eta}\,,$ (49)

where the topological charge $Q_{\rm CS}=\frac{\alpha_{S}}{16\pi}\int
d^{3}x\,K^{0}$ and $K^{\mu}$ is the Chern-Simons current. A neat
interpretation Lappi and McLerran (2006); Chen  _et al._ (2015a) of this
result is that the YM equations at $\tau=0^{+}$ can be expressed as
$\nabla\cdot E=\rho_{\rm el.}$ and $\nabla\cdot B=\rho_{\rm mag.}$, where
$\rho_{\rm el.}$, $\rho_{\rm mag.}$ are respectively electric and magnetic
charges densities333333These induced charge densities are proportional to the
commutators $\delta^{ij}[A_{1,{\rm cl.}}^{i},A_{1,{\rm cl.}}^{j}]$ and
$\epsilon^{ij}[A_{1,{\rm cl.}}^{i},A_{1,{\rm cl.}}^{j}]$ respectively. on the
gluon shockwaves after the collision.

As sketched in Fig. 8, the induced electric and magnetic charges generate a
“stringy” Glasma flux tube Dumitru  _et al._ (2008a) of chromo-electromagnetic
fields that is uniform in rapidity stretching between the fragmentation
regions of the nuclei and are color screened Krasnitz  _et al._ (2003c) on
transverse distance scales $\geq 1/Q_{S}$.

One can straightforwardly compute the energy densities and pressures in the
Glasma from the different components of the stress-energy tensor343434Note
that
$T^{\mu\nu}=-g^{\mu\alpha}g^{\nu\beta}g^{\gamma\delta}F_{\alpha\gamma}F_{\beta\delta}+\frac{1}{4}g^{\mu\nu}g^{\alpha\gamma}g^{\beta\delta}F_{\alpha\beta}F_{\gamma\delta}.$
(50) . We obtain $\mathcal{E}=2\mathcal{P}_{T}+\mathcal{P}_{L}$ where,

$\displaystyle\mathcal{P}_{T}\equiv\frac{1}{2}\left(T^{xx}+T^{yy}\right)={\rm
Tr}\left(F_{xy}+E_{\eta}^{2}\right)\,\,$
$\displaystyle\mathcal{P}_{L}\equiv\tau^{2}T^{\eta\eta}=\frac{1}{\tau^{2}}{\rm
Tr}\left(F_{\eta i}^{2}+E_{i}^{2}\right)-{\rm
Tr}\left(F_{xy}+E_{\eta}^{2}\right)\,.$ (51)

At the earliest times after the collision $\tau=0^{+}$, as noted, only the
longitudinal $E_{\eta}$ and $B_{\eta}=F_{xy}$ fields are non-zero. The above
equation then immediately gives $\mathcal{P}_{T}=\mathcal{E}$ and
$\mathcal{P}_{L}=-\mathcal{E}$. Thus at the earliest times, the pressure in
the Glasma is purely transverse; after initial transverse dynamics, the
longitudinal pressure $\mathcal{P}_{L}\rightarrow 0$ from below by $\tau\sim
1/Q_{S}$. Since the Glasma at LO is conformal, the energy density satisfies
$\mathcal{E}=2\mathcal{P}_{T}$ at this time.

Figure 8: Glasma flux tubes: Boost invariant LO Glasma configurations of
transverse size $1/Q_{S}$ at $\tau=0^{+}$ with parallel $E_{\eta}$ and
$B_{\eta}$, corresponding to finite Chern-Simons charge. Such configurations
decay rapidly and are unstable to quantum fluctuations. Figure from Dumitru
_et al._ (2008a).

Stringy models capture essential features of confining dynamics in QCD Bali
(2001). In high energy collisions, they have a long history and capture the
bulk features of the spectrum of multi-particle production Andersson  _et al._
(1979); Artru (1983); they underlie event generators such as PYTHIA Andersson
_et al._ (1983). These models however screen color at distance scales
$1/\Lambda_{\rm QCD}$ and only carry electric flux and no magnetic flux;
particle production is assumed to arise from the Schwinger mechanism Andersson
_et al._ (1979). It is remarkable nevertheless to observe that similar stringy
solutions emerge from the more fundamental framework of classical YM
equations.

Motivated by this stringy picture, we expect the number of gluons per unit
rapidity equals the number of flux tubes ($S_{\perp}/(1/Q_{S}^{2})$) times the
gluon occupancy in a flux tube
($2(N_{c}^{2}-1)/{{\bar{\alpha}}_{S}}/(2\pi)^{3}$) multiplied by a non-
perturbative coefficient of ${\mathcal{O}}(1)$. Extracting the number density
from the correlator of gauge fields at $\tau\sim 1/Q_{S}$ Krasnitz and
Venugopalan (2001), one indeed finds that353535Here and henceforth, for
simplicity of notation, the path integral over gauge fields (moot at LO), and
over sources, $\langle\langle\rangle\rangle$ is implicit.

${dN_{\rm LO}\over
dY}=c_{N}\frac{2\,(N_{c}^{2}-1)}{(2\pi)^{3}}\,\frac{Q_{S}^{2}S_{\perp}}{{\bar{\alpha}}_{S}}\,,$
(52)

where $S_{\perp}$ is the transverse area of the collision,
${\bar{\alpha}}_{S}=\alpha_{S}N_{c}/\pi$ and $c_{N}$ is a gluon liberation
coefficient Mueller (2000) estimated from the numerical simulations to be
$c_{N}=1.1$ with $10\%$ accuracy Lappi (2008).

The YM simulations can also be extended to compute two particle correlations
in the Glasma Lappi  _et al._ (2010):

$\frac{d^{2}N_{\rm LO}^{\rm
conn.}}{dY_{1}d^{2}p_{\perp}dY_{2}d^{2}k_{\perp}}=\frac{\kappa_{2}}{(N_{c}^{2}-1)Q_{S}^{2}S_{\perp}}\,{dN_{\rm
LO}\over dY_{1}d^{2}p_{\perp}}\,\,{dN_{\rm LO}\over dY_{2}d^{2}k_{\perp}}\,,$
(53)

where $\kappa_{2}$ is a non-perturbative constant363636The results have a weak
dependent on the ratio $m/Q_{S}$, where $m$ is an infrared lattice regulator..
Again, the numerical simulations bear out the Glasma flux tube interpretation:
the likelihood that two particles are correlated is suppressed by the number
of flux tubes, and non-factorizable color connected graphs by
${\mathcal{O}}(1/N_{c}^{2})$. Perturbative arguments suggest that this picture
can be extended to $n$-particle cumulants and that the $n$-particle
multiplicity distribution that generates these cumulants is a negative
binomial distribution Gelis  _et al._ (2009). For $n$-particle multiplicities,
this expectation is confirmed by non-perturbative numerical simulations
Schenke  _et al._ (2012b).

####  IV.3.3 The IP-Glasma model

In the discussion thus far, color charge fluctuations on the scale $1/Q_{S}$
provide the only structure in the colliding gluon shockwaves. However nucleon
distributions in nuclei are not uniformly smooth and can fluctuate from event
to event. These fluctuations in nucleon positions are extremely important to
understand key features of the data such as the azimuthal moments $v_{n}$ of
the flow distributions at low momenta Alver and Roland (2010); Alver  _et al._
(2010). Another important ingredient in the realistic modeling of heavy-ion
collisions is the dependence of the saturation scale in the nuclei on $x$ (or
equivalently, $\sqrt{s}$), which describes the variations of particle
multiplicites in energy and rapidity at RHIC and the LHC.

We will outline here the IP-Glasma model Schenke  _et al._ (2012a, b, 2014a,
2014b), and improvements thereof, which incorporates the fluctuations in the
nucleon positions to construct event-by-event lumpy color charge distributions
and corresponding gluon field configurations in the LO Glasma framework. As we
will also discuss, the energy dependence of these configurations at a given Y
or $\sqrt{s}$ is determined by the saturation scales in the two nuclei.

An essential input is the dipole cross-section of the proton. The model we
consider here is the IP-Sat saturation model Bartels  _et al._ (2002);
Kowalski and Teaney (2003) which, as discussed in Sec. III.4, is an impact
parameter dependent generalization of the MV model. As noted, high precision
combined data from the H1 and ZEUS collaborations Aaron  _et al._ (2010);
Abramowicz  _et al._ (2013) are used to constrain the parameters of the
model373737A “b-CGC” model, incorporating the geometrical scaling shown in Eq.
(31) also provides good agreement with HERA data Rezaeian and Schmidt (2013).
Next-to-leading log BK computations too give excellent fits to inclusive data
Ducloué  _et al._ (2020) but do not include the $b_{\perp}$-dependence
essential for sub-nucleon structure. In the IP-Sat model, $b_{\perp}$
distributions are constrained by HERA exclusive vector meson data. and
excellent fits are obtained Rezaeian  _et al._ (2013).

The dipole cross-section for each nucleus at a given $x$ is constructed by
taking the product of the S-matrices corresponding to the dipole cross-
sections of overlapping nucleons at a given spatial location
${\mathbf{x}_{\perp}}$. It can be expressed as Kowalski  _et al._ (2008)

$\frac{1}{2}\frac{\mathrm{d}\sigma^{\textrm{A}}_{\textrm{dip}}}{{\mathrm{d}}^{2}{\mathbf{x}_{\perp}}}=\left[1-e^{-\frac{\pi^{2}}{2N_{c}}{\mathbf{r}_{\perp}}^{2}\alpha_{S}(Q^{2})xG(x,Q^{2})\sum_{i=1}^{A}T_{p}({\mathbf{x}_{\perp}}-{\mathbf{x}_{T}}^{i})}\right]\,,$
(54)

where $T_{p}$ stands for the Gaussian thickness function for each of the $A$
nucleons in each nucleus and $Q^{2}=4/\mathbf{r}_{\perp}^{2}+Q_{0}^{2}$, with
$Q_{0}$ fixed by the HERA inclusive data. The gluon distribution $xG(x,Q^{2})$
is parametrized at the initial scale $Q_{0}^{2}$ and then evolved up to the
scale $Q^{2}$ using LO DGLAP-evolution. We define the nuclear saturation scale
$Q_{S}=1/\sqrt{\mathbf{r}_{\perp,s}^{2}}$, at the $r_{\perp}=r_{\perp,s}$ for
which the argument of the exponential in Eq. (54) equals one-half. To obtain
the spatial dependence of $Q_{S}$, one self-consistently solves
$x=0.5\,Q_{S}({\mathbf{x}_{\perp}},x)/\sqrt{s}$ for every
${\mathbf{x}_{\perp}}$.

The result of this procedure is a lumpy distribution of
$Q_{S}^{2}({\mathbf{x}_{\perp}},x)$ denoting the sub-nucleon structure of the
nucleus. Since the IP-Sat model is a simple generalization of the MV model,
one can extract the variance of the color charge density
$g^{2}\mu_{A}^{2}({\mathbf{x}_{\perp}})$ at each $x$ from
$Q_{S}^{2}({\mathbf{x}_{\perp}},x)$ Lappi (2008). One then samples random
color charges $\rho^{a}({\mathbf{x}_{\perp}})$ on a transverse lattice,

$\langle\rho_{k}^{a}({\mathbf{x}_{\perp}})\rho_{l}^{b}({\mathbf{y}_{\perp}})\rangle=\delta^{ab}\delta^{kl}\delta^{2}({\mathbf{x}_{\perp}}-{\mathbf{y}_{\perp}})\frac{g^{2}\mu_{A}^{2}({\mathbf{x}_{\perp}})}{N_{y}}\,,$
(55)

where the indices $k,l=1,2,\dots,N_{y}$ label the $N_{y}$ points of
representing the width of the nucleus in $x^{-}$. The path ordered Wilson line