# The McDipper: A novel saturation-based 3+1D initial state model for Heavy
Ion Collisions

Oscar Garcia-Montero  Fakultät für Physik, Universität Bielefeld, D-33615
Bielefeld, Germany  Hannah Elfner  GSI Helmholtzzentrum für
Schwerionenforschung, Planckstr. 1, 64291 Darmstadt, Germany  Institute for
Theoretical Physics, Goethe University, Max-von-Laue-Strasse 1, 60438
Frankfurt am Main, Germany  Frankfurt Institute for Advanced Studies, Ruth-
Moufang-Strasse 1, 60438 Frankfurt am Main, Germany  Helmholtz Research
Academy Hesse for FAIR (HFHF), GSI Helmholtz Center, Campus Frankfurt, Max-
von-Laue-Straße 12, 60438 Frankfurt am Main, Germany  Sören Schlichting
Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany

(August 28, 2024)

###### Abstract

We present a new 3D resolved model for the initial state of ultrarelativistic
heavy-ion collisions, based on the $k_{\perp}$-factorized Color Glass
Condensate hybrid approach. The McDipper framework responds to the need for a
rapidity-resolved initial-state Monte Carlo event generator which can deposit
the relevant conserved charges (energy, charge and baryon densities) both in
the midrapidity and forward/backward regions of the collision. This event-by-
event generator computes the gluon and (anti-) quark phase-space densities
using the IP-Sat model, from where the relevant conserved charges can be
computed directly. In the present work we have included the leading order
contributions to the light flavor parton densities. As a feature, the model
can be systematically improved in the future by adding next-to-leading order
calculations (in the CGC hybrid framework), and extended to lower energies by
including sub-eikonal corrections the channels included. We present relevant
observables, such as the eccentricities and flow decorrelation, as tests of
this new approach.

##  I Introduction

It is widely accepted that during heavy ion collisions (HICs) a strongly
interacting medium of quarks and gluons is produced, the Quark-Gluon Plasma
(QGP) Shuryak (2017); Arslandok  _et al._ (2023). The QGP medium expands into
the vacuum and after a violent initial pre-equilibrium stage its subsequent
space-time evolution can be effectively described by relativistic viscous
hydrodynamics. During the expansion of this fluid, spatial anisotropies of the
initial state are converted into momentum space final-state particle
multiplicity correlations. These correlations and anisotropies inherited by
the final state can be found along both transverse and longitudinal
directions, with respect to the collisional axis.

It is a known result of this area of study that particle correlations in the
transverse plane reveal strong azimuthal angle modulation, quantified by the
so-called flow harmonics, $v_{n}$ Heinz and Snellings (2013). However,
recently the scrutiny of the longitudinal structure of the $v_{n}$ and also
higher correlations show a non-trivial rapidity dependence Petersen  _et al._
(2011); Nie (2019); Khachatryan  _et al._ (2015); Aaboud  _et al._ (2018).
These interesting results call for an improvement of the hydrodynamical models
used, which often focus on the midrapidity region, and thus rely on a boost-
invariant prescription. The evolution of the medium in a 3+1D fashion has been
successfully applied by such hydrodynamical codes as MUSIC Schenke  _et al._
(2010) and VHLLE Karpenko  _et al._ (2014). Once the reader has chosen the
model of their preference, what is left is to specify the initial condition
for the 3+1D, where they can chose from an array of models, including extended
MC-Glauber modeling Moreland  _et al._ (2015); Shen and Alzhrani (2020),
deposition via string-hadron transport Weil  _et al._ (2016); Schäfer  _et
al._ (2022), and color string dynamics Bozek and Broniowski (2016), amongst
others Ehehalt and Cassing (1995); Geiss  _et al._ (1998); Bratkovskaya  _et
al._ (2011).

In this work, we present a new framework to compute event-by-event (EbE)
longitudinally resolved initial conditions, based on the general ideas of
saturation physics. So far initial state models based on saturation physics,
such as IP-Glasma Schenke  _et al._ (2012a, b), EKRT Eskola  _et al._ (2000);
Niemi  _et al._ (2016) or MC-KLN Drescher and Nara (2007) have primarily been
designed to describe the initial state around mid-rapidity in 2+1D hydro
simulations. Despite various attempts to extend saturation based models to
3+1D Schenke and Schlichting (2016); Schlichting and Singh (2021); Schenke
_et al._ (2022); Ipp  _et al._ (2021, 2022), so far all of these calculations
exploit the fact that in the high energy limit the cross-section for gluon
production increases rapidly, quickly overtaking that of stopping of valence
quarks. However, even in the high energy limit, the so-called fragmentation
region in the forward and backward rapidity directions is populated by the
stopping of -near- collinear partons containing a high portion of the energy
and charge of the colliding nuclei. Evidently, in the forward and backward
rapidity regions the net-densities of the conserved electric ($Q$) and baryon
($B$) charges are non-vanishing due to the stopping of valence partons. Hence,
to properly describe the 3+1D initial state in heavy-ion collisions, one needs
to account not only for the (gluon-dominated) energy deposition at mid-
rapidity , but also one for the (valence-quark dominated) energy and charge
deposition in the fragmentation region. Evidently, a proper treatment of the
fragmentation region also becomes important, when trying to connect saturation
physics based initial state models developed at high-energy, with initial
state models at lower energies based e.g. on hadronic transport approaches
Schäfer  _et al._ (2022).

The central objective of this paper is to exploit first principles
calculations of quark and gluon scattering processes in high-energy QCD to
develop an initial state model for the deposition of energy and the conserved
$(u,d,s)$ charges in heavy-ion collisions. Specifically, we will use leading
order cross-section calculations in the ${\bf k}_{T}$-factorized dilute-dense
limit of the Color Glass Condensate effective theory of high-energy QCD Gribov
_et al._ (1983); Kovchegov and Tuchin (2002) to deduce the initial spectra of
quark Lappi and Mantysaari (2013); Dumitru and Jalilian-Marian (2002a, b) and
gluon production Gribov  _et al._ (1983); Kovchegov and Tuchin (2002);
Gyulassy and McLerran (1997); Blaizot  _et al._ (2010); Dumitru and McLerran
(2002) and infer the densities of energy and conserved $(u,d,s)$ charges by
taking appropriate moments of the spectra. By including event-by-event nucleon
position fluctuations as well as realistic parametrizations of the valence
parton distributions Buckley  _et al._ (2015) and the dipole scattering
amplitude (IPSat) Kowalski  _et al._ (2006); Kowalski and Teaney (2003), we
develop an efficient Monte Carlo generator – McDipper for the initial energy
and charge density profiles. Within this first study, we further benchmark the
model at the hand of multiplicity measurements and provide additional results
for the longitudinal structure of initial charge density profiles as well as
the longitudinal decorrelation of initial state eccentricities.

This paper is organised as follows, in sec. II we present an extensive outline
of the McDipper overview, where in sec. II.1 we detail the master formulas of
the framework, namely single gluon and quark production from the
$k_{\perp}$-factorized CGC formalism, as well as the building blocks for the
framework. In sec. II.2 we introduce two different sets of input models, the
GBW and the IP-Sat model, which will be used to benchmark the McDipper
framework. In sec. II.3 we briefly explain the Monte Carlo sampling used to
generate the event-by-event initial conditions. Subsequently, sec. II.4
contains a description of the code, namely an explanation of how the
aforementioned parts and models are used to produce rapidity resolved EbE HICs
initial states. The final subsection in the section corresponds to the fitting
procedure to the free-parameter of the framework, the so called $K$-factor.

In sec. III we give an overview of the generalities of the framework, a study
on the different dependencies on the energy and charge deposition. We explore
the thickness, energy and space-time rapidity dependencies, for both
saturation models presented in section II.2. This is performed system-
independently, meaning that we analyse charge deposition before the event
generation. In sec. IV we explore the EbE observables of the code, where we
estimate multiplicity, finding good agreement with data. We also present
results for rapidity resolved transverse energy and charge deposition, as well
as the long-distance initial correlations, the spatial eccentricities. We
discuss also the rapidity correlations which arise with the $x$-dependence of
the saturation models, where we present computations for flow-decorrelations.
We end by summarizing our findings and presenting the future possible research
avenues for this framework in sec. V. Additionally, in Appendix A we review
the collinear limit of the single gluon production formula, valid in the far-
fragmentation regions.

##  II Model Description

The central idea of our initial state model is to calculate the energy and
charge deposition in high-energy heavy-ion collisions within the dilute-dense
approximation of the CGC Gribov  _et al._ (1983); Kovchegov and Tuchin (2002);
Kowalski  _et al._ (2006); Kowalski and Teaney (2003). Based on this formalism
the single inclusive gluon
$(\mathrm{d}N_{g}/{\mathrm{d}^{2}\mathbf{x}\mathrm{d}^{2}\mathbf{p}\mathrm{d}y})$
and (net-) quark distributions
$\mathrm{d}N_{\bar{q}-q}/{\mathrm{d}^{2}\mathbf{x}\mathrm{d}^{2}\mathbf{p}\mathrm{d}y}$
as a function of (momentum) rapidity $y$ transverse momentum $\mathbf{p}$ and
transverse position $\mathbf{x}$ are computed by evaluating the leading order
cross-sections for the transverse momentum dependent dipole gluon
distributions from the IP-Sat model Bartels  _et al._ (2002); Kowalski and
Teaney (2003); Rezaeian  _et al._ (2013) and collinear parton distributions
(PDFs) from the LHAPDF library Buckley  _et al._ (2015). Since at leading
order in the high-energy limit the phase-space distributions are proportional
to $\delta(y-\eta_{s})$ Greif  _et al._ (2017), where $\eta_{s}$ is the space-
time rapidity, these distributions also provide the structure of initial state
in space-time rapidity upon the straightforward identification of
$y=\eta_{s}$. By taking moments of the single inclusive gluon and (net-) quark
distributions, we can then directly compute the initial 3+1D profiles of
energy ($e$) and $u,d,s$ charge distributions $(n_{u/d/s})$. By including also
fluctuating nucleon positions according to the Monte-Carlo Glauber model
Miller  _et al._ (2007); Alver and Roland (2010), our publicly available
implementation McDipper Oscar Garcia-Montero (2023) thereby provides event-by-
event initial conditions for all the conserved charges in the light flavour
sector of QCD.

Since our model is based on high-energy QCD cross-section calculations, it has
the advantages that a) the number of free parameters in the model is extremely
small and b) it can in principle be improved to higher order accuracy as next-
to-leading order cross-section calculations become available. Since most of
the input information can be extracted alternatively e.g. from Deep Inelastic
Scattering (DIS) experiments, the current version of the model only has a
single free parameter, namely – as we will discuss shortly – a so-called
$K$-factor, entering the normalization of the gluon production cross-section.

Figure 1: Diagrammatic representation for the leading order processes included
in the McDipper framework. (_Left_) Single gluon production of a gluon with
transverse moment $\mathbf{p}$ and rapidity $y$ from the interaction of two
CGCs with gluon distributions $\Phi_{1}(x_{1},\mathbf{q})$ and
$\Phi_{2}(x_{2},\mathbf{k})$. (_Right_) Deposition of a single quark with
transverse momentum $\mathbf{p}$ an rapidity $y$ from the stopping of a
collinear quark from the incoming projectile with momentum fraction $x$ and
virtuality $\mathbf{p}^{2}$.

###  II.1 Single inclusive gluon and quark production

Now that we have laid out the basic ideas underlying the McDipper we proceed
to provide the detailed expressions for the single inclusive gluon production
and valence quark stopping cross-sections from the literature along with the
description of the calculation of the associated contributions to the energy
and charge densities. In order to provide a compact overview the leading order
diagrams included in McDipper are depicted in Fig. 1.

Based on the $k_{\perp}$-factorized CGC formalism the single inclusive
spectrum of gluons $\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}$ produced
in a high-energy collision as a function of transverse momentum $\mathbf{p}$,
momentum rapidity $y$ and transverse position $\mathbf{x}$, can be expressed
as the convolution of the un-integrated gluon distributions
$\Phi(x,\mathbf{x},\mathbf{q})$ of the projectile and target as Dumitru and
McLerran (2002); Lappi and Schlichting (2018)

$\begin{split}\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}=&\frac{g^{2}\,N_{c}}{4\pi^{5}(N_{c}^{2}-1)\,\mathbf{p}^{2}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\leavevmode\nobreak\
\\\ &\times\leavevmode\nobreak\
\Phi_{1}(x_{1},\mathbf{x},\mathbf{q})\leavevmode\nobreak\
\Phi_{2}(x_{2},\mathbf{x},\mathbf{k})\\\ &\times\leavevmode\nobreak\
(2\pi)^{2}\delta^{(2)}(\mathbf{q}+\mathbf{k}-\mathbf{p})\\\ \end{split}$ (1)

where $x_{1/2}=\frac{\mathbf{p}}{\sqrt{s_{NN}}}e^{\pm y}$ denotes the light-
cone momentum fraction in the projectile and target, $g$ is the strong
coupling constant, and $N_{c}=3$ denotes the number of colors. We note that
the above formula, which exhibits $\mathbf{k}_{\bot}$-factorization, can be
rigorously derived in the case of dilute-dense or dense-dilute treatment of
the projectile-target system. However, it has also been demonstrated
explicitly that Eq. (1) provides a very good description of the energy
deposition even in the case of dense-dense collisions Blaizot  _et al._
(2010); Schlichting and Skokov (2020), which justify its use in the present
context.

Within the CGC framework of small-x physics, the unintegrated gluon
distribution function (uGDF) of the participant nucleus $i=(1,2)$ can be
expressed in terms of the correlation function of light-like Wilson lines
Blaizot  _et al._ (2004a); Gelis and Peshier (2002a, b). Specifically,

$\Phi_{i}(x_{i},\mathbf{x},\mathbf{q})=\frac{4\pi^{2}(N^{2}_{c}-1)}{g^{2}N_{c}}\,\mathbf{k}^{2}\,D_{\rm
adj}(x,\mathbf{x},\mathbf{q})\,$ (2)

where the adjoint dipole $D_{\rm adj}(x,\mathbf{x},\mathbf{q})$ is defined in
momentum space as

$D_{\rm
adj}(x,\mathbf{x},\mathbf{q})=\frac{1}{N_{c}^{2}-1}\int_{\mathbf{s}}\tr_{\rm
adj}[U^{\rm adj}_{\mathbf{x}+\mathbf{s}/2}U^{\rm
adj,\dagger}_{\mathbf{x}-\mathbf{s}/2}]\leavevmode\nobreak\
e^{i\mathbf{q}\cdot\mathbf{s}}\;.$ (3)

Here $U^{\rm adj}_{\mathbf{x}+\mathbf{s}/2}$ and $U^{\rm
adj,\dagger}_{\mathbf{x}-\mathbf{s}/2}$ denote adjoint light-like Wilson lines
at transverse positions $\mathbf{x}\pm\mathbf{s}/2$, such that the transverse
momentum scale, $\mathbf{q}$, is the Fourier transform of the relative size of
the dipole, $\mathbf{s}$, with fixed average position $\mathbf{x}$. By
substituting the above definitions into Eq. (1) one obtains the following
expression for single gluon production explicitly in terms of the (adjoint)
dipole amplitude,

$\begin{split}\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}=&\frac{N_{c}^{2}-1}{16\pi^{4}\alpha_{s}N_{c}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\leavevmode\nobreak\
\frac{\mathbf{q}^{2}\mathbf{k}^{2}}{\mathbf{p}^{2}}\\\
&\times\leavevmode\nobreak\ D_{\rm
adj}(x_{1},\mathbf{x},\mathbf{q})\leavevmode\nobreak\ D_{\rm
adj}(x_{2},\mathbf{x},\mathbf{k})\\\ &\times\leavevmode\nobreak\
(2\pi)^{2}\delta(\mathbf{q}+\mathbf{k}-\mathbf{p})\\\ \end{split}$ (4)

While at leading order the single inclusive production of gluons is due to
radiative processes, the leading order contribution to single inclusive quark
production in the high-energy limit originates from the stopping of collinear
quarks inside the projectile/target nuclei due to multi-scattering off the
constituents of the other nucleus. By accounting for the stopping of quarks
and anti-quarks from both nuclei, the single inclusive distribution
$\frac{dN_{q_{f}}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}$ for an (anti-)quark of
flavor $f(\bar{f})$ can be expressed as Dumitru and Jalilian-Marian (2002b);
Dumitru  _et al._ (2006)

$\begin{split}\frac{dN_{q_{f}}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}&=\frac{x_{1}q^{A}_{f}(x_{1},\mathbf{p}^{2},\mathbf{x})\leavevmode\nobreak\
D_{\rm fun}(x_{2},\mathbf{x},\mathbf{p})}{(2\pi)^{2}}\\\
&+\frac{x_{2}q^{A}_{f}(x_{2},\mathbf{p}^{2},\mathbf{x})\leavevmode\nobreak\
D_{\rm fun}(x_{1},\mathbf{x},\mathbf{p})}{(2\pi)^{2}}\,.\end{split}$ (5)

where $q^{A}_{f}(x_{1/2},Q^{2},\mathbf{x})$ are the collinear quark
distributions sinde the two nuclei (see below) and the kinematics of $x_{1,2}$
is the same as for single inclusive gluon production. Intuitively, this
expression can be understood as the deflection of a (collinear) quark inside
the projectile due to the transverse momentum transfer accumulated in multiple
interactions with gluons from the target nucleus. By including all
possibilities of multiple scattering, the quark receives a transverse position
dependent color rotation expressed in terms of a light-like (fundamental)
Wilson line in the amplitude and complex conjugate amplitude. It is then
straightforward to show that the probability of producing an in-medium quark
(non-collinear) as in Eq. (5), is just the product of the probabilities of
finding a collinear quark of flavor $f$ in the incoming projectile with the
probability of a single quark to be deflected by scattering off the target
Mäntysaari (2015), as described by the fundamental dipole distribution

$D_{\rm fun}(x,\mathbf{x},\mathbf{q})=\frac{1}{N_{c}}\int_{\mathbf{s}}\tr_{\rm
fun}[U_{\mathbf{x}+\mathbf{s}/2}U^{\dagger}_{\mathbf{x}-\mathbf{s}/2}]\leavevmode\nobreak\
e^{\rm i\mathbf{q}\cdot\mathbf{s}}\\\ $ (6)

By assuming local Gaussian correlations of color charges inside the nucleus,
the fundamental and adjoint dipole gluon distribution in Eq. (4) and (5) can
be related to each other (in coordinate space) as

$D_{\rm adj}(x,\mathbf{x},\mathbf{s})=[D_{\rm
fun}(x,\mathbf{x},\mathbf{s})]^{C_{A}/C_{F}}\;,$ (7)

and we will employ this relation throughout this work. With regards to the
collinear quark/anti-quark distributions, we simply assume that these can be
described in terms of uncorrelated partons, which are distributed in
transverse space according to the local density of protons ($p$) and neutrons
($n$). The corresponding distributions are then given by

$\begin{split}u^{A}(x,Q^{2},\mathbf{x})&=u_{p}(x,Q^{2})\;T_{p}(\mathbf{x})+u_{n}(x,Q^{2})\;T_{n}(\mathbf{x})\,,\\\
d^{A}(x,Q^{2},\mathbf{x})&=d_{p}(x,Q^{2})\;T_{p}(\mathbf{x})+d_{n}(x,Q^{2})\;T_{n}(\mathbf{x})\,,\\\
s^{A}(x,Q^{2},\mathbf{x})&=T_{p}(\mathbf{x})\,s_{p}(x,Q^{2})+T_{n}(\mathbf{x})\,s_{n}(x,Q^{2})\,.\\\
\end{split}$ (8)

for $u,d,s$ flavor quarks/anti-quarks, where by
$T_{p/n}(x)=\sum_{i=1}^{A}\leavevmode\nobreak\
t(\mathbf{x}-\mathbf{x}_{i})\,\delta_{i,p/n}$ we denote the thickness function
of protons and neutrons inside the nucleus, where

$t(\mathbf{x})=\frac{1}{2\pi
B_{G}}\exp\left[-\frac{\mathbf{x}^{2}}{2\,B_{G}}\right]\,,$ (9)

where the nucleon size, $B_{G}=0.156\mathrm{fm}^{2}$, is determined through
fits to the HERA data Kowalski and Teaney (2003); Rezaeian  _et al._ (2013).
The thickness function of an individual nucleon is then normalized as $\int
d^{2}\mathbf{x}\leavevmode\nobreak\ t(\mathbf{x})=1$ such that $\int
d^{2}\mathbf{x}\,[T_{p}(\mathbf{x})+T_{n}(\mathbf{x})]=A$ for each nucleus. In
this work we will only use a gaussian profile for the individual nucleons.

Since the parton distributions of the neutron are not well constrained, we
will assume isospin symmetry, which means that the neutron PDFs are equivalent
to the proton PDFs after applying $u\leftrightarrow d$, i.e.
$u_{p/n}=d_{n/p}$, such that Eq. (LABEL:eq:PDFS-pn) simplifies to

$\begin{split}u^{A}(x,Q^{2},\mathbf{x})&=u(x,Q^{2})\;T_{p}(\mathbf{x})+d(x,Q^{2})\;T_{n}(\mathbf{x})\,,\\\
d^{A}(x,Q^{2},\mathbf{x})&=d(x,Q^{2})\;T_{p}(\mathbf{x})+u(x,Q^{2})\;T_{n}(\mathbf{x})\,,\\\
s^{A}(x,Q^{2},\mathbf{x})&=\left(T_{p}(\mathbf{x})+T_{n}(\mathbf{x})\right)\,s(x,Q^{2})\,.\\\
\end{split}$ (10)

where we dropped the sub-script $(p)$ in all distributions.

While at first sight, the expression in Eqns. (1) and (5) appear to be rather
distinct, it is worth noting that in the forward limit, the single gluon
production formula yields a result equivalent in form to Eq. (5). In this
limit, the gluon distribution of the projectile is evaluated at large $x$,
where the projectile becomes dilute, and the uGDF can be expanded to leading
twist Blaizot  _et al._ (2004b, a). By integrating over the small intrinsic
transverse momentum of the projectile gluon (up to the transverse monentum
transfer from the target), one the uGDF then reduces to a collinear gluon PDF
yielding an expression similar to Eq. (5) (see Appendix A). We finally note
that – in the same spirit – the renormalization scale $Q^{2}$ of the quark
distributions in eq. (5) is set equal to the transverse momentum transfer from
the target nucleus ($Q^{2}=\mathbf{p}^{2}$), as this choice ensures that the
transverse momentum of the quark remains smaller than the transverse momentum
transfer $\mathbf{p}$ from the target, justifying the collinear factorization
treatment.

Based on the single inclusive spectra of quarks and gluons in Eqns. (1) and
(5), it is then straightforward to compute the local energy density and charge
distributions. Specifically, the total energy density is determined as the
first moment of the single particle distributions, namely 111We note that in
our numerical implementation, we add an infrared regulator by virtue of the
replacement $1/\mathbf{p}^{2}\to 1/(\mathbf{p}^{2}+m^{2})$ to circumvent
numerical problems when computing the contribution around an integrable
singularity. However, since the integral is UV dominated and infrared finite
even in the absence of a regulator, this regulator has no significant effect
on the results.

$(e\tau)_{0}=\int d^{2}\mathbf{p}\leavevmode\nobreak\
|\mathbf{p}|\leavevmode\nobreak\
\left[K_{g}\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}+\sum_{f,\bar{f}}\frac{dN_{q_{f}}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}\right]_{y=\eta_{s}}\;,$
(11)

while the conserved charges, $n_{i}\quad\text{with}i=(u,d,s)$, are computed
using the zeroth moment of the distribution, namely

$(n_{f}\tau)_{0}=\int d^{2}\mathbf{p}\leavevmode\nobreak\
\left[\frac{dN_{q_{f}}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}-\frac{dN_{\bar{q}_{f}}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}\right]_{y=\eta_{s}}$
(12)

We note that the above expressions for single inclusive quark and gluon
production provide the leading order process in the high-energy limit. In the
future, this treatment could be systematically improved by including higher
order corrections, such as e.g. the contributions of quark-antiquark pair
production by gluon fusion Gelis and Venugopalan (2004); Blaizot  _et al._
(2004a). Naturally, a systematic order-by-order improvement of the model also
applies to the inclusive gluon production processes, which will include
radiative corrections at next-to-leading order Iancu  _et al._ (2016);
Altinoluk  _et al._ (2018); Dumitru  _et al._ (2006).

###  II.2 Saturation models

Now that we have established the formalism to compute the energy and conserved
charge deposition, we continue to discuss the relevant input from saturation
physics. Within this exploratory study, we will consider two different
saturation models, namely the GBW Model and the IP-Sat Model, noting that in
both cases, the relevant parameters have been determined previously from the
analysis of DIS data.

####  II.2.1 GBW Model

We first consider a simple impact parameter dependent generalization of the
Golec-Biernat Wusthoff (GBW) model Golec-Biernat and Wusthoff (1999); Kowalski
and Teaney (2003); Kowalski  _et al._ (2006), where the fundamental dipole
gluon distribution is given by

$\displaystyle D_{\rm
fun}(\mathbf{x},\mathbf{q})=\frac{4\pi}{Q_{F}^{2}(x,\mathbf{x})}\exp\left[-\frac{\mathbf{k}^{2}}{Q_{F}^{2}(x,\mathbf{x})}\right]\,.$
(13)

Clearly, the key advantage of this model lies in its simplicity. Due to the
particularly simple Gaussian form of the dipole gluon distribution, it is
possible to obtain analytic expressions for the the dipole gluon distribution
both in transverse position and momentum space. Similarly, the adjoint dipole
distribution also exhibits the same Gaussian shape, with the difference that
the saturation scale is replaced by the adjoint saturation scale
$Q^{2}_{F}\rightarrow Q^{2}_{A}$ with $Q_{A}^{2}=(C_{A}/C_{F})\,Q_{F}^{2}$.
Based on fits to DIS data Kowalski and Teaney (2003); Rezaeian  _et al._
(2013), the kinematic dependence of the saturation scale is given
phenomenologically by the relation

$Q_{F}^{2}(x,\mathbf{x})=Q_{p,0}^{2}\,x^{-\lambda}\,(1-x)^{\delta}(T_{p}(\mathbf{x})+T_{n}(\mathbf{x}))\sigma_{0}$
(14)

where $Q_{p,0}^{2}=0.152\,$GeV2, $\lambda=0.215$, $\sigma_{0}\equiv 2\pi
B_{G}$ denotes the effective transverse area of the nucleon, such that
$(T_{p}(\mathbf{x})+T_{n}(\mathbf{x}))\sigma_{0}$ essentially counts the
density of nucleons per unit transverse area and $\delta=1$ is introduced to
regulate the large $x$ behavior.

####  II.2.2 IP-Sat Model

We will also consider the IP-Sat model Kowalski  _et al._ (2006); Kowalski and
Teaney (2003) which provides a more realistic expression for the interaction
of a color dipole with a nucleus, as it includes the perturbative tail of
large transverse momentum transfers. In the IP-Sat model, the fundamental
dipole is given by

$D_{\rm
fun}(x,\mathbf{x},\mathbf{s})=\exp\left[-\frac{\pi^{2}\mathbf{s}^{2}}{2\,Nc}\alpha_{S}(\mu^{2})\,xg(x,\mu^{2})\,T(\mathbf{x})\right]$
(15)

where $\mu^{2}=\mu^{2}_{0}+C/\mathbf{s}^{2}$, and $g(x,\mu^{2})$ is the parton
distribution function, which is initialized at a scale $\mu_{0}$ via the
parametrization

$xg(x,\mu^{2}_{0})=a_{g}\,x^{-\lambda_{g}}\,(1-x)^{5.6}\,.$ (16)

This initial condition is evolved along the kinematic range of $\mu^{2}$ using
LO DGLAP evolution. In this work we include the results using the parameter
set extracted in ref. Rezaeian  _et al._ (2013) which fixes the mass of the
charm quark mass to $m_{c}=1.27$ GeV, called in this framework Set 1. In the
code we also include the necessary files to generate data based on a second
parameter set, for which the mass of the charm is fixed to $m_{c}=1.4$ GeV, or
Set 2.

Based on the above IP-SAT parametrization, the adjoint dipole can be extracted
in coordinate space according to Eq. (7) and the corresponding momentum space
gluon distributions required for the calculation of the energy and charge
deposition are then obtained by a numerical Fourier transform with respect to
the relative transverse position $\mathbf{s}$.

###  II.3 Event-by-event sampling

Based on the above formalism, we can obtain the average energy and charge
deposition profiles for a given configuration of the nucleon positions.
Evidently, it would be possible to include further sub-nucleonic fluctuations
e.g. of the spatial and momentum distribution of valence partons, however this
is beyond the scope of the present work and we will restrict ourselves to
event-by-event fluctuations of the nucleon positions. Our procedure follows
common practice, but for the sake of completeness, we will document it below.

Starting point of the event-generation is the sampling of the positions of the
individual nucleons in each nucleus from a Wood-Saxon’s distribution, for
which we include deformation effects for the nuclei which present such
properties d’Enterria and Loizides (2021).

Now that we know the positions of each nucleon, two incoming nucleons are set
to interact via rejection sampling of the nucleon overlap probability
distribution McLerran and Tribedy (2016),

$\frac{\mathrm{d}P(b)}{\mathrm{d}^{2}\mathbf{b}}=\frac{1-\exp\left(-\sigma_{g}\,T_{NN}(\mathbf{b})\right)}{\int\mathrm{d}^{2}\mathbf{b}\left[1-\exp\left(-\sigma_{g}\,T_{NN}(\mathbf{b})\right)\right]}$
(17)

where $\sigma_{g}$ represents the total effective partonic cross section at a
given collision energy, and $T_{NN}$ stands for the nucleon-nucleon overlap,
which can be computed

$T_{NN}(\mathbf{b})=\int\mathrm{d}^{2}\mathbf{s}\,\,t\left(\mathbf{s}+\frac{\mathbf{b}}{2}\right)\,t\left(\mathbf{s}-\frac{\mathbf{b}}{2}\right)\,.$
(18)

Specifically, for the Gaussian nucleon thickness in Eq. (9), the normalization
factor of the overlap distribution has the following analytical solution,

$\begin{split}\sigma_{\rm eff}^{\rm
inel}(\tilde{\sigma})&=\int\mathrm{d}^{2}b\left[1-\exp\left(-\sigma_{g}T_{\mathrm{NN}}(b)\right)\right]\\\
&=4\pi^{2}\,B_{G}\left[\log\left(\frac{\sigma_{g}}{4\pi
B_{G}}\right)+\Gamma\left(0,\frac{\sigma_{g}}{4\pi
B_{G}}\right)+\gamma_{E}\right]\end{split}$ (19)

where $\Gamma(w,z)=\int_{z}^{\infty}dt\,t^{w-1}\,e^{-t}$ is the incomplete
Gamma function, and $\gamma_{E}$ is the Euler-Mascheroni constant. By matching
$\sigma_{\rm eff,inel}(\sigma_{g})=\sigma_{inel}(\sqrt{s_{\mathrm{NN}}})$ to
the inleastic nucleon-nucleon cross-section, for which we employ the
parametrization in ref.McLerran and Tribedy (2016),

Figure 2: Energy (_upper panels_) and electric and baryon charge deposition
(_lower panels_) as a function of space-rapidity, $\eta_{s}$ for three
different collisional energies, , computed using the IP-Sat model with the
CT18NNLO PDF set. The results here shown are presented for two different sets
of fixed thickness, $T_{1,p}=T_{2,p}=2\,{\rm fm}^{-2}$ (_left panels_) and
$T_{1,p}=1\,{\rm fm}^{-2}$, $T_{2,p}=4\,{\rm fm}^{-2}$ (_right panels_). The
neutron density is fixed to zero for both projectiles for this representation,
meaning $T_{1,n}=T_{2,n}=0$.

$\begin{split}\sigma_{\mathrm{inel}}&=25.2+0.05\log(\sqrt{s_{\mathrm{NN}}})+0.56\log^{2}(\sqrt{s_{\mathrm{NN}}})\\\
&+45.2(\sqrt{s_{\mathrm{NN}}})^{-0.9}+33.8(\sqrt{s_{\mathrm{NN}}})^{-1.1}\text{mb}\end{split}$
(20)

we can then determine the $\sigma_{g}$ at each center-of-mass energy222What it
is useful from eq. (19) is that one can find general dimensionless solutions
if one expresses it as
$\tilde{\sigma}_{inel}\equiv\frac{\sigma_{inel}}{2\pi\,B_{G}}=\pi\,\left[\log\left(\tilde{\sigma}_{g}\right)+\Gamma\left(0,\tilde{\sigma}_{g}\right)+\gamma_{E}\right]$
(21) where we have used that $\tilde{\sigma}=\sigma/(2\pi\,B_{G})$. This is
very convenient, as it will allow the code to have a permanent solution,
parametrized in a fit function, regardless of the $B_{G}$ parameter. The
solution of (21) is very well behaved, and can be fitted easily with
polynomial function. We use a fifth order polynomial to describe the function
accurately (to an error $<0.05\%$ in the relevant ranges,
$50-13000\;\text{GeV}$ )..

Nucleons that undergo at least one inelastic interaction with another nucleon
are counted as participants in the collision, and contribute to the overall
thickness function $T_{p/n}(\mathbf{x})$ of the respective nucleus which
determines the energy and charge deposition. Nucleons that do not interact
inelastically with any other nucleon are commonly referred to as spectators,
and simply disregarded in the calculation of energy and charge deposition.

###  II.4 Structure of the Code

Now that we have described the model, we briefly comment on the numerical
implementation McDipper Oscar Garcia-Montero (2023). Clearly, the bottleneck
of the numerical implementation of Eqns. (11) and (12) is the numerical
integration required for the calculation of the energy and charge densities
for every point in the three dimensional ($\mathbf{x}$, $\eta_{s}$) space.
However, since for a given center-of-mass energy $\sqrt{s}$ and PDF set the
transverse profiles only depend on the space-time rapidity $\eta$ and the
nuclear thickness $T_{1/2}$ of the colliding nuclei, the relevant information
can be pre-tabulated, allowing for an efficient generation of event-by-event
profiles of the conserved charges.

Below we illustrate this at the example of quark production, and refer to the
documentation Oscar Garcia-Montero (2023) of McDipper for a more concise
introduction to the code itself. Single inclusive quark production of a quark
with flavor $f$ depends linearly on the nuclear thickness $T_{i}(\mathbf{x})$
of the projectile (the nucleus containing the collinear quark) and non-
linearly on the thickness $T_{j}(\mathbf{x})$ of the target (the nucleus
containing small-x gluons), such that e.g. the $u$ charge density per unit
rapidity can be expressed as (c.f. Eq. (12)

$\begin{split}(n_{u}\tau)_{0}\leavevmode\nobreak\ =&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
T_{1,p}(\mathbf{x})\left(\,n^{(0)}_{u,12}(T_{2}(\mathbf{x}))-\,n^{(0)}_{\bar{u},12}(T_{2}(\mathbf{x}))\right)\\\
&+T_{1,n}(\mathbf{x})\left(\,n^{(0)}_{d,12}(T_{2}(\mathbf{x}))-\,n^{(0)}_{\bar{d},12}(T_{2}(\mathbf{x}))\right)\\\
&+T_{2,p}(\mathbf{x})\left(\,n^{(0)}_{u,21}(T_{1}(\mathbf{x}))-\,n^{(0)}_{\bar{u},21}(T_{1}(\mathbf{x}))\right)\\\
&+T_{2,n}(\mathbf{x})\left(\,n^{(0)}_{d,21}(T_{1}(\mathbf{x}))-\,n^{(0)}_{\bar{d},21}(T_{1}(\mathbf{x}))\right)\end{split}$
(22)

where $T_{i/j,p/n}(\mathbf{x})$ is the density of proton/neutrons in nucleus
$i/j$ and

$n^{(k)}_{f,ij}(T_{j}(\mathbf{x}))=\int\mathrm{d}^{2}\mathbf{p}\,|\mathbf{p}|^{k}\,\frac{x_{i}q_{f}(x_{i},\mathbf{p}^{2})\leavevmode\nobreak\
D_{\rm fun}(x_{j},T_{j}(\mathbf{x}),\mathbf{p})}{(2\pi)^{2}}$ (23)

are simple moments of the single-inclusive quark/anti-quark distribution,
which depend on the transverse position $\mathbf{x}$ only through the nuclear
thickness $T_{j}(\mathbf{x})$. By pre-tabulating the expressions in Eq. (23)
for a given center-of-mass energy $\sqrt{s}$ as a function of rapidity $\eta$
and thickness $T_{j}(\mathbf{x})$, it is then possible to directly evaluate
Eq. (22) without the need for any additional numerical integration.

Similarly, one can also calculate and pre-tabulate the deposited energy per
unit rapidity $(e\tau)_{0}(\mathbf{x})$, as a three dimensional function of
the space-time rapidity $\eta$ and the nuclear thickness
$T_{1}(\mathbf{x}),T_{2}(\mathbf{x})$ of the two nuclei for a given center-of-
mass energy $\sqrt{s}$. Even though the calculation of these tables typically
requires 2-3 CPU hours333The running time depends on the collisional energy.
Lower energies require relatively shorter computation times, since the more
rapidity cells are kinematically restricted or prohibited, which means that
the $x$ values in the integrand are outside of the allowed ranges, and
therefore less points in the grid need to be accessed. -for a 3D grid with 101
points in all $T_{1},T_{2}$ and $\eta_{s}$ directions-, the generation of an
individual event is then very efficient, and only requires 20-25 CPU seconds,
when computed on a grid with 101 points in rapidity and 256x256 points in the
transverse plane.

Figure 3: Energy (_upper panels_) and electric and baryon charge deposition
(_lower panels_) as a function of space-rapidity, $\eta_{s}$, for
$\sqrt{s_{\rm NN}}=5.02$ TeV, computed using both the IP-Sat and GBW model
with the CT18NNLO PDF set. The results here shown are presented for two
different sets of fixed thickness, $T_{1,p}=T_{2,p}=2\,{\rm fm}^{-2}$ (_left
panels_) and $T_{1,p}=1\,{\rm fm}^{-2}$, $T_{2,p}=4\,{\rm fm}^{-2}$ (_right
panels_). The neutron density is fixed to zero for both projectiles for this
representation, meaning $T_{1,n}=T_{2,n}=0$.

###  II.5 Determination of the $K$-factor

We already explained above, that the there is only one free parameter in the
model $K_{g}$, which enters the calculation of the energy deposition in Eq.
(11) and supposedly accounts for higher-order corrections in the leading order
expression to single inclusive gluon production in Eq. (1). We note that the
introduction of such a parameter is a well known procedure, and has been
employed in many previous works using the hybrid formalism Ma  _et al._
(2018a, b); Benić  _et al._ (2019). In this work we will determine $K_{g}$ by
matching our results to the experimentally measured values of transverse
energy per unit rapidity at mid-rapidity in min.bias p+p and p+Pb collisions,
i.e.

$\left\langle\frac{\mathrm{d}E_{\perp}}{\mathrm{d}y}\right\rangle_{\rm
exp}=\left\langle\frac{\mathrm{d}E_{q,\perp}}{\mathrm{d}y}\right\rangle+K_{g}\,\left\langle\frac{\mathrm{d}E_{g,\perp}}{\mathrm{d}y}\right\rangle\,,$
(24)

Experimental values for the transverse energy per unit rapidity are determined
from the ALICE data as in ref. ALI (2022),

$\begin{split}\left\langle\frac{\mathrm{d}E_{\perp}}{\mathrm{d}y}\right\rangle&=\langle
m_{\perp}\rangle\frac{1}{f_{\rm tot}}\frac{\mathrm{d}N_{\rm
ch}}{\mathrm{d}y}\\\ &\approx\frac{\langle m\rangle}{f_{\rm
tot}}\sqrt{1+a^{2}}J^{-1}(a,\eta)\frac{\mathrm{d}N_{\rm
ch}}{\mathrm{d}\eta}\end{split}$ (25)

where $\langle m_{\perp}\rangle$ and $\langle m\rangle=0.215\pm 0.001$ are the
average transverse mass and rest mass of all measured particles, respectively.
The prefactor $f_{\rm tot}=0.55\pm 0.01$ stands for the ratio of charged
particles to the total number produced. The parameter $a=\langle
p_{\perp}\rangle/\langle m\rangle$ is used to approximate $\langle
m_{\perp}\rangle$. Additionally, we need the Jacobian
$J(a,\eta)=\sqrt{1+1/(a\cosh{\eta})^{2}}$ to convert the pseudorapidity
dependent measured results into rapidity.

The transverse energy deposited is then matched to the mid-rapidity value for
p-p and p-Pb minimum bias collisions, using the relation (24). The
experimental errors are propagated accordingly to the $K_{g}$ factor. Using
this prefactor fitting scheme, we obtain $K_{g,\rm GBW}=1.25$ and $K_{g,\rm
IP-Sat}=1.85$.

##  III General features of the model

Before we proceed to present results for observables relating to the
transverse and longitudinal structure of our initial state model, we find it
important to further illustrate some of its general features.

We first showcase the rapidity dependence of energy and charge deposition in
Fig. 2, where the top panels show the initial energy $(e\tau)_{0}$ per unit
rapidity, carried by quark $(q)$ and gluon $(g)$ degrees of freedom, while the
bottom panels show the deposited amount of baryon $(B)$ and electric $(Q)$
charge $(\tau n)_{0}$ per unit rapidity. Different columns show the space-time
rapidity dependence of the deposition profile, for different configurations of
the nuclear thickness in the colliding nuclei, for a symmetric configuration
$T_{1,p}=T_{2,p}=2fm^{-2}$ (left) and an asymmetric configuration
$T_{1,p}=1fm^{-2}\;,\leavevmode\nobreak\ T_{2,p}=4fm^{-2}$ (right). Different
colored curves in each panel correspond to two different center of mass
energies, $\sqrt{s}_{NN}=200{\rm GeV}$ (blue) and $\sqrt{s}_{NN}=2.76{\rm
TeV}$ (red). Before we comment on the detailed features, it is important to
emphasize, that all the dependencies on space-time rapidity, center-of-mass
energy and nuclear thickness are entirely determined by the $x$ and $T$
dependence of the uGDS and PDFs, and if not stated otherwise we employ the IP-
Sat and the CT18NNLO PDF sets for the following figures.

When considering symmetric thickness profiles, shown in the left column of
Fig. 2, we find that the energy deposition at mid-rapidity is dominated by the
gluonic contribution $e_{g}$. Nevertheless both at RHIC and LHC energies the
quark contribution to the energy deposition $e_{q}$ is non-negligible. By
moving out to larger rapidities, the relative contribution of quarks to the
energy deposition increases; at the same time valence charge is deposited at
large rapidities resulting in pronounced peaks in the net baryon (B) and net
electric charge (Q) distribution. By comparing the results for different
center of mass energies, one finds that the gluon dominance at mid-rapidity
increases with increasing center of mass energy, which results in a rapid rise
of the overall energy deposition. Conversely, the overall amount of net-charge
deposition exhibits a much weaker dependence on the center of mass energy, as
the peaks mainly shift out to larger rapidity for increasing center of mass
energy.

Figure 4:  Energy (_upper panels_) and electric and baryon charge deposition
(_lower panels_) as a function of the proton thicknes, $T_{2,p}$, for fixed
$T_{1,p}=2\,{\rm fm}^{-2}$ and collisional energy $\sqrt{s_{\rm NN}}=5.02$
TeV, computed using both the IP-Sat and GBW model with the CT18NNLO PDF set.
The results here shown are presented for two different sets of fixed space-
time rapidity, $\eta_{s}=0$ (_left panels_) and $\eta_{s}=3$ (_right panels_).

By looking at the results in the right column of Fig. 2, one finds that the
energy and charge deposition for asymmetric thickness profiles, features
highly asymmetric rapidity distributions. Clearly, the peaks in the energy and
charge deposition move away from mid rapidity, indicating that geometric
fluctutations in the transverse plane can induce non-trivial modulations of
the rapidity profiles on an event-by-event basis. While the energy deposition
due to quarks only shows a relatively weak forward-backward asymmetry, this
effect is much more pronounced for the energy deposition of gluons and the
charge deposition due to valence quark stopping.

By comparing the results for the IP-Sat and GBW model, as depicted examplarily
in Fig. 3 for a center of mass energy $\sqrt{s_{NN}}=2.76$ TeV, we can gauge
the sensitivity to the underlying UGDs. While the total energy and charge
deposition turns out to be rather similar between the two saturation models,
moderate differences emerge in the rapidity dependence of the energy and
charge deposition, which we attribute to the rather different $\mathbf{p}$
dependence of the UGDs in the IP-Sat and GBW model.

Next we investigate the dependence of the energy and charge deposition on the
nuclear thickness, as shown in Fig. 4, where we fix the nuclear thickness
$T_{1}^{p}=2{\rm fm}^{-2}$ of protons of one of the nuclei and vary the
nuclear thickness of protons $T_{2}^{p}$ of the other nucleus. Neutron
thicknesses of both nuclei $T_{n}^{1/2}$ are set to zero in this comparison.
Irrespective of the underlying saturation model (GBW or IP-Sat), the $T_{1}$
and $T_{2}$ dependence energy and charge deposition is governed by the $T$
dependence of the valance parton flux and the saturation scales of both
nuclei. Interestingly, the quark and gluon contributions to the energy
deposition exhibit a rather different $T$ dependence. In the case of the gluon
contribution $e_{g}$, one observes a rapid rise for $Q_{s}^{(2)}<Q_{s}^{(1)}$
followed by a saturation of growth for $Q_{s}^{(2)}>Q_{s}^{(1)}$. Conversely,
the energy deposition due to quark stopping $e_{q}$ exhibits a steady rise
across the entire range of $T$ values depicted in Fig. 4. It is particularly
important to point out, that due to the asymmetry in the saturation scales,
the energy deposition profiles at different rapidities and center of mass
energies, are different both in regards to the $T_{1}$ and $T_{2}$ dependence
as well as in the overall normalization. Notably, this behavior can be
intuitively understood in the GBW model, if the $x$ dependence of the UGDs in
Eq. (1) is approximated by self-consistently evaluating
$x_{1/2}=Q_{s}^{1/2}(x)/\sqrt{s_{NN}}e^{\pm y}$ instead of accounting for full
$\mathbf{p}$ dependence of the UGDs. Since the $\mathbf{p}$ dependence of $x$
is ignored in this approximation, one can integrate analytically as in ref.
Borghini  _et al._ (2022) to obtain an expression for the energy density. Such
an exression will depend on the saturation scales parametrically as follows

Figure 5: Comparison of the energy density (_top_), electric charge density
$Q$ (_middle_) and baryon charge density $B$ (_bottom_) deposition when using
different PDFs. Here are pictured the conserved charges when using the
CT18NLO, CT18NLO and CT14NNLO sets. As the reader can easily see, while there
are significant changes in the quark energy density, $Q$ and $B$ do not
change, just as expected.

$\displaystyle
e_{g}\propto\frac{Q_{s,1}^{2}Q_{s,2}^{2}\big{(}Q_{s,1}^{4}+\frac{7}{2}Q_{s,1}^{2}Q_{s,1}^{2}+Q_{s,2}^{4}\big{)}}{(Q_{s,1}^{2}+Q_{s,2}^{2})^{5/2}}\,,$
(26)

where according in the small-$x$ limit,
$x_{1/2}=Q_{s}^{1/2}(x)/\sqrt{s_{NN}}e^{\pm y}$ can be self-consistently
solved to fin the analytical expression,

$\displaystyle Q_{s,1/2}^{2}\propto
T_{1/2}^{\frac{2}{2+\lambda}}e^{\mp\frac{2\lambda}{2+\lambda}\eta}\Big{(}\sqrt{{s}_{NN}}\Big{)}^{\frac{2\lambda}{2+\lambda}}.$
(27)

We depict this result, Eq. 26, as a gray line in Fig. 4. Even though such a
simple parametrization can approximately describe the behavior over a large
range of nuclear thickness, it fails to describe important details that emerge
from the non-trivial $x$-dependence of the underlying parton distributions.

Before we turn to phenomenological results, we finally investigate the
sensitivity of the energy and charge deposition to the underlying PDF set. We
provide a compact summary of our results in Fig. 5, where we compare the
rapidity dependence of energy and charge deposition for the CT18NNLO,CT18NLO
and CT14NNLO PDF sets. Since the large $x$ structure of the valence quark PDFs
is well constrained, the charge deposition profiles are virtually unchanged.
Slight differences emerge for the energy deposition, that probes the sea quark
PDFs, in particular at larger rapidities where the small $x$ is explored.
Nevertheless, the profiles are remarkably stable with regards to the PDF
uncertainties.

Figure 6: Charged particle yield as a function of centrality for $\sqrt{s_{\rm
NN}}=200$ GeV (_top_), $\sqrt{s_{\rm NN}}=2.76$ GeV (_center_) and
$\sqrt{s_{\rm NN}}=5.02$ GeV (_Bottom_). Experimental data shown in the figure
for Au-Au at $\sqrt{s_{\rm NN}}=200$ GeV is from ref. Abelev  _et al._ (2009),
Pb-Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV from ref. Abbas  _et al._
(2013); Aamodt  _et al._ (2011), and $\sqrt{s_{\rm NN}}=5.02$ TeV from ref.
Adam  _et al._ (2017, 2016). The bands represent the uncertainties in $\eta/s$
and $C_{\infty}$ for the estimates here presented, where they correspond to
variations of $C_{\infty}=0.8-1.15$ and $\eta/s=0.08-0.24$.

##  IV Observables

Now that we have established the basic features of our initial state model, we
continue to constrast the model with experimental results at RHIC and LHC.
Clearly, the results presented below are only meant to provide a general idea
of the behavior of the McDipper model as we anticipate more detailed
phenomenological studies to be performed based on full 3+1D Initial State +
Hydrodynamics + Hadronic Cascade simulations in the future.

Figure 7:  Eccentricities in the IP-Sat (full symbols) and GBW (empty symbols)
models as a function of centrality, $\eta_{s}$. We show here the $n=2,3,4$
harmonics (blue squares, red circles and yellow diamonds, repectively) for Pb-
Pb collisions at $2.76$ TeV at midrapidity $\eta_{s}=0$. Included is a
comparison to ${\rm T}_{\rm R}$ENTo Moreland  _et al._ (2015); Borghini  _et
al._ (2022) and IP-Glasma Schenke  _et al._ (2020) models.

Before we explore the full 3D structure of the model it is important to cross-
check that observables at mid-rapidity can be reasonably well described within
we model. We first investigate the centrality dependence of the charged
particle multiplicity $\mathrm{d}N_{\rm ch}/\mathrm{d}y$, which following
Giacalone  _et al._ (2019) can be estimated directly from the initial energy
density profile as

$\begin{split}\left\langle\frac{\mathrm{d}N_{\rm
ch}}{\mathrm{d}y}\right\rangle=&\frac{4}{3}\frac{N_{\rm
ch}}{S}\,C_{\infty}^{3/4}\,\left(4\pi\,\frac{\eta}{s}\right)^{1/3}\left(\frac{\pi^{2}}{30}\nu_{\rm
eff}\right)^{1/3}\\\ &\times\int\mathrm{d}^{2}\mathbf{x}\left[\tau
e(y,\mathbf{x})\right]^{2/3}_{0}\end{split}$ (28)

where $\nu_{eff}=40$, $S/N_{\rm ch}=7.5$. We present the results in Fig. 6,
where we compare the centrality dependence of $\mathrm{d}N_{\rm
ch}/\mathrm{d}y$ to experimental measurement at RHIC Abelev  _et al._ (2009)
and LHC energies Adam  _et al._ (2017, 2016); Abbas  _et al._ (2013); Aamodt
_et al._ (2011). Bands signify the uncertainties in $\eta/s$ and $C_{\infty}$
in these estimates and correspond to variations of $C_{\infty}=0.8-1.15$ and
$\eta/s=0.08-0.24$. We find that within these uncertainties, the model
describes the centrality dependence of the data rather well; in particular
also the energy dependence of the charged particle multiplicity is well
described within the McDipper.

Next we turn to behaviour of the event-by-event eccentricities at mid-
rapidity, which provide knowledge about collective flow and geometric
fluctuations. We present our results for the cumulants
$\varepsilon_{2}\\{2\\},\varepsilon_{3}\\{2\\}$ and $\varepsilon_{4}\\{2\\}$
in Fig. 7, where we plot $\varepsilon_{n}\\{2\\}$ as a function of centrality.
The eccentricities are computed using the definition

$\varepsilon_{n}(y)=\frac{\int\mathrm{d}^{2}r_{\perp}|\mathbf{r}_{\perp}|^{n}\,e^{\mathrm{i}n(\phi-\Phi_{n})}\left[e(y,\mathbf{r}_{\perp})\tau\right]_{0}}{\int\mathrm{d}^{2}r_{\perp}\left[e(y,\mathbf{r}_{\perp})\tau\right]_{0}}$
(29)

with
$\varepsilon_{n}\\{2\\}=\sqrt{\left\langle|\varepsilon_{n}|^{2}\right\rangle}$,
see ref. Schenke  _et al._ (2022). Different panels show the results for
$Au+Au$ collisions at RHIC (top) and Pb+Pb collisions at LHC (center,bottom)
and we also compare the results of the McDipper to calculations in the ${\rm
T}_{\rm R}$ENTo Moreland  _et al._ (2015); Borghini  _et al._ (2022)444We use
parameters extracted from a Bayesian analysis from ref. Liyanage  _et al._
(2023), namely $p=0.038$, $w=0.985$ in the ${\rm T}_{\rm R}$ENTo
parametrization. For the parameter $N$, we used $N=6.39$ for Au-Au at
$\sqrt{s_{\rm NN}}=200$ GeV, $N=20.013$ for Pb-Pb at $\sqrt{s_{\rm NN}}=2760$
GeV and $N=25.96$ for Pb-Pb at $\sqrt{s_{\rm NN}}=5020$ GeV. and IP-Glasma
Schenke  _et al._ (2020) models. We find excellent agreement with the IP-
Glasma model for central collisions $<50\%$ centrality and small differences
for larger centralities. We find that, in comparison to the ${\rm T}_{\rm
R}$ENTo model, the eccentricities obtained in our model are rather large, as
has already been observed in earlier comparisions of the MC-Glauber and MC-KLN
models, which are similar in spirit. Given the fact that the two rather
different models TrENTO and IP-Glasma can both describe the experimental
results rather well, it seems likely – with suitable adjustment of
hydrodynamic evolution parameters – experimental results for collective flow
can also be described by the McDipper.

Figure 8: Energy (a) and charge (b) profiles for different rapidity windows in
a single 5.02 TeV Pb+Pb event in the $0-5\%$ centrality class. The energy
profiles are given for gluon energy deposition (a,upper pannels), and quark
energy density (a, lower pannels). The rapidity windows are chosen for
$\eta=0$ (left pannels) and $\eta=4$ (left pannels). In (b) the deposited
electric (b, upper pannels) and baryon (b, lower pannels) charges are shown
for the same rapidity windows.

Now that we have characterized the most important properties our initial state
model at mid rapidity, we continue with the exploration of the three
dimensional structure of the model. Before we investigate any concrete
observables, it proves insightful to illustrate the features of the model at
the example of the energy and charge profile in a single event. Different
panels in the left part of the figure show the energy deposition due to
quarks(bottom) and gluons (top) at mid-rapidity (left) and forward rapidity
(right); the right part of the figure shows the corresponding electric (top)
and baryon (bottom) charge deposition profiles at mid-rapidity and forward
rapidity. Generally, the geometric fluctuations of the energy and charge
profiles are dictated by the positions of the nucleons inside the colliding
nuclei. However, due to the different dependencies on the nuclear thickness of
the two nuclei (c.f. Sec. II.3), these fluctuations leave a different imprint
on the energy and charge profiles at different rapidities. One observes for
example, that the relative fluctuations of the energy density profile are
larger at mid-rapidity than at forward rapidity. Similarly, the geometric
fluctuations of the charge density are also different than those of the energy
density, and on an event-by-event basis none of the charge profiles is
symmetric w.r.t to the forward and backward rapidity direction. Since all of
these features naturally emerge within our framework, it will be interesting
to further explore the phenomenological consequences of these findings in the
future.

Figure 9: Transverse deposition of energy (_top pannels_), electric charge
(_central pannels_) and baryon charge (_lower pannels_). Results are presented
for Au-Au collisions at $\sqrt{s_{\rm NN}}=200$ GeV (_left pannels_) and Pb-Pb
collisions at $\sqrt{s_{\rm NN}}=200$ GeV (_right pannels_) for different
centrality classes. The computation is performed using the IP-Sat model with
the CT18NNLO PDF set.

Next we consider the rapidity dependence of the energy and charge deposition
profiles in Fig. 9, which we obtain by averaging over several events and
integrating over transverse space. Different columns in Fig. 9 show the
results for RHIC and LHC energies, while different curves in each panel show
the results for different centralities. While the energy profiles are relative
flat around mid-rapdity, they do not exhibit a clear mid-rapidity plateau,
which as discussed in Sec. III can primarily be attributed to the rapidity
dependence of the gluon energy deposition. Charge deposition profiles for the
net-electric and net-baryon charge are peaked at forward and backward
rapidities, with the peaks moving out to larger rapidities with increasing
center of mass energy. By comparing the profiles at RHIC and LHC energies, one
also finds a stronger suppression of net-baryon number and net-electric charge
at mid-rapidity for LHC energies, which is consistent with measurements of the
hadrochemistry. However, one should caution that over the course of the space-
time evolution of a heavy-ion collision the initial charge deposition profiles
are subject to diffusion, and can thus not be directly compared with
experimental data.

Despite the fairly pronounced rapidity dependence of the energy density and
charge profiles, we find that the rapidity dependence of the spatial
eccentricities is rather weak. This can be observed from Fig. 10, in which the
three panels show the rapidity dependence of the event-by-event eccentricities
$\varepsilon_{2}\\{2\\}$,$\varepsilon_{2}\\{3\\}$, and
$\varepsilon_{2}\\{4\\}$, for three selected centrality classes, $0-5\%$,
$10-20\%$ and $30-40\%$. Except at very large rapidities in the forward and
backward direction, the eccentricities are almost flat across the entire range
of rapidities. Notably, in this regime, the results of the GBW and IP-Sat
saturation models can also differ significantly from each other, which can be
attributed to the fact that the two models exhibit a rather different behavior
at large $x$, which is probed in the far forward and backward regions.

Figure 10: Eccentricities in the IP-Sat (full symbols) and GBW (empty symbols)
models as a function of space-time rapidity, $\eta_{s}$. We show here the
$n=2,3,4$ harmonics (in descending order) for Pb-Pb collisions at $2.76$ TeV
and three centrality classes, $0-5\%$, $10-20\%$ and $30-40\%$.

Ultimately, we also investigate the longitudinal decorrelation of the event
geometry, which following the experimental procedure Nie (2019); Khachatryan
_et al._ (2015) can be quantified by the forward/backward ratio Petersen  _et
al._ (2011); Pang  _et al._ (2018, 2016)

$r_{n}(\eta_{a},\eta_{b})=\frac{\left\langle\epsilon_{n}(-\eta_{a})\epsilon^{*}_{n}(\eta_{b})\right\rangle}{\left\langle\epsilon_{n}(\eta_{a})\epsilon^{*}_{n}(\eta_{b})\right\rangle}$
(30)

Our results are compactly summarized in Figs. 11 and 12, where we present the
decorrelation ratios $r_{2}$ (left), $r_{3}$ (middle) and $r_{4}$ (right) for
$5-10\%$ (top) and $30-40\%$ centrality (bottom) at RHIC and LHC energy.
Besides the results from the McDipper, we also show experimental measurements
from STAR Nie (2019) (preliminary) and CMS Khachatryan  _et al._ (2015) for
the ratios $r_{2}$ and $r_{3}$ based on the final state $v_{n}$s. We first
note that in comparison to experimental data, our model generally tends to
underestimate the de-correlation of the event geometry with varying
$\eta_{a}$. One obvious explanation is that the model is missing additional
sources of fluctuations; in particular the current model treats the valence
content of the nuclei in terms of inclusive (average) PDFs and does not take
into account event-by-event fluctuations of the energy-momentum fractions
carried by valence partons. Naturally including such fluctuations, e.g. along
the lines of Shen and Schenke (2022), will induce additional sources of
geometric fluctuations in the model and thereby increase the slope of the
rapidity decorrelation. We also observe large discrepancies for the
longitudinal decorrelation between the IP-Sat and GBW model, which we found
can be attributed to the rather different large $x$-behavior in these
models555The approach $x\rightarrow 1$ can be quantified by taking the second
moment of the dipole function, namely
$\langle\mathbf{k}^{2}\rangle_{i}=\int\mathrm{d}^{2}\mathbf{k}\mathbf{k}^{2}\,D_{\rm
i}(x,\mathbf{x},\mathbf{q})/\int\mathrm{d}^{2}\mathbf{k}\,D_{\rm
i}(x,\mathbf{x},\mathbf{q})$, where the representation is given by $i=(\rm
adj,fun)$ for the fundamental case. In the gaussian GBW case, eq. (15), this
corresponds trivially to $Q^{2}_{i}\sim(1-x)^{\delta}$, with $\delta=1$ as
defined. Nevertheless, the IP-Sat corresponds to a more nuance case, as the
effective extracted $\delta$ presents dependence on the the thickness
function, with $\delta\sim 4.6-5.6$. , which is probed at the reference
rapidity $\eta_{b}$. In fact when comparing the results across such large
rapidity intervals, the underlying physics changes dramatically, as e.g. in
our model the energy densities in the forward regions are no longer dominated
by gluon radiation, but rather by quark stopping. Generally speaking, the
simultaneous description of the mid-rapidity and fragmentation region provides
a challenge for any theoretical description that is not QCD. Even though the
choice of a large reference rapidity $\eta_{b}$ certainly has advantages from
an experimental point of view, it would be very interesting and highly
informative to perform such measurements within a more restricted rapidity
range, where theoretical models developed for the physics around mid-rapidity
remain applicable.

Figure 11: Flow decorrelation of initial spatial eccentricities, $r_{n}$, in
Au-Au collisions $\sqrt{s}=200$ GeV for the GBW and IP-Sat models. In this
figure we present the results for the harmonics $n=2,3,4$ (_(_ left), _(_
center) and _(_ right) panels, respectively) and a reference rapidity of
$\eta_{b}=2.5-4.7$ Nie (2019) Figure 12: Comparison of the decorrelation of
initial spatial eccentricities,$r_{n}$, of Pb-Pb collisions $\sqrt{s}=2.76$
TeV for the GBW and IP-Sat models. In this figure we present the results for
the harmonics $n=2,3,4$ and a reference rapidity of $\eta_{b}=4.4-5.0$. The
black squares correspond to experimental data from the CMS collaboration
Khachatryan  _et al._ (2015).

##  V Conclusions

We presented a novel approach to calculate energy and charge deposition in
high-energy heavy-ion collisions, based on the leading order cross-section
calculations in the Color Glass Condensate effective theory of high-energy
QCD. Since the underlying quark and gluon distributions can be independently
constrained from DIS measurements, the resulting model – McDipper – only
features a single free parameter $K_{g}$, that is supposed to account for
higher order corrections to the leading order perturbative cross-sections.

Based on this framework, the model makes concrete predictions for the
dependence of energy and charge deposition on the nuclear thickness, center of
mass energy and space-time rapidity, which naturally emerge from the
$x$,$Q^{2}$ and $T$ dependence of the underlying (nuclear) parton
distributions. We therefore believe that the model is not only useful for
immediate phenomenological applications, but may also aid the development of
parametric models of the longitudinal and transverse structure of the initial
state of high-energy heavy-ion collisions.

Our first comparisons with experimental measurements show promising results,
except for the longitudinal decorrelation of the transverse geometry which is
underestimated in the current implementation of the model. However, as we
discussed this may hint at the necessity to include additional sources of
event-by-event fluctuations, in particular of the valence quark distributions,
which provides an interesting direction for future investigations.

Based on our analysis, we also noticed that the physics in the far forward and
backward region, which is dominated by nuclear fragmentation, is distinctly
different than the physics around mid-rapidity, which is dominated by new
particle production, making it hard to describe and interpret measurements
that correlate the different rapidity regions. Even though the development of
a comprehensive model remains an interesting theoretical challenge, it may
also be beneficial to study e.g. de-correlation observables within a more
restricted rapidity range, where the underlying physics remains the same.

We finally note that the model, as it is firmly grounded in QCD, is in
principle systematically improvable by including higher-order perturbative
corrections Gelis and Venugopalan (2004) and sub-leading terms in the eikonal
expansion Chirilli (2019). Even though this is a rather challenging task, we
believe that this is still a promising avenue, as with the operation of a
future Electron-Ion Collider, a more and more precisely knowledge on the
underlying (nuclear) parton distributions will emerge.

## Acknowledgements

OGM would like to thank Travis Dore, Philip Plaschke and Stephan Ochsenfeld
for discussions. We would like to thank Bjoern Schenke and Wenbin Zhao for
providing the IP-Glasma data shown in this work. Additionally we would like to
thank Clemens Werthmann for providing us with the ${\rm T}_{\rm R}$ENTo data.
This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the CRC-TR 211 ‘Strong-interaction matter under
extreme conditions’ — project number 315477589 — TRR 211. OGM and SS
acknowledge also support by the German Bundesministerium für Bildung und
Forschung (BMBF) through Grant No. 05P21PBCAA.

##  Appendix A Collinear limit

In this section we will show a brief derivation of the collinear limit for the
gluon production formula, eq. (1). This formula has also been derived in the
hybrid formalism Dumitru and Jalilian-Marian (2002b), and it is valid for
gluon forward/backward kinematics, where a large-$x$ gluon scatters of a dense
target. For us this means particularly the large $|\eta_{s}|$ regions.

We start by integrating the momentum conservation delta function,

$\begin{split}\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}=&\frac{g^{2}\,N_{c}}{4\pi^{5}(N_{c}^{2}-1)\,\mathbf{p}^{2}}\int\frac{\mathrm{d}^{2}\mathbf{q}}{(2\pi)^{2}}\\\
&\times\leavevmode\nobreak\
\Phi_{1}(x_{1},\mathbf{x},\mathbf{q})\leavevmode\nobreak\
\Phi_{2}(x_{2},\mathbf{x},\mathbf{p}-\mathbf{q})\end{split}$ (31)

For simplicity we will focus specifically in the forward region, where
$\eta_{s}>0$ and large. In this case, the right moving nucleon has $x_{1}\sim
1$ while the left moving nucleon possesses the opposite behaviour, $x_{2}\ll
1$. Since the characteristic scale, the saturation scale, of these
distribution increases with decreasing x, at forward rapidities $\Phi_{1}$
will peak at small values of $k_{\perp}$. Since $\Phi_{2}$ will be dominated
by large modes we can take the limit in which $|\mathbf{p}|>|\mathbf{q}|$ and
expand around it. For the purposes of this work, we can keep only the zeroth
order,
$\Phi_{2}(x_{2},\mathbf{x},\mathbf{p}-\mathbf{q})\rightarrow\Phi_{2}(x_{2},\mathbf{x},\mathbf{p})$,
while higher corrections come as gradients in $\mathbf{q}$. In this way we get
a simplified formula, namely

$\begin{split}\frac{dN_{g}}{d^{2}\mathbf{x}d^{2}\mathbf{p}dy}=&\frac{1}{(2\pi)^{2}}x\,g_{A}(x_{1},\mathbf{x},\mathbf{p})\,D_{\rm
adj}(x_{2},\mathbf{x},\mathbf{p})\end{split}$ (32)

where we have defined the nuclear PDF as

$\begin{split}xg_{A}(x,\mathbf{x},\mathbf{p}^{2})=&\frac{N_{c}^{2}-1}{16\pi^{4}\alpha_{s}N_{c}}\leavevmode\nobreak\
\int^{|\mathbf{q}|<|\mathbf{p}|}\mathrm{d}^{2}\mathbf{q}\\\
&\,\,\times\mathbf{q}^{2}\leavevmode\nobreak\ D_{\rm
adj}(x,\mathbf{x},\mathbf{q})\end{split}$ (33)

Notice that our definition of the gluon distribution is defined with
dependence to the transverse position in the nucleus, which makes it in fact a
parton density. For the case of the total nuclear PDF, we start by integrating
with respect to the transverse position $\mathbf{x}$,

$\begin{split}xg_{A}(x,\mathbf{p}^{2})=&\frac{N_{c}^{2}-1}{16\pi^{4}\alpha_{s}N_{c}}\leavevmode\nobreak\
\int\mathrm{d}^{2}\mathbf{x}\int\mathrm{d}^{2}\mathbf{q}\\\
&\,\,\times\mathbf{q}^{2}\leavevmode\nobreak\ D_{\rm
adj}(x,\mathbf{x},\mathbf{q})\,\\\
=&\frac{N_{c}^{2}-1}{16\pi^{4}\alpha_{s}N_{c}}\leavevmode\nobreak\
\int\mathrm{d}^{2}\mathbf{x}\int\mathrm{d}^{2}\mathbf{s}\int\mathrm{d}^{2}\mathbf{q}\,\,\mathbf{q}^{2}\\\
&\,\,\times e^{\mathrm{i}\mathbf{q}\cdot\mathbf{s}}\leavevmode\nobreak\
\frac{1}{N_{c}^{2}-1}\tr_{\rm adj}[U^{\rm adj}_{\mathbf{x}+\mathbf{s}/2}U^{\rm
adj,\dagger}_{\mathbf{x}-\mathbf{s}/2}]\,,\end{split}$ (34)

where we have used the definition of the dipole, eq. 3, in the last line. By
defining $\mathbf{x}_{1}=\mathbf{x}+\mathbf{s}/2$ and
$\mathbf{x}_{2}=\mathbf{x}-\mathbf{s}/2$, we can transform our formula to

$\begin{split}xg_{A}(x,\mathbf{p}^{2})=&\frac{1}{16\pi^{4}\alpha_{s}N_{c}}\leavevmode\nobreak\
\int\mathrm{d}^{2}\mathbf{x}_{1}\int\mathrm{d}^{2}\mathbf{x}_{2}\int\mathrm{d}^{2}\mathbf{q}\\\
&\,e^{\mathrm{i}\mathbf{q}\cdot(\mathbf{x}_{1}-\mathbf{x}_{2})}\,\tr_{\rm
adj}[(\partial_{i}U^{\rm adj}_{\mathbf{x}_{1}})(\partial_{i}U^{\rm
adj,\dagger}_{\mathbf{x}_{2}})]\,.\end{split}$ (35)

In the last line is obtain via partial integration. Using the definition of
the Wilson line

$U^{\rm
adj}_{\mathbf{x}}=\exp\left[-\mathrm{i}g\int_{-\infty}^{\infty}\mathrm{d}x^{+}\,A^{-}_{a}(x^{+},\mathbf{x})\,T^{a}\right]\,,$
(36)

where $T^{a}$ are the generators of the $SU(3)$ in the adjoint representation,
one can express the trace in eq. (38) in the dilute limit to be

$\begin{split}\tr_{\rm adj}[(\partial_{i}U^{\rm
adj}_{\mathbf{x}_{1}})(\partial_{i}U^{\rm
adj,\dagger}_{\mathbf{x}_{2}})]=g^{2}\,N_{c}\int_{-\infty}^{\infty}\mathrm{d}x_{1}^{+}\\\
\times\int_{-\infty}^{\infty}\mathrm{d}x_{2}^{+}\left\langle(\partial_{i}A^{-}_{a}(x_{1}^{+},\mathbf{x}_{1}))(\partial_{i}A^{-}_{a}(x_{2}^{+},\mathbf{x}_{2}))\right\rangle\end{split}$
(37)

Mixing all together, we get the expression

$\begin{split}xg_{A}(x,\mathbf{p}^{2})=&\frac{2}{(2\pi)^{3}}\leavevmode\nobreak\
\int\mathrm{d}^{2}\mathbf{q}\,\mathbf{q}^{2}\int\mathrm{d}^{2}\mathbf{x}_{1}\int\mathrm{d}^{2}\mathbf{x}_{2}\\\
&\times\int_{-\infty}^{\infty}\mathrm{d}x_{1}^{+}\int_{-\infty}^{\infty}\mathrm{d}x_{2}^{+}\,e^{\mathrm{i}\mathbf{q}\cdot(\mathbf{x}_{1}-\mathbf{x}_{2})}\\\
&\times\left\langle
A^{-}_{a}(x_{1}^{+},\mathbf{x}_{1})A^{-}_{a}(x_{2}^{+},\mathbf{x}_{2})\right\rangle\end{split}$
(38)

This result is consistent to the definition of the gluon transverse momentum
distribution (TMD) in ref. Dumitru  _et al._ (2019), which can be obtained
from the definitions of the TMD found in the literature, see Petreska (2018);
Mulders and Rodrigues (2001); Marquet  _et al._ (2016); Kotko  _et al._ (2015)
and references therein.

## References

  * Shuryak (2017) E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017), arXiv:1412.8393 [hep-ph] .
  * Arslandok  _et al._ (2023) M. Arslandok _et al._ , (2023), arXiv:2303.17254 [nucl-ex] .
  * Heinz and Snellings (2013) U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th] .
  * Petersen  _et al._ (2011) H. Petersen, V. Bhattacharya, S. A. Bass, and C. Greiner, Phys. Rev. C 84, 054908 (2011), arXiv:1105.0340 [nucl-th] .
  * Nie (2019) M. Nie (STAR), Nucl. Phys. A 982, 403 (2019).
  * Khachatryan  _et al._ (2015) V. Khachatryan _et al._ (CMS), Phys. Rev. C 92, 034911 (2015), arXiv:1503.01692 [nucl-ex] .
  * Aaboud  _et al._ (2018) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C 78, 142 (2018), arXiv:1709.02301 [nucl-ex] .
  * Schenke  _et al._ (2010) B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010), arXiv:1004.1408 [hep-ph] .
  * Karpenko  _et al._ (2014) I. Karpenko, P. Huovinen, and M. Bleicher, Comput. Phys. Commun. 185, 3016 (2014), arXiv:1312.4160 [nucl-th] .
  * Moreland  _et al._ (2015) J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Shen and Alzhrani (2020) C. Shen and S. Alzhrani, Phys. Rev. C 102, 014909 (2020), arXiv:2003.05852 [nucl-th] .
  * Weil  _et al._ (2016) J. Weil _et al._ (SMASH), Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Schäfer  _et al._ (2022) A. Schäfer, I. Karpenko, X.-Y. Wu, J. Hammelmann, and H. Elfner (SMASH), Eur. Phys. J. A 58, 230 (2022), arXiv:2112.08724 [hep-ph] .
  * Bozek and Broniowski (2016) P. Bozek and W. Broniowski, Phys. Lett. B 752, 206 (2016), arXiv:1506.02817 [nucl-th] .
  * Ehehalt and Cassing (1995) W. Ehehalt and W. Cassing, (1995), arXiv:hep-ph/9507274 .
  * Geiss  _et al._ (1998) J. Geiss, W. Cassing, and C. Greiner, Nucl. Phys. A 644, 107 (1998), arXiv:nucl-th/9805012 .
  * Bratkovskaya  _et al._ (2011) E. L. Bratkovskaya, W. Cassing, V. P. Konchakovski, and O. Linnyk, Nucl. Phys. A 856, 162 (2011), arXiv:1101.5793 [nucl-th] .
  * Schenke  _et al._ (2012a) B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. C 86, 034908 (2012a), arXiv:1206.6805 [hep-ph] .
  * Schenke  _et al._ (2012b) B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. Lett. 108, 252301 (2012b), arXiv:1202.6646 [nucl-th] .
  * Eskola  _et al._ (2000) K. J. Eskola, K. Kajantie, P. V. Ruuskanen, and K. Tuominen, Nucl. Phys. B 570, 379 (2000), arXiv:hep-ph/9909456 .
  * Niemi  _et al._ (2016) H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph] .
  * Drescher and Nara (2007) H. J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007), arXiv:nucl-th/0611017 .
  * Schenke and Schlichting (2016) B. Schenke and S. Schlichting, Phys. Rev. C 94, 044907 (2016), arXiv:1605.07158 [hep-ph] .
  * Schlichting and Singh (2021) S. Schlichting and P. Singh, Phys. Rev. D 103, 014003 (2021), arXiv:2010.11172 [hep-ph] .
  * Schenke  _et al._ (2022) B. Schenke, S. Schlichting, and P. Singh, Phys. Rev. D 105, 094023 (2022), arXiv:2201.08864 [nucl-th] .
  * Ipp  _et al._ (2021) A. Ipp, D. I. Müller, S. Schlichting, and P. Singh, Phys. Rev. D 104, 114040 (2021), arXiv:2109.05028 [hep-ph] .
  * Ipp  _et al._ (2022) A. Ipp, M. Leuthner, D. I. Müller, S. Schlichting, and P. Singh, EPJ Web Conf. 274, 05017 (2022), arXiv:2212.09363 [hep-ph] .
  * Gribov  _et al._ (1983) L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rept. 100, 1 (1983).
  * Kovchegov and Tuchin (2002) Y. V. Kovchegov and K. Tuchin, Phys. Rev. D 65, 074026 (2002), arXiv:hep-ph/0111362 .
  * Lappi and Mantysaari (2013) T. Lappi and H. Mantysaari, Nucl. Phys. A 908, 51 (2013), arXiv:1209.2853 [hep-ph] .
  * Dumitru and Jalilian-Marian (2002a) A. Dumitru and J. Jalilian-Marian, Phys. Lett. B 547, 15 (2002a), arXiv:hep-ph/0111357 .
  * Dumitru and Jalilian-Marian (2002b) A. Dumitru and J. Jalilian-Marian, Phys. Rev. Lett. 89, 022301 (2002b), arXiv:hep-ph/0204028 .
  * Gyulassy and McLerran (1997) M. Gyulassy and L. D. McLerran, Phys. Rev. C 56, 2219 (1997), arXiv:nucl-th/9704034 .
  * Blaizot  _et al._ (2010) J. P. Blaizot, T. Lappi, and Y. Mehtar-Tani, Nucl. Phys. A 846, 63 (2010), arXiv:1005.0955 [hep-ph] .
  * Dumitru and McLerran (2002) A. Dumitru and L. D. McLerran, Nucl. Phys. A 700, 492 (2002), arXiv:hep-ph/0105268 .
  * Buckley  _et al._ (2015) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, Eur. Phys. J. C 75, 132 (2015), arXiv:1412.7420 [hep-ph] .
  * Kowalski  _et al._ (2006) H. Kowalski, L. Motyka, and G. Watt, Phys. Rev. D 74, 074016 (2006), arXiv:hep-ph/0606272 .
  * Kowalski and Teaney (2003) H. Kowalski and D. Teaney, Phys. Rev. D 68, 114005 (2003), arXiv:hep-ph/0304189 .
  * Bartels  _et al._ (2002) J. Bartels, K. J. Golec-Biernat, and H. Kowalski, Phys. Rev. D 66, 014001 (2002), arXiv:hep-ph/0203258 .
  * Rezaeian  _et al._ (2013) A. H. Rezaeian, M. Siddikov, M. Van de Klundert, and R. Venugopalan, Phys. Rev. D 87, 034002 (2013), arXiv:1212.2974 [hep-ph] .
  * Greif  _et al._ (2017) M. Greif, C. Greiner, B. Schenke, S. Schlichting, and Z. Xu, Phys. Rev. D 96, 091504 (2017), arXiv:1708.02076 [hep-ph] .
  * Miller  _et al._ (2007) M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007), arXiv:nucl-ex/0701025 .
  * Alver and Roland (2010) B. Alver and G. Roland, Phys. Rev. C 81, 054905 (2010), [Erratum: Phys.Rev.C 82, 039903 (2010)], arXiv:1003.0194 [nucl-th] .
  * Oscar Garcia-Montero (2023) Oscar Garcia-Montero, “McDipper,” https://github.com/Non-Equilibrium-QCD/McDipper (2023), accessed: 2023-08-21.
  * Lappi and Schlichting (2018) T. Lappi and S. Schlichting, Phys. Rev. D 97, 034034 (2018), arXiv:1708.08625 [hep-ph] .
  * Schlichting and Skokov (2020) S. Schlichting and V. Skokov, Phys. Lett. B 806, 135511 (2020), arXiv:1910.12496 [hep-ph] .
  * Blaizot  _et al._ (2004a) J. P. Blaizot, F. Gelis, and R. Venugopalan, Nucl. Phys. A 743, 57 (2004a), arXiv:hep-ph/0402257 .
  * Gelis and Peshier (2002a) F. Gelis and A. Peshier, Nucl. Phys. A 697, 879 (2002a), arXiv:hep-ph/0107142 .
  * Gelis and Peshier (2002b) F. Gelis and A. Peshier, Nucl. Phys. A 707, 175 (2002b), arXiv:hep-ph/0111227 .
  * Dumitru  _et al._ (2006) A. Dumitru, A. Hayashigaki, and J. Jalilian-Marian, Nucl. Phys. A 765, 464 (2006), arXiv:hep-ph/0506308 .
  * Mäntysaari (2015) H. Mäntysaari,  _Scattering off the Color Glass Condensate_ , Ph.D. thesis, Jyvaskyla U. (2015), arXiv:1506.07313 [hep-ph] .
  * Blaizot  _et al._ (2004b) J. P. Blaizot, F. Gelis, and R. Venugopalan, Nucl. Phys. A 743, 13 (2004b), arXiv:hep-ph/0402256 .
  * Gelis and Venugopalan (2004) F. Gelis and R. Venugopalan, Phys. Rev. D 69, 014019 (2004), arXiv:hep-ph/0310090 .
  * Iancu  _et al._ (2016) E. Iancu, A. H. Mueller, and D. N. Triantafyllopoulos, JHEP 12, 041 (2016), arXiv:1608.05293 [hep-ph] .
  * Altinoluk  _et al._ (2018) T. Altinoluk, N. Armesto, A. Kovner, M. Lublinsky, and E. Petreska, JHEP 04, 063 (2018), arXiv:1802.01398 [hep-ph] .
  * Golec-Biernat and Wusthoff (1999) K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D 60, 114023 (1999), arXiv:hep-ph/9903358 .
  * d’Enterria and Loizides (2021) D. d’Enterria and C. Loizides, Ann. Rev. Nucl. Part. Sci. 71, 315 (2021), arXiv:2011.14909 [hep-ph] .
  * McLerran and Tribedy (2016) L. McLerran and P. Tribedy, Nucl. Phys. A 945, 216 (2016), arXiv:1508.03292 [hep-ph] .
  * Ma  _et al._ (2018a) Y.-Q. Ma, R. Venugopalan, K. Watanabe, and H.-F. Zhang, Phys. Rev. C 97, 014909 (2018a), arXiv:1707.07266 [hep-ph] .
  * Ma  _et al._ (2018b) Y.-Q. Ma, P. Tribedy, R. Venugopalan, and K. Watanabe, Phys. Rev. D 98, 074025 (2018b), arXiv:1803.11093 [hep-ph] .
  * Benić  _et al._ (2019) S. Benić, K. Fukushima, O. Garcia-Montero, and R. Venugopalan, Phys. Lett. B 791, 11 (2019), arXiv:1807.03806 [hep-ph] .
  * ALI (2022) (2022), arXiv:2204.10210 [nucl-ex] .
  * Borghini  _et al._ (2022) N. Borghini, M. Borrell, N. Feld, H. Roch, S. Schlichting, and C. Werthmann, (2022), arXiv:2209.01176 [hep-ph] .
  * Abelev  _et al._ (2009) B. I. Abelev _et al._ (STAR), Phys. Rev. C 79, 034909 (2009), arXiv:0808.2041 [nucl-ex] .
  * Abbas  _et al._ (2013) E. Abbas _et al._ (ALICE), Phys. Lett. B 726, 610 (2013), arXiv:1304.0347 [nucl-ex] .
  * Aamodt  _et al._ (2011) K. Aamodt _et al._ (ALICE), Phys. Rev. Lett. 106, 032301 (2011), arXiv:1012.1657 [nucl-ex] .
  * Adam  _et al._ (2017) J. Adam _et al._ (ALICE), Phys. Lett. B 772, 567 (2017), arXiv:1612.08966 [nucl-ex] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), Phys. Rev. Lett. 116, 222302 (2016), arXiv:1512.06104 [nucl-ex] .
  * Schenke  _et al._ (2020) B. Schenke, C. Shen, and P. Tribedy, Phys. Rev. C 102, 044905 (2020), arXiv:2005.14682 [nucl-th] .
  * Giacalone  _et al._ (2019) G. Giacalone, A. Mazeliauskas, and S. Schlichting, Phys. Rev. Lett. 123, 262301 (2019), arXiv:1908.02866 [hep-ph] .
  * Liyanage  _et al._ (2023) D. Liyanage, O. Sürer, M. Plumlee, S. M. Wild, and U. Heinz, (2023), arXiv:2302.14184 [nucl-th] .
  * Pang  _et al._ (2018) L.-G. Pang, H. Petersen, and X.-N. Wang, Phys. Rev. C 97, 064918 (2018), arXiv:1802.04449 [nucl-th] .
  * Pang  _et al._ (2016) L.-G. Pang, H. Petersen, G.-Y. Qin, V. Roy, and X.-N. Wang, Eur. Phys. J. A 52, 97 (2016), arXiv:1511.04131 [nucl-th] .
  * Shen and Schenke (2022) C. Shen and B. Schenke, Phys. Rev. C 105, 064905 (2022), arXiv:2203.04685 [nucl-th] .
  * Chirilli (2019) G. A. Chirilli, JHEP 01, 118 (2019), arXiv:1807.11435 [hep-ph] .
  * Dumitru  _et al._ (2019) A. Dumitru, V. Skokov, and T. Ullrich, Phys. Rev. C 99, 015204 (2019), arXiv:1809.02615 [hep-ph] .
  * Petreska (2018) E. Petreska, Int. J. Mod. Phys. E 27, 1830003 (2018), arXiv:1804.04981 [hep-ph] .
  * Mulders and Rodrigues (2001) P. J. Mulders and J. Rodrigues, Phys. Rev. D 63, 094021 (2001), arXiv:hep-ph/0009343 .
  * Marquet  _et al._ (2016) C. Marquet, E. Petreska, and C. Roiesnel, JHEP 10, 065 (2016), arXiv:1608.02577 [hep-ph] .
  * Kotko  _et al._ (2015) P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta, and A. van Hameren, JHEP 09, 106 (2015), arXiv:1503.03421 [hep-ph] .