# Transverse mass scaling of dilepton radiation off a quark-gluon plasma

Maurice Coquet<EMAIL_ADDRESS>Xiaojian Du  Jean-Yves Ollitrault
Sören Schlichting  Michael Winn  Université Paris-Saclay, Centre d’Etudes de
Saclay (CEA), IRFU, Département de Physique Nucléaire (DPhN), Saclay, France
Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany
Université Paris Saclay, CNRS, CEA, Institut de physique théorique, 91191 Gif-
sur-Yvette, France

###### Abstract

The spectrum of dileptons produced by the quark-gluon plasma in an
ultrarelativistic nucleus-nucleus collision depends only, to a good
approximation, on the transverse mass $M_{t}$ of the dilepton. This scaling is
exact as long as transverse flow is negligible, and the system is in local
thermal equilibrium. We implement a state-of-the-art modelization of kinetic
and chemical equilibration in the early stages of the evolution to study the
modifications of the spectrum. Violations of $M_{t}$ scaling resulting from
these effects are evaluated as a function of the shear viscosity to entropy
ratio ($\eta/s$) that controls the equilibration time. We determine the
dependence of the spectrum on system size, centrality, rapidity, and collision
energy. We show that the quark-gluon plasma produces more dileptons than the
Drell-Yan process up to invariant masses of order $M\sim 4$ GeV at LHC
energies. Due to different kinematics, for a given $M_{t}$, the dependence of
the dilepton yield on $M$ is opposite for the two processes, so that
experiment alone can in principle determine which process dominates.

††journal: NPA

##  1 Introduction

An ultrarelativistic nucleus-nucleus collision produces strongly-interacting
matter which rapidly thermalizes into a hot quark-gluon plasma (QGP) [1]. This
plasma expands freely into the vacuum and eventually cools down into a gas of
hadrons. Electron-positron and muon-antimuon pairs, referred to as dileptons,
are created throughout the history of the QGP by quark-antiquark annihilation.
Once produced, they reach the detector without any further interaction, so
that they probe the entire space-time dynamics, including the early stages of
the collision. In particular, they carry unique information about the
thermalization of the QGP [2]. Dileptons produced by the QGP can be separated
from those produced later on in the hadronic phase using the invariant mass,
$M$, of the pair as a selection criterion. Specifically, the contribution of
the QGP dominates for $M\gtrsim 1.2$ GeV [3, 4].

This QGP dilepton production can be studied not only as a function of the
invariant mass $M$, but also as a function of the momentum of the dilepton. It
has long been known [5] that at a given rapidity $y$, the spectrum
$dN^{l^{+}l^{-}}/d^{4}K$, where $K$ is the 4-momentum of the dilepton, should
depend only on the transverse mass $M_{t}\equiv\sqrt{M^{2}+k_{t}^{2}}$, where
$k_{t}$ is the transverse momentum of the dilepton, provided that the QGP is
in local equilibrium, and that the production occurs early enough that
transverse flow is negligible. We improve over this seminal work by
implementing a state-of-the-art treatment of pre-equilibrium dynamics [6],
still neglecting transverse flow. We study the effect of kinetic and chemical
equilibration on the $M_{t}$ spectrum, and the deviations from $M_{t}$ scaling
that they induce.

In Sec. 2, we explain when and why $M_{t}$ scaling is expected, and we
rederive the expression obtained by McLerran and Toimela for the $M_{t}$
distribution [5]. In Sec. 3, we discuss qualitatively, using dimensional
analysis, the effects of pre-equilibrium dynamics. Quantitative results are
presented in Sec. 4. The setup of our calculation is the same as in our
previous work [6], in which we only calculated the mass spectrum of dileptons,
integrated over momentum. We first briefly recall this setup, and we show our
results for the QGP dilepton spectrum in Pb+Pb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV. We show how the shear viscosity over entropy ratio $\eta/s$ at
early times, which governs the thermalization of the QGP, can be extracted
from the spectrum. In Sec. 5, we evaluate the dilepton spectrum resulting from
the Drell-Yan process, i.e, the annihilation of quarks and antiquarks
belonging to incoming nuclei [7], and we compare this background with the QGP
spectrum.

##  2 Transverse mass scaling

It has long been known that in a hadronic gas in thermal equilibrium, the
transverse momentum spectra of all identified hadrons fall on the same curve
when plotted as a function of their transverse mass $M_{t}$ [8, 9, 10, 11].
This follows from the fact that the phase-space distribution of particles
within the gas is a Boltzmann factor $dN/d^{3}pd^{3}x=\exp(-E/T)$, where $E$
is the energy of the particle, $T$ the temperature, and we have neglected the
small effects of quantum statistics. Writing $E=M_{t}\cosh y$, where $y$ is
the rapidity, and integrating over the rapidity, the resulting distribution
only depends on $M_{t}$.

This argument does not immediately apply to dileptons because they cease to
interact as soon as they are produced. Therefore, dileptons produced by an
equilibrated QGP are not themselves in thermal equilibrium. However,
transverse mass scaling still holds, for reasons which we now explain.

Dilepton production occurs through the production of a virtual photon, which
then decays into a lepton-antilepton pair. What one calls the 4-momentum of
the dilepton, $K$, is actually the 4-momentum of the virtual photon. To
leading order in perturbation theory, the production of the virtual photon
occurs through a $2\to 1$ process: The annihilation of a quark, with
4-momentum $P_{1}$, and an antiquark, with 4-momentum $P_{2}$, into a virtual
photon with 4-momentum $K=P_{1}+P_{2}$. We denote the phase space
distributions of quarks and antiquarks by $f_{q}(x,{\bf p}_{1})$ and
$f_{\bar{q}}(x,{\bf p}_{2})$, respectively, where $x$ denotes space-time
coordinates. The rate of dilepton production is obtained by integrating the
transition rate over all possible values of ${\bf p}_{1}$ and ${\bf p}_{2}$,
taking energy-momentum conservation into account:

$\frac{dN^{l^{+}l^{-}}}{d^{4}xd^{4}K}=\int{\frac{d^{3}p_{1}}{(2\pi)^{3}2p_{1}}\frac{d^{3}p_{2}}{(2\pi)^{3}2p_{2}}f_{q}(x,{\bf
p}_{1})f_{\bar{q}}(x,{\bf p}_{2})|{\cal
A}|^{2}(2\pi)^{4}\delta^{(4)}(P_{1}+P_{2}-K)},$ (1)

where ${\cal A}$ is the Lorentz-invariant amplitude of the process, which is
not modified by the thermal medium to leading order in perturbation theory. We
have neglected quark masses so that the invariant phase-space element is just
$d^{3}p_{i}/p_{i}$, with $i=1,2$, and $p_{i}\equiv|{\bf p}_{i}|$. In a
baryonless thermal fluid at rest, $f_{q}(x,{\bf p})=f_{\bar{q}}(x,{\bf
p})=\exp(-p/T(x))$ (assuming Boltzmann statistics), where $T(x)$ is the
temperature at space-time point $x$. Conservation of energy then implies
$f_{q}(x,{\bf p}_{1})f_{\bar{q}}(x,{\bf p}_{2})=\exp(-k_{0}/T(x))$, which can
be moved outside the integral:

$\frac{dN^{l^{+}l^{-}}}{d^{4}xd^{4}K}=e^{-k_{0}/T(x)}\int{\frac{d^{3}p_{1}}{(2\pi)^{3}2p_{1}}\frac{d^{3}p_{2}}{(2\pi)^{3}2p_{2}}|{\cal
A}|^{2}(2\pi)^{4}\delta^{(4)}(P_{1}+P_{2}-K)}.$ (2)

The pre-factor on the right-hand side is the Boltzmann factor corresponding to
the dilepton, and the rest is a non-trivial kinematic integral involving the
scattering amplitude, which could in principle depend on the four-momentum
$K$. We show that it is in fact independent of $K$.

First, one notes that the integrand is a Lorentz scalar. If one neglects quark
and lepton masses, the only Lorentz-invariant scale is the invariant mass
$M=(K^{\mu}K_{\mu})^{1/2}$ of the dilepton, hence the integral can only depend
on $M$. This dependence can be obtained through dimensional analysis. The
left-hand side of Eq. (2) is dimensionless in natural units $\hbar=c=1$,
therefore, the integral in the right-hand side is also dimensionless. This
implies that it is actually independent of $M$.

For a fluid at rest, we have demonstrated that

$\frac{dN^{l^{+}l^{-}}}{d^{4}xd^{4}K}=C\exp\left(-\frac{k_{0}}{T(x)}\right),$
(3)

where $C$ is a dimensionless constant. Now, since
$dN^{l^{+}l^{-}}/d^{4}xd^{4}K$ is a Lorentz scalar, the result for a moving
fluid is identical, provided that one replaces $k_{0}$ with the dilepton
energy in the rest frame of the fluid.

The dilepton spectrum is obtained by integrating the production rate over the
space-time coordinates $x^{\mu}$. We assume that the QGP is invariant under
longitudinal boosts [12]. Then, its space-time volume can be rewritten as
$d^{4}x=d^{2}{\bf x}_{\perp}\tau d\tau dy_{f}$, where ${\bf x}_{\perp}$ is the
transverse position, $\tau\equiv\sqrt{t^{2}-z^{2}}$ the proper time and
$y_{f}={\rm artanh}(z/t)$ the fluid rapidity. Finally, we neglect transverse
flow. Then, the dilepton energy in the fluid rest frame is
$M_{t}\cosh(y-y_{f})$. The dilepton spectrum is:

$\frac{dN^{l^{+}l^{-}}}{d^{4}K}=C\int d{\bf x}_{\perp}\int_{0}^{\infty}\tau
d\tau\int_{-\infty}^{+\infty}dy_{f}\exp\left(-\frac{M_{t}\cosh(y-y_{f})}{T({\bf
x}_{\perp},\tau)}\right).$ (4)

It is independent of $y$, as a consequence of the assumed longitudinal boost
invariance. It depends on $k_{t}$ and $M$ only through $M_{t}$, which is the
property of transverse mass scaling. As can be seen from the above argument,
$M_{t}$ scaling is a robust property of the leading-order dilepton production,
which is independent of the detailed space-time dynamics. As long as the
system is longitudinally boost invariant and the transverse expansion can be
neglected it simply follows from symmetry and dimensional analysis. One
expects it to be broken by next-to-leading order perturbative corrections,
which are smaller than the leading-order contribution [13], and by non-
perturbative dynamics [14, 15]. For leading-order production, an explicit
calculation gives the expression of the proportionality constant $C$ in Eq.
(3) [6, 16]:

$C=\frac{N_{c}\alpha^{2}}{12\pi^{4}}\sum_{f}q_{f}^{2},$ (5)

where $N_{c}=3$ is the number of quark colors, $\alpha$ is the fine structure
constant, $\sum_{f}$ denotes the summation over quark flavors, $q_{f}$ is the
quark electric charge, $\frac{2}{3}$ for $u$ and $-\frac{1}{3}$ for $d$ and
$s$.

We now derive the explicit form of the $M_{t}$ spectrum assuming that the
equation of state of the QGP is conformal [17], which is approximately true at
high temperatures, and implies that $\tau T({\bf x}_{\perp},\tau)^{3}$ is
independent of $\tau$ [12]. Then, the integral over $\tau$ in Eq. (4) can
easily be done analytically using the change of variables $\tau=x^{3}$. The
integral over the fluid rapidity $y_{f}$ can also be done analytically. We
further simplify the problem (although this simplification is not essential)
by assuming that the temperature profile is uniform within a transverse area
$A_{\perp}$, that is, $T({\bf x}_{\perp},\tau)$ is independent of ${\bf
x}_{\perp}$. We obtain:

$\left(\frac{dN^{l^{+}l^{-}}}{d^{4}K}\right)_{\rm
ideal}=\frac{32N_{c}\alpha^{2}\sum_{f}q_{f}^{2}}{\pi^{4}}\frac{A_{\perp}(\tau
T^{3})^{2}}{M_{t}^{6}},$ (6)

which is the McLerran-Toimela $M_{t}^{-6}$ spectrum [5]. The subscript ideal
refers to the fact that local equilibrium holds at all times, which in turn
implies that the expansion is ruled by ideal hydrodynamics.

Note that the $M_{t}$ spectrum is a power law. This seems to contradict the
expectation from lower energies that the $M_{t}$ spectrum should be
exponential, with the inverse slope measuring the effective temperature probed
by dileptons [18, 19, 20, 21]. The contradiction is only apparent. Larger
values of $M_{t}$ are produced at earlier times, when the temperature is
higher. Therefore, the effective temperature depends on $M_{t}$:

$T_{\rm
eff}(M_{t})\equiv-\left[\frac{d}{dM_{t}}\ln\left(\frac{dN^{l^{+}l^{-}}}{d^{4}K}\right)\right]^{-1}=\frac{M_{t}}{6}.$
(7)

It is the integration over time which converts the exponential spectrum into a
power law [5].

The constant $\tau T^{3}$ in Eq. (6) is proportional to the charged
multiplicity per unit rapidity [22], and inversely proportional to $A_{\perp}$
[6]. Our estimates for Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV near
mid-rapidity in the $0-5\%$ centrality range are:

$\displaystyle A_{\perp}$ $\displaystyle=$ $\displaystyle
104\leavevmode\nobreak\ {\rm fm}^{2}=2670\leavevmode\nobreak\ {\rm GeV}^{-2}$
(8) $\displaystyle\tau T^{3}$ $\displaystyle=$ $\displaystyle
10.3\leavevmode\nobreak\ {\rm fm}^{-2}=0.40\leavevmode\nobreak\ {\rm
GeV}^{2}.$ (9)

These values will be used in Sec. 4, where we show that our numerical results
smoothly converge to Eq. (6) when the viscosity over entropy ratio $\eta/s$,
which controls the deviations from equilibrium, goes to zero.

Finally, the spectrum (6) can be integrated over $M_{t}$ for fixed $M$. The
resulting invariant mass spectrum is proportional to $M^{-3}$ [5]:

$\displaystyle\left(\frac{dN^{l^{+}l^{-}}}{dMdy}\right)_{\rm ideal}$
$\displaystyle=$ $\displaystyle 2\pi
M\int_{M}^{\infty}M_{t}dM_{t}\left(\frac{dN^{l^{+}l^{-}}}{d^{4}K}\right)_{\rm
ideal}$ (10) $\displaystyle=$
$\displaystyle\frac{16N_{c}\alpha^{2}\sum_{f}q_{f}^{2}}{\pi^{3}}\frac{A_{\perp}(\tau
T^{3})^{2}}{M^{3}}.$ (11)

##  3 Pre-equilibrium dynamics: qualitative discussion

We now discuss qualitatively the effects of pre-equilibrium dynamics. Several
effects must be taken into account:

  * 1.

The time dependence of the temperature is modified due to the anisotropy of
the momentum distribution of quarks and gluons.

  * 2.

The quark momentum distribution entering the production rate (1) is
anisotropic.

  * 3.

Quarks are underpopulated relative to gluons.

The first two effects correspond to kinetic equilibration, while the third
corresponds to chemical equilibration. At the end of this section, we also
discuss the qualitative effects of transverse flow, which is not included in
our numerical results.

For this qualitative discussion, we model the departure from local thermal
equilibrium by replacing ideal hydrodynamics with Navier-Stokes viscous
hydrodynamics. The relative order of magnitude of viscous corrections can then
be derived on the basis of dimensional analysis. The largest term in the
energy-momentum tensor of an ideal fluid is proportional to $\epsilon+P=Ts$
[23]. The correction involving the shear viscosity $\eta$ is a gradient [17].
In the early stages of the collision, due to the fast longitudinal expansion,
the largest gradient is the time derivative, which is of order $1/\tau$ for
dimensional reasons. Hence, the viscous term is of order $\eta/\tau$, while
the ideal term is of order $Ts$. The ratio of the two is the inverse Reynolds
number, which depends on $\tau$:

${Re}^{-1}(\tau)\equiv\frac{\eta}{s}\frac{1}{\tau T(\tau)}.$ (12)

Now, dileptons with a transverse mass $M_{t}$ are dominantly produced when the
temperature $T(\tau)$ is of the order of $M_{t}$. This occurs at a time $\tau$
proportional to $(\tau T^{3})/M_{t}^{3}$, where we recall that $\tau T^{3}$ is
approximately constant. Inserting these orders of magnitude of $T(\tau)$ and
$\tau$ into Eq. (12), we obtain the order of magnitude of the relevant
Reynolds number, which now depends on $M_{t}$:

${Re}^{-1}(M_{t})\equiv\frac{\eta}{s}\frac{M_{t}^{2}}{\tau T^{3}}.$ (13)

The relative correction to the dilepton yield due to pre-equilibrium dynamics
is of the order of ${Re}^{-1}(M_{t})$. As we shall see in Sec. 4, this
dimensional reasoning is confirmed by numerical calculations.

We now discuss, still at the qualitative level, the breaking of $M_{t}$
scaling which is expected when the plasma is not in local equilibrium. The
first effect is that the quark momentum distribution is no longer isotropic.
Due to the fast longitudinal expansion, longitudinal momenta in the comoving
frame are typically much smaller than tranverse momenta [2]. In order to
evaluate the qualitative effect of this momentum anisotropy on dilepton
emission, we consider the extreme case where quark distributions are purely
transverse, still assuming, for simplicity, that they are Boltzmann
distributions:

$f_{q,\bar{q}}(x,{\bf p})\propto\delta(p_{z})\exp(-p/T(x)),$ (14)

where the proportionality factor has dimension of energy. Inserting this
expression into Eq. (1), and dropping the constant proportionality factors,
one obtains

$\frac{dN^{l^{+}l^{-}}}{d^{4}xd^{4}K}\propto
e^{-k^{0}/T(x)}\delta(k_{z})\int{\frac{d^{2}p_{1}}{p_{1}}\frac{d^{2}p_{2}}{p_{2}}|{\cal
A}|^{2}\delta^{(3)}(P_{1}+P_{2}-K)},$ (15)

where the integration only runs over the transverse momenta, and we have
factored out $\delta(k_{z})$, so that the Dirac constraint inside the integral
is now in 2+1 dimensions (transverse momentum and energy). The integrand is
invariant under Lorentz transformations in 2+1 dimensions, which again implies
that the integral can only depend on the invariant mass of the dilepton, $M$.
Dimensional analysis of Eq. (1) shows that the scattering amplitude ${\cal A}$
is dimensionless, such that the integral in Eq. (15) has the mass dimension
$-1$, and is therefore proportional to 1/M. Explicitly, the integral evaluates
to $2\pi|{\cal A}|^{2}/M$. The factor $\delta(k_{z})$ can be rewritten as
$(1/M_{t})\delta(y)$, where $y$ is the rapidity of the dilepton. In a
reference frame where the fluid has rapidity $y_{f}$, this becomes
$(1/M_{t})\delta(y-y_{f})$. Finally, since $k_{z}=0$, the energy of the
dilepton coincides with its transverse mass, and one obtains:

$\frac{dN^{l^{+}l^{-}}}{d^{4}xd^{4}K}\propto\exp\left(-\frac{M_{t}}{T(x)}\right)\frac{1}{M_{t}}\delta(y-y_{f})\frac{1}{M}.$
(16)

The dilepton spectrum is obtained by integrating over the space-time history
of the fluid, as in Eq. (4). One obtains

$\frac{dN^{l^{+}l^{-}}}{d^{4}K}\propto\frac{1}{M_{t}M}\int d{\bf
x}_{\perp}\int_{0}^{\infty}\tau d\tau\exp\left(-\frac{M_{t}}{T({\bf
x}_{\perp},\tau)}\right).$ (17)

$M_{t}$ scaling is broken by the factor $1/M$, which results from the reduced
dimensionality of the phase-space integral in Eq. (15). For a given $M_{t}$,
$dN^{l^{+}l^{-}}/d^{4}K$ is smaller for larger values of $M$. As we shall see
in Sec. 4, this hierarchy is borne out by numerical calculations.

The other effect of pre-equilibrium dynamics is that the relative abundances
of quarks and antiquarks are smaller than thermal abundances in the early
stages of the collision. The collision between the incoming nuclei creates
mostly gluons [24]. Quark-antiquark pairs are then gradually produced by
collisions between gluons [25, 33]. Since dileptons are produced by quark-
antiquark annihilation, quark suppression implies a suppression of dilepton
production [26].

Studies of QCD thermalization in kinetic theory [33] suggest that kinetic and
chemical equilibration are governed by a single time scale, $\tilde{w}=\tau
T/(4\pi\eta/s)$ , which can be seen as the age of the system $\tau$ in units
of the equilibrium relaxation time $\sim(4\pi\eta/s)/T$. Hence, within our
effective description the shear viscosity over entropy ratio $\eta/s$ controls
not only the kinetic equilibration but also the approach to chemical
equilibration. Therefore, the suppression of dilepton production due to quark
suppression follows the above dimensional analysis, and should scale like
$Re^{-1}(M_{t})$ in Eq. (13). It depends only on $M_{t}$, so that quark
suppression by itself should not break $M_{t}$ scaling. However, it occurs in
the early stages where the largest breaking of $M_{t}$ scaling is expected.
Therefore, one expects the breaking to be milder when quark suppression is
included. We will check this in Sec. 4.

The last effect which breaks $M_{t}$ scaling is transverse flow. The
transverse fluid velocity is proportional to $\tau$ at early times [23, 27,
28]. Therefore, transverse flow becomes more and more important as time goes
by. Since the time of dilepton production decreases with $M_{t}$ like
$M_{t}^{-3}$, one expects that effects of transverse flow become negligible if
$M_{t}$ is large enough [19]. The qualitative effect of transverse flow on
$M_{t}$ scaling is the following [29]: For a given $M_{t}$, the transverse
boost enhances dilepton production for larger $k_{t}$ or, equivalently,
smaller values of $M$. Note that this effect is qualitatively similar to the
effect of pre-equilibrium dynamics discussed above. It has been seen
experimentally by the NA60 Collaboration [30]. We do not model transverse
flow, therefore, we cannot assess quantitatively the breaking of $M_{t}$
scaling resulting from it. However, we will estimate in Sec. 4 the range of
$M_{t}$ for which transverse flow is likely to be important.

##  4 Pre-equilibrium dynamics: quantitative results

We now present quantitative estimates of QGP dilepton production in Pb+Pb
collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. The calculation is essentially the
same as in Ref. [6], therefore we only recall the essential steps. Compared to
the calculation of Sec. 2, the main difference lies in the quark momentum
distribution $f_{q,\bar{q}}({\bf p})$ in Eq. (1). The momentum anisotropy is
modeled by carrying out the following replacement [2] in the Boltzmann111The
only difference with Ref. [6] is that we use Boltzmann distributions instead
of Fermi-Dirac distributions in Eq. (1). The advantage of this simplification
is that our results converge to the McLerran-Toimela spectrum (6) as
$\eta/s\to 0$, which is a useful benchmark. We have checked that the dilepton
yields decrease only by a few percent if one uses Fermi-Dirac instead of
Boltzmann. distribution:

$|{\bf p}|\rightarrow\sqrt{p_{t}^{2}+\xi^{2}p_{z}^{2}},$ (18)

where $\xi>1$ is the anisotropy parameter. Note that this ansatz implicitly
assumes that the tail of the momentum distribution is exponential. Therefore,
our modelization does not take into account the possibility that the falloff
at large momentum is slower than exponential, corresponding to the presence of
an increased number of high-momentum partons in the early stages, usually
referred to as “minijets” [31]. Note that, on the other hand, some choices of
initial conditions inspired by the color glass picture imply a falloff at
large momentum which is faster than exponential [32].

We take quark suppression into account by multiplying the Boltzmann
distribution by a global “quark suppression” factor $q_{s}$, which is smaller
than unity. The anisotropy parameter $\xi$ and the quark suppression factor
are computed as a function of time using QCD kinetic theory [33]. More
precisely, we use QCD kinetic theory to evaluate the pressure anisotropy and
the fraction of energy density carried by quarks. We then match the anisotropy
parameter and the quark suppression factor to these results. Note that we
could have used QCD kinetic theory to calculate directly the quark
distribution. The reason why we choose not to do so is the following. There is
by now strong theoretical evidence that the evolution of the pressure
anisotropy is fairly universal [34, 35] and does not depend on the details of
the microscopic dynamics [36]. We therefore believe that the results we obtain
in this indirect way, through a minimal distortion of the Boltzmann
distribution, provide an efficient, transparent and robust way to investigate
the dilepton spectrum. We recall that both out-of-equilibrium parameters,
namely the anisotropy parameter $\xi$ and quark suppression $q_{s}$ depend on
the scaling variable $\tilde{w}=\tau T/(4\pi\eta/s)$ as described in detail in
[6]. Hence the only free parameter in the calculation is the viscosity over
entropy ratio $\eta/s$, which is assumed to be constant and controls the
proper-time dependence of quark production and kinetic equilibration, e.g. for
larger values of $\eta/s$, the quark distribution will approach isotropy and
chemical equilibrium later on ($\tau\sim 3$fm/c for $\eta/s=0.32$) than for
lower values ($\tau\sim 1$fm/c for $\eta/s=0.16$).

The other difference with the calculation of Sec. 2 is that the temperature
decreases more slowly than $\tau^{-1/3}$ at early times, due to the smaller
longitudinal pressure. One recovers the $\tau^{-1/3}$ dependence at late
times. The temperature is determined by matching the value of $\tau T^{3}$ at
late times to the observed multiplicity. That is, the value of $\tau T^{3}$ at
late times is the same as in the ideal case.

Our numerical calculations are carried out for Pb+Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV in the 0-5% centrality range, and the
corresponding normalizations are given by Eq. (8). We have carried out
calculations with and without quark suppression, for four different values of
$\eta/s$: $0.04$, $0.08$ (not shown), $0.16$, and $0.32$. The expected value
for QCD, in the temperature range spanned by the early evolution, typically
lies between $0.16$ and $0.32$ [37]. Smaller values $0.04$ and
$0.08\simeq\frac{1}{4\pi}$ [38] have also been implemented, in order to check
numerically that our results converge smoothly to McLerran-Toimela spectrum
(6) in the limit $\eta/s\to 0$.

Figure 1: (Color online) Ratio of the dilepton spectrum, evaluated with
$\eta/s=0.04$, to the McLerran-Toimela spectrum (6), as a function of the
transverse mass of the dilepton $M_{t}$, for several values of the invariant
mass $M$. Calculations are for $0-5\%$ central Pb+Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV. Closed symbols correspond to our numerical
calculations with all pre-equilibrium effects taken into account. Open symbols
correspond to calculations where chemical equilibrium is assumed at all times,
i.e., quark suppression is not taken into account. The thin lines are the
global fit of our results using Eq. (19) (see text).

The results with $\eta/s=0.04$ are displayed in Fig. 1 for five equally-spaced
values of the invariant mass $M$.222The lowest value $M=1$ GeV is shown only
for the sake of illustration, as hadronic production, which we do not
consider, is significant for $M<1.2$ GeV. We have divided the spectrum
calculated numerically with the analytic result for ideal hydrodynamics,
Eq.(6). The ratio is smaller than unity, which confirms the expectation that
pre-equilibrium effects inhibit dilepton emission. It is naturally smaller
when quark suppression is taken into account, as can be seen by comparing
closed symbols with open symbols. The ratio is very close to unity for small
$M_{t}$. The deviation from unity increases as a function of $M_{t}$ as
expected from the larger deviations from equilibrium at the time of production
in Eq. (13). For a fixed $M_{t}$, the yield decreases as the invariant mass
$M$ increases, in line with the expectation from Eq. (17).

The dependence of the dilepton yield on the parameters $\eta/s$, $M_{t}$ and
$M$ is well captured by the following formula:

$\frac{dN^{l^{+}l^{-}}}{d^{4}K}\simeq\left(\frac{dN^{l^{+}l^{-}}}{d^{4}K}\right)_{\rm
ideal}\frac{\left(1+a\frac{\displaystyle\eta}{\displaystyle
s}M_{t}^{2}/n\right)^{-n}}{\sqrt{1+b\frac{\displaystyle\eta}{\displaystyle
s}M^{2}}}$ (19)

where the first term in the right-hand side is the McLerran-Toimela spectrum
(6), and $a$, $b$, $n$ are adjustable parameters. The parameter $a$ quantifies
the dependence of the suppression on $M_{t}$, according to Eq. (13). The
parameter $b$ quantifies the breaking of $M_{t}$ scaling due to pre-
equilibrium dynamics. The functional form (19) guarantees that the deviations
from ideal hydrodynamics are linear in $\eta/s$ in the limit
$\eta/s\rightarrow 0$, as implied by the dimensional analysis in Sec. 3. The
fact that this functional form gives a satisfactory fit of our numerical
results for a wide range of values of $\eta/s$ is a clear indication that our
dilepton spectrum converges smoothly to the McLerran-Toimela spectrum in the
limit $\eta/s\to 0$. The parameter $n$ specifies the dependence of pre-
equilibrium effects on the Reynolds number, in the non-linear regime where
these effects are large. Note that the mass spectrum $dN^{l^{+}l^{-}}/dM$
obtained by integrating Eq. (19) over $M_{t}$ is a much better approximation
of our numerical results than Eq. (28) of [6].

The parametrization (19) implies that the spectrum is proportional to $1/M$ in
the limit of large $\eta/s$, in agreement with Eq. (17). However, the
calculation leading to Eq. (17) is not strictly equivalent to our quantitative
calculation for the following reason. The hypothesis leading to Eq. (17) is
Eq. (14), namely, that the momentum distribution has zero width in $p_{z}$ and
is exponential in $p_{t}$. In the quantitative calculation, the width of the
$p_{z}$ distribution also goes to 0 in the limit of large $\eta/s$, but it is
not strictly exponential in $p_{t}$ (it is a Bessel function
$K_{0}(p_{t}/T)$). Nevertheless, the dependence of the dilepton spectrum on
$M$ ends up being essentially the same in both cases.

Figure 2: (Color online) Full symbols: Expected dilepton invariant yield per
event in $0-5\%$ central Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, in
the central rapidity window $|y|<1$, as a function of the transverse mass, for
several values of the invariant mass $M$, and two values of the shear
viscosity over entropy ratio $\eta/s$. The thick line is the McLerran-Toimela
spectrum (6). The thin lines are the global fit of our results using Eq. (19).
Open symbols: Dilepton yield from the Drell-Yan process. We only plot the
central value of the NLO+NLL calculation (see Sec. 5). Uncertainties are shown
in Fig. 3.

For each setup of our calculation, i.e., with or without quark suppression
taken into account, we have carried out a global fit of all our results for
$1.5<M_{t}<7$ GeV using Eq. (19). The best-fit values with quark suppression
are $a=0.61$ GeV-2, $b=1.6$ GeV-2, $n=3.1$. Without quark suppression, they
are $a=0.07$ GeV-2, $b=2.4$ GeV-2, $n=1.2$.333The error on $n$ is large for
the results without quark suppression. If one fixes the value of $n$ to the
same value as with quark suppression, the fit is almost as good, and the
parameters $a$ and $b$ are not significantly modified, so that the smaller
value of $n$ returned by the fit seems of little significance. As expected
from the discussion of Sec. 3, the breaking of $M_{t}$ scaling is larger
without quark suppression, resulting in a larger value of $b$. The parameter
$a$ is an order of magnitude smaller without quark suppression, which means
that 90% of the pre-equilibrium effects on the $M_{t}$ spectrum come from
chemical equilibration. Note that from mere dimensional analysis, by comparing
Eq. (19) with Eq. (13), one expects $a\sim b\sim 1/(\tau T^{3})\sim 2.5$
GeV-2, where the numerical estimate is given by Eq. (8). The values of $b$
returned by the fit are comparable, while those of $a$ are significantly
smaller, in particular when quark suppression is not implemented.

Fig. 2 displays the dilepton yield per event for two values of the viscosity
which roughly span the expected range in QCD [37]. By comparing with the
McLerran-Toimela spectrum, one sees that pre-equilibrium dynamics suppresses
dilepton production by at least a factor 10 for $M_{t}>5$ GeV. The dilepton
yield is still mostly determined by $M_{t}$, and the breaking of $M_{t}$
scaling is a modest effect.

We now evaluate the robustness of our results with respect to transverse flow,
which we have neglected. Transverse flow develops gradually over a time of the
order of the nuclear radius. It becomes important for $\tau\gtrsim 5$ fm/c. If
the fraction of the dileptons produced after $5$ fm/c is small, the dilepton
yield is likely to have little sensitivity to transverse flow. We have
calculated this fraction numerically and found that it only depends on
$M_{t}$. It is roughly 25% for $M_{t}=2$ GeV, but only 4% for $M_{t}=3$ GeV.
We conclude that for $M_{t}\lesssim 2$ GeV, sizable corrections from
transverse flow are to be expected. At RHIC, it has been argued on the basis
of simple dimensional arguments that these corrections are small in the
intermediate mass region [19]. However, hydrodynamic calculations have shown
that they are visible up to $M_{t}=2.5$ GeV [29]. Effects of transverse flow
are larger at LHC than at RHIC. We intend to study them in a future
publication.

The dilepton spectrum is often characterized by its effective temperature
$T_{\rm eff}$, defined as the inverse slope of the $M_{t}$ spectrum (Eq. (7)).
It is interesting to note that the spectrum depends on $M$ only through a
global factor in Eq. (19), so that $T_{\rm eff}$ still solely depends on
$M_{t}$, as long as transverse flow can be neglected.444The rise and fall of
$T_{\rm eff}$ as a function of $M$ observed by NA60 [18, 39] in the low-mass
region can be ascribed to transverse flow. In the limit of small $\eta/s$, Eq.
(19) gives:

$T_{\rm
eff}(M_{t})\simeq\frac{M_{t}}{6+2a\frac{\displaystyle\eta}{\displaystyle
s}M_{t}^{2}},$ (20)

where $a\simeq 0.61$ GeV-2. This equation shows that pre-equilibrium dynamics
decreases the effective temperature, and that the shear viscosity over entropy
ratio at early times $\eta/s$ can be extracted from the inverse slope. Note
that the inverse slope obtained in Ref. [4] using the Parton-Hadron String
Dynamics (PHSD) model exceeds the McLerran-Toimela value (7) already at RHIC
energies. We believe that this is due to the presence of hard particles of
jets/mini-jets in the PHSD Monte Carlo simulations. This source of large
invariant mass dileptons warrants further investigation.

We now discuss the centrality and system-size dependence of QGP dilepton
production. This dependence is encapsulated in the transverse area,
$A_{\perp}$, and in the value of $\tau T^{3}$ at late times. Both quantities
satisfy simple scaling laws as a function of the charged hadron multiplicity
per unit pseudorapidity, $dN_{\rm ch}/d\eta\propto A_{\perp}\tau T^{3}$. The
observation that the mean transverse momentum of hadrons $\langle
p_{t}\rangle$ depends weakly on centrality and system size [40] implies that
$A_{\perp}$ varies approximately like $(dN_{\rm ch}/d\eta)^{2/3}$ [41] as a
function of centrality and system size for fixed rapidity and collision
energy. This in turn implies that $\tau T^{3}$, which is proportional to
$(dN_{\rm ch}/d\eta)/A_{\perp}$ [6], scales like $(dN_{\rm ch}/d\eta)^{1/3}$.
Eq. (6) then shows that the McLerran-Toimela spectrum varies like $(dN_{\rm
ch}/d\eta)^{4/3}$. In other words, the dilepton yield scales like the space-
time volume, while the hadron yield scales like the volume at freeze-out. The
time component explains the extra factor $(dN_{\rm ch}/d\eta)^{1/3}$.

Eq. (13) then shows that the relative modification of the dilepton yield due
to pre-equilibrium dynamics varies with system size and centrality like
$(dN_{\rm ch}/d\eta)^{-1/3}$. This implies that the parameters $a$ and $b$ in
Eq. (19) are also proportional to $(dN_{\rm ch}/d\eta)^{-1/3}$. However, note
that local event-to-event fluctuations of the initial density [42], which we
neglect, will break this simple scaling.

Similar dimensional arguments can be used to predict the dependence of QGP
dilepton production on the collision energy $\sqrt{s_{\rm NN}}$. For a given
collision system, the transverse area $A_{\perp}$ is approximately independent
of $\sqrt{s_{\rm NN}}$, while the hadron multiplicity $dN_{\rm ch}/d\eta$
increases with $\sqrt{s_{\rm NN}}$ [43]. Therefore, $\tau T^{3}$ scales with
energy like $dN_{\rm ch}/d\eta$. This implies that the McLerran-Toimela
spectrum (6) is proportional to $(dN_{\rm ch}/d\eta)^{2}$. On the other hand,
the coefficients $a$ and $b$, which govern the modifications due to pre-
equilibrium effects, are proportional to the inverse Reynolds number (13),
i.e., to $(dN_{\rm ch}/d\eta)^{-1}$.

Finally, the dependence on rapidity $y$ follows the same scaling rules as the
dependence on collision energy, up to the replacement of $dN_{\rm ch}/d\eta$
with $dN_{\rm ch}/dy$.555Our calculation setup assumes longitudinal boost
invariance, but the results can still be applied if the multiplicity depends
on rapidity, since the longitudinal pressure gradient has a negligible effect
at LHC energy [23]. The McLerran-Toimela spectrum is proportional to $(dN_{\rm
ch}/dy)^{2}$, while $a$ and $b$ are proportional to $(dN_{\rm ch}/dy)^{-1}$.
The dilepton spectrum is therefore maximum at mid-rapidity, where $dN_{\rm
ch}/dy$ is maximum [44].

##  5 Background from the Drell-Yan process

Figure 3: (Color online) Comparison between the dilepton yield from production
in the QGP and from the Drell-Yan process for two fixed values of the
invariant mass $M=3.5$ GeV (top) and $M=4.5$ GeV (bottom). Lines: QGP
production with pre-equilibrium dynamics and quark suppression taken into
account, for $\eta/s=0.16$ (upper line) and $\eta/s=0.32$ (lower line). Dark
shaded band: Drell Yan process, with uncertainty from the parton distribution
function. Light shaded band: uncertainty on Drell-Yan from the renormalization
and factorization scales. Both bands are obtained by taking the envelope of
the results obtained by varying the model parameters.

The main backgrounds to QGP dilepton production in the intermediate mass
region, besides the large $J/\psi$ peak around $M\simeq 3.1$ GeV, are
semileptonic decays of heavy quark hadrons, and Drell-Yan production in the
initial state. The number of dileptons from charm hadron decays is expected to
exceed the QGP dilepton yield for $\sqrt{s_{\rm NN}}>40$ GeV [4]. Charm decays
are indeed observed to be the dominant source of dileptons at LHC energies
[45]. This source of background can be rejected based on the finite lifetime
of the heavy quark hadrons since the leptons from this source do not originate
from the primary vertex. Despite the small lifetime of charm hadron ground
states of $c\tau\approx$ 50-410 $\mu$m, the partial rejection of these decay
leptons in heavy-ion collisions is feasible and will strongly improve with the
future detector projects LHCb Upgrade 2 [46] and ALICE 3 [47]. We consider in
this phenomenological publication only the irreducible background from the
Drell-Yan process.

We compute the production of dileptons by the Drell-Yan process using the
Drell-Yan Turbo software package [48]. The cross section is evaluated at Next-
to-Leading Order (NLO) in the strong coupling constant $\alpha_{s}$. In
addition, the calculation employs a resummation at small transverse momentum
at next-to-leading logarithm (NLL). For the non-perturbative contribution to
the form factor, the same Gaussian form is chosen as in [48]. We evaluate the
uncertainty on the Drell Yan spectrum in the following way: We vary the
renormalization and factorization scales by a factor two independently
resulting in 8 variations with respect to the default choice. We take into
account the uncertainty on the parton distribution function in the EPPS
parametrization [49].

The central value of our calculation is plotted in Fig. 2 for two values of
the invariant mass $M$. The slope of the Drell-Yan $M_{t}$ spectrum is roughly
similar to that of the QGP spectrum, albeit slightly flatter. The main
difference between the two spectra is the normalization, which depends on the
invariant mass $M$. Drell-Yan production is enhanced for larger values of $M$
at a given $M_{t}$, contrary to QGP production. The physical explanation is
that the momenta of the incoming quark and antiquark responsible for Drell-Yan
production are mostly longitudinal, so that smaller values of the transverse
momentum $k_{t}$ (corresponding to larger values of $M$ at a given $M_{t}$)
are preferred. The kinematics of early QGP production is opposite, in the
sense that longitudinal momenta in the QGP are typically smaller than
transverse momenta. This means that the breaking of $M_{t}$ scaling alone
provide a handle to distinguish between QGP and Drell-Yan production.

Drell-Yan gradually takes over QGP production as $M$ increases. In order to
evaluate where the transition occurs, we compare both spectra in Fig. 3 for
two values of $M$ above the $J/\psi$ peak, and for the two values of $\eta/s$
used in Fig. 2. The uncertainty on the Drell-Yan spectrum from the
renormalization and factorization scales is displayed as a light shaded band.
The uncertainty from parton distribution function is displayed as a dark
shaded band. One sees that the scale is the dominant source of uncertainty.
The top panel shows that QGP production is likely to dominate over Drell-Yan
for $M$ up to $3.5$ GeV, at least for the lowest values of $M_{t}$. Looking at
the bottom panel, it seems unlikely that QGP dileptons can be isolated above
$M=4.5$ GeV. The precise value of $M$ above which Drell-Yan dominates over QGP
production depends on the value of $\eta/s$ at early times. Drell-Yan
calculations can be also performed down to lower invariant masses. However, at
low invariant masses the pQCD calculation exhibits increasingly large scale
uncertainties. Furthermore, the non-perturbative set-up of the used code
(NLO+NLL) was developed based on high invariant masses (Z-mass) and we hence
do not want to stretch into a regime where the chosen approach may not be the
most appropriate one that one could apply. Nevertheless, based on our
estimates, the production from the QGP and the preequilibrium phase is
expected to outshine the Drell-Yan contribution at smaller invariant masses.

The calculations displayed in the figures are carried out at mid-rapidity. The
rapidity dependence of the Drell-Yan is milder than that of the QGP spectrum
and goes in the opposite direction: Its minimum is at midrapidity. On the
other hand, the background from semileptonic decays of heavy quarks is easier
to eliminate at larger rapidities because the secondary vertex is farther from
the collision point.

##  6 Conclusions

We have calculated the spectrum of dileptons produced by the quark-gluon
plasma in ultrarelativistic heavy-ion collisions for invariant masses larger
than $1.2$ GeV, with input from state-of-the-art QCD kinetic theory to model
the kinetic and chemical equilibration of the QGP at early times. The
invariant spectrum $dN^{l^{+}l^{-}}/d^{4}K$ depends mostly on the transverse
mass $M_{t}$. The underpopulation of quarks at early times results in a
steeper $M_{t}$ spectrum. The viscosity over entropy ratio, which determines
the equilibration time of the QGP, can be inferred by measuring the slope of
the spectrum.

The anisotropy of the momentum distribution at early times breaks transverse
mass scaling, by suppressing the production of higher invariant masses $M$.
Interestingly, the trend is opposite for dileptons produced by the Drell-Yan
process, which is enhanced for larger $M$. Therefore, one can distinguish
experimentally QGP production from Drell-Yan production by studying the
variation of the dilepton yield as a function of $M$ at fixed $M_{t}$.

We have introduced a simple parametrization (19), from which one can infer the
dependence of QGP dilepton production on centrality, system size, rapidity,
and collision energy. Our modelization can be improved in several ways. For
the larger values of $M_{t}$, the contribution of minijets [4, 31] warrants
additional study. For the smaller values of $M_{t}$, one should take
transverse flow into account, as its effect on dilepton production is likely
to be significant for transverse masses $M_{t}\lesssim 2$ GeV. It has been
studied in detail in the low-mass region [50]. In the intermediate mass
region, detailed studies of transverse flow have been carried out [51, 52].
However, the effect of transverse flow on the slope on the $M_{t}$ spectrum
has only been studied at RHIC energy [29], and deserves further studies at
higher energies.

## Acknowledgments

This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the CRC-TR 211 “Strong-interaction matter
under extreme conditions” project number 315477589 – TRR 211 and in part in
the framework of the GLUODYNAMICS project funded by the “P2IO LabEx
(ANR-10-LABX-0038)” in the framework “Investissements d’Avenir”
(ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR),
France.

## References

  * [1] W. Busza, K. Rajagopal and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68 (2018), 339-376 doi:10.1146/annurev-nucl-101917-020852 [arXiv:1802.04801 [hep-ph]].
  * [2] M. Martinez and M. Strickland, Phys. Rev. C 78 (2008), 034917 doi:10.1103/PhysRevC.78.034917 [arXiv:0805.4552 [hep-ph]].
  * [3] R. Rapp, Adv. High Energy Phys. 2013 (2013), 148253 doi:10.1155/2013/148253 [arXiv:1304.2309 [hep-ph]].
  * [4] T. Song, W. Cassing, P. Moreau and E. Bratkovskaya, Phys. Rev. C 97 (2018) no.6, 064907 doi:10.1103/PhysRevC.97.064907 [arXiv:1803.02698 [nucl-th]].
  * [5] L. D. McLerran and T. Toimela, Phys. Rev. D 31 (1985), 545 doi:10.1103/PhysRevD.31.545
  * [6] M. Coquet, X. Du, J. Y. Ollitrault, S. Schlichting and M. Winn, Phys. Lett. B 821 (2021), 136626 doi:10.1016/j.physletb.2021.136626 [arXiv:2104.07622 [nucl-th]].
  * [7] S. D. Drell and T. M. Yan, Phys. Rev. Lett. 25 (1970), 316-320 [erratum: Phys. Rev. Lett. 25 (1970), 902] doi:10.1103/PhysRevLett.25.316
  * [8] R. Hagedorn, Nuovo Cim. Suppl. 3 (1965), 147-186 CERN-TH-520.
  * [9] J. Bartke et al. [Aachen-Berlin-Bonn-CERN-Cracow-Heidelberg-Warsaw], Nucl. Phys. B 120 (1977), 14-22 doi:10.1016/0550-3213(77)90092-X
  * [10] F. Becattini and U. W. Heinz, Z. Phys. C 76 (1997), 269-286 [erratum: Z. Phys. C 76 (1997), 578] doi:10.1007/s002880050551 [arXiv:hep-ph/9702274 [hep-ph]].
  * [11] L. Altenkämper, F. Bock, C. Loizides and N. Schmidt, Phys. Rev. C 96 (2017) no.6, 064907 doi:10.1103/PhysRevC.96.064907 [arXiv:1710.01933 [hep-ph]].
  * [12] J. D. Bjorken, Phys. Rev. D 27 (1983), 140-151 doi:10.1103/PhysRevD.27.140
  * [13] J. Ghiglieri and G. D. Moore, JHEP 12 (2014), 029 doi:10.1007/JHEP12(2014)029 [arXiv:1410.4203 [hep-ph]].
  * [14] H. T. Ding, O. Kaczmarek and F. Meyer, Phys. Rev. D 94 (2016) no.3, 034504 doi:10.1103/PhysRevD.94.034504 [arXiv:1604.06712 [hep-lat]].
  * [15] G. Jackson and M. Laine, JHEP 11 (2019), 144 doi:10.1007/JHEP11(2019)144 [arXiv:1910.09567 [hep-ph]].
  * [16] M. Laine, PoS CPOD2014 (2015), 065 doi:10.22323/1.217.0065 [arXiv:1502.05796 [hep-ph]].
  * [17] R. Baier, P. Romatschke, D. T. Son, A. O. Starinets and M. A. Stephanov, JHEP 04 (2008), 100 doi:10.1088/1126-6708/2008/04/100 [arXiv:0712.2451 [hep-th]].
  * [18] R. Arnaldi et al. [NA60], Eur. Phys. J. C 61 (2009), 711-720 doi:10.1140/epjc/s10052-009-0878-5 [arXiv:0812.3053 [nucl-ex]].
  * [19] R. Rapp and H. van Hees, Phys. Lett. B 753 (2016), 586-590 doi:10.1016/j.physletb.2015.12.065 [arXiv:1411.4612 [hep-ph]].
  * [20] R. A. Tripolt, Nucl. Phys. A 1005 (2021), 121755 doi:10.1016/j.nuclphysa.2020.121755 [arXiv:2001.11232 [hep-ph]].
  * [21] J. Adamczewski-Musch et al. [HADES], Nature Phys. 15 (2019) no.10, 1040-1045 doi:10.1038/s41567-019-0583-8
  * [22] K. Kajantie, M. Kataja, L. D. McLerran and P. V. Ruuskanen, Phys. Rev. D 34 (1986), 811 doi:10.1103/PhysRevD.34.811
  * [23] J. Y. Ollitrault, Eur. J. Phys. 29 (2008), 275-302 doi:10.1088/0143-0807/29/2/010 [arXiv:0708.2433 [nucl-th]].
  * [24] E. V. Shuryak, Phys. Rev. Lett. 68 (1992), 3270-3272 doi:10.1103/PhysRevLett.68.3270
  * [25] A. Kurkela and A. Mazeliauskas, Phys. Rev. D 99 (2019) no.5, 054018 doi:10.1103/PhysRevD.99.054018 [arXiv:1811.03068 [hep-ph]].
  * [26] E. V. Shuryak and L. Xiong, Phys. Rev. Lett. 70 (1993), 2241-2244 doi:10.1103/PhysRevLett.70.2241 [arXiv:hep-ph/9301218 [hep-ph]].
  * [27] J. Vredevoogd and S. Pratt, Phys. Rev. C 79 (2009), 044915 doi:10.1103/PhysRevC.79.044915 [arXiv:0810.4325 [nucl-th]].
  * [28] A. Kurkela, A. Mazeliauskas, J. F. Paquet, S. Schlichting and D. Teaney, Phys. Rev. Lett. 122 (2019) no.12, 122302 doi:10.1103/PhysRevLett.122.122302 [arXiv:1805.01604 [hep-ph]].
  * [29] J. Deng, Q. Wang, N. Xu and P. Zhuang, Phys. Lett. B 701 (2011), 581-586 doi:10.1016/j.physletb.2011.06.027 [arXiv:1009.3091 [nucl-th]].
  * [30] R. Arnaldi et al. [NA60], Phys. Rev. Lett. 100 (2008), 022302 doi:10.1103/PhysRevLett.100.022302 [arXiv:0711.1816 [nucl-ex]].
  * [31] R. Paatelainen, K. J. Eskola, H. Niemi and K. Tuominen, Phys. Lett. B 731 (2014), 126-130 doi:10.1016/j.physletb.2014.02.018 [arXiv:1310.3105 [hep-ph]].
  * [32] J. Churchill, L. Yan, S. Jeon and C. Gale, Phys. Rev. C 103 (2021) no.2, 024904 doi:10.1103/PhysRevC.103.024904 [arXiv:2008.02902 [hep-ph]].
  * [33] X. Du and S. Schlichting, Phys. Rev. D 104 (2021) no.5, 054011 doi:10.1103/PhysRevD.104.054011 [arXiv:2012.09079 [hep-ph]].
  * [34] M. P. Heller and M. Spalinski, Phys. Rev. Lett. 115 (2015) no.7, 072501 doi:10.1103/PhysRevLett.115.072501 [arXiv:1503.07514 [hep-th]].
  * [35] M. P. Heller, A. Kurkela, M. Spaliński and V. Svensson, Phys. Rev. D 97 (2018) no.9, 091503 doi:10.1103/PhysRevD.97.091503 [arXiv:1609.04803 [nucl-th]].
  * [36] G. Giacalone, A. Mazeliauskas and S. Schlichting, Phys. Rev. Lett. 123 (2019) no.26, 262301 doi:10.1103/PhysRevLett.123.262301 [arXiv:1908.02866 [hep-ph]].
  * [37] N. Christiansen, M. Haas, J. M. Pawlowski and N. Strodthoff, Phys. Rev. Lett. 115 (2015) no.11, 112002 doi:10.1103/PhysRevLett.115.112002 [arXiv:1411.7986 [hep-ph]].
  * [38] G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 87 (2001), 081601 doi:10.1103/PhysRevLett.87.081601 [arXiv:hep-th/0104066 [hep-th]].
  * [39] R. Arnaldi et al. [NA60], Eur. Phys. J. C 59 (2009), 607-623 doi:10.1140/epjc/s10052-008-0857-2 [arXiv:0810.3204 [nucl-ex]].
  * [40] S. Acharya et al. [ALICE], Phys. Lett. B 788 (2019), 166-179 doi:10.1016/j.physletb.2018.10.052 [arXiv:1805.04399 [nucl-ex]].
  * [41] F. G. Gardim, G. Giacalone and J. Y. Ollitrault, Phys. Lett. B 809 (2020), 135749 doi:10.1016/j.physletb.2020.135749 [arXiv:1909.11609 [nucl-th]].
  * [42] C. E. Aguiar, Y. Hama, T. Kodama and T. Osada, Nucl. Phys. A 698 (2002), 639-642 doi:10.1016/S0375-9474(01)01447-6 [arXiv:hep-ph/0106266 [hep-ph]].
  * [43] J. Adam et al. [ALICE], Phys. Rev. Lett. 116 (2016) no.22, 222302 doi:10.1103/PhysRevLett.116.222302 [arXiv:1512.06104 [nucl-ex]].
  * [44] J. Adam et al. [ALICE], Phys. Lett. B 772 (2017), 567-577 doi:10.1016/j.physletb.2017.07.017 [arXiv:1612.08966 [nucl-ex]].
  * [45] S. Acharya et al. [ALICE], Phys. Rev. C 99 (2019) no.2, 024002 doi:10.1103/PhysRevC.99.024002 [arXiv:1807.00923 [nucl-ex]].
  * [46] R. Aaij et al. [LHCb], [arXiv:1808.08865 [hep-ex]].
  * [47] D. Adamová, G. Aglieri Rinella, M. Agnello, Z. Ahammed, D. Aleksandrov, A. Alici, A. Alkin, T. Alt, I. Altsybeev and D. Andreou, et al. [arXiv:1902.01211 [physics.ins-det]].
  * [48] S. Camarda, M. Boonekamp, G. Bozzi, S. Catani, L. Cieri, J. Cuth, G. Ferrera, D. de Florian, A. Glazov and M. Grazzini, et al. Eur. Phys. J. C 80 (2020) no.3, 251 [erratum: Eur. Phys. J. C 80 (2020) no.5, 440] doi:10.1140/epjc/s10052-020-7757-5 [arXiv:1910.07049 [hep-ph]].
  * [49] K. J. Eskola, P. Paakkinen, H. Paukkunen and C. A. Salgado, Eur. Phys. J. C 77 (2017) no.3, 163 doi:10.1140/epjc/s10052-017-4725-9 [arXiv:1612.05741 [hep-ph]].
  * [50] G. Vujanovic, C. Young, B. Schenke, R. Rapp, S. Jeon and C. Gale, Phys. Rev. C 89 (2014) no.3, 034904 doi:10.1103/PhysRevC.89.034904 [arXiv:1312.0676 [nucl-th]].
  * [51] R. Ryblewski and M. Strickland, Phys. Rev. D 92 (2015) no.2, 025026 doi:10.1103/PhysRevD.92.025026 [arXiv:1501.03418 [nucl-th]].
  * [52] B. S. Kasmaei and M. Strickland, Phys. Rev. D 99 (2019) no.3, 034015 doi:10.1103/PhysRevD.99.034015 [arXiv:1811.07486 [hep-ph]].