# Study of the azimuthal asymmetry in heavy ion collisions combining initial
state momentum orientation and final state collective effects

L.S. Moriggi1<EMAIL_ADDRESS>1 Universidade Estadual do Centro-
Oeste (UNICENTRO), Campus Cedeteg, Guarapuava 85015-430, Brazil  É.S. Rocha2
<EMAIL_ADDRESS>M.V.T. Machado2<EMAIL_ADDRESS>2 High Energy Physics
Phenomenology Group, GFPAE. Institute of Physics, Federal University of Rio
Grande do Sul (UFRGS)
Caixa Postal 15051, CEP 91501-970, Porto Alegre, RS, Brazil

###### Abstract

In the present work we investigate the source of azimuthal asymmetry for
nuclear collision using a model that contemplates particles produced in the
initial hard collisions and the collective effects described by a Blast-Wave
like expansion. The latter is described by the relaxation time approximation
of the Boltzmann transport equation. The parameters regarding collective flow
and asymmetry are fitted by the experimental data from $p_{T}$ spectrum and
$v_{2}$ for PbPb and XeXe collisions at different centrality classes. As a by-
product the ratio of final elliptic flow with the initial anisotropy,
$v_{2}/\epsilon_{2}$, and the average transverse momentum are predicted.

##  I Introduction

The available relativistic heavy ion collision experiments allow us to
investigate the underlying dynamics of quarks and gluons at very high energy
density [1]. The thermalized deconfined parton system, i.e. the Quark-Gluon
Plasma (QGP) [2, 3], generated in these reactions has important signatures
like the collective flow, parton energy loss, quarkonium production
suppression and many others [4, 5]. In particular, the elliptic flow [6, 7, 8,
9] ($v_{2}$, the second harmonic coefficient of the azimuthal Fourier
decomposition of the momentum distribution) measures the non-uniformity of the
flow in all directions as viewed along the beam-line [10, 11, 12, 13]. Namely,
it characterizes the azimuthal momentum space anisotropy of particle emission
from non-central heavy-ion collisions in the plane transverse to the beam
direction. As the anisotropy is largest in the first instance of the system
evolution, the $v_{2}$-coefficient is quite sensible to the early stages of
the collisions. Thus, elliptic flow encodes the residual asymmetry of the
particle density in momentum space referring to the reaction plane subsequent
to the hadronization.

The main motivation of present work is to address the simultaneous description
of the nuclear modification factor, $R_{AA}$, and the elliptic flow in a
consistent way. It is already known that parton energy loss models that have
been tuned to describe $R_{AA}$ in identified hadron production underestimate
the elliptic flow $v_{2}$ at intermediate transverse momentum [14, 15, 16, 17,
18, 19, 20]. This has been named in literature as the $R_{AA}\otimes v_{2}$
puzzle. Both observables are nicely described by the QGP´s hydrodynamic
expansion as a strongly coupled fluid at low transverse momentum
($p_{T}\lesssim 2$ GeV) [21, 22, 23, 24, 11, 25, 26, 27, 28, 29]. At high
momentum, typically $p_{T}\gtrsim 10$ GeV, they can be correctly described in
terms of jet quenching due to hard parton energy loss in their propagation
across the hot QGP medium [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45, 46]. On the other hand, it is a challenge to describe $R_{AA}$ and
$v_{2}$ consistently in the intermediate transverse momentum region
characterized by the soft and hard physics confluence. In our case, the
parton/gluon saturation physics embedded into the theoretical approach allows
to use weak coupling methods in this interface region. Recently, studies using
a combination of the state-of-art on transport models, hadronization and
hydrodynamic evolution have been carried out. For instance, in Ref. [47] both
quark coalescence and a hadronic afterburner in the COLBT-hydro model are
implemented. That investigation combines event-by-event hydrodynamics, jet
quenching as well as hadron cascade and simultaneously describe $R_{AA}$ and
$v_{2}$ in the full range of transverse momentum. Similar approaches have been
proposed, see [48, 49], which will shed light into further investigations of
the topic. In present study the question is to what extent the elliptic flow
can be determined by the initial state effects and how these effects can be
detached from the final state ones.

Concerning identified particle spectra in transverse momenta ($p_{T}$), the
hadron production can be described within the $k_{T}$-factorization formalism
(including the primordial parton transverse momenta) and one considers that
the cold matter nuclear effects are generated in the hard interaction of the
nuclei at initial states of the corresponding collision. On the other hand,
afterwards those systems undergo a hydrodynamics evolution to freeze-out which
alters the corresponding $p_{T}$-spectra. In the context of the relaxation
time approximation (RTA) of the Boltzmann transport equation (BTE) [50] the
$p_{T}$ spectrum can be described by performing a temporal separation between
hadrons produced in initial state hard collision and those produced in the
equilibrium situation [51, 52, 53]. In Refs. [54, 55] we considered this
approach to describe the spectra of light hadrons in lead-lead (PbPb)
collisions at the Large Hadron Collider (LHC) and at the Relativistic Heavy
Ion Collider (RHIC) as well. In particular, the approach has been used to
describe particle production in small systems at RHIC [55]. The nuclear
modification factors for pAl, pAu, dAu and HeAu have been successfully
reproduced as a function of $p_{T}$ for different centralities. An important
result was that the thermal parametrization can be considerably modified by
taking into consideration the nuclear effects embedded in the gluon
distribution function of the target.

In this work we investigate the possible sources of azimuthal asymmetry for
relativistic heavy ion collisions using a theoretical model that contemplates
particles produced in the initial hard collisions and the collective effects
described by a Blast-Wave like expansion. The interface between the hard
process described within the QCD $k_{T}$-factorization formalism and the final
state collective effects is studied in details. In the hard part of spectrum a
contribution to the elliptic flow is included associated to the azimuthal
orientation of the momentum of the produced particles in relation to the
reaction plane. This effect is introduced in the QCD color dipole amplitude
from which the unintegrated gluon distribution is obtained. The charged hadron
production in proton-proton cross section has been described successfully
using this approach in Ref. [56]. In addition, azimuthal anisotropy is also
incorporated in the BTE within the RTA approximation. There, the transverse
rapidity variable $\rho$ has been modified in order to include an anisotropy
dependence in the flow. In this last context we are following closely the
Refs. [57, 58, 59, 60]. The nuclear modification factor $R_{AA}$ and $v_{2}$
are described simultaneously for PbPb and XeXe collisions at the LHC.

This paper is organized as follows. In Sec. II the theoretical model for the
charged hadron production is presented in the context of the
$k_{T}$-factorization approach including azimuthal asymmetry and nuclear
shadowing. This initial distribution is then embedded in the formalism of
hydrodynamical blast wave model (BTE-RTA approach). Expressions for both the
nuclear modification factor, $R_{AA}$, and elliptic flow, $v_{2}$ are
provided. In Sec. III the comparison of the results to the LHC data are done
and the interplay between initial and collective azimuthal asymmetries is
investigated. Finally, we summarize our conclusions in the last section IV.

##  II Theoretical framework and main predictions

We shall study the anisotropic flow, which is one of the key observables to
analyse the transport properties of the QGP. It is determined by the flow
harmonic coefficients $v_{n}$ obtained from the Fourier decomposition of the
azimuthal distribution of the produced particles in the following way [6, 7,
8, 9]:

$E\frac{d^{3}N}{d^{3}\vec{p}}=\frac{d^{2}N}{2\pi
p_{T}dp_{T}dy}\left(1+2\sum\limits_{n=1}^{\infty}v_{n}\mathrm{cos}(n\phi)\right),$
(1)

where $\phi$ is the azimuthal angle with respect to the reaction plane.
Accordingly, $E$ is the particle energy with $p_{T}$ and $y$ being its
transverse momentum value and rapidity.

The corresponding anisotropy for an identified hadron $h$ in a inclusive
nucleus-nucleus collision, $AA\rightarrow hX$, is characterized by the set of
Fourier flow coefficients $v_{n}$, defined as:

$\displaystyle
v_{n}(p_{T})\equiv\frac{\int_{0}^{2\pi}d\phi\cos(n\phi)\frac{d^{3}\sigma(AA\rightarrow
hX)}{dyd^{2}\vec{p}_{T}}}{\int_{0}^{2\pi}d\phi\frac{d^{3}\sigma(AA\rightarrow
hX)}{dyd^{2}\vec{p}_{T}}},$ (2)

with $y$ and $\vec{p}_{T}$ being the rapidity and transverse momentum vector
of the produced hadron, respectively. The integration is done over the
azimuthal angle $\phi$ in coordinate space.

Here, the RTA approximation of the BTE will be considered which is an
effective model where the collisional term has the form
$\mathrm{C}[f]=-(f-f_{eq})/t_{r}$. The Boltzmann local equilibrium
distribution, $f_{eq}$, for the distribution of particles $f$ is typified by a
freeze-out temperature $T_{eq}$ and $t_{r}$ is the relaxation time. The latter
corresponds to the time for a non-equilibrium system to reach the equilibrium
whereas $t_{f}$ is the freeze-out time parameter. Given $\mathrm{C}[f]$, the
BTE is then solved with the following initial conditions, $f(t=0)=f_{in}$ and
$f(t=t_{f})=f_{fin}$. Therefore, the final distribution $f_{fin}$ is evaluated
as a function of the ratio $t_{f}/t_{r}$ [52],

$\displaystyle
f_{fin}=f_{eq}+\left(f_{in}-f_{eq}\right)e^{-\frac{t_{f}}{t_{r}}},$ (3)

with $f_{fin}\rightarrow f_{eq}$ if the system is given enough time or
$t_{f}>t_{r}$.

Let us now focus on the determination of the initial distribution of particles
in the context of the high energy factorization approach. In this case, the
inclusive cross section for producing identified particles is given in terms
of the convolution of the unintegrated gluon distribution (UGD) for both
target and projectile and the gluon-gluon sub-process cross section. The
distribution $f_{in}$ is proportional to the production multiplicity,
$d^{3}N(AA\rightarrow hX)/dyd^{2}\vec{p}_{T}$. The UGDs containing initial
state nuclear effects can be computed within the QCD color dipole framework.
The main advantage is that the input dipole-target scattering amplitude can be
constrained from experimental data in the small-$x$ regime. The dipole-nuclei
$S$-matrix at transverse distance $r$ coordinate space for a given impact
parameter $b$ can be computed from the dipole-proton cross-section,
$\sigma_{p}(x,r)$, considering the Glauber-Mueller multiple-scattering
formalism [61, 62],

$\displaystyle S_{A}(x,r,b)=e^{-\frac{1}{2}T_{A}(b)\sigma_{p}(x,r)},$ (4)

where $T_{A}(b)$ is the nuclear thickness function and $b$ is the distance to
the nuclei center. The nuclear $S$-matrix is a function of the gluon
longitudinal momentum fraction, $x$, the transverse size of the QCD color
dipole $r$ and $b$. We are considering the proton as a homogeneous target in
the impact parameter, with transverse area $\pi R_{p}^{2}=\sigma_{0}/2$, so
that $\sigma_{p}(x,r)=2\int d^{2}b(1-S_{p}(x,r,b))=\sigma_{0}(1-S_{p}(x,r))$.
The last quantity is the dipole-proton cross section. The scattering matrix
for the proton case, $S_{p}$, can be obtained from the Fourier transform of
the UGD as a function of transverse momentum. In this work we are considering
the UGD parametrization proposed by Moriggi, Peccini, and Machado fitted to
the $p_{T}$ spectrum at different energies in proton-proton collisions [56].
The corresponding analytical dipole-proton cross section is given by,

$\displaystyle\sigma_{p}(\tau_{r})=\sigma_{0}\left(1-\frac{2(\frac{\tau_{r}}{2})^{\xi}K_{\xi}(\tau_{r})}{\Gamma(\xi)}\right),$
(5)

where $\tau_{r}=Q_{s}(x)r$ is the scaling variable as a function of $r$ and
$K_{\xi}(\tau_{r})$ the Bessel function of the second kind. The saturation
scale and the power parameter $\xi$ are fitted considering the geometric
scaling in the momentum spectrum $\tau_{Q}=Q^{2}/Q_{s}^{2}(x)$, so that
$Q^{2}=p_{T}^{2}$,

$\displaystyle\xi$ $\displaystyle=$ $\displaystyle 1+a\tau_{Q}^{b},$ (6)
$\displaystyle Q_{s}^{2}(x)$ $\displaystyle=$
$\displaystyle\left(\frac{x_{0}}{x}\right)^{0.33}.$ (7)

The parameters $a$, $b$, and $x_{0}$ are given by Moriggi et al in Ref. [56].
The main nuclear effect in Eq. (4) is to introduce shadowing and the Cronin
peak[54]. For $r\rightarrow 0$ the QCD dipole-nucleus cross section is
$\sigma_{A}(x,r)\approx A\sigma_{p}(x,r)$, which implies in a scaling of
binary collisions $N_{coll}$ for $p_{T}\gg Q_{s}(x)$ and consequently the
nuclear modification factor $R_{AA}\rightarrow 1$.

Figure 1: Multiplicity of $\pi^{\pm}$ produced in PbPb and XeXe collisions at
$\sqrt{s}=5.02$ TeV and $\sqrt{s}=5.44$ TeV, respectively, compared with
experimental data from ALICE collaboration [63, 64] for different centrality
classes. Each curve has been multiplied by a factor of $2^{i}$ for a better
displaying.

An azimuthal asymmetry can be generated if there is a dipole orientation
described by the transverse size $\vec{r}$ regarding the impact parameter
$\vec{b}$ (the angle between them is denoted by $\phi_{rb}$). Namely, the
scattering amplitude will depend upon the angle between $\vec{r}$ and
$\vec{b}$ which should introduce a dependence on angle $\phi_{p}$ between
$\vec{p}_{T}$ and $\vec{B}$ (the impact parameter of the nucleus-nucleus
collision). The evaluation of the amplitude taking into account the dipole
orientation is not an easy task (see Refs. [65, 66, 67, 68, 69] for more
detailed theoretical analyzes). The numerical solution of the Balitsky-
Kovchegov (BK) equation including dipole orientation has been studied in Ref.
[70, 71]. In what follows we consider a way to introduce such asymmetry in a
more fundamental level as proposed by Iancu and Rezaeian [65]. It was shown
that by taking into account a non homogeneous color charge sources
distribution for the target in the context of the Color Glass condensate
Condensate (CGC) formalism one produces $S_{A}=\exp[-{\cal{N}}_{2g}]$ with
${\cal{N}}_{2g}(r,b,\phi_{rb})={\cal{N}}_{0}(r,b)+{\cal{N}}_{\theta}(r,b)\cos(2\phi_{rb})$.
Expressions are provided in [65] for the auxiliary functions
${\cal{N}}_{0,\theta}$ for the case of a Gaussian distribution.

In this paper, we chose a simpler way to incorporate this effect from a
phenomenological perspective, by angular modulation, given by the substitution

$r^{2}\rightarrow r^{2}\left[1+a_{r}\cos(2\phi_{rb})\right].$ (8)

The parameter $a_{r}$ will be fitted from the experimental data of each
centrality class, and it measures the amount of asymmetry needed to describe
$v_{2}$. Such modification should give a non-zero $v_{2}$ coefficient in the
momentum spectrum of produced gluons. The cross section for inclusive gluon
production with transverse momentum $p_{T}$ in the $k_{T}$-factorization
formalism, can be described with the dipole scattering matrix $S_{A}(x,r,b)$
in the position space and the corresponding initial particle distribution is
given by:

$\displaystyle f_{in}(\phi_{p})$ $\displaystyle=$
$\displaystyle\frac{1}{p_{T}^{2}}\frac{2C_{F}}{(2\pi)^{4}\alpha_{s}}\int
d^{2}bd^{2}re^{i\vec{p}_{T}\cdot\vec{r}}\nabla_{r}^{2}S_{A}(x_{1},r,b)$ (9)
$\displaystyle\times$ $\displaystyle\nabla_{r}^{2}S_{A}(x_{2},r,b^{\prime}),$

such that $\vec{b}^{\prime}=\vec{b}-\vec{B}$ and $x_{1,2}=p_{T}e^{\pm
y}/\sqrt{s}$ are the momentum fractions carried by the gluons for each nuclei
and $\nabla_{r}^{2}$ is the Laplacian with respect to $r$ coordinate. We can
decompose the distribution $p_{T}$ into its harmonic components using the
identity

$e^{ip_{T}r\cos(\phi_{p}-\phi_{r})}=\sum_{n=-\infty}^{n=\infty}i^{n}J_{n}(p_{T}r)e^{in(\phi_{p}-\phi_{r})},$
(10)

where $J_{n}(x)$ is the Bessel function of first kind of order $n$ and
$\phi_{r}$ is the angle between $\vec{r}$ and $\vec{B}$. The second
distribution harmonic is given by

$\displaystyle\int_{0}^{2\pi}f_{in}(\phi_{p})\cos(2\phi_{p})d\phi_{p}=\frac{-2\pi}{p_{T}^{2}}\int
bdbrdrd\phi_{b}d\phi_{r}\,{\cal{I}}_{2},$
$\displaystyle{\cal{I}}_{2}=J_{2}(p_{T}r)\cos(2\phi_{r})\nabla_{r}^{2}S_{A}(r,b)\nabla_{r}^{2}S_{A}(r,b^{\prime}),$
(11)

where $\phi_{b}$ is the angle between $\vec{b}$ and $\vec{B}$. In addition,
the integral over $\phi_{p}$ is expressed as,

$\displaystyle\int_{0}^{2\pi}f_{in}(\phi_{p})d\phi_{p}=\frac{2\pi}{p_{T}^{2}}\int
bdbrdrd\phi_{b}d\phi_{r}\,{\cal{I}}_{0},$
$\displaystyle{\cal{I}}_{0}=J_{0}(p_{T}r)\nabla_{r}^{2}S_{A}(x_{1},r,b)\nabla_{r}^{2}S_{A}(x_{2},r,b^{\prime}).$
(12)

The gluon jet decay with mass $m_{j}$ models the hadron production with
transverse momentum $p_{T_{h}}$, where the hadron carries the average momentum
fraction $\left<z\right>$. The initial produced hadron spectrum can be written
as

$\displaystyle
f_{in}(p_{T_{h}},\phi_{p})=\frac{K}{\left<z\right>^{2}}f_{in}\left(p_{T}^{2}=p_{T_{h}}^{2}/\left<z\right>^{2}+m_{j}^{2},\phi_{p}\right).$
(13)

Here the parameters $K$, $\left<z\right>$, and $m_{j}$ were determined for
pion production in $pp$ collisions by Moriggi, Peccini, and Machado for
different energy collisions [56].

Using the Boltzmann-Gibbs-Blast-Wave (BGBW) model for produced hadrons in
thermal equilibrium, we could model the collective flow and subsequent
hydrodynamic expansion [72]. This phenomenological model takes into account
the main characteristics of hydrodynamic evolution. Thus we consider a
velocity profile $\rho_{0}=tanh^{-1}(\beta_{r})$ (with
$\beta_{r}=\beta_{s}(\xi)^{n}$), determined by the surface expansion velocity
$\beta_{s}$. Here, $\beta_{r}$ is the radial flow with $\xi=r/R_{0}$ being the
ratio between the radial distance of the transverse plane and the fireball
radius $R_{0}$. It is assumed a linear velocity profile, i.e. $n=1$. The
azimuthal asymmetry can be parameterized by modulation in the velocity profile
$\rho=\rho_{0}+\rho_{a}\cos(2\phi_{s})$, where $\phi_{s}$ is the azimuthal
angle in relation to reaction plane, as proposed by Huovinen et al [59]. The
quantity $\rho_{a}$ is the anisotropy parameter in the flow. We observe the
need to incorporate an extra parameter $s_{2}$ to describe the small
$p_{T_{h}}$ region, whose purpose is introduce an azimuthal variation of the
source density as described by STAR Collaboration at RHIC [73]. Hence, the
equilibrium distribution will be given by

$\displaystyle f_{eq}(\phi_{p})\propto
m_{T_{h}}\int_{0}^{R_{0}}rdr\int_{0}^{2\pi}d\phi_{s}K_{1}\left(\frac{m_{T_{h}}}{T_{eq}}\cosh(\rho(\phi_{s}))\right)$
$\displaystyle\exp\left(\frac{p_{T_{h}}}{T_{eq}}\sinh(\rho(\phi_{s}))\cos(\phi_{s}-\phi_{p})\right)\left[1+2s_{2}\cos(2\phi_{s})\right],$

where $m_{T_{h}}=\sqrt{p_{T_{h}}^{2}+m_{h}^{2}}$ is the transverse mass of
produced hadron and $K_{1}$ is the modified Bessel function of second kind.
The model above assumes that elliptic flow from collective origin is generated
by a blending of an azimuthal velocity alteration and a spatially anisotropic
freeze-out hyper-surface. In the expression above Bjorken correlation in
rapidity, $\eta=y$, is assumed where $\eta$ is the space-time rapidity.

The equilibrium temperature $T_{eq}$ (in units of GeV) and the dimensionless
parameters $\beta_{s}$, $\rho_{a}$, $s_{2}$ are dependent on the collision
impact parameter, $\vec{B}$, since the final momentum asymmetry must be
dependent on the initial geometry of the nuclear overlapping area. Such
quantities are fitted from experimental data for each centrality class
regarding $v_{2}$ and $\frac{dN_{AA}}{dp_{T_{h}}dy}$. As already mentioned,
the final $p_{T}$ spectrum can be obtained from RTA-BTE approach [52],

$f_{fin}(\phi_{p})=f_{eq}(\phi_{p})+\left[f_{in}(\phi_{p})-f_{eq}(\phi_{p})\right]e^{-t_{f}/t_{r}},$
(15)

where $t_{r}$ and $t_{f}$ are respectively the relaxation and freeze-out time.
The initial hard distribution $f_{in}$, given by the Eqs. (9) and (13),
evolves until reach the equilibrium distribution given by Eq. (II), where the
ratio $t_{f}/t_{r}$ is fitted for each centrality class. The second harmonic
of $f_{fin}(\phi_{p})$ gives the final elliptic flow coefficient $v_{2}$,
which can be written as,

$v_{2}=\frac{\int_{0}^{2\pi}f_{fin}(\phi_{p})\cos(2\phi_{p})d\phi_{p}}{\int_{0}^{2\pi}f_{fin}(\phi_{p})d\phi_{p}}.$
(16)

It is worth pointing out that not only the azimuthal asymmetry of the initial
distribution has its origin in the dipole orientation, that depends on the
impact parameter, but also generates a momentum asymmetry of Fourier Transform
given by the integration shown in Eq. (9). Meanwhile, in the equilibrium
equation given by the Eq. (II), the asymmetry can arise from the geometry of
the initial collision, even with a vanishing second harmonic of the initial
distribution. An important point is that in our approach nuclear effects
(nuclear shadowing) are already present in the nuclear UGD. It was shown in
Refs. [54, 55] that the nuclear shadowing effect changes significantly the
spectrum at low-$p_{T}$ and consequently alters the fitted parameters in the
equilibrium distribution.

The eccentricity $\varepsilon_{2}$ of the initial collision can be obtained by
the integration of Eq. (13) with respect to transverse momentum $p_{T_{h}}$.
We shall write this result as

$\varepsilon_{2}=\frac{\left<b_{y}^{2}-b_{x}^{2}\right>}{\left<b_{y}^{2}+b_{x}^{2}\right>},$
(17)

where $b_{x}$ is the component of $\vec{b}$ in the $\vec{B}$ direction.

In the present analysis, we consider only the spectrum of produced pions due
the fact that its distribution is well described up to $p_{T_{h}}\lesssim 10$
GeV for LHC energies in $pp$ collisions [56]. We need a good description of
the region $p_{T_{h}}\gtrsim 5$ GeV to find appropriate values of the yield,
Eq. (18), from a suitable parameter fitting. In this region the thermal
distribution of equilibrium, given by Eq. (II), is very small. The opposite
scenario results in an unrealistic increase of $T_{eq}$ and $\beta_{s}$ values
to fit the spectrum in the large $p_{T_{h}}$ region. After these
consideration, we can write down the hadron productions in nuclear collision
as

$\frac{dN_{AA}}{p_{T_{h}}dp_{T_{h}}dyd\phi_{p}}=e^{-t_{f}/t_{r}}f_{in}(\phi_{p})+(1-e^{-t_{f}/t_{r}})f_{eq}(\phi_{p}),$
(18)

and the nuclear modification factor will be given by,

$R_{AA}=\frac{\frac{d^{3}N_{AA}}{d\vec{p}^{3}}}{\langle
T_{AA}\rangle\frac{d^{3}\sigma_{pp}}{d\vec{p}^{3}}},$ (19)

where $\langle T_{AA}\rangle$ is the average value of the nuclear overlapping
function.

| centrality(%) | $T_{eq}\,(\mathrm{GeV})$ | $\beta_{s}$ | $\rho_{a}$ | $s_{2}$ | $t_{f}/t_{r}$ | $\chi^{2}/dof$
---|---|---|---|---|---|---|---
PbPb | 05-10 | 0.229 | 0.504 | 0.0386 | 0.0590 | 2.05 | 0.990
| 10-20 | 0.148 | 0.787 | 0.0678 | 0.0670 | 1.86 | 1.126
| 20-30 | 0.122 | 0.856 | 0.103 | 0.0835 | 1.64 | 1.389
| 30-40 | 0.116 | 0.867 | 0.128 | 0.110 | 1.38 | 1.488
| 40-50 | 0.100 | 0.900 | 0.167 | 0.132 | 1.17 | 1.683
| 50-60 | 0.090 | 0.911 | 0.202 | 0.164 | 0.89 | 1.122
| 60-70 | 0.091 | 0.910 | 0.210 | 0.254 | 0.36 | 0.766
XeXe | 05-10 | 0.190 | 0.671 | 0.0351 | 0.057 | 2.03 | 0.556
| 10-20 | 0.121 | 0.862 | 0.0782 | 0.065 | 1.65 | 0.232
| 20-30 | 0.131 | 0.829 | 0.0941 | 0.106 | 1.29 | 0.167
| 30-40 | 0.116 | 0.868 | 0.129 | 0.129 | 1.02 | 0.105
| 40-50 | 0.125 | 0.843 | 0.113 | 0.187 | 0.68 | 0.250
| 50-60 | 0.151 | 0.767 | 0.0771 | 0.324 | 0.31 | 0.098
Table 1: Kinetic freeze-out parameters in each centrality class for production
of charged pions in PbPb and XeXe collisions at $\sqrt{s}=5.02$ TeV and
$\sqrt{s}=5.44$ TeV at the LHC, respectively. Parameters
$\beta_{s},\,\rho_{a},\,s_{2}$ and $t_{f}/t_{r}$ are dimensionless. Only
central values are shown and the errors are typically of 20%.

##  III Results and discussion

Centrality | 5-10 | 10-20 | 20-30 | 30-40 | 40-50 | 50-60 | 60-70
---|---|---|---|---|---|---|---
$PbPb$ | $0.668\pm 0.068$ | $0.261\pm 0.015$ | $0.151\pm 0.086$ | $0.130\pm 0.075$ | $0.113\pm 0.073$ | $0.114\pm 0.01$ | $0.114\pm 0.015$
$XeXe$ | $0.70\pm 0.42$ | $0.190\pm 0.054$ | $0.1809\pm 0.033$ | $0.128\pm 0.031$ | $0.171\pm 0.038$ | $0.189\pm 0.057$ |
Table 2: Results of fitting the dimensionless parameter $a_{r}$ at different
centralities for PbPb and XeXe collisions. Figure 2: Nuclear modification
factors for PbPb collisions compared with data of ALICE collaboration [63] for
different centrality classes. The results for XeXe collisions are are also
presented.

We now proceed to the analysis of the available data at the LHC for the
nuclear modification factor and elliptic flow by using the theoretical
approach presented in previous section. The $p_{T_{h}}$ spectrum and $v_{2}$
coefficient were fitted in the interval $0.1<p_{T_{h}}<8$ GeV for data of
ALICE collaboration [74, 63, 75, 64], regarding pion production in PbPb and
XeXe collisions at $\sqrt{s}=5.02$ TeV and $\sqrt{s}=5.44$ TeV, respectively.
Since that initial distribution given by Eq. (9) is not capable of generating
enough asymmetry to describe $v_{2}$ in these very central collisions, we
excluded the centrality class of 0-5% in our analysis. The BGBW distribution
parameters and the ratio $t_{f}/t_{r}$ are shown in the Table 1. The Woods-
Saxon nuclear density was considered to compute the nuclear thickness function
$T_{A}(b)$, using parameters from De Vries et al [76]. It is worth mentioning
that in the present analysis we did not consider the deformation of Xe nuclei.

Thus, this set of parameters are in good agreement with the expected behavior
where $T_{eq}$ decreases with the centrality, whereas the quantities that
define azimuthal asymmetry $\rho_{a}$ and $s_{2}$ increase with the
centrality. In order to discuss the freeze-out parameters, for PbPb collisions
it is shown that $\langle\beta_{r}\rangle=\left(\frac{2}{2+n}\right)\beta_{s}$
($\langle\beta_{r}\rangle=(2/3)\beta_{s}$, for $n=1$) decreases with
centrality, reaching $\langle\beta_{r}\rangle=0.336\pm 0.024$ in $5-10\%$
central collisions, while the equilibrium temperature increases going from
$T_{eq}=(0.091\pm 0.015)$ GeV to $T_{eq}=(0.229\pm 0.084)$ GeV. As considered
in a previous analysis [54], the radial velocity $\beta_{s}$ is anti-
correlated with $T_{eq}$ and therefore diminishes with the centrality.
Although the values of $\chi^{2}/dof$ shows sizable deviation from the unit
for the PbPb case, the resulting parameters have values physically consistent.
Namely, the transverse momentum spectra data cover greater values of
$p_{T_{h}}$ and they help to restrict $T_{eq}$ and $\beta_{s}$, generating
smaller errors in the evaluation of these parameters. It should be stressed
that the value of the exponent of the expansion velocity profile, $n$, is
equal unity in the present case. The situation is different in the fits of
spectra with a blast-wave function done by ALICE Collaboration in Refs. [77,
63] where $n$ is about 0.74 in central collisions and it increases up to 2.52
in more peripheral collisions. Here, a variation of $n$ was not necessary to
reproduce the large-$p_{T}$ tail of the spectra as it is not thermal over the
full range of transverse momentum.

In general, the centrality dependence of the thermal parameters $T_{eq}$ and
$\beta_{s}$ appearing in Tab. 1 has the opposite behavior compared to the BGBW
approaches [57, 58, 59]. The main reason is that the shadowing/anti-shadowing
effect in the small-$p_{T}$ region alters the corresponding particle
distribution and modifies the parameter associated to equilibrium temperature.
Let us take for instance a Boltzmann distribution $\sim\exp(-p_{T}/T)$: the
temperature quantifies how much disperse the $p_{T}$ spectrum is. On the other
hand, the initial state nuclear effects also produce a broadening of that
distribution. In more peripheral collisions the hard (initial) distribution
starts to become more prominent and its effect is almost sufficient to
describe particle spectrum. In this sense, the temperature can be smaller in
this case compared to more central collisions.

The dimensionless parameter $a_{r}$ associated to azimuthal modulation in the
dipole cross section is more relevant in the large $p_{T}$ region. In the
$p_{T_{h}}\sim 5$ GeV region the distribution $f_{eq}$ is very small and the
final asymmetry is basically related to $f_{in}$ and therefore the parameter
$a_{r}$ can be determined separately. Further, we defined $v_{2in}$ associated
only with $f_{in}$, which was fitted within the region of $p_{T_{h}}>5$ GeV.
The results are shown in the Table 2. We can see that in more central
collisions (5-10%) the asymmetry generated by the initial distribution is very
small and $a_{r}$ should be large to fit the data. For bigger values of
centrality, $a_{r}$ is almost constant, with $a_{r}\approx 0.18$. At this
stage of the phenomenological study it is not clear how to explain the
dependence of $a_{r}$ on the centrality. Based on the work of Ref. [65], the
nuclear color dipole amplitude with dipole orientation can be written as
$N_{2g}(r,B,\phi_{rb})={\cal{N}}_{0}^{A}(r,B)\left[1+\kappa(B)\cos(2\phi_{rb}\right]$,
with $\kappa(B)\sim[T_{A}^{\prime\prime}(B)-T_{A}^{\prime}(B)/B]/T_{A}(B)$.
Here, $T_{A}^{\prime\prime}$ and $T_{A}^{\prime}$ are the second and first
derivative of the thickness function, respectively. Using by simplification of
a Gaussian thickness $T_{A}\propto\exp(-B^{2}/R_{A}^{2})$, with $R_{A}$ being
the nuclear radius, one obtains $\kappa(B)\sim(2B/R_{A}^{2})^{2}$ meaning that
the asymmetry is more intense for peripheral collisions. This feature is also
confirmed by the numerical solution of BK evolution equation when dipole
orientation is included [70, 71].

Figure 3: Elliptic flow coefficient $v_{2}(p_{T_{h}})$ for XeXe and PbPb
collisions compared to data from ALICE collaboration [75, 74]. The results for
PbPb have been multiplied by factor 2 for a better visualization.

We present in the Figure 1 the final spectrum, given by the Eq. (18), in
comparison with experimental data of ALICE collaboration for XeXe and PbPb
collisions. The spectrum of produced pions in $pp$ collisions is also shown.
It is important to note that we need a good agreement with the $pp$ comparison
case to have a suitable description of nuclear modification factor $R_{AA}$.
If this is not the case, we can still have a good fit of multiplicity but not
a suitable description of $R_{AA}$. In the Figure 2, the resulting nuclear
modification factor from the fitting is compared with data from ALICE
collaboration [63] of charged pion production. The dashed lines are the
fitting results of the spectrum for XeXe collision evaluated with the Eq.
(19). For central collisions, $R_{AA}$ is basically equal for PbPb and XeXe.
In our model, this is due mainly to the fact that in both cases the ratio
$t_{f}/t_{r}$ is essentially the same, implying that both systems evolve in a
similar time. In more peripheral collisions, the ratio $t_{f}/t_{r}$ for XeXe
is remarkably smaller, which generates less suppression of initial spectrum.

Figure 4: Detailed analysis for the $f_{in}$ and $f_{eq}$ contributions for
the final asymmetry $v_{2}$ (left plot) and nuclear modification $R_{AA}$
(right plot) for PbPb collision at centrality (50-60)%.

The differential elliptic flow $v_{2}(p_{T_{h}})$ evaluated by Eq. (16) is
shown in the Figure 3. The results are contrasted with data for XeXe and PbPb
collisions in several centralities. We can see that our model has better
agreement with $v_{2}$ for XeXe case, where $\chi^{2}$ is smaller than the
PbPb one. At a given value of $p_{T_{h}}$, the more peripheral collisions have
the largest value of elliptic flow, and $v_{2}$ decreases for more central
collisions. For all the centralities, in the small $p_{T_{h}}$ range, the
transverse momentum dependence of the flow for charged pions is approximately
linear.

The interface between the initial and final distribution is illustrated by
Figure 4, where the nuclear modification factor $R_{AA}$ (plot on the right)
and $v_{2}$ (plot on the left) as function of $p_{T_{h}}$ are shown. In these
results, the dashed lines represent the $f_{in}$ (curves in color red) and
$f_{eq}$ (curve in color black) contributions and the solid lines the
summation. To be clear, Fig. 4 shows the relative contribution to the
$p_{T}$-spectrum coming from initial and equilibrium distributions appearing
in Eq. (18). Namely, it is shown the contributions driven by
$f_{in}(p_{T})\times e^{-t_{f}/t_{r}}$ and $f_{eq}\times(1-e^{t_{f}/t_{r}})$
terms. The nuclear modification factor $R_{AA}$ is given by the sum of these
two contributions. Both $R_{AA}$ and $v_{2}$ grows up to their highest values
for $p_{T_{h}}\sim 2$ GeV, where the $f_{eq}$ contribution is maximum. After
this region, $R_{AA}$ and $v_{2}$ decreases in value, indicating a
predominance of the contribution of produced particles at initial stages of
the collision. In our model, the growth of $v_{2}$ and the Cronin peak has a
common origin as the interface between these two distributions. It also
indicates a limit to our approach, since that for $p_{T_{h}}\gtrsim 10$ GeV
$v_{2}$ decreases, while $R_{AA}$ increases. In this picture, models that
consider energy loss or color transparency can be considered. Please, see [18,
66, 47, 78] for further discussion on this topic. We can see that the
azimuthal asymmetry of the momentum generated by initial collision is
basically given by partons with large $p_{T}$, whereas the asymmetry in the
equilibrium occurs mainly in small $p_{T}$.

Figure 5: The $\phi_{p}$ dependence of the initial, final and equilibrium
distributions. Contributions for $p_{T}=2$ GeV and 50-60$\%$ centrality class
in PbPb collisions are presented. Distributions $f_{in}$ and $f_{eq}$ are
respectively labeled by the upper and lower dot-dashed curves whereas
$f_{fin}$ is given by the solid curve.

In order to illustrate the contributions arising from $f_{in}$ and $f_{eq}$ as
a function of the angle between $\vec{p}_{T}$ and impact parameter $\vec{B}$,
$\phi_{p}$, in Fig. 5 the corresponding distributions are shown. Accordingly,
as in Fig. 4 the distributions $f_{in}$ (upper dot-dashed curve), $f_{eq}$
(lower dot-dashed curve) and $f_{fin}$ (solid line curve) are presented as a
function of $\phi_{p}$ for PbPb collisions at 50-60$\%$ centrality class.
Transverse momentum $p_{T}=2$ GeV has been considered in which $f_{in}$ and
$f_{eq}$ give similar order of magnitude contributions to the final
distribution. The values of the final distribution, $f_{fin}$, are in the
middle of the equilibrium $f_{eq}$ and initial $f_{in}$ cases.

In the remainder we investigated the behaviour of $v_{2}$ as a function of
$p_{T_{h}}$ for different centralities and its dependence on the anisotropies
generated at initial states and from the generalization of the blast wave
model. Before going into the conclusions we discuss additional observables
sensitive to the ingredients considered in our approach. As examples, the
ratio $v_{2}/\varepsilon_{2}$ and the average transverse momentum of particles
are discussed.

Figure 6: Ratio $v_{2}/\varepsilon_{2}$ as function of $f_{in}/A_{T}$ for PbPb
and XeXe at different centralities. The solid line represents the linear fit
following Eq. 20. Plot in the bottom shows the ratio data/model.

In the work of Voloshin and Poskanze [79], is proposed that the ratio of final
elliptic flow with the initial anisotropy, $v_{2}/\varepsilon_{2}$, can
indicate the equilibrium level of produced system. In the limit of low density
one expects

$\frac{v_{2}}{\varepsilon_{2}}\propto\frac{1}{A_{T}}f_{in},$ (20)

where $A_{T}$ is the transverse area given by the region of nuclear
overlapping area and $f_{in}$ is the integrated initial distribution in
$p_{T}$. On the other hand, in the hydrodynamic limit, representing the
complete thermalization of the final system, we should have
$v_{2}\propto\varepsilon_{2}$. Therefore, scaling deviation given by the Eq.
(20) can indicate different mechanisms of particle production. The centrality
dependency of $v_{2}$ is considered for each centrality by integrating $v_{2}$
in $p_{T_{h}}$,

$v_{2}(B)=\frac{\int dp_{T_{h}}p_{T_{h}}v_{2}(p_{T_{h}})f_{fin}}{\int
dp_{T_{h}}p_{T_{h}}f_{fin}}.$ (21)

The elliptic flow replicates the space and momentum correlation developed
because of early stage pressure gradient. On the other hand,
$A_{T}^{-1}f_{in}$ can be associated to a measure of the transverse particle
density. Therefore, the plot $v_{2}/\varepsilon_{2}\times A_{T}^{-1}f_{in}$
could be viewed as an analogous of pressure versus energy density plot. The
Figure 6 the variation of eccentricity scaled elliptic flow with such
transverse density. It shows the ratio $v_{2}/\varepsilon_{2}$ as function of
$f_{in}/A_{T}$, with $f_{in}$ integrated over $p_{T_{h}}$ for different
centralities for both PbPb (open circles) and XeXe collisions (open squares).
The expected linear scaling (solid line in figure) is a good approximation in
more peripheral collisions. For large multiplicities, there is a substantial
deviation of the fitted line. In our model, such effect results of the
increasing in $f_{eq}$ contribution for the spectrum in more central
collisions.

Figure 7: (Left) Initial and final mean transverse momentum $\langle
p_{T_{h}}\rangle$ for different centralities. (Right) Initial and final ratio
$\frac{dN_{AA}/dy}{N_{p}}$ (right) as a function of centrality.

Finally, the initial and final mean transverse momentum is presented in Figure
7 (left plot) for different centralities. Initial $\langle p_{T_{h}}\rangle$
are represented by dot-dashed line for PbPb and dashed line for XeXe, whereas
final $\left<p_{T_{h}}\right>$ is represented by open circles for PbPb and
open squares for XeXe. For more central collisions, there is a large
increasing of $\left<p_{T}\right>$ due the temperature increasing of the
distribution $f_{eq}$. However, for more peripheral collisions, the final
transverse momentum comes closer of initial transverse momentum. The
multiplicity by wounded nucleon is also shown. While the final distribution
scales with $N_{p}$, $f_{in}$ grows faster than $N_{p}$ in more central
collisions.

We did not discuss universality aspects related to the elliptic flow in the
context of parton saturation physics used in our theoretical formalism for
$f_{in}$. Interestingly, it was demonstrated in Refs. [80, 81] that the
available data on the $v_{2}$ of charged particles for collisions at RHIC and
LHC for different centralities present a scaling law. Namely, the elliptic
flow normalized by the product of the inverse of the Knudsen number,
$\mathrm{Kn}^{-1}=Q_{s,A}L$, and eccentricity $\epsilon_{1}$ satisfy
geometrical scaling. Here, $Q_{s,A}=Q_{s,A}(x,B;A)$ is the nuclear saturation
scale and $L$ is the length related to the size of the collision area. The
normalized $v_{2}$ is a function only the $\tau_{A}=p_{T_{h}}^{2}/Q_{s,A}^{2}$
variable,

$\displaystyle\frac{v_{2}(p_{T_{h}})}{\epsilon_{1}/\mathrm{Kn}}=F(\tau_{A})=a\tau_{A}^{b},$
(22)

with the parameters $a$ and $b\sim 1/3$ fitted from data [80, 81] and
$0\leq\tau_{A}\leq 1$. It is a relative hard task to derive such a scaling law
from the geometrical scaling phenomenon. It is argued in [82] that the scaling
law can be traced back to the parton energy lost in nucleus-nucleus
collisions.

##  IV Summary and conclusions

Let us now organize the main results and draw some conclusions. In this work
we presented a phenomenological model that incorporates the main
characteristics of produced particle spectrum in heavy ions collisions. For
evaluate the inclusive gluon production, it has been considered hard initial
collision, nuclear shadowing of partonic distribution, and azimuthal asymmetry
in this model. Also the process makes use of dipole scattering matrix in the
position space, where an initial asymmetry has been added. However, collective
effects are described by BGBW distribution, leading to an increasing of the
momentum asymmetry for small $p_{T_{h}}$. One advantage of the present
parameterization is that it allow us to describe an interface between
different mechanisms of hadron production providing a good description of both
soft and hard contributions. In particular, we argue that the behavior of
$v_{2}$ and $R_{AA}$ has a common origin, i.e., the increasing in the
contribution of produced particles in thermal equilibrium and the subsequent
predominance of the distribution of produced particle at initial hard
peripheral collisions. The nuclear modification factor for pion production in
PbPb and XeXe collisions and respective elliptic flow $v_{2}(p_{T})$ have been
simultaneously described. The fitted parameters $T_{eq}$, $\beta_{s}$,
$\rho_{a}$ and $s_{2}$ are associated to the modified BGBW approach whereas
the parameter $a_{r}$ is connected to the azimuthal asymmetry in the early
stages of the collision. The quality of fit is good and the role played by
$f_{in}$ and $f_{eq}$ in description of $R_{AA}$ and $v_{2}$ has been
discussed. Finally, the ratio $v_{2}/\varepsilon_{2}$ as function of
$f{in}/A_{T}$ is computed and analyzed in terms of linear scaling between
these two quantities. The average transverse momentum, $\langle p_{T}\rangle$,
has been also investigated.

## Acknowledgments

This research was funded by the Brazilian National Council for Scientific and
Technological Development (CNPq) under the contract number 303075/2022-8.

## References

  * Arslandok  _et al._ [2023] M. Arslandok _et al._ , arXiv e-prints , arXiv:2303.17254 (2023), arXiv:2303.17254 [nucl-ex].
  * Shuryak [1978] E. V. Shuryak, Phys. Lett. B 78, 150 (1978).
  * Shuryak [1980] E. V. Shuryak, Phys. Rept. 61, 71 (1980).
  * Das  _et al._ [2022] S. K. Das _et al._ , Int. J. Mod. Phys. E 31, 12 (2022), arXiv:2208.13440 [nucl-th].
  * Ritter and Stock [2014] H. G. Ritter and R. Stock, J. Phys. G 41, 124002 (2014), arXiv:1408.4296 [nucl-ex].
  * Ollitrault [1992] J.-Y. Ollitrault, Phys. Rev. D 46, 229 (1992).
  * Voloshin [2009] S. A. Voloshin, Nucl. Phys. A 827, 377C (2009).
  * Voloshin and Zhang [1996] S. Voloshin and Y. Zhang, Z. Phys. C 70, 665 (1996).
  * Poskanzer and Voloshin [1998] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58, 1671 (1998).
  * Snellings [2011] R. Snellings, New J. Phys. 13, 055008 (2011), arXiv:1102.3010 [nucl-ex].
  * Heinz and Snellings [2013] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013), arXiv:1301.2826 [nucl-th].
  * Aggarwal [2021] M. M. Aggarwal, in  _Advances in Nuclear Physics_ , edited by R. K. Puri, J. Aichelin, S. Gautam, and R. Kumar (Springer Singapore, Singapore, 2021) pp. 161–188.
  * Ollitrault [2023] J.-Y. Ollitrault, (2023), arXiv:2308.11674 [nucl-ex].
  * Molnar and Sun [2013] D. Molnar and D. Sun, (2013), arXiv:1305.1046 [nucl-th].
  * Noronha-Hostler  _et al._ [2016] J. Noronha-Hostler, B. Betz, J. Noronha, and M. Gyulassy, Phys. Rev. Lett. 116, 252301 (2016), arXiv:1602.03788 [nucl-th].
  * Zhang and Liao [2013] X. Zhang and J. Liao, (2013), arXiv:1311.5463 [nucl-th].
  * Liao and Shuryak [2009] J. Liao and E. Shuryak, Phys. Rev. Lett. 102, 202302 (2009), arXiv:0810.4116 [nucl-th].
  * Kopeliovich  _et al._ [2012] B. Z. Kopeliovich, J. Nemchik, I. K. Potashnikova, and I. Schmidt, Phys. Rev. C 86, 054904 (2012), arXiv:1208.4951 [hep-ph].
  * Cao  _et al._ [2017] S. Cao, L.-G. Pang, T. Luo, Y. He, G.-Y. Qin, and X.-N. Wang, Nucl. Part. Phys. Proc. 289-290, 217 (2017).
  * Andres  _et al._ [2020] C. Andres, N. Armesto, H. Niemi, R. Paatelainen, and C. A. Salgado, Phys. Lett. B 803, 135318 (2020), arXiv:1902.03231 [hep-ph].
  * Romatschke and Romatschke [2007] P. Romatschke and U. Romatschke, Phys. Rev. Lett. 99, 172301 (2007), arXiv:0706.1522 [nucl-th].
  * Song  _et al._ [2011] H. Song, S. A. Bass, U. Heinz, T. Hirano, and C. Shen, Phys. Rev. Lett. 106, 192301 (2011), [Erratum: Phys.Rev.Lett. 109, 139904 (2012)], arXiv:1011.2783 [nucl-th].
  * Gale  _et al._ [2013a] C. Gale, S. Jeon, B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev. Lett. 110, 012302 (2013a), arXiv:1209.6330 [nucl-th].
  * Huovinen [2013] P. Huovinen, Int. J. Mod. Phys. E 22, 1330029 (2013), arXiv:1311.1849 [nucl-th].
  * Gale  _et al._ [2013b] C. Gale, S. Jeon, and B. Schenke, Int. J. Mod. Phys. A 28, 1340011 (2013b), arXiv:1301.5893 [nucl-th].
  * Niemi  _et al._ [2016] H. Niemi, K. J. Eskola, and R. Paatelainen, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph].
  * Bernhard  _et al._ [2016] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu, and U. Heinz, Phys. Rev. C 94, 024907 (2016), arXiv:1605.03954 [nucl-th].
  * McDonald  _et al._ [2017] S. McDonald, C. Shen, F. Fillion-Gourdeau, S. Jeon, and C. Gale, Phys. Rev. C 95, 064913 (2017), arXiv:1609.02958 [hep-ph].
  * Zhao  _et al._ [2017] W. Zhao, H.-j. Xu, and H. Song, Eur. Phys. J. C 77, 645 (2017), arXiv:1703.10792 [nucl-th].
  * Wang and Gyulassy [1992] X.-N. Wang and M. Gyulassy, Phys. Rev. Lett. 68, 1480 (1992).
  * Wang [2000] X.-N. Wang, Phys. Rev. C 61, 064910 (2000), arXiv:nucl-th/9812021.
  * Vitev and Gyulassy [2002] I. Vitev and M. Gyulassy, Phys. Rev. Lett. 89, 252301 (2002), arXiv:hep-ph/0209161.
  * Wang and Wang [2002] E. Wang and X.-N. Wang, Phys. Rev. Lett. 89, 162301 (2002), arXiv:hep-ph/0202105.
  * Eskola  _et al._ [2005] K. J. Eskola, H. Honkanen, C. A. Salgado, and U. A. Wiedemann, Nucl. Phys. A 747, 511 (2005), arXiv:hep-ph/0406319.
  * Qin  _et al._ [2008] G.-Y. Qin, J. Ruppert, C. Gale, S. Jeon, G. D. Moore, and M. G. Mustafa, Phys. Rev. Lett. 100, 072301 (2008), arXiv:0710.0605 [hep-ph].
  * Schenke  _et al._ [2009] B. Schenke, C. Gale, and S. Jeon, Phys. Rev. C 80, 054913 (2009), arXiv:0909.2037 [hep-ph].
  * Chen  _et al._ [2011] X.-F. Chen, T. Hirano, E. Wang, X.-N. Wang, and H. Zhang, Phys. Rev. C 84, 034902 (2011), arXiv:1102.5614 [nucl-th].
  * Majumder and Shen [2012] A. Majumder and C. Shen, Phys. Rev. Lett. 109, 202301 (2012), arXiv:1103.0809 [hep-ph].
  * Buzzatti and Gyulassy [2012] A. Buzzatti and M. Gyulassy, Phys. Rev. Lett. 108, 022301 (2012), arXiv:1106.3061 [hep-ph].
  * Zapp  _et al._ [2013] K. C. Zapp, F. Krauss, and U. A. Wiedemann, JHEP 03, 080 (2013), arXiv:1212.1599 [hep-ph].
  * Gyulassy  _et al._ [2001] M. Gyulassy, I. Vitev, and X. N. Wang, Phys. Rev. Lett. 86, 2537 (2001), arXiv:nucl-th/0012092.
  * Kumar  _et al._ [2017] A. Kumar, E. Bianchi, J. Elledge, A. Majumder, G.-Y. Qin, and C. Shen, Nucl. Phys. A 967, 536 (2017), arXiv:1706.07547 [nucl-th].
  * Zigic  _et al._ [2020] D. Zigic, B. Ilic, M. Djordjevic, and M. Djordjevic, Phys. Rev. C 101, 064909 (2020), arXiv:1908.11866 [hep-ph].
  * Mehtar-Tani  _et al._ [2013] Y. Mehtar-Tani, J. G. Milhano, and K. Tywoniuk, Int. J. Mod. Phys. A 28, 1340013 (2013), arXiv:1302.2579 [hep-ph].
  * Qin and Wang [2015] G.-Y. Qin and X.-N. Wang, Int. J. Mod. Phys. E 24, 1530014 (2015), arXiv:1511.00790 [hep-ph].
  * Blaizot and Mehtar-Tani [2015] J.-P. Blaizot and Y. Mehtar-Tani, Int. J. Mod. Phys. E 24, 1530012 (2015), arXiv:1503.05958 [hep-ph].
  * Zhao  _et al._ [2022] W. Zhao, W. Ke, W. Chen, T. Luo, and X.-N. Wang, Phys. Rev. Lett. 128, 022302 (2022), arXiv:2103.14657 [hep-ph].
  * Barreto  _et al._ [2022] L. Barreto, F. M. Canedo, M. G. Munhoz, J. Noronha, and J. Noronha-Hostler, (2022), arXiv:2208.02061 [nucl-th].
  * Stojku  _et al._ [2022] S. Stojku, J. Auvinen, M. Djordjevic, P. Huovinen, and M. Djordjevic, Phys. Rev. C 105, L021901 (2022), arXiv:2008.08987 [nucl-th].
  * Florkowski and Ryblewski [2016] W. Florkowski and R. Ryblewski, Phys. Rev. C 93, 064903 (2016), arXiv:1603.01704 [nucl-th].
  * Tripathy  _et al._ [2016] S. Tripathy, T. Bhattacharyya, P. Garg, P. Kumar, R. Sahoo, and J. Cleymans, The European Physical Journal A 52 (2016), 10.1140/epja/i2016-16289-4.
  * Tripathy  _et al._ [2017] S. Tripathy, A. Khuntia, S. K. Tiwari, and R. Sahoo, Eur. Phys. J. A 53, 99 (2017), arXiv:1703.02416 [nucl-th].
  * Qiao  _et al._ [2020] L. Qiao, G. Che, J. Gu, H. Zheng, and W. Zhang, Journal of Physics G: Nuclear and Particle Physics 47, 075101 (2020).
  * Moriggi  _et al._ [2021] L. S. Moriggi, G. M. Peccini, and M. V. T. Machado, Phys. Rev. D 103, 034025 (2021), arXiv:2012.05388 [hep-ph].
  * Moriggi and Machado [2022] L. Moriggi and M. V. T. Machado, MDPI Physics 4, 787 (2022), arXiv:2207.07794 [hep-ph].
  * Moriggi  _et al._ [2020] L. S. Moriggi, G. M. Peccini, and M. V. T. Machado, Phys. Rev. D 102, 034016 (2020), arXiv:2005.07760 [hep-ph].
  * Tripathy  _et al._ [2018] S. Tripathy, S. K. Tiwari, M. Younus, and R. Sahoo, Eur. Phys. J. A 54, 38 (2018), arXiv:1709.06354 [hep-ph].
  * Younus  _et al._ [2020] M. Younus, S. Tripathy, S. K. Tiwari, and R. Sahoo, Adv. High Energy Phys. 2020, 4728649 (2020), arXiv:1803.01578 [hep-ph].
  * Huovinen  _et al._ [2001] P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen, and S. A. Voloshin, Phys. Lett. B 503, 58 (2001), arXiv:hep-ph/0101136.
  * Akhil and Tiwari [2023] A. Akhil and S. K. Tiwari, (2023), arXiv:2309.06128 [hep-ph].
  * Glauber [1955] R. Glauber, Phys. Rev. 100, 242 (1955).
  * Mueller [1990] A. H. Mueller, Nucl. Phys. B 335, 115 (1990).
  * Acharya  _et al._ [2020] S. Acharya _et al._ (ALICE), Phys. Rev. C 101, 044907 (2020), arXiv:1910.07678 [nucl-ex].
  * Acharya  _et al._ [2021a] S. Acharya _et al._ (ALICE), Eur. Phys. J. C 81, 584 (2021a), arXiv:2101.03100 [nucl-ex].
  * Iancu and Rezaeian [2017] E. Iancu and A. H. Rezaeian, Phys. Rev. D 95, 094003 (2017), arXiv:1702.03943 [hep-ph].
  * Kopeliovich  _et al._ [2021] B. Z. Kopeliovich, M. Krelina, and J. Nemchik, Phys. Rev. D 103, 094027 (2021), arXiv:2102.06106 [hep-ph].
  * Altinoluk  _et al._ [2016] T. Altinoluk, N. Armesto, G. Beuf, and A. H. Rezaeian, Physics Letters B 758, 373 (2016).
  * Mäntysaari  _et al._ [2019] H. Mäntysaari, N. Mueller, and B. Schenke, Phys. Rev. D 99, 074004 (2019).
  * Salazar and Schenke [2019] F. Salazar and B. Schenke, Phys. Rev. D 100, 034007 (2019).
  * Golec-Biernat and Stasto [2003] K. J. Golec-Biernat and A. M. Stasto, Nucl. Phys. B 668, 345 (2003), arXiv:hep-ph/0306279.
  * Berger and Staśto [2011] J. Berger and A. M. Staśto, Phys. Rev. D 83, 034015 (2011).
  * Schnedermann  _et al._ [1993] E. Schnedermann, J. Sollfrank, and U. Heinz, Phys. Rev. C 48, 2462 (1993).
  * Adler  _et al._ [2001] C. Adler _et al._ (STAR), Phys. Rev. Lett. 87, 182301 (2001), arXiv:nucl-ex/0107003.
  * Acharya  _et al._ [2018] S. Acharya _et al._ (ALICE), JHEP 09, 006 (2018), arXiv:1805.04390 [nucl-ex].
  * Acharya  _et al._ [2021b] S. Acharya _et al._ (ALICE), JHEP 10, 152 (2021b), arXiv:2107.10592 [nucl-ex].
  * De Vries  _et al._ [1987] H. De Vries, C. De Jager, and C. De Vries, Atomic Data and Nuclear Data Tables 36, 495 (1987).
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044910 (2013), arXiv:1303.0737 [hep-ex].
  * Nemchik  _et al._ [2015] J. Nemchik, R. Pasechnik, and I. Potashnikova, Eur. Phys. J. C 75, 95 (2015), arXiv:1407.2781 [hep-ph].
  * Voloshin and Poskanzer [2000] S. A. Voloshin and A. M. Poskanzer, Phys. Lett. B 474, 27 (2000), arXiv:nucl-th/9906075.
  * Andrés  _et al._ [2015] C. Andrés, J. Dias de Deus, A. Moscoso, C. Pajares, and C. A. Salgado, Phys. Rev. C 92, 034901 (2015), arXiv:1405.2177 [hep-ph].
  * Andrés Casas [2017] C. Andrés Casas,  _Phenomenological studies of initial state effects and jet quenching in High-Energy Nuclear Collisions at LHC_ , Ph.D. thesis, Santiago de Compostela U., IGFAE (2017).
  * Andrés  _et al._ [2017] C. Andrés, M. Braun, and C. Pajares, Eur. Phys. J. A 53, 41 (2017), arXiv:1609.03927 [hep-ph].