# A new Microscopic Model for $J/\psi$ Production in Heavy Ion Collisions

Denys Yen Arrebato Villar, Joerg Aichelin, Pol Bernard Gossiaux  SUBATECH,
Université de Nantes, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France

###### Abstract

We present a new model for the creation of $J/\psi$ mesons in
ultrarelativistic heavy ion collisions, which allows to follow the individual
heavy quarks from their creation until the detector through the Quark Gluon
Plasma (QGP), which is formed in these collisions and described by the EPOS2
event generator. The $c$ and $\bar{c}$ quarks interact via a potential, based
on results of lattice gauge calculations. The annihilation and creation of
$J/\psi$ is described by a density matrix approach whose time evolution is
studied in the expanding system. The comparison with PbPb data at
$\sqrt{s}$=5.02 TeV shows that this model can describe simultaneously the
nuclear modification factor $R_{AA}$ and the elliptic flow $v_{2}$ of the
$J/\psi$ .

###### pacs:

12.38Mh

##  I Introduction

There is overwhelming evidence that in ultrarelativistic heavy ion collisions
a plasma of quarks and gluons (QGP) is created, which evolves in time and
disintegrates at the end of its lifetime into hadrons. The multiplicity of
light and strange hadrons is well described by statistical model calculations
Andronic et al. (2018). The consequence of this observation is that light and
strange hadrons cannot provide direct information about the time evolution of
the QGP from its creation to hadronization. To study this time evolution and
to get insight into the early phase of the heavy ion collision one has to
focus on probes which do not come to equilibrium with the expanding QGP. They
include electromagnetic probes, jets as well as hadrons which contain heavy
quarks.

Among these probes, especially the hidden heavy flavour meson $J/\psi$,
composed of a $c$ and a $\bar{c}$, quark has recently gained a lot of
interest. This is due to two experimental results, which came as a surprise.

  * •

The nuclear modification factor
$R_{AA}=\frac{d\sigma_{AA}/dp_{T}}{N_{coll}d\sigma_{pp}/dp_{T}}$, where
$N_{coll}$ is the number of initial binary collisions in the AA system, stays
almost constant as a function of the centrality in heavy ion collisions at LHC
Adam et al. (2017) whereas it decreases strongly at RHIC energies Adare et al.
(2011).

  * •

In elementary collisions $J/\psi$ , as well as charm quarks, are produced
isotropically in azimuthal direction and show therefore no elliptic flow. The
$J/\psi$’s observed in experiment show a strong elliptic flow which follows
the systematics of the $v_{2}$ observed for light hadrons Acharya et al.
(2020a). This contradicts the idea that the $J/\psi$ traverses the QGP as a
color-neutral, weakly interacting object.

In this paper we study how these observations can be understood and what we
can learn from the $J/\psi$ about the properties of the QGP, created in heavy
ion collisions.

The idea to use $J/\psi$’s for such studies goes back to the seminal paper of
Matsui and Satz Matsui and Satz (1986) who argued that in strongly interacting
thermal matter the color charges of the $c$ and $\bar{c}$ are screened by
color charges of the medium to the extent that the $J/\psi$ melts if the
density of these charges becomes too large. Later this melting has been
confirmed by lattice gauge calculations Mocsy and Petreczky (2007); Digal et
al. (2001) but the exact dissociation temperature, $T_{\rm diss}$, is still
unknown.

The Wilson loop allows to determine the free energy between the c and
$\bar{c}$ as a function of their distance. The lattice gauge results for the
Wilson loop as a function of the temperature allowed to develop a static
$c\bar{c}$ potential which can be employed in a Schrödinger equation and
allows for studying how the ground state energy of the $c\bar{c}$ pair
develops as a function of the temperature of the QGPBala et al. (2021). Also
these calculations came also to the conclusion that there is a limiting
temperature above which the $J/\psi$ becomes unstable. For a recent review we
refer to Ref.Rothkopf (2020).

If produced in ultrarelativistic heavy ion collisions the situation is more
complex than in a static medium. Shortly after the initial binary collisions
of the nucleons of projectile and target a high temperature QGP is formed in
which a $J/\psi$ cannot survive. It can only be produced when the temperature
of the expanding system gets lower than $T_{\rm diss}$. Therefore, the c and
$\bar{c}$ of the final $J/\psi$ and those $c$ and $\bar{c}$ quarks, which are
finally part of open heavy flavour hadrons, traverse initially the same QPG.
Therefore the knowledge, which one has acquired in the last years about open
heavy flavour mesons, is of use also for the study of $J/\psi$ .

Open heavy flavour hadrons, produced in heavy ion collisions, have been
extensively studied in the last years, experimentally and theoretically.
Recently the theoretical models, which differ in details, have been compared
Xu et al. (2019); Cao et al. (2019); Beraudo et al. (2018). This comparison
suggests that the initial c and $\bar{c}$ quarks are created in elementary
baryon-baryon collisions at the beginning of the heavy ion reaction and that
the initial transverse momentum distribution is well described by FONLL (first
order next to leading log) calculations Cacciari et al. (2008, 2012). The
heavy quarks interact subsequently with the QGP constituents, the quarks and
gluons, in elementary collisions, which are described by pQCD (perturbative
QCD) Born diagrams, and pick up by these collisions a finite elliptic flow.
Finally, they convert into heavy hadrons when the QGP hadronizes. The last
process is usually described by a combination of coalescence and
fragmentation. Thereafter hadronic rescattering modifies still the momentum
distribution, as shown in Ref.Song et al. (2015).

When the local temperature of the QGP gets lower than $T_{\rm diss}$ of the
$J/\psi$ , the $J/\psi$ can be formed but also destroyed by a elastic or a
inelastic collision of one of its constituents with a QGP parton. The
difference between formation and collisional decay determines the $J/\psi$
spectrum at the end of the QGP phase. A modification of the $J/\psi$
distribution by hadronic rescattering is possible, Ref. Linnyk et al. (2008),
but beyond the scope of the present article.

Recently it has been shown that also in the statistical hadronization model
the relative multiplicity between charmed hadrons in central collisions at LHC
energies can be understood assuming that all the charmed hadrons are formed at
chemical freeze out, when also the light hadrons are produced Andronic et al.
(2021). This questions whether multiplicities in central collisions are
sufficient to distinguish between a creation at the static freeze out and a
creation during the expansion of the QGP. To explore this one has to employ
other observables.

Transport models have been advanced to study the dynamical production of
$J/\psi$ Du and Rapp (2015); Zhou et al. (2014). The model of Du et al.Du and
Rapp (2015) describes the $J/\psi$ production in central PbPb collisions at
$\sqrt{s}=2.76$ by a kinetic rate equation applied in an expanding fireball.
The rate for $J/\psi$ production in the fireball is inspired from lattice
results for the free energy and only dependent on the temperature of the
system. This production in the fireball is supplemented by substantial
regeneration contribution in the expanding hadron gas after the fireball has
been disintegrated into hadrons. Being in equilibrium with the expanding
fireball, the c and $\bar{c}$ acquire a finite elliptical flow when the
geometrical anisotropy in coordinate space is converted into an anisotropy in
momentum space. The absolute value of the elliptic flow is, however,
underestimated. Introducing in addition some space-momentum correlations
between the $c$ and $\bar{c}$ the value of the elliptic flow can be increased
He et al. (2022).

Zhou et al.Zhou et al. (2014) have advanced a dynamical model for dissociation
and regeneration of $J/\psi$ when the c and $\bar{c}$ pass through the QGP,
which is modeled by hydrodynamics. Dissociation and regeneration are
calculated via the $\sigma_{gJ/\psi}$ cross section, assuming that the charm
quarks are in equilibrium in the QPG.

There are also two more recent and interesting approaches, which have not yet
yield quantitative predictions. The one is the treatment of the $J/\psi$
production under the aspect of an open quantum system Akamatsu and Miura
(2022); Delorme et al. (2022) whose time evolution is given by the Lindblad
equation. The other is the description of the $J/\psi$ by the time evolution
of a reduced density matrix Blaizot and Escobedo (2018). Both treat the $c$
and $\bar{c}$ pairs as quantum systems, a challenging as well as complex task.

In this paper we advance a model for the $J/\psi$ production which follows the
c and $\bar{c}$s from the initial creation until hadronization. By this we
avoid one body transport approaches like Boltzmann or Fokker-Planck equations.
These equations are not appropriate to study two body correlations, which are
at the origin of the $J/\psi$ formation. While travelling through the QGP the
heavy quarks have energy and momentum conserving collisions with the
constituents of the QGP and interact among themselves by a potential derived
from lattice QCD. The Lagrangian, which we employ, includes relativistic
corrections up to the order $\gamma-1$ where $\gamma=1/\sqrt{1-\beta^{2}}$.
Below $T_{\rm diss}$ the $J/\psi$ s are described by the Wigner density in
relative coordinates with a root mean square (rms) radius, which depends on
the temperature of the QGP. Above $T_{\rm diss}$ a $J/\psi$ cannot be
produced. The rate of production and dissociation is obtained by solving the
von Neumann equation for the two body $c\bar{c}$ system in the expanding
medium, following a formalism which has been developed by Remler et al. Remler
(1981); Gyulassy et al. (1983); Aichelin and Remler (1987) for the production
of deuterons in heavy ion collisions. It has also been employed in the study
of quarkonia production in pp collisions within the PHSD approach Song et al.
(2017).

The paper is organized as follows: In section II we present our model. We
introduce the density matrix formalism introduced by Remler and study the rate
of $J/\psi$ production for time independent $Q\bar{Q}$ Wigner densities. This
is followed by a description of the interaction of heavy quarks with the QGP
partons. Finally we discuss the non relativistic $Q\bar{Q}$ Wigner density and
its relativistic extension. In section III we extend our formalism to the case
that the $Q\bar{Q}$ Wigner density gets time dependent. Section IV is devoted
to the potential interaction between $Q$ and $\bar{Q}$. In section V we report
about the initial distribution of the heavy (anti)quarks. In section VI we
present numerical details of our approach and study the consequences of the
different ingredients on the observables. In section VII we compare our
results with experimental data before we draw our conclusions in section VIII.

##  II The Model

We start out with an outline of the approach of Remler, which we employ,
adapted to the problem of heavy quarks: In the initial collisions between
projectile and target nucleons heavy (anti)quarks $Q(\bar{Q})$ are created,
which we assume to be uncorrelated in momentum space. Their transverse
momentum distribution is given by FONLL calculations Cacciari et al. (2008,
2012). When the QGP is created the heavy quarks enter the QGP which is, as the
initial interaction points, modelled by EPOS2. The expansion of the QGP is
described by ideal hydrodynamics Drescher et al. (2001); Werner et al. (2010).
While traversing the QGP the heavy quarks interact with the plasma
constituents Gossiaux and Aichelin (2008); Gossiaux et al. (2009). At the same
time $Q\bar{Q}$ pairs interact among themselves via a chromoelectical
potential. When the QGP has cooled down locally to $T_{\rm diss}^{\Phi}$, the
dissociation temperature of a heavy $Q\bar{Q}$ meson of type $\Phi$, these
mesons can be created but also destroyed. We employ the Remler formalism to
describe their creation and annihilation. These processes cease when the heavy
quarks hadronize to open heavy flavour hadrons. The Remler formalism gives us
the final momentum distribution of the quarkonia. Final state interactions
with the hadronic matter are presently not taken into account.

###  II.1 The Remler density matrix formalism

All information about a N-particle system is encoded in the N-body density
operator, $\rho_{N}(t)$, of the system. Among the N particles there may be one
or several $c\bar{c}$ pairs. Because the relative motion of heavy quarks in
bound heavy quark systems is small compared to the heavy quark mass, we use
here non relativistic kinematics and discuss the extension to towards a
relativistic treatment later.

The density operator obeys the von Neumann equation Gyulassy et al. (1983):

$\partial\rho_{N}/\partial t=-\frac{i}{\hbar}[H,\rho_{N}]$ (1)

where $H$ is the Hamiltonian of the full system

$H=\Sigma_{i}K_{i}+\Sigma_{i>j}V_{ij}.$ (2)

$K_{i}$ is the kinetic energy operator of the particle $i$ and $V_{ij}$ is the
interaction between the particles $i$ and $j$. Quarkonia, like a $J/\psi$, are
two body objects described by the two-body density operator
$\rho^{\Phi}=|\Phi><\Phi|$. $\Phi$ is the wave function of the eigenstate
$\Phi$ of the two body $Q\bar{Q}$ system. Thus

$P^{\Phi}(t)={\rm Tr}[\rho^{\Phi}\rho_{N}(t)],$ (3)

where the trace is taken over all N-body coordinates (which includes the $Q$
and $\bar{Q}$ degrees of freedom), measures the probability of finding the $Q$
and the $\bar{Q}$ at time t in the eigenstate $|\Phi>$. We are in particular
interested in the value of $P^{\Phi}(t\to+\infty)$, as it corresponds to
experimental measurements. From the viewpoint of heavy quarks, our standard
EPOS2+MC@HQ is quite similar to an intranuclear cascade model to which the
Remler algorithm was originally applied: In heavy ion reactions the QGP
expands until hadronization. Without potential the distance between the $Q$
and $\bar{Q}$ quarks increases and at the end of the QGP expansion it is large
with respect to the radius of the eigen state $\Phi$. Therefore
$P^{\Phi}(t\to\infty)$ tends to zero. To circumvent this issue, we resort to
the method of Remler’s original work: We express the probability to observe
$Q\bar{Q}$ pairs in the eigen state $\Phi$ at a time $\tilde{t}$ as the
integral of the rate of decay and formation of pairs in the eigen state
$\Phi$, $\Gamma^{\Phi}(t)$:

$\displaystyle P^{\Phi}(\tilde{t})$ $\displaystyle=$ $\displaystyle
P^{\Phi}(0)+\int_{0}^{\tilde{t}}\Gamma^{\Phi}(t)dt$ (4)

with the rate $\Gamma^{\Phi}$ defined as

$\Gamma^{\Phi}(t)=\frac{dP^{\Phi}}{dt}=\frac{d}{dt}{\rm
Tr}[\rho^{\Phi}\rho_{N}(t)].$ (5)

In our numerical scheme, we have introduced an attractive potential acting
between $Q$ quarks and $\bar{Q}$ antiquarks, see section IV. Consequently for
some pairs the relative distance between the $Q$ and $\bar{Q}$ remains finite
when $t\to\infty$. However, in our semi-classical modelling, which is best
suited when many momentum exchanges occur but less reliable to describe the
long time dynamics of this quantity, the formulation based on the rate, Eq. 4,
is more accurate than the direct projection Eq. 3. In a full quantum evolution
of $\rho_{N}$ both methods would give identical results.

Proceeding with the time derivative inside the ${\rm Tr[\cdots]}$, assuming
that $\rho^{\Phi}$ is time independent and using the von Neumann equation (1),
one gets

$\displaystyle P^{\Phi}(\tilde{t})$ $\displaystyle=$ $\displaystyle
P^{\Phi}(0)+\int_{0}^{\tilde{t}}Tr[\rho^{\Phi},\frac{\partial\rho_{N}}{\partial
t}]dt$ (6) $\displaystyle=$ $\displaystyle
P^{\Phi}(0)-\frac{i}{\hbar}\int_{0}^{\tilde{t}}Tr[\rho^{\Phi},[H,\rho_{N}]]dt.$

We first focus on the case that among the N particles we find only a single
$Q\bar{Q}$ pair. We assign to this $Q\bar{Q}$ pair the indices 1 and 2 and
decompose the total Hamiltonian:

$H=H_{1,2}+H_{N-2}+U_{1,2}$ (7)

where

$H_{1,2}=K_{1}+K_{2}+V_{1,2}$ (8)

is the two particle Hamiltonian of the $Q\bar{Q}$ pair,
$H_{N-2}=\Sigma_{i}K_{i}+\Sigma_{j>i\geq 3}V_{ji}$ is the Hamiltonian of the
remaining N-2 body system and $U_{1,2}$ is the interaction of the heavy quarks
1 and 2 with the rest of the system

$U_{1,2}=\Sigma_{j}V_{1,j}+\Sigma_{j}V_{2,j}.$ (9)

We replace in Eq. 6 the full Hamiltonian of the system by this decomposition
and profit from the relations

$[\rho^{\Phi},H_{1,2}]=0$ (10)

because $|\Phi\rangle$ is an eigenstate of $H_{1,2}$ and

$[\rho^{\Phi},H_{N-2}]=0$ (11)

because $H_{N-2}$ does act only on the remaining N-2 particles due to the
cyclic property of the trace. Therefore, we can write

$\frac{dP^{\Phi}(t)}{dt}=\Gamma^{\Phi}(t)=\frac{-i}{\hbar}Tr[\rho^{\Phi}[U_{1,2},\rho_{N}(t)]].$
(12)

This is the starting point of our approach. With eq.12 we calculate the
probability $P^{\Phi}(\tilde{t})$ that a $Q$ and a $\bar{Q}$ are in a bound
state $\Phi$ at $t=\tilde{t}$ by integrating the rate from $t=0$ to
$t=\tilde{t}$. $P^{\Phi}(\tilde{T}\to\infty)$ is then the probability that at
the end of the heavy ion reaction a meson of type $\Phi$ is observed. To make
calculations possible we have to know $\rho_{N}(t)$. A full quantum treatment
of the evolution of $\rho_{N}$ or of the equivalent N-body Wigner density,
$W_{N}$, defined as

$\displaystyle
W_{N}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)=\frac{1}{h^{3N}}\int
d^{3}y_{1}...d^{3}y_{N}(e^{\mathbf{p}_{1}\cdot\mathbf{y}_{1}}...e^{\mathbf{p}_{N}\cdot\mathbf{y}_{N}})$
$\displaystyle\langle\mathbf{r}_{1}+\frac{\mathbf{y}_{1}}{2},...,\mathbf{r}_{N}+\frac{\mathbf{y}_{N}}{2}|\rho(t)|\mathbf{r}_{1}-\frac{\mathbf{y}_{1}}{2},...,\mathbf{r}_{N}-\frac{\mathbf{y}_{N}}{2}\rangle,$
(13)

where $\mathbf{r}_{i}$ and $\mathbf{p}_{i}$ are the coordinates and momentum
of the particles in the Wigner representation, is out of reach but in the past
it turned out that many observables in heavy ion collisions can be well
described if one replaces the N-body Wigner density by an average over
classical N-body phase space densities $W^{c}_{N}$

$W_{N}\approx\langle W^{c}_{N}\rangle$ (14)

with

$W_{N}^{c}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)=\prod_{i}^{N}\delta(\mathbf{r_{i}}-\mathbf{r}_{i0}(t))\delta(\mathbf{p_{i}}-\mathbf{p}_{i0}(t)).$
(15)

$W_{N}^{c}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)$ as well as
$W_{N}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)$ are normalized to 1.

$\displaystyle\int\prod_{i=1}^{N}d^{3}r_{i}d^{3}p_{i}W_{N}^{c}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)=1$
$\displaystyle\int\prod_{i=1}^{N}d^{3}r_{i}d^{3}p_{i}W_{N}(\\{\mathbf{r_{i}}\\},\\{\mathbf{p_{i}}\\},t)=1.$
(16)

###  II.2 The rate for time independent $Q\bar{Q}$ Wigner densities

In this section we assume that, as in the original Remler formalism,
$W^{\Phi}$, the Wigner density of density matrix of the eigenstates of the
$Q\bar{Q}$ Hamiltonian, $|\Phi\rangle\langle\Phi|$, is time independent. The
extension to a time dependent $W^{\Phi}(t)$ will be discussed in section III.

Employing Wigner densities we can rewrite the rate, Eq.5. We assume again that
among the N particles there is only one heavy quark Q which carries the index
1 and one heavy antiquark $\bar{Q}$ with the index 2. Then we find

$\displaystyle\frac{dP^{\Phi}(t)}{dt}$ $\displaystyle=$
$\displaystyle\Gamma^{\Phi}(t)=h^{3}\frac{d}{dt}\int\prod_{j}^{N}d^{3}{r}_{j}d^{3}{p}_{j}W^{\Phi}_{12}W^{c}_{N}(t)$
(17) $\displaystyle=$ $\displaystyle
h^{3}\int\prod_{j}^{N}d^{3}\mathbf{r}_{j}d^{3}\mathbf{p}_{j}\
W^{\Phi}_{12}\frac{\partial}{\partial t}W^{c}_{N}(t).$

where
$W^{\Phi}_{12}=W^{\Phi}(\mathbf{r_{1}}-\mathbf{r_{2}},\mathbf{p_{1}}-\mathbf{p_{2}})$.
The form of $W^{\Phi}$, the Wigner density of the quarkonium $\Phi$, will be
discussed in section II.4.

The interaction between the $N$ partons is of short range (as compared to the
mean free path). This means that we consider that the QGP partons and the
heavy quarks move on straight line trajectories between the collisions whose
strength is given by cross sections.

We can number the collisions between a given couple of scattering partners $i$
and $j$ by $n$, up to $n^{\rm max}_{ij}$. We define as $t_{ij}(n)$ the time at
which the $n$-th collision, involving the partons $i$ and $j$, takes place.
This allows to calculate the momentum of particle $i$ at time $t$ as

$p_{i}(t)=p_{i}(0)+\Sigma_{j\neq i}\Sigma_{n=1}\Theta(t-t_{ij}(n))\Delta
p_{ij}(n)$ (18)

where $\Delta p_{ij}(n)$ is the momentum transfer in the $n-$th collision and
where the sum on $n$ runs from 1 $\to n^{\rm max}_{ij}$. This notation we use
implicitly from now on. $\Delta p_{ij}(n)$ is equal $-\Delta p_{ji}(n)$. With
this choice of time dependent momenta in the Wigner density $W_{N}^{c}$
(Eq.15) we can calculate the time evolution of the N-body Wigner density,
(Eq.17)

$\displaystyle\frac{\partial}{\partial t}W_{N}^{c}(t)$ $\displaystyle=$
$\displaystyle\Sigma_{i}v_{i}\cdot\partial_{r_{i}}W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t)$
$\displaystyle+$ $\displaystyle\Sigma_{j\geq i}\Sigma_{n}\delta(t-t_{ij}(n))$
$\displaystyle\cdot$
$\displaystyle(W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t+\epsilon)-W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t-\epsilon)).$

The first term arises from the straight line motion of the particles between
the collisions while the second is due to the impulse received at the time
$t_{ij}(n)$ when the $n$-th collision between particle $i$ and $j$ takes
place. The $\delta(t-t_{ij}(n))$ assures that a momentum transfer takes place
exactly at the time of collisions. We can separate the change of $W_{N}$ due
to collisions and straight line motion by writing

$\partial\rho_{N}(t)/\partial
t=-\frac{i}{\hbar}\Sigma_{j}[K_{j},\rho_{N}(t)]-\frac{i}{\hbar}\Sigma_{k>j}[V_{jk},\rho_{N}(t)].$
(20)

From the comparison between the equations (II.2) and (20) we find that

$\displaystyle-\frac{i}{\hbar}\Sigma_{j}[K_{j},\rho_{N}(t)]\equiv\langle\Sigma_{i}v_{i}\cdot\partial_{r_{i}}W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t)\rangle$
(21)

and (renaming indices)

$\displaystyle-\frac{i}{\hbar}\Sigma_{k>j}[V_{jk},\rho_{N}(t)]\equiv\langle\Sigma_{k>j}\Sigma_{n}\delta(t-t_{jk}(n))$
$\displaystyle\cdot$
$\displaystyle(W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t+\epsilon)-W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t-\epsilon))\rangle.$

Strictly speaking we assume that the equivalence holds for each term of the
sum separately. This means that like in cascade calculations the interaction
range is small as compared to the mean free path. Substituting in Eq. 12 the
square bracket by the r.h.s of Eq. II.2 we find

$\displaystyle\Gamma^{\Phi}(t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2}\sum_{j\geq
3}^{N}\sum_{n}\delta(t-t_{ij}(n))\int\prod^{N}_{k=1}d^{3}\mathbf{r}_{k}d^{3}\mathbf{p}_{k}$
$\displaystyle\cdot$ $\displaystyle
h^{3}W^{\Phi}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{p}_{1},\mathbf{p}_{2})$
$\displaystyle\cdot$
$\displaystyle[W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\};t+\epsilon)-W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\};t-\epsilon)],$

where $j$ sums over the light quarks and gluons of the QGP. Hence, in the
Remler formalism, collisions of the $Q$ and $\bar{Q}$ with the QGP medium
determine the rate of creation and destruction of quarkonia.

Figure 1: Visualization of Eq.II.2. A quark (gluon) from the QGP collides with
a heavy quark Q causing a momentum change of the heavy quark. The probability
that the heavy quark formed a quarkonia $\Phi$ with the $\bar{Q}$ before the
collision ($W-=W^{\Phi}W_{N}^{c}(t-\epsilon)$) hence differs from the
probability that it forms a $\Phi$ after the collision
($W+=W^{\Phi}W_{N}^{c}(t+\epsilon)$). The difference $W+\ $-$\ W-$ is
therefore the change of the $J/\psi$ multiplicity due to this collision.

Fig. 1 visualizes Eq.II.2. A parton from the QGP collides with a heavy quark Q
leading to a momentum change of Q. The Q may form with each of the heavy
antiquarks $\bar{Q}$ a quarkonium. We calculate $W^{\Phi}$ before ($W_{-}$)
and after ($W_{+}$) the collision for the $Q\bar{Q}$ pairs which the $Q$ can
form.

It is useful to explain this equation a bit more. At time $t=t_{1j}(n)$ a
heavy quark $Q$, to which we assign the coordinates
$\mathbf{r}_{1},\mathbf{p}_{1}$, has a collision with a QGP parton $j\geq$ 3\.
To the heavy antiquark we assign the coordinates
$\mathbf{r}_{2},\mathbf{p}_{2}$, and define the two body Wigner density of the
$Q\bar{Q}$ pair as

$W_{2}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{p}_{1},\mathbf{p}_{2},t)=\prod_{k=3}^{N}\int
d^{3}r_{k}d^{3}p_{k}W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t).$ (24)

We can calculate the contribution of this collision to the yield of the state
$\Phi$. We define the relative and center of mass coordinates of the
$Q\bar{Q}$ pair $\mathbf{q}=\frac{\mathbf{p}_{1}-\mathbf{p}_{2}}{2}\
(\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2})$ and
$\mathbf{P}=\mathbf{p}_{1}+\mathbf{p}_{2}\
(\mathbf{R}=\frac{\mathbf{r}_{1}+\mathbf{r}_{2}}{2})$. The contribution of
this Q-parton $n^{\rm th}$ collision to the $\Phi$ production rate can then be
expressed as

$\displaystyle\Gamma^{\Phi}_{1,2;j}(n;t)=h^{3}\delta(t-t_{1j}(n))\int
d^{3}Pd^{3}Rd^{3}rd^{3}qW{{}^{\Phi}}(\mathbf{r},\mathbf{q})$ (25)
$\displaystyle\cdot$
$\displaystyle\left(W_{2}(\mathbf{R},\mathbf{r},\mathbf{P},\mathbf{q},t+\epsilon)-(W_{2}(\mathbf{R},\mathbf{r},\mathbf{P},\mathbf{q},t-\epsilon)\right).$

This allows for expressing the total rate in a form suitable for Monte Carlo
implementations

$\Gamma_{1,2}^{\Phi}(t)\equiv\langle\sum_{i=1,2}\sum_{j\geq
3}\sum_{n}\Gamma_{i,3-i;j}^{\Phi}(n;t)\rangle,$ (26)

where $\Gamma_{2,1;j}^{\Phi}(n;y)$ is the equivalent quantity as
$\Gamma_{1,2;j}^{\Phi}$ when parton $j$ collides with the $\bar{Q}$.

Let us now consider the general situation with $N_{Q}$ $Q$-quarks as well as
$N_{\bar{Q}}$ $\bar{Q}$-quarks in the N-body system, and let us assign index
$i\in[1,N_{Q}]$ for $Q$ and $j\in[N_{Q}+1,N_{Q}+N_{\bar{Q}}]$ for $\bar{Q}$.
The total rate of quarkonia formation then writes

$\displaystyle\Gamma^{\Phi}(t)$ $\displaystyle=$
$\displaystyle\sum_{i}\sum_{j}\sum_{k>N_{Q}+N_{\bar{Q}}}\sum_{n}\left[\delta(t-t_{ik}(n))+\delta(t-t_{jk}(n))\right]$
$\displaystyle\int\prod^{N}_{l=1}d^{3}\mathbf{r}_{l}d^{3}\mathbf{p}_{l}\,h^{3}W^{\Phi}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{p}_{i},\mathbf{p}_{j})$
$\displaystyle\cdot$
$\displaystyle[W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\};t+\epsilon)-W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\};t-\epsilon)],$

where collisions between heavy quarks are neglected as they are rare. We have
to sum over all possible $Q\bar{Q}$ pairs because they can all lead to the
formation of a $\Phi$ meson after the scattering of either the $Q$
($\delta(t-t_{ik}(n))$) or the $\bar{Q}$ ($\delta(t-t_{jk}(n))$) with light
particles. One can then generalize Eq. 24 to

$W_{2}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{p}_{i},\mathbf{p}_{j},t)=\prod_{\underset{\underset{\scriptstyle
l\neq j}{\scriptstyle l\neq i}}{l=1}}^{N}\int
d^{3}r_{l}d^{3}p_{l}W_{N}^{c}(\\{\mathbf{r}\\},\\{\mathbf{p}\\},t).$ (28)

and Eq. 25 to

$\displaystyle\Gamma^{\Phi}_{ij;k}(n;t)=h^{3}\delta(t-t_{ik}(n))\int
d^{3}Pd^{3}Rd^{3}rd^{3}qW{{}^{\Phi}}(\mathbf{r},\mathbf{q})$ (29)
$\displaystyle\cdot$
$\displaystyle(W_{2}(\mathbf{R},\mathbf{r},\mathbf{P},\mathbf{q},t+\epsilon)-(W_{2}(\mathbf{R},\mathbf{r},\mathbf{P},\mathbf{q},t-\epsilon)$

as well as $\Gamma^{\Phi}_{2,1;k}(n;t)$ to $\Gamma^{\Phi}_{ji;k}(n;t)$. One
thus obtains

$\Gamma^{\Phi}(t)\equiv\langle\sum_{i}\sum_{j}\sum_{k\geq
N_{Q}+N_{\bar{Q}}}\sum_{n}(\Gamma_{ij;k}^{\Phi}(n;t)+\Gamma_{ji;k}^{\Phi}(n;t))\rangle,$
(30)

where $i$ runs from 1 to $N_{Q}$ and $j$ from $N_{Q}+1$ to
$N_{Q}+N_{\bar{Q}}$, allowing to take into account all possible $Q\bar{Q}$
pairs, independent of whether the entrained heavy quarks come from the same
vertex or from different vertices.

###  II.3 Heavy Quark - parton interactions

Our approach for the $J/\psi$ production, derived in the last section, is
based on the collisions of heavy quarks with partons from the QGP. The study
of these collisions was already presented quite a while ago Gossiaux and
Aichelin (2008); Gossiaux et al. (2009); Nahrgang et al. (2016) to investigate
the production of open heavy flavour mesons. In this study we calculate the
interaction rate for heavy quark - parton interactions and determine whether a
collision takes place by a Monte Carlo procedure. If a collision is taking
place we determine randomly, from the local equilibrium distribution, the
momentum of the parton. The interaction of the gluons and quarks with the
heavy quarks is then described by Born type matrix elements. These matrix
elements have two inputs: The running coupling constant and the infrared
regulator. The running coupling constant remains finite at zero momentum
transfer and agrees with the analysis of $\tau$ decays and $e^{+}e^{-}$
scattering Gossiaux and Aichelin (2008). The infrared regulator has been
chosen to make the result independent of the scale which separates the hard
thermal loop dominated low energy behaviour and the high momentum transfer
region which is described by Born terms. This approach has been successfully
used to describe the open heavy flavor meson production in ultrarelativistic
heavy ion collisions Gossiaux and Aichelin (2009). For the calculations
presented here we limit us to elastic collisions and employ a $K$ scaling
factor of 1.5 for the collision probabilities of heavy quarks with light
partons. As shown such a scaling factor compensates for radiative collisions,
which are not considered here Nahrgang et al. (2015).

###  II.4 Wigner density of quarkonia

To make use of Eq. II.2 we need to know the Wigner density of the quarkonia.
For quarkonia, which are created in the vacuum in pp collision, the Wigner
density has been discussed and employed in ref Song et al. (2017). Here we
follow this approach. Considering first quarkonia in non relativistic motion,
the center of mass ($\mathbf{R},\mathbf{P}$) motion of the state is given by a
plane wave. Due to the large mass of the heavy quark we assume that the
relative wave functions of the different eigenstates $\Phi$, $|\Phi\rangle$,
of the $Q\bar{Q}$ pair in vacuum, as well as its Wigner density,
$W^{\Phi}(\mathbf{r},\mathbf{p})$, can be calculated by solving the
Schrödinger equation for the relative motion in a Cornell potential Das et al.
(2016). The calculation of the $J/\psi$ production in heavy ion collisions
becomes much easier if we replace $W^{\Phi}(\mathbf{r},\mathbf{p})$ for the
s-wave states by a Gaussian Wigner density

$W^{\Phi}_{1s}(\mathbf{r},\mathbf{p})=\frac{8g}{h^{3}}e^{-\frac{r^{2}}{\sigma_{1s}^{2}}-p^{2}\frac{\sigma_{1s}^{2}}{\hbar^{2}}}$
(31) $\frac{3}{2}\sigma_{1s}^{2}=\langle r^{2}_{1s}\rangle,$

where $\sqrt{\langle r^{2}_{1s}\rangle}\approx 0.4$ fm for $J/\psi$ .
$g=\frac{3}{4}$ is the spin factor of a vector meson. The factor
$\frac{8}{h^{3}}$ is due to the normalization of the Wigner density. To
simplify the calculation we assign initially to each heavy quark-antiquark
pair whether it is in a color singlet (probability = 1/9) or in a color octet
state (probability = 8/9) and stick then to this assignment. Thus we do not
follow the color flow. This is foreseen as a future project.

###  II.5 Operational summary

We come back now to the relation between the probability and the rate.
Following Eq.4, the total probability that a $Q\bar{Q}$ pair, which has the
coordinates $\\{1,2\\}$, forms a quarkonium state at time $t$ is given by a
time integration of the rate (Eq.12)

$P^{\Phi}_{1,2}(\tilde{t})=P^{\rm
prim}_{1,2}(0)+\int^{\tilde{t}}_{t_{0}}\Gamma_{1,2}(t)dt$ (32)

with $\tilde{t}<t_{\rm hadr}$, the time when the QGP fully hadronizes. $P^{\rm
prim}_{1,2}$ is the probability that at the moment of their creation the
$Q\bar{Q}$ pair forms a quarkonium (see Song et al. (2017), as well as section
V), $t_{0}$ is the time when the QGP is formed111Assuming for the interval
$[0,t_{0}]$ no contribution to quarkonia production.. The rate (Eq. 12) is
treated in the Monte Carlo adopting the Remler method, leading to Eq.26. The
time integral of the rate thus accumulates the change of the probability that
the $Q\bar{Q}$ pairs form a quarkonium caused by all collisions with plasma
partons suffered either by the $Q$ or the $\bar{Q}$ during the time evolution
until the time $\tilde{t}$.

In the Remler formalism, equation (32), which refers only to a single
$Q\bar{Q}$ pair, naturally extends to all $Q\bar{Q}$ pairs that are in the QGP
at a given time $\tilde{t}$. Using Eq.29, the later can be expressed by
summing over all possible pair combination at a given time. As exhibited in
Eq.30, the total rate of quarkonium formation, the sum of the rate due to the
scattering of the heavy quark and that of the heavy antiquark, at a given time
can be expressed as

$\Gamma(t)=\sum^{N_{Q}}_{i=1}\sum^{N_{Q}+N_{\bar{Q}}}_{j=N_{Q}+1}\left(\Gamma_{i,j}(t)+\Gamma_{j,i}(t)\right).$
(33)

In practice, this sum over the rates is performed in the numerical program
according to Eq. 30.

We would like to stress that in our approach the $\Phi$ mesons are not
represented by pseudo particles, produced and destroyed by $2\to 2$, like
$c\bar{c}\leftrightarrow{J/\psi}+g$, or $3\to 2$, like
$Xc\bar{c}\leftrightarrow{J/\psi}X$, processes, as done in standard cascade
approaches. Instead we sum coherently the contributions to the rate of the
different $Q\bar{Q}$ pairs, which offers the advantage to add coherently all
possible contributions, what is not possible in standard MC approach based on
pseudo-particles. The non-trivial effect of adding the diagonal and off
diagonal components for the primordial contribution –$\sum_{i}\sum_{j}P^{\rm
prim}_{i,j}(0)$ – has already been discussed in Song et al. (2017).

Finally, it should be noted that the Monte Carlo implementation of the rate
can be formulated locally – see Eq. 29 – as a sum of a gain and a loss term.
If one bins the phase space along any variable (f.i. transverse momentum ${\bf
P}_{T}$), one can thus reformulate the Monte Carlo process as a depletion of
some ${\bf P}_{T}$ bin and the population of a ${\bf P}_{T}+\Delta{\bf P}_{T}$
bin, where $\Delta{\bf P}_{T}$ is the transverse momentum transferred of the
light parton to the heavy quark, its scattering partner. This opens the
possibility to evaluate differential $\Phi$-spectra by book-keeping these gain
and loss terms.

###  II.6 Generalization for relativistic quarkonia

Up to now we formulated the Wigner density in a nonrelativistic approach. As
shown in the appendix the corresponding relativisitc Wigner density can be
written as

$\displaystyle W_{i}^{\Phi}(y,\mathbf{u}_{T};\mathbf{r}^{\rm
cm},\mathbf{q}^{\rm cm})$ $\displaystyle=$
$\displaystyle\frac{\delta(y-y_{\Phi})}{(2\pi)^{3}}\delta^{2}(\mathbf{u}_{T,\Phi}-\mathbf{u}_{T})\times$
(34) $\displaystyle W_{i,\rm NR}(\mathbf{r}^{\rm cm},\mathbf{q}^{\rm cm}).$

In this expression, $u_{T}$ is the transverse component of the 4-velocity

$\mathbf{u}_{T,\Phi}=\frac{\mathbf{P}_{T}}{m_{\Phi}}$ (35)

where $\mathbf{P}_{T}$ is the total transverse momentum of the $Q\bar{Q}$
center of mass, $y_{\Phi}$ is the rapidity of the quarkonium, while
$\mathbf{r}^{\rm cm}$ and $\mathbf{q}^{\rm cm}$ are the (relative) coordinates
in the center of mass frame. $\Phi$ and the index NR indicates that the Wigner
density of the relative coordinates are evaluated in a non relativistic
framework (see Eq.(31)). This is justified because $Q\bar{Q}$ pairs with a
large relative momentum do not form quarkonia. Here it is important to mention
that only for those states, for which we impose a well-defined center of mass
4-velocity and a well defined relative momentum with respect to the center of
mass, we have been able to successfully derive a prescription which allows to
evaluate our Wigner density at any system of reference as a function of the
Wigner density in the center of mass frame. The later condition comes from the
fact that, even if we can always define a total momentum for the center of
mass, the mass of the quarkonium state $m_{\Phi}$ depends on the relative
momentum of the pair $\mathbf{q}$ as shown in the equation (78) by the on-
shell condition. This implies that one cannot, in our construction, impose
both, a fixed total momentum AND a fixed velocity for the 2-body state. To
overcome this problem one has to solve the Bethe Salpeter equation what is
beyond the scope of the present approach.

The result obtained in the equation (34) allows us to study the formation of
$J/\psi$ in the center of mass of the $c\bar{c}$ pair and at the same time be
able to evaluate the Wigner density at any time in any other system of
reference. The latter is rather important, due to the fact our multi-particle
dynamics requires to adopt a common computational frame, as can be seen from
the definition of the global rate, see Eq. (II.2). The standard computational
frame is the center of mass frame of the heavy ion collisions. Benefiting from
the boost invariance of the phase space, it is nevertheless possible to define
the equivalent Wigner density in this center of mass frame, called lab frame
to distinguish it from the center of mass frame of the $c\bar{c}$ pair, (see
equivalence between Eq. (113) and Eq. (114)) by expressing $\mathbf{q}^{\rm
cm}$ as a function of $\mathbf{q}^{\rm lab}$ as well as $\mathbf{r}^{\rm cm}$
as a function of $\mathbf{r}^{\rm lab}$ (while taking ${x^{\rm lab}}^{0}=0$ in
the Lorentz transform). This leads to a Wigner density in the lab frame:

$\displaystyle W_{i}(y,\mathbf{u}_{T},\mathbf{q}^{\rm lab},\mathbf{p}^{\rm
lab})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\delta(y-y_{\Phi})\delta^{2}(\mathbf{u}_{T,\Phi}-\mathbf{u}_{T})$
(36) $\displaystyle\cdot$ $\displaystyle W_{i,\rm NR}(\mathbf{r}^{\rm
cm}(\mathbf{r}^{\rm lab}),\mathbf{q}^{\rm cm}(\mathbf{q}^{\rm lab}))$

##  III Appropriate Basis for the Quarkonium States in the QGP

The Remler formalism was originally developed for two body systems for which
the vacuum eigenstates provide the appropriate basis. In this case the density
operator
$\rho^{\Phi}(\mathbf{r}_{1},\mathbf{r^{\prime}}_{1},\mathbf{r}_{2},\mathbf{r^{\prime}}_{2})$
in Eq. (10) corresponds to the two-body vacuum density operator. Lattice
results Mocsy and Petreczky (2007) show that the potential between the $Q$ and
$\bar{Q}$ changes with temperature and at high temperatures the quarkonia
melt. To cope with these results we introduce a temperature dependent
potential between the $Q$ and $\bar{Q}$, taken from Lafferty and Rothkopf
(2020). This renders the two body Hamiltonian temperature dependent and the
eigenstates of the relative motion of the quarkonia need to be chosen
accordingly in order to fulfil Eq.(10). We assume that also at finite
temperature the $J/\psi$ wave function can be approximated by a Gaussian. To
obtain the temperature dependence of the Gaussian width we solve the two-body
Schrödinger equation with a Lafferty-Rothkopf potential Lafferty and Rothkopf
(2020) and determine the rms radius of the $J/\psi$ wave function Katz (2015);
Katz and Gossiaux (2014). The rms radius of the $J/\psi$ wave function is
related to the Gaussian width by $\sigma^{2}(T)=\frac{2}{3}\langle
r^{2}(T)\rangle$. This calculation shows as well that the $J/\psi$ melts at
$T_{\rm diss}=0.4\,{\rm GeV}$. For $T\to 0$ the temperature-dependent
potential becomes the potential in vaccuum and therefore we recover Eq.31. The
dependence of the Gaussian width $\sigma$ on the local temperature $T$ is
displayed in Fig. 2.

Figure 2: Width of the Gaussian Wigner density, $\sigma(T)$, as function of
the temperature, obtained by solving the Schrödinger equation with a potential
taken from Lafferty and Rothkopf (2020) in the interval 0.15 GeV $<T<$ 0.4
GeV. Above $T$= 0.4 GeV = $T_{\rm diss}$, the $J/\psi$ is unstable, below
$T=0.15$ a QGP does not exist.

In an expanding QGP, the temperature changes rapidly as a function of time.
Therefore the temperature dependence of the width is equivalent to a time
dependence.

Introducing a temperature dependent potential creates an additional term in
Eq. (6) when replacing
$\rho_{\Phi}(\mathbf{r}_{1},\mathbf{r^{\prime}}_{1},\mathbf{r}_{2},\mathbf{r^{\prime}}_{2}$)
by
$\rho_{\Phi}(\mathbf{r}_{1},\mathbf{r^{\prime}}_{1},\mathbf{r}_{2},\mathbf{r^{\prime}}_{2},T(t))$.
This additional rate term we call local rate.

$\displaystyle\Gamma_{\rm eff}$ $\displaystyle=$
$\displaystyle\Gamma+\Gamma_{\rm loc}$ $\displaystyle=$ $\displaystyle
Tr[\rho^{\Phi}(\mathbf{r},\mathbf{r}^{{}^{\prime}},T(t))\dot{\rho}_{N}(t)]+Tr[\dot{\rho}^{\Phi}(\mathbf{r},\mathbf{r}^{{}^{\prime}},T(t))\rho_{N}(t)].$

The first term is the rate from the Remler formalism Eq. (6) and the second is
the new $\Gamma_{\rm loc}$.

Performing the trace integral, we obtain

$\Gamma_{\rm loc}=\sum\int
d^{3}\mathbf{r}d^{3}\mathbf{r}^{{}^{\prime}}\dot{\rho}^{\Phi}(\mathbf{r},\mathbf{r}^{{}^{\prime}},T(t))\rho_{Q\bar{Q}}(\mathbf{r},\mathbf{r}^{{}^{\prime}},t)$
(38)

where $\rho_{Q\overline{Q}}(\mathbf{r},\mathbf{r}^{{}^{\prime}},t)$is the
density operator of the N-body system integrated over the positions of the
remaining $N-2$ particles which are not part of the pair, while $\rho^{\Phi}$
is the density operator of the bound quarkonium states. The sum runs over all
possible $Q\bar{Q}$ pairs. Converting $\rho^{\Phi}$ into the corresponding
Wigner density

$\rho^{\Phi}(\mathbf{r}+\frac{\mathbf{r^{\prime}}}{2},\mathbf{r}-\frac{\mathbf{r}^{\prime}}{2})=\int
d^{3}pe^{-i\frac{\mathbf{p}\cdot\mathbf{r^{\prime}}}{\hbar}}W^{\Phi}(\mathbf{r},\mathbf{p})$
(39)

we obtain

$\Gamma_{\rm loc}=(2\pi\hbar)^{3}\int d^{3}rd^{3}p\
\dot{W}^{\Phi}(\mathbf{r},\mathbf{p},T(t))W_{Q\bar{Q}}(\mathbf{r},\mathbf{p},t).$
(40)

For classical phase space densities, Eq.(15), we arrive at

$\displaystyle\Gamma_{\rm loc}$ $\displaystyle=$ $\displaystyle
8\dot{\sigma}(T(t))\partial_{\sigma}e^{-(\frac{\mathbf{r}^{2}}{\sigma^{2}}+\frac{\sigma^{2}\mathbf{p}^{2}}{\hbar})}$
$\displaystyle=$
$\displaystyle\frac{16}{\pi^{3}}\dot{\sigma}(T(t))\left(\frac{\mathbf{r}^{2}}{\sigma^{3}(T)}-\frac{\sigma(T)\mathbf{p}^{2}}{\hbar^{2}}\right)e^{-(\frac{\mathbf{r}^{2}}{\sigma^{2}}+\frac{\sigma^{2}\mathbf{p}^{2}}{\hbar})}.$

where $\dot{\sigma}(T(t))=\dot{T}(t)\sigma^{\prime}(T)$. This local rate is
non zero if the temperature changes with time and therefore the temperature
dependent Gaussian width becomes time dependent. Including the local rate, the
probability that a quarkonium state is formed at time $t$ from a single
$Q\bar{Q}$ pair now reads

$P^{\Phi}_{Q\bar{Q}}(t)=P^{\rm init}(t^{Q,\bar{Q}}_{\rm
init})+\int^{t}_{t^{Q,\bar{Q}}_{\rm init}}(\Gamma_{{\rm
coll},Q\bar{Q}}(t^{{}^{\prime}})+\Gamma_{{\rm
loc},Q\bar{Q}}(t^{{}^{\prime}}))dt^{{}^{\prime}}.$ (42)

$\Gamma_{\rm coll}$ is the contribution of the collisions to the rate and
$P^{\rm init}$ is the initial probability for the production of a quarkonium
once the local $Q\bar{Q}$ temperature falls below $T_{\rm diss}$. Contrary to
Eq.(32), obtained from the projection on the vacuum basis, $t^{Q,\bar{Q}}_{\rm
init}$, the time when the $Q\bar{Q}$ passes $T_{\rm diss}$, is dependent on
the environment of the $Q\bar{Q}$ pair.

##  IV In-medium $Q\bar{Q}$ propagation

Besides the usual short distance collisions between heavy quarks and light QGP
partons, we have updated our MC@HQ scheme by implementing a long distance
$Q\bar{Q}$ interaction. This is detailed in the present section.

###  IV.1 Global strategy

First, let us recall that the implementation of the Remler method, described
in section II, requires trajectory calculations in Minkowski time steps
$\Delta t$. The EPOS2 event generator, which we use to simulate the time
evolution of the QGP and to obtain the initial interaction points at which the
$c$ and $\bar{c}$ quarks are produced, uses, like almost all hydrodynamical
calculations, Milne coordinates (Bjorken time $\tau^{\rm
Bjorken}=\sqrt{t^{2}-z^{2}}$ as well as the space time rapidity
$\eta=\frac{1}{2}\ln\frac{t+z}{t-z}$). The coordinate and momentum space
variables are updated employing a finite time step $\Delta\tau^{\rm Bjorken}$.
Therefore, before being able to apply the Remler algorithm, we have to
construct the trajectories of the $Q$ and $\bar{Q}$ quarks in Minkowski space
from the Milne coordinates. If $\Delta\tau^{\rm Bjorken}$ is too large and, as
a consequence, the differences in momentum of the heavy quarks between two
Bjorken time steps becomes sizable, these trajectories become kinky.

Solving the relativistic dynamics of an ensemble or particles in mutual
potential-interaction is a more involved problem and nowadays still a
challenge. A first option consists in considering retarded potentials,
generalizing the Liénard-Wiechert potential from classical electrodynamics.
This approach suffers from the non-conservation of energy and angular momentum
associated to radiative field emission. Such non-conserving features can be
cured at the price of adding the retarded and advanced propagator Wheeler and
Feynman (1949) which, however, leads to advanced interactions from the future
to the past. A second choice is the constrained Hamilton dynamics which
reduces the 8N dimensional phase space to a 6+1 dimensional by imposing time
and energy constraints Sorge et al. (1989); Marty and Aichelin (2013). This
is, however, only possible if the potential has a form, which is invariant
under a Lorentz transformation.

In the absence of an exact way of solving the multi-body dynamics, we have
adopted the following strategy: We proceed to the evolution of each quark $Q$
from $\tau^{\rm Bjorken}$ and $\tau^{\rm Bjorken}+\Delta\tau^{\rm Bjorken}$ by
considering the closest $\bar{Q}$ partners. For each such pair, we first
transform the $Q\bar{Q}$ coordinates into the CM system of the pair, where the
evolution can be performed exactly – see next subsection – until the quark $Q$
reaches times $\tau^{\rm Bjorken}+\Delta\tau^{\rm Bjorken}$. Such an evolution
leads to a variation $\delta x_{Q}$ with respect to the free propagation.

The total evolution of the $Q$ under consideration is thus defined as the sum
of the various $\delta x_{Q}$ from the interaction with the different
$\bar{Q}$ superposed to the free motion. Such a linearized algorithm leads to
acceptable results when $\Delta\tau^{\rm Bjorken}$ is small with respect to
the revolution time of the $Q\bar{Q}$ pairs.

If the distance between the quarks is large, the potential does not affect the
trajectories. Only the trajectories of neighboring $Q\bar{Q}$ pairs are
therefore concerned and usually one finds for each heavy quark $Q$ not more
than one $\bar{Q}$ (or two, in the early stage), which is sufficiently close
that the potential has an influence on the trajectory.

The drawback of such a sequence of Lorentz transformations is that small
shifts in space-time coordinates are introduced at each step due to Lorentz
transformations222Equal time $t(Q)=t(\bar{Q})$ in the CM does not generally
correspond to a unique time in the computational frame.. As those shifts are
proportional to the ${Q\bar{Q}}$ CM velocity, this approach is not suited for
the calculation of high $p_{T}$ pairs.

###  IV.2 The $Q\bar{Q}$ potential interaction

The simplest relativistic modification to the movement of a particle in a
central potential field, $V(r)$, is described by the Lagrangian Lemmon and
Mondragon (2010).

$\mathcal{L}=-\gamma^{-1}mc^{2}-V(r)$ (43)

with $\gamma^{-1}=\sqrt{1-v^{2}/c^{2}}$ and $m$ being the particle mass. With

$\frac{\partial\mathcal{L}}{\partial v_{i}}=p_{i}=\gamma mv_{i}$ (44)

and

$\displaystyle H=p_{i}v_{i}-\mathcal{L}$ $\displaystyle=$ $\displaystyle\gamma
mv_{i}v_{i}+\frac{mc^{2}}{\gamma}+V(r)=\gamma mc^{2}+V(r)$ (45)
$\displaystyle=$ $\displaystyle\sqrt{m^{2}c^{4}+p^{2}c^{2}}+V(r)=E.$

$E$, the energy, is a constant of motion. If we employ spherical coordinates
$\dot{r}^{2}=\dot{r}^{2}_{r}+r^{2}\dot{\theta}^{2}$ we find the corresponding
momenta from the Euler - Lagrange equations (we employ now c=1)

$\displaystyle p_{\theta}=\frac{\partial{\mathcal{L}}}{\partial\dot{\theta}}$
$\displaystyle=$ $\displaystyle\gamma mr^{2}\dot{\theta}$ $\displaystyle
p_{r}=\frac{\partial{\mathcal{L}}}{\partial\dot{r}}$ $\displaystyle=$
$\displaystyle\gamma m\dot{r}$ (46)

Expressing the Hamiltonian in terms of $p_{r}$ and $p_{\theta}$

$H=\sqrt{m^{2}+p_{r}^{2}+\frac{p^{2}_{\theta}}{r^{2}}}+V(r)$ (47)

we obtain the equations of motion from the Hamilton’s equation

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial
p_{r}}=\frac{p_{r}}{\sqrt{m^{2}+p_{r}^{2}+\frac{p^{2}_{\theta}}{r^{2}}}}$
$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial
p_{\theta}}=\frac{p^{2}_{\theta}r^{2}}{\sqrt{m^{2}+p_{r}^{2}+\frac{p^{2}_{\theta}}{r^{2}}}}$
$\displaystyle\dot{p}_{r}$ $\displaystyle=$ $\displaystyle-\frac{\partial
H}{\partial
r}=\frac{p_{\theta}^{2}r}{\sqrt{m^{2}+p_{r}^{2}+\frac{p^{2}_{\theta}}{r^{2}}}}-\frac{\partial
V}{\partial r}$ $\displaystyle=$
$\displaystyle\frac{p_{\theta}\dot{\theta}}{r}-\frac{\partial V}{\partial r}$
$\displaystyle\dot{p}_{\theta}$ $\displaystyle=$ $\displaystyle-\frac{\partial
H}{\partial\theta}=0\to p_{\theta}={\rm const}=L$ (48)

The last equation states that the generalized angular momentum

$L=\gamma mr^{2}\dot{\theta}$ (49)

is conserved in this ansatz. Thus one only needs to solve the radial equations
of motion on $r$ and $p_{r}$ numerically and can then integrate the
differential on $\dot{\theta}$ .

For two heavy quarks in their center of mass system we can formulate a
spinless Hamiltonian, generalizing Eq. 45

$H_{2}=\sqrt{m_{1}^{2}+p_{1}^{2}}+\sqrt{m_{2}^{2}+p_{2}^{2}}-V(r_{12})=E$ (50)

where ${\mathbf{p}}_{1}=-{\mathbf{p}}_{2}$. $m_{1}$ and $m_{2}$ are the masses
of the heavy quarks and ${\mathbf{r}}_{12}$ is their relative distance. For
the case we consider here, $m_{1}=m_{2}$, the 2-body dynamics is directly
mapped on the single one, discussed above, by taking $V(r_{12})=V(2\times
r/2)$ in Eq. 45. $V(r_{12})$ is here taken as the Lafferty-Rothkopf potential
Lafferty and Rothkopf (2020), which depends on the temperature.

###  IV.3 Calculation of $\Gamma_{\rm coll}$

At the boundaries of each time interval $[t,t+\Delta t$]we compare the
momentum change of each $c(\bar{c})$. If the $c\bar{c}$ potential is not
active, means switched off, and the heavy quark has changed its momentum in
this time interval we know that this heavy quark had a collision with a QGP
parton and we calculate the Wigner density of this quark with the antiquarks,
which are in a hydro cell with $T<T_{\rm diss}$, to determine $\Delta
W=W(t+\Delta t)-W(t)$. $\Delta W$ is then the contribution of this $c\bar{c}$
pair to the $J/\psi$ multiplicity and the sum of $\Delta W$ over all
$\bar{c}(c)$ is the contribution of this collision to the $J/\psi$
multiplicity.

If the potential is active, means switched on, the heavy quark momentum in the
$\Delta t$ interval can also be changed by the potential. This change is small
but increases with the size of $\Delta t$. To avoid to calculate the momentum
change due to the potential as a collision we add a contribution to
$\Gamma_{\rm col}$ only if

$\eta=\frac{|\mathbf{p}^{c,\bar{c}}_{n}-\mathbf{p}^{c,\bar{c}}_{n-1}|}{|\mathbf{p}^{c,\bar{c}}_{n}|}$
(51)

is above $\eta^{thres}$. $\mathbf{p}^{c,\bar{c}}_{n}$ and
$\mathbf{p}^{c,\bar{c}}_{n-1}$ are the momenta of the $c$ or $\bar{c}$ quark
at the time step $n$ and $n-1$, respectively. We find that $\eta^{thres}=0.05$
for $\Delta t$= 0.025fm/c excludes effectively the unwanted contributions from
the potential interaction.

###  IV.4 Consequences of the $Q\bar{Q}$ Interaction

Fig. 3 shows the influence of $Q\bar{Q}$ potential on the time evolution of
the $c\bar{c}$ pairs. We display there the number of $c\bar{c}$ pairs whose
constituents have an invariant distance, the relative distance between the $c$
and $\bar{c}$ quarks, measured in their CM, of $r\leq 1$ fm as a function of
the Minkowski time.

Figure 3:  Number of $c\bar{c}$ pairs whose relative distance in the CM is
lower than 1 fm as function of time under three conditions: without any
$c\bar{c}$ interaction (blue line), with the interaction screened by the
medium (dashed orange line) and with the vacuum interaction (dotted black
line).

We show this quantity for three scenarios: a) without any $c\bar{c}$
interaction (blue line), b) with the interaction screened by the medium
Lafferty and Rothkopf (2020) (dashed orange line) and c) with the vacuum
interaction (dotted black line). Without potential the number of $c\bar{c}$
pairs that stay close indeed decreases strongly with time. We see that a
medium screened $Q\bar{Q}$ potential keeps the $c\bar{c}$ pairs longer
together, especially at the final stage of the evolution ($t\approx 4-8$
fm/c). This is quite important because in a fast expanding medium the heavy
quarks tend to move away from each other. Adding an interaction potential
enhances the recombination (regeneration) process, especially at the latest
stage of QGP evolution. In the final stage $t\gtrsim 8$ fm/c, the $c$ and
$\bar{c}$ quarks escape from the QGP as mesons and are not counted anymore in
Fig. 3.

For comparison we display in Fig.3 also the (unphysical) case of an unscreened
potential.

The influence of the $Q\bar{Q}$ potential, both for the collision rate and for
the local rate, will be studied in section VI.

##  V $J/\psi$ Production in pp and Initial $c\bar{c}$ State Calibration

Although there are some well established models and formalisms available to
deal with quarkonia formation in individual pp collisions, we have chosen, for
consistency reasons, to evaluate this production in the same Wigner density
coalescence approach as the one used for AA collisions.

For this purpose, we start from the double differential spectrum eq. (114) or
its Monte Carlo equivalent eq. (C), considering that one single $Q\bar{Q}$
pair is produced, so normalizing
$\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm lab}d^{3}x^{\rm lab}_{r}}$
to unity. In our present approach, we neglect momentum correlations between
the initial $Q$ and $\bar{Q}$ quarks, although they can become significant at
large $p_{T}$. Accordingly, we express more conveniently

$\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm lab}d^{3}x^{\rm
lab}_{r}}=\mathcal{J}\,e_{Q\bar{Q}}^{\rm
lab}\,\frac{d^{3}\tilde{N}_{Q\bar{Q}}}{d^{3}x_{r}^{\rm
lab}}\times\frac{d^{3}N_{Q}}{d^{3}p_{Q}^{\rm
lab}}\times\frac{d^{3}N_{\bar{Q}}}{d^{3}p_{\bar{Q}}^{\rm lab}}$ (52)

where $\frac{d^{3}\tilde{N}_{Q\bar{Q}}}{d^{3}x_{r}^{\rm lab}}$ is the
normalized distribution to observe a distance ${\bf x}_{r}^{\rm lab}$ between
the $Q$ and $\bar{Q}$ quark. $\mathcal{J}=s+u_{T}^{2}\frac{\partial
s}{\partial u_{T}^{2}}$ is the Jacobian of the variable transformation. As the
longitudinal $z$-space is Lorentz contracted, we moreover assume that

$\frac{d^{3}\tilde{N}_{Q\bar{Q}}}{d^{3}x_{r}^{\rm
lab}}\propto\delta(z)\,e^{-\left(\frac{x_{rT}^{\rm
lab}}{\sigma_{r}}\right)^{2}}$ (53)

where $x_{rT}^{\rm lab}$ is the transverse initial distance between $Q$ and
$\bar{Q}$ quark. A last assumption is to consider factorisation of the
individual $Q$ and $\bar{Q}$ production along longitudinal and transverse
direction according to

$\frac{d^{3}N_{Q}}{d^{3}p_{Q}^{\rm
lab}}=\frac{d^{2}\tilde{N}_{Q}}{dp_{T,Q}^{2}}\frac{d\tilde{N}_{Q}}{dp_{L,Q}}=\frac{1}{e_{Q}}\frac{d^{2}\tilde{N}_{Q}}{dp_{T,Q}^{2}}\frac{d\tilde{N}_{Q}}{dy_{Q}},$
(54)

where both distributions are normalized to unity and where the transverse
spectrum is taken from the FONLL approach. From these hypothesis, it is
possible to establish

$\frac{dN_{\Phi}}{dy}\approx
P_{\Phi}\times\frac{d\tilde{N}_{Q}}{dy}\frac{d\tilde{N}_{\bar{Q}}}{dy},$ (55)

where all rapidity distributions are evaluated at the same rapidity and where
$P_{\Phi}$ represents a kind of conditional probability for a quarkonia $\Phi$
to be formed in the $Q\bar{Q}$ Wigner density coalescence, see eq.(3), where
P(t) is in this case time independent. It can be evaluated semi-analytically
once the width $\sigma_{1s}$ in eq.(31), $\sigma_{r}$ and the $p_{T}$ spectrum
of the individual quarks are specified. It is important to realize that in our
Wigner density coalescence model, the quarkonium production at a given
rapidity scales quadratically with the local abundance of heavy quarks. Next,
one obtains an equivalent relation for the production cross section in pp by
noticing a) that the conditional probability to produce a $Q\bar{Q}$ pair in
such collision is $\frac{\sigma_{Q}}{\sigma_{\rm tot}}$ (${\sigma_{\rm tot}}$
being the total pp cross section) and b) that the normalized distribution
$\frac{d\tilde{N}_{\bar{Q}}}{dy}=\frac{d\sigma_{Q}/dy}{\sigma_{Q}}$. The
average distribution per pp collision is therefore

$\frac{dN_{\Phi}}{dy}\approx\frac{P_{\Phi}}{\sigma_{\rm
tot}\sigma_{Q}}\times\frac{d\sigma_{Q}}{dy}\frac{d\sigma_{\bar{Q}}}{dy}$ (56)

leading to

$\frac{d\sigma_{\Phi}}{dy}\approx\frac{P_{\Phi}}{\sigma_{Q}}\times\left(\frac{d\sigma_{Q}}{dy}\right)^{2}\Leftrightarrow\frac{\frac{d\sigma_{\Phi}}{dy}}{\frac{d\sigma_{Q}}{dy}}=P_{\Phi}\times\frac{\frac{d\sigma_{Q}}{dy}}{\sigma_{Q}},$
(57)

where the last factor can be seen as an inverse rapidity width. The last
relation allows to calibrate the model using experimental results for $c$
quarks. In our case, we only adjust the $\sigma_{r}$ parameter because
$\sigma_{1s}=\sigma_{J/\psi}$ is constrained by the vacuum wave function to
$\sigma_{1s}=0.35\,{\rm fm}$ (see eq. 31). For $y_{\rm CM}\approx 0$, one
takes $\frac{\frac{d\sigma_{c}}{dy}}{\sigma_{c}}=0.125$, a value in agreement
with NLO calculations 333
http://www.lpthe.jussieu.fr/$\sim$cacciari/fonll/fonllform.html, while
$\frac{d\sigma_{J/\psi}}{dy}$ and $\frac{d\sigma_{c}}{dy}$ were respectively
taken as 6 $\mu$b and 1.165 mb following Acharya et al. (2019) and Acharya et
al. (2022). Tuning $\sigma_{r}$ in our MC code, we obtain the corresponding
$P_{J/\psi}$ for $\sigma_{r}=0.25\,{\rm fm}$, which is a reasonable value
according to the $m_{c}$ scale. Once the parameters are fixed, the $p_{T}$
distribution of $J/\psi$ production in pp can be calculated directly. It will
be discussed in section VI.1.

##  VI Results

As discussed in previous sections, the production and disintegration of
$J/\psi$s is a complex process. Therefore we start out with a short overview
and some definitions.

Whereas in pp collisions the $c$ and $\bar{c}$ in the $J/\psi$ come almost
exclusively from the same interaction vertex, in heavy ion reactions, when
several $c\bar{c}$ pairs are produced, this is not necessarily the case. In
our analysis we call those $J/\psi$ , which contain a $c\bar{c}$ pair from the
same vertex, diagonal $J/\psi$ , the others are called off-diagonal. In
central heavy ion collisions the system forms a QGP. If its temperature is
higher than $T_{\rm diss}$ the $J/\psi$ are not stable and a bound state
cannot be formed. At these high temperatures only unbound $c$ and $\bar{c}$
exist, which interact, however, among each other and with the QGP
constituents. We call the distribution of $c$ and $\bar{c}$ at the moment of
their production in initial hard collisions primordial distribution. The
distribution at $T_{\rm diss}$, when $J/\psi$ formation starts, is named
initial distribution. The distribution of $c$ and $\bar{c}$ at that moment
differs considerable from their primordial distribution, due to collisions of
the $c$ and $\bar{c}$ with the QGP constituents, due to the potential
interaction between the $c$ and $\bar{c}$ and due to free streaming of the $c$
and $\bar{c}$ quarks. Therefore, when applying the same $J/\psi$ Wigner
density to the primordial and to the initial distribution of the heavy quark
pairs, to determine the $J/\psi$ yield, we expect large differences.

Below $T_{\rm diss}$ the $J/\psi$ rate has two contributions: The collision
rate, which is a consequence of the $c$ or $\bar{c}$ collision with the QGP
partons, described by the Remler formalism, and the local rate, a consequence
of the change of the width of the $J/\psi$ Wigner density with temperature and
hence with time. The rates are non-zero until the QGP hadronizes or until the
$c$ or $\bar{c}$ leaves the QGP. During the hadronization of the QGP no
further $J/\psi$ will be produced. We neglect here also hadronic rescattering
of the $J/\psi$ .

From a more theoretical viewpoint it is known (see for instance Blaizot and
Escobedo (2018)) that compact white objects do not interact with the QGP. The
interaction rate of those objects in a QGP increases quadratically with their
size $r_{\rm sing}$, until $r_{\rm sing}\sim l_{\rm corr}$, the correlation
length of gluon thermal fields. From this value on, both, $Q$ and $\bar{Q}$,
interact independently with the QGP. At the same time an interaction of the
$J/\psi$ with the QGP gluons is only possible if their wave length is smaller
than $r_{\rm sing}$. Otherwise a gluon does not see the color of the color
neutral $c\bar{c}$ dipole.

One may consider that the traditional models, based on a non-interacting
initial singlet component, and our approach explore the both extreme facets of
the more involved reality. Presently only the finite value of the elliptic
flow, observed in experiments and discussed in section VII.3, presents strong
evidence that the $J/\psi$ or its constituents interact with the partons of
the QGP.

###  VI.1 Scaled Proton-Proton Production

One of the key observables in the study of $J/\psi$ production in heavy ion
collisions is the nuclear modification factor

$R_{AA}^{{J/\psi}}(p_{T})=\frac{\frac{d\sigma^{J/\psi}_{AA}}{dp_{T}}}{N_{coll}\frac{d\sigma^{J/\psi}_{pp}}{dp_{T}}}.$
(58)

$N_{coll}$ is the number of the initial hard pp collisions. The $R_{AA}$
calculation requires the knowledge of
$N_{coll}\frac{\sigma^{J/\psi}_{pp}}{dp_{T}}$. As said, we do not consider
shadowing in this paper. Therefore we do not have to simulate pp calculations
separately. Instead, to obtain
$N_{coll}\frac{d\sigma^{{J/\psi}}_{pp}}{dp_{T}dy}$, we can use the initial
$p_{T}$ and $y$ distribution of the $c$ and $\bar{c}$ quarks in AA collisions,
which is for diagonal pairs - up to a constant - identical for pp and AA.
Admitting only diagonal pairs we calculate the scaled pp production as

$N_{coll}\frac{dN^{J/\psi}_{pp}}{dp_{T}dy}=\frac{dN^{{J/\psi},{\rm
diag}}_{AA}}{dp_{T}dy}$ (59)

where the rhs is evaluated in the initial stage of the evolution with the help
of Eq. (C), with $u_{T,J/\psi}=\frac{p_{T}}{M_{J/\psi}}$, selecting $c$ and
$\bar{c}$ coming from the same vertex. We recall that in Eq. (C) $\bf{r}^{\rm
cm}$ and $\bf{q}^{\rm cm}$ are the relative distance in coordinate (momentum)
space of the $Q$ and $\bar{Q}$ in the system defined by $\\{y_{\Phi},{\bf
u}_{T,\Phi}\\}$, the rapidity and the transverse 4-velocity of the quarkonium,
while $W_{\rm NR}$ is the Wigner density defined in $Q\bar{Q}$ CM system .

We can compare our primordial diagonal A+A distribution with the
experimentally measured pp cross section on prompt $J/\psi$ production

$\frac{d\sigma^{J/\psi}_{pp}}{dp_{T}dy}\cdot T_{{AA}}=\frac{dN^{{J/\psi},{\rm
diag}}_{AA}}{dp_{T}dy}.$ (60)

For the [0-20%] centrality class, we display in Fig. 4 the midrapidity $p_{T}$
distribution of $J/\psi$ in central PbPb collisions at $\sqrt{s}=5.02\ {\rm
TeV}$. The result is compared with the prompt pp data of the ALICE
collaboration Acharya et al. (2021a), scaled by Eq.60 with an associated
nuclear overlap function $T_{\rm PbPb}$ of 20.55 ${\rm mb}^{-1}$, compatible
with the one extracted from EPOS. The data points are marked as full squares,
our MC results are given by a blue line.

Figure 4: Comparison of our model prediction for the primordial $J/\Psi$
distribution with the prompt experimental pp cross section from the ALICE
collaboration Abelev et al. (2012), scaled by $T_{PbPb}=20.55$ $mb^{-1}$.

We see that our approach reproduces quite nicely the experimental data.
Deviations are seen at large $p_{T}$. This may be due to the fact that, in
lack of having a better approach available, the $c$ and $\bar{c}$ are created
uncorrelated, in $p_{T}$ as well as in the opening angle. It is also important
to emphasize (as was mentioned in chapter 4) that our model calculates the
direct $J/\psi$ production. Therefore, one has to be cautious to compare our
results with the experimental prompt data since decay from excited quarkonium
states may contribute to the spectra. For forward rapidities the ALICE results
Acharya et al. (2021b) show that this contribution is about 15%.

In Fig. 5, we display the same analysis for the forward rapidity ($2.5\leq
y\leq 4$) data at $\sqrt{s}=2.76\ {\rm TeV}$. Since the non-prompt $J/\psi$
fraction increases from roughly 0.08 to 0.2 in this $p_{T}$ interval Aaij et
al. (2011) we expect deviations between our results for prompt $J/\psi$ and
the experimental inclusive results at higher $p_{T}$ values. We display as
well the prompt cross section, measured by the LHCb collaboration Aaij et al.
(2021), for $\sqrt{s}=5.02$ TeV.

We can conclude from this comparison that our formalism reproduces the
$J/\psi$ production in elementary pp collisions. Thus we confirm the results
obtained in Song et al. (2017).

Figure 5: Comparison of our model prediction for the primordial $J/\Psi$
distribution in PbPb with the inclusive experimental pp cross section from the
ALICE collaboration Abelev et al. (2012) as well as with the prompt pp cross
section, measured by the LHCb collaboration Aaij et al. (2021) (for
$\sqrt{s}=5.02$ TeV), scaled by $T_{PbPb}=20.55$ $mb^{-1}$.

###  VI.2 Heavy Ion Collisions

####  VI.2.1 Initial Distribution

In heavy ion collisions the $J/\psi$s, which would be produced in the initial
hard NN collisions in absence of a medium, are not observed finally because
they dissolve into $c$ and $\bar{c}$ quarks when they pass the (high
temperature) QPG. $c$ and $\bar{c}$ can only form a stable $J/\psi$ when the
QGP temperature falls below $T^{{J/\psi}}_{\rm diss}$. In practice we apply
the following description: When a $c$-quark arrives for the first time in a
region with $T<{T^{{J/\psi}}_{\rm diss}}$ we calculate its probability to form
a $J/\psi$ with all $\bar{c}$, which are already satisfying this condition:

$N_{k}^{{J/\psi},{\rm init}}=\sum^{n_{\bar{c}}}_{l=1}W_{{\rm
NR}}^{J/\psi}(r^{cm}(k,l),q^{cm}(k,l)).$ (61)

$n_{\bar{c}}$ is the number of active charm anti-quarks (means from a region
of the QGP with $T<{T^{{J/\psi}}_{\rm diss}}$). $k$ is the index of the
$c$-quark, which has passed the dissociation temperature at time t. l is the
index of a $\bar{c}$ quark which is active. $r^{cm}(k,l)(q^{cm}(k,l))$ stand
for the relative distance in coordinate (momentum) space of the kl pair in the
pair center of mass system. The sum of all these contributions for all
$c$-quarks (and analogously those for all $\bar{c}$ quarks) is the initial
$J/\psi$ distribution

We can also define the initial rapidity distribution of the $J/\psi$ s, which
contain a $c$ or $\bar{c}$ quark, which passed $T_{\rm diss}$ between t and dt

$\displaystyle\frac{dN^{\rm init}}{dy}(t,t+{\rm d}t)=\sum^{N_{\rm
first}(t,t+{\rm d}t)}_{k=1}\sum^{n_{\bar{c}}(n_{c})}_{l=1}$ (62)
$\displaystyle\times$ $\displaystyle\int
d^{2}u_{T}W(y,u_{T},r^{cm}(k,l),q^{cm}(k,l))$

where $N_{\rm first}(t,t+{\rm d}t)$ stands for the number of $c$ or $\bar{c}$
quarks, which passed the temperature threshold $T_{\rm diss}$ between t and
$t+{\rm d}t$

The time evolution of the initial $J/\psi$ production, $\frac{dN^{\rm
init}}{dy}$, is shown in Fig.6. We display this quantity, normalized to
$(\frac{dN^{c}(t)}{dy})^{2}$ at midrapidity for different centrality intervals
and for the reaction PbPb at $\sqrt{s}=5.02\ {\rm TeV}$ (black dashed line for
[0-20%], red dashed-dotted line for [20-40%], olive full line for [30-50%]).
In a rate equation the creation of $J/\psi$ would be proportional to this
normalization in a system of fixed volume and we would expect a time
independent production. Deviations show the influence of the expansion.

Figure 6: Time evolution of the $J/\psi$ initial production for different
centralities and for two scenarios: without collisions of $c$ and $\bar{c}$
with QGP partons (top) and including these collisions (bottom). Figure 7: Time
integrated initial $p_{T}$ distribution of $J/\psi$ in midrapidity [0-20%]
central PbPb collisions at $\sqrt{s}=5.02$ TeV, obtained for the four
different evolution conditions combining $c\bar{c}$ potential interaction
($ON$ and $OFF$) and elastic collisions ($ON$ and $OFF$) with QGP partons.

We see a wide distribution of times when a $c$ or $\bar{c}$ quarks can form a
$J/\psi$ for the first time, in mid-central as well as in central collisions.
In central collisions the QGP is larger and therefore present for a longer
time. Therefore the distribution is shifted to later times. These
distributions are rather independent on whether the potential interaction
between the $c$ and $\bar{c}$ is active or not but depend quite strongly on
whether we admit collisions of the heavy quarks with the QGP partons. As we
will see below, these collisions lower the $p_{T}$ momentum of the heavy
quarks and therefore decelerate the expansion. In addition, heavy quarks with
lower momenta have a higher chance to form a $J/\psi$ .

We come now to the $p_{T}$ distribution of the initial $J/\psi$s. It is
displayed for $|y|\leq 0.9$ and for [0-20%] central PbPb collisions at
$\sqrt{s}=5.02\ {\rm TeV}$ in Fig. 7. We display 4 different scenarios to show
the consequences of the collisions of $c$ and $\bar{c}$ with the QCD
constituents and of the presence of the potential between the $c$ and
$\bar{c}$ . The black dashed line shows our result if neither collisions occur
nor the $c\bar{c}$ potential is active. For the blue long dashed line we
switched on the $c\bar{c}$ potential. The orange dashed dotted and the red
full line show the results without and with $c\bar{c}$ potential if collisions
between the heavy quarks and QGP partons are admitted.

The potential has little influence on the inital $p_{T}$ spectrum but
collisions shift the $J/\psi$ distribution to lower $p_{T}$ values, which
corresponds to a $c$($\bar{c}$) quenching between the time when the
$c$($\bar{c}$) enters the QGP and the ”initial” time (when $T=T_{\rm diss})$.

####  VI.2.2 Impact of the potential between the $c$ and $\bar{c}$ quarks

As discussed in section 3, the $c$ and $\bar{c}$ quarks interact via a
potential interaction whose parameters are taken from Lafferty and Rothkopf
(2020). We calculate the heavy quark trajectories in Minkowski space using the
equations of motion of Eq. 48 after boosting the pair into their cm frame. To
show the influence of this potential we display in Fig. 8,

Figure 8: Number of $c\bar{c}$ pairs that fulfil the condition of $r^{\rm
cm}_{\rm rel}\cdot p^{\rm cm}_{\rm rel}<\hbar$ as a function of
$p_{T}^{J/\psi}$, integrated over time.

as function of the total transverse momentum $p_{T}^{J/\psi}$ of the
$c\bar{c}$ pair, the number of those pairs, which have been in phase space
closer than $r^{\rm cm}_{\rm rel}\cdot p^{\rm cm}_{\rm rel}=\hbar$. This is
studied for central reactions of PbPb at $\sqrt{s}=5.02\ {\rm TeV}$ and at
midrapidity, $|y|\leq 0.9$. The color coding is the same as in Fig.7. When
collisions of the $c$ and $\bar{c}$ with the QGP partons are admitted, the
number of close pairs increases strongly at low $p_{T}$ due to the momentum
change of the quarks in these collisions (see Fig. 7). At high $p_{T}$ we see
less pairs if the collisions are active because the number of $c$ and
$\bar{c}$ quarks there is strongly decreased due to the momentum shift. If we
allow furthermore for a potential interaction between the $c$ and $\bar{c}$
quarks we further enhance the number of close pairs by roughly a factor of 2.
$c$ and $\bar{c}$ pairs with a small distance in phase space have a higher
probability to form a $J/\psi$ . So we can conclude that both, potential as
well as collisions, contribute to an enhancement of the $J/\psi$ production
with respect to free streaming. Both effects nearly factorize.

####  VI.2.3 Impact of the local rate

In an expanding QGP the temperature decreases as a function of time. Therefore
the temperature dependent Wigner density for $J/\psi$ , which we employ, is
time dependent. This leads in the Remler formalism, extended to temperature
dependent eigenstates (see Eq.(III)), to a local production rate for $J/\psi$
,

$\Gamma_{\rm loc}=\frac{dT}{dt}\frac{d\sigma}{dT}\frac{dW}{d\sigma}$ (63)

where $\sigma$ is the width of the Gaussian Wigner density $W_{{\rm
NR}}^{J/\psi}(r^{cm}(k,l),q^{cm}(k,l))$ and T is the temperature of the QGP
region in which the $J/\psi$ is located at the time $t$. The $c$ and $\bar{c}$
may be in regions of slightly different temperature although the most relevant
$c\bar{c}$ contributions to $\Gamma$ are those for which $T_{c}\approx
T_{\bar{c}}$. In our calculation we take the average value

$T_{c\bar{c}}=\frac{1}{2}(T_{c}+T_{\bar{c}}).$ (64)

In the numerical program we use a fixed time step $\Delta t$. Therefore we
replace

$\frac{d\sigma(t)}{dt}\approx\frac{\Delta\sigma(t)}{\Delta t}$ (65)

where $\Delta t$ is the computational Minkowski time step measured in the
center of mass of the pair. This $c\bar{c}$ pair contributes to the final
$J/\psi$ multiplicity with

$N^{\rm loc}_{\rm pair}=\int^{\infty}_{t^{\rm first}}\Gamma^{c\bar{c}}_{\rm
loc}(t^{{}^{\prime}})dt^{{}^{\prime}}$ (66)

where $t^{\rm first}$ is the time in which the later of the two ($c$ or
$\bar{c}$) passes $T_{\rm diss}$. For the multiplicity we sum over all
$c\bar{c}$ pairs.

Figure 9: Integral (time integrated) local rate as a function of $p_{T}$ for
$J/\psi$ production in PbPb in (0-20%) central PbPb collisions at $\sqrt{s}$ =
5.02 TeV for four different evolution scenarios: $c\bar{c}$ potential
interaction ($ON$ and $OFF$) and elastic collisions with QGP partons ($ON$ and
$OFF$). When neither the potential nor the collisions are active the rate is
negative (and thus not displayed in this semilog plot). The color coding is
the same as in Fig.7

.

In Fig. 9 we display the $p_{T}^{J/\psi}$ dependence of the integrated local
rate for [0-20%] central PbPb collisions at $\sqrt{s}=5.02\ {\rm TeV}$ and for
the different scenarios discussed in Fig. 7. The color coding is the same as
in Fig.7. The Wigner density gets larger with decreasing distance between the
$c$ and $\bar{c}$ . The potential interaction brings $c$ and $\bar{c}$ quarks
closer together (see Fig.8) and therefore it is expected that the local rate
contribution gets larger when the potential interaction between the heavy
quarks is active. This is indeed seen in the calculations. We see as well that
the collisions of the $c$ and $\bar{c}$ quarks with the QGP partons influences
strongly the local rate. These collisions shift the transverse momentum
distribution of the heavy quarks towards lower values (Fig. 7) and therefore
also the relative momentum between the $c$ and $\bar{c}$ , which enters the
Wigner density, gets smaller. The potential interaction does not change the
form of the integrated local yield (which does not depend on the centrality),
only the absolute value. In the scenario without potential and without
collisions the local rate is negative.

####  VI.2.4 Impact of the collision rate

Whenever a $c$ or $\bar{c}$ quark collides at time t with a parton of the QGP,
the heavy quark changes its momentum and therefore for all $c\bar{c}$ pairs,
which include that scattered heavy quark, the overlap with the $J/\psi$ Wigner
density is different before and after the collision. These collisions lead
only to the production of $J/\psi$ if the QGP around the heavy quark has a
temperature $T<T_{\rm diss}$. In this case the difference $\Delta
W=W(t+\epsilon)-W(t-\epsilon)$ is the contribution of this collision to the
integrated collision rate of $J/\psi$s. We count these collisions if
$\eta>0.05$ (see Eq.(51)).

Fig. 10 shows how the integrated collision rate depends on the final center of
mass momentum of the $J/\psi$ . We display this rate for central [0-20%] PbPb
collisions at $\sqrt{s}=5.02\ {\rm TeV}$ for all $J/\psi$ which have at the
end of the reaction a rapidity $|y|<0.9$. We use the same coding as for Fig.
7. We see, first of all, that the integrated rate is zero or negative. This
means that collisions reduce the multiplicity of $J/\psi$s. For small values
of $p_{T}$ the momentum transfer is limited by kinematics and therefore the
momentum of the heavy quark changes little. For very large $J/\psi$ momenta
the cross section leads to a small angle scattering what also limits the
possible momentum transfer. In both cases $\Delta W$ is small and so the
contribution to the integrated rate is close to zero.

The integrated collision rate is obviously zero if we admit neither collisions
nor potential. If we allow only for potential interactions we see only a
spurious rate caused by the potential (see discussion in section IV.3). It is
much smaller than the integrated rate if collisions are active. Without
potential interaction the integrated collisional rate is most negative. The
$Q\bar{Q}$ potential keeps the $c$ and $\bar{c}$ closer together hence the
overlap with the $J/\psi$ Wigner density gets larger with the consequence of a
less negative integrated collisional yield.

Figure 10:  $(-1)\times\frac{dN_{\rm coll}}{dp_{T}}$, integrated over the
collision rate, for a collision coefficient $\eta_{cut}=5\%$ and for the four
different evolution scenarios combining $c\bar{c}$ potential interaction ($ON$
and $OFF$) and elastic collisions ($ON$ and $OFF$) with QGP partons.

####  VI.2.5 Total Multiplicity

Adding the local and collisional rates as well as the initial production we
can study the net-multiplicity, Eq. (42), the difference between production
and disintegration, of the $J/\psi$ par unit rapidity,
$\frac{dN^{J/\psi}}{dy}$. In Fig.11 we display the net-multiplicity as a
function of time for three different centralities [0-20%], black dashed line,
[20-40%], brown long dashed line, and [40-60%], magenta dashed dotted line. On
top we display this quantity if the potential interaction between $c$ and
$\bar{c}$ is not active, in the bottom figure it is included. We see that the
production starts early, when the first $c$ and $\bar{c}$ pass $T_{\rm diss}$
and increases strongly until around 1 fm/c. Without potential at later times
there are more $c$ and $\bar{c}$ dissociating than produced as a consequence
of the negative $\Delta W$, displayed in Fig. 10. For an active $c\bar{c}$
potential the $c$ and $\bar{c}$ stay closer together (Fig. 8) with the
consequence of a steadily increasing the yield until 10 fm/c, in agreement
with the less negative integrated collision rate, see Fig.10 . This increase
is rather independent of the centrality of the collision.

Thus our results do not support an instantaneous coalescence mechanism during
hadronization, which is the basis of statistical models at fixed volume and
which is also frequently applied to describe the total or part of the $J/\psi$
yield in heavy ion collisions. This is one of the key messages of our
approach, that we expect to hold irrespective of the calculation details.

Figure 11:  Time evolution of the $J/\psi$ production $\frac{dN}{dy}$,
integrated over $p_{T}$, as a function of time for different centrality bins
and for two different scenarios: $c\bar{c}$ interaction potential $ON$ (top)
and $OFF$ (bottom) for $\eta_{cut}$=5%, Minkowski and Bjorken time steps 0.25
fm/c and 0.1 fm/c, respectively.

##  VII Comparison with Experiment

In this section we compare the production of $J/\psi$ , obtained from Eq. 42,
with the corresponding experimental results and make predictions where data
have not been published yet.

###  VII.1 $p_{T}$-Spectrum of $J/\psi$

The calculated midrapidity $p_{T}$ spectra, for $|y|<0.9$, of $J/\psi$
produced in PbPb collisions at $\sqrt{s}=5.07\ {\rm TeV}$ is displayed in Fig.
12 for three different centralities [0-20%], [20-40%] and [40-60%].

Figure 12: Final $J/\psi$ $p_{T}$ spectrum, $\frac{dN^{J/\psi}}{dp_{T}dy}$, at
midrapidity obtained for different centrality bins. The results were obtained
with the standard values of the parameters: $\eta_{cut}=5\%$, Minkowski time
step value 0.25 fm/c and Bjorken time step 0.1 fm/c.

Here collisions with the QGP partons as well as the $Q\bar{Q}$ potential
interactions are included. For all centralities the form of the curve is
rather similar. We observe that at high $p_{T}$ the slope of the curves are
rather similar. In central collisions the maximum is shifted towards lower
$p_{T}$ values due to the increasing number of collisions of the heavy quarks.

Figure 13: Final $p_{T}$-spectra at mid-rapidity for [0-20%] central PbPb
collisions, for the four scenarios (switching on/off interaction potential and
in-medium elastic collision) and for our standard value of the collision
coefficient $\eta_{\rm cut}=5\%$, and Minkowski and Bjorken time steps of 0.25
fm/c and 0.1 fm/c, respectively. The full squares mark the ALICE experimental
data Acharya et al. (2020b).

In Fig. 13 we investigate in detail how the different ingredients of our model
influence the final $p_{T}$ distribution at midrapidity for [0-20%] central
PbPb collisions at $\sqrt{s}=5.07\ {\rm TeV}$. We present the results for the
four possible combinations if we activate/deactivate collisions and potential.
The color code corresponds to that in Fig. 7. We observe that if collisions
are active the maximum of the $J/\psi$ distribution is shifted to lower values
of $p_{T}$ as a consequence of the shift of the heavy quarks spectra under
this condition. The $c\bar{c}$ potential enhances the yield without changing
the high $p_{T}$ slope. In Fig. 13 we compare as well our results with the
data of the ALICE collaboration Acharya et al. (2020b), which are shown as
black points. We see that our results are close to the experimental data, if
both, potential and collisions, are active,

###  VII.2 $J/\psi$ nuclear modification factor

One of the most interesting observables is the nuclear modification factor,
$R_{AA}$, Eq. 58. Its deviation from unity shows how the $p_{T}$ spectra are
modified by nuclear effects. With the proton reference spectrum discussed in
subsection 5.A we present the results of our model in comparison with the
ALICE data Acharya et al. (2020b) in Fig. 14 for PbPb collisions at
$\sqrt{s}=5.02\ {\rm TeV}$. The ALICE data present the results for inclusive
$J/\psi$ whereas we calculate only the directly produced $J/\psi$ . B meson
decay as well as the decay of excited states, like the $\psi^{\prime}$ and
$\chi$, contribute to the experimental $J/\psi$ distribution and therefore the
comparison between our results and the data has to be made with caution. Fig.
14 shows $R_{AA}(p_{T})$, on top for central ([0-20%]) on bottom for mid
central [20-40%] collisions. Our results are presented as red lines, the ALICE
data as black points. In the top figure we display as well the $R_{AA}$
values, which we obtain for the different combinations of switching on/off the
$c\bar{c}$ potential and the collision of the heavy quarks with the QGP
partons. We see that our results, as the experimental data, show an
enhancement at small $p_{T}$. For large $p_{T}$ $R_{AA}$ decreases in the
calculations as in the data. We recall that in our relativistic trajectory
calculation we include only the first order (in $\gamma$-1) correction.

Figure 14: Comparison between our model prediction for the nuclear
modification factor $R_{AA}$ for two different centrality ranges, $0-20\%$
(top) and $20-40\%$ (bottom), and the corresponding inclusive experimental
data from the ALICE collaboration Acharya et al. (2020b).

It is remarkable that we obtain in AA collisions a $R_{AA}<1$. As explained in
the last section, in AA collisions we have off-diagonal contributions which
dominate the primordial multiplicity (see Song et al. (2017)), so naively one
would expect a strong enhancement. The reason for the suppression is that in
AA collisions the $J/\psi$ are created later (when $T<T_{\rm diss}$), where
the average distance between $c$ and $\bar{c}$ is larger and therefore the
overlap with the Wigner density is smaller. Despite collisions and $c\bar{c}$
potential, which enhance the yield, this overcompensates the primordial
enhancement.

This is demonstrated in Fig. 15 where we compare the initial and the
primordial $J/\psi$ production as a function of $p_{T}$. We display the ratio,
called $R_{AA}^{prim}$, of $J/\psi$s obtained from the $c\bar{c}$ which have
passed $T_{\rm diss}$ and the diagonal $J/\psi$s produced in the initial hard
collisions. The ratio is displayed for the four different scenarios with the
same color coding as in Fig. 7. The ratio is, despite of the off-diagonal
contribution, smaller than one for almost all $p_{T}$ values. Collisions
enhance $R_{AA}^{prim}$ at small $p_{T}$ and lower the ratio at large $p_{T}$,
whereas the potential has little influence on the ratio.

Figure 15: $R_{AA}$ as a function of the transverse momentum at the initial
time, when the charm quarks are below $T_{\rm diss}$, at midrapidity,
$|y|<0.9$ for [0-20%] central PbPb collisions at $\sqrt{s}$ = 5.02 TeV. We
display this result for the 4 possible combinations of the set up. For details
see text.

The importance of the off diagonal contribution to the $J/\psi$ yield is
demonstrated in Fig.16. It shows for the same reaction and for central
collisions the total yield (short dashed black line) and the diagonal
contribution (dashed blue line).

Figure 16: Comparison of our model prediction for the final yield (short
dashed black line) at midrapidity, $|y|<0.9$, for [0-20%] central PbPb
collisions at $\sqrt{s}$ = 5.02 TeV. We display also our model prediction for
the diagonal contribution (dashed blue line). For this calculation medium
elastic collision and $c\bar{c}$ potential are included and we employ our
standard parameter values: collision coefficient $\eta_{\rm cut}=5\%$ and
Minkowski and Bjorken time steps 0.25 fm/c and 0.1 fm/c, respectively.

For low $p_{T}$ the diagonal part of the yield is up to one order of magnitude
smaller than the total yield, so in most of the $J/\psi$s the two heavy quarks
come from different vertices. At low $p_{T}$ this is a consequence of the
observation that low $p_{T}$ heavy quarks come to thermal equilibrium with the
QGP Nahrgang et al. (2014) at least in azimuthal direction and therefore the
ratio of diagonal to off-diagonal contribution is determined by statistics.
With increasing $p_{T}$ the

Figure 17: Nuclear modification factor as a function of the average number of
participants in PbPb obtained for active interaction potential and medium
elastic collisions, employing Bjorken and Minkwoski time steps of 0.1 fm/c and
0.25 fm/c, respectively, and $\eta_{cut}=5\%$. We compare the results with
ALICE data (black squares Acharya et al. (2020b).

off-diagonal contribution decreases but even at the largest $p_{T}$,
investigated in this study, it does not fall below 50% of the total yield.

Fig. 17 shows the centrality dependence of $R_{AA}$ for central PbPb
collisions at midrapidity and at $\sqrt{s}=5.02\ {\rm TeV}$. The ALICE data
Acharya et al. (2020b) are presented as black points, the result of our
calculation as a red line. We observe in theory as well as in experiment an
enhancement of $R_{AA}$ for central collisions, where the number of produced
$c$ and $\bar{c}$ is large and therefore recombination is more probable, as
well as a decrease with decreasing centrality. In peripheral reactions, which
resemble pp collisions, $R_{AA}$ returns to one.

###  VII.3  Elliptic Flow of $J/\psi$

The azimuthal distribution of the transverse momentum distribution can be
expanded in a Fourier series

$\frac{d^{2}N}{d^{2}p_{T}}=\frac{1}{2\pi
p_{T}}\frac{dN}{dp_{T}}\left(1+2\sum^{\infty}_{n=1}v_{n}\cos(n(\phi-\Psi_{RP}))\right),$
(67)

where $\phi$ is the azimuthal angle of the $J/\psi$ and $\Psi_{RP}$ is the
angle of the reaction plane. The elliptic flow, $v_{2}$, the second
coefficient of the expansion, can be expressed (if the xz plane is the
reaction plane) as

$v_{2}=\langle\frac{p^{2}_{x}-p^{2}_{y}}{p^{2}_{T}}\rangle.$ (68)

$v_{2}$ is another key observable in heavy quark physics. The eccentricity of
the almond shaped interaction region in coordinate space is converted, during
the hydrodynamical expansion of the QGP, into an eccentricity in momentum
space and hence into a finite $v_{2}$ value. The production of the initial
heavy quarks in hard collisions is azimuthally isotropic and hence initially
$v_{2}=0$. The $v_{2}$ value, observed for final $J/\psi$s, is therefore a
measure of their interaction (or that of their predecessors, the $c$ and
$\bar{c}$ quarks) with the QGP because only in collisions with the medium they
can acquire a finite $v_{2}$. This is true for moderate $p_{T}$. At higher
$p_{T}$ values the dependence of the path length in the medium on the
azimuthal angle starts to play the leading role.

Figure 18: Elliptic flow, $v_{2}$, in the [30-50%] centrality range between
for PbPb at $\sqrt{s}$ = 5.02 TeV and for mid-rapidity ($|y|\leq 0.9$. We
compare our results, employing the standard parameters, with the data from the
ALICE collaboration Acharya et al. (2017).

In Fig. 18 we compare our results for the standard parametrization and
including collisions and potential (red line) for $|y|<0.9$ with midrapidity
data from the ALICE collaboration Acharya et al. (2020a) (black squares) for
the reaction PbPb at $\sqrt{s}=5.02\ {\rm TeV}$ and for the [30-50%]
centrality interval. Our calculation reproduces the experimental data but the
error bars are large.

Figure 19: Comparison between our model predictions for the elliptic flow,
$v_{2}$, at mid-rapidity ($|y|\leq 0.9$) with experimental data from the ALICE
collaboration for the rapidity ($2.5<y<4.0$) Acharya et al. (2020a). On top
(bottom) we show the centrality ranges [20-40%] ([30-50%]).

In Fig. 19 we compare our results for $|y|<0.9$ and for different centrality
bins with the results of the ALICE collaboration for $2.5<y<4$. The top figure
shows the results for [20-40%], the bottom figure those for [30-50%]
centrality. This bin corresponds to the centrality bin shown in Fig. 18.
Comparing both figures, we see that $v_{2}$ at mid and forward rapidity is
rather similar. Therefore it is meaningful to compare the midrapidy
calculations with data at forward rapidity. We include in this figure, as blue
dotted line, $v_{2}$ of the diagonal $J/\psi$ , those where the $c$ and
$\bar{c}$ come from the same elementary vertex. We see that they have a
considerably smaller $v_{2}$, as also seen in the recent calculation by He et
al. He et al. (2022).

We can conclude that our approach reproduces quite well the unexpected large
$v_{2}$ values, observed in experiment, if we include the non-diagonal
$J/\psi$ . $J/\psi$s with a $c$ and $\bar{c}$ from the same elementary vertex
have a considerably lower $v_{2}$.

##  VIII Conclusions

In this work, we presented a new theoretical approach to understand the
experimental data on $J/\psi$ production in pp and AA collisions at LHC
energies. We follow the individual $c$ and $\bar{c}$ quarks from their
creation in initial hard collisions until their final observation in the
hidden heavy flavor mesons. Our treatment is thus in the same spirit as recent
open quantum system approaches, which explicitly study the time evolution of
the quarkonia density matrix under the influence of an approximated density
matrix of the whole system.

For pp collisions, as in ref. Song et al. (2017), the production of $J/\psi$s
is described by a sudden approximation which has been found to give a good
description of the experimental findings, not only for $J/\psi$ and
$\psi^{\prime}$ but also for $\chi_{c}$. In this approach the initial
distribution of $c$ and $\bar{c}$ is projected on the Wigner densities of the
quarkonia states.

This primordial distribution of $J/\psi$ and its excited states are of no
relevance for the $J/\psi$ production in central heavy ion collisions. There a
QGP is produced. Lattice gauge calculations reveal that if its temperature is
above $T=T_{diss}$, $J/\psi$s are not stable and that below $T_{diss}$ but
finite T, the $J/\psi$ wave function is quite different from the vacuum wave
function due to the interactions of the $c$ and $\bar{c}$ with the QGP
environment. This renders it more complicated to model the dynamics of
$J/\psi$sin heavy ion collisions.

To cope with this observation we

  * •

employ a mutual potential interaction of the $c\bar{c}$ pairs with a
potential, which is adjusted to the lattice data. It is active between all
$c\bar{c}$ pairs in singlet states.

  * •

modify below $T_{diss}$ the Wigner density of the $J/\psi$ in the QGP medium
by introducing a temperature dependent width in the Gaussian parametrization
of the $J/\psi$ Wigner density. This width reproduces the $J/\psi$ in medium
radius given by a potential model based on lattice data.

For the collisions of the $c$ or $\bar{c}$ with the QGP constituents, the
quarks and gluons, we use the model which has already been successfully
applied for the studies of open charm mesons ref.Gossiaux and Aichelin (2008);
Gossiaux et al. (2009). They change the Wigner density of the $c$ or $\bar{c}$
quark which was involved in that collisions and therefore the convolution of
the $c\bar{c}$ Wigner density with the Wigner density of the $J/\psi$ changes
as well.

In central heavy ion collisions one would expect a strong enhancement of the
$J/\psi$ multiplicity as compared to the multiplicity in a single pp collision
multiplied by the number of these collisions in the heavy ion reaction. The
reason is that in heavy ion collisions $c$ and $\bar{c}$ from different
elementary vertices can form a $J/\psi$ . This enhancement is, however,
(over)compensated by the late production of $J/\psi$s. $J/\psi$s can only be
produced when the plasma temperature has fallen below $T_{diss}$. There the
average distance between the $c$ and $\bar{c}$ is considerable larger and
therefore the convolution with the Wigner density is smaller.

We observe furthermore that the collisions between QGP partons as well as the
potential interaction between $c$ and $\bar{c}$ quarks enhance the $J/\psi$
yield as compared to a free streaming scenario. The collisions shift the
$p_{T}$ distribution of heavy quarks towards lower values, the potential keeps
$c\bar{c}$ pairs closer together. Both processes shift therefore the two-body
Wigner density of the $c\bar{c}$ pairs to regions where the $J/\psi$ Wigner
density is large.

We find in our approach a quite satisfying agreement with the experimental
data, simultaneously for $R_{AA}$ and for $v_{2}$. The latter is created by
the $c$ and $\bar{c}$ , which have not yet formed a $J/\psi$ , due to their
collisions with the QGP partons.

We have employed in our approach a simple model for the color degrees of
freedom and we have concentrated on $J/\psi$ mesons. It should be noted that
the gluo-dissociation mechanism $g+\Phi\to Q+\bar{Q}$ as well as its detailed
balance counterpart – expected to become significant for deeply bound states,
thus around $T_{c}$ – were not included in our dynamical treatment. To improve
on these aspects will be the subject of an upcoming publication. The new EPOS4
approach will also allow to treat correctly the correlations between the
initially formed $c\bar{c}$ pairs.

## Acknowledgement

The authors thank Taesoo Song for providing the program for the calculation of
$J/\psi$ in pp collisions as well as Taesoo Song and Elena Bratkovskaya for
inspiring discussions. This study is part of a project that has received
funding from the European Union’s Horizon 2020 research and innovation program
under grant agreement STRONG – 2020 - No 824093. We are also pleased to
acknowledge the support from the Region Pays de la Loire, under contract n°
2015-08473.

## Appendix

The HQ are created in the early phase of the HIC through hard scattering
processes and can have relativistic energies in the computational frame, the
nucleus-nucleus center of mass frame. It is crucial that the Wigner function,
used in our coalescence approach, is adequate to describe HQ under these
conditions.

To obtain such a Wigner density we need first to build an orthogonal set of
states that represents the bound states of $Q\bar{Q}$ and to ensure the
correct normalization and relativistic invariance. From that we can then
proceed to derive the Wigner function associated with the $Q\bar{Q}$ states,
expressed in this orthogonal basis.

##  Appendix A Finding Orthogonal Two-Particle Basis for non-relativistic
states

Given a longitudinal direction and a transverse plane, a quantum state for a
two-particle non relativistic system (two scalar particles) can be written as

$|\Phi\rangle=\int\frac{dp_{L,1}d^{2}\vec{p}_{T,1}}{2E_{1}}\frac{dp_{L,2}d^{2}\vec{p}_{T,2}}{2E_{2}}f(\vec{p}_{1},\vec{p}_{2})|1\rangle|2\rangle$
(69)

where the single particle states are normalized to

$\langle
1^{\prime}|1\rangle=2E_{1}\delta(\vec{p}_{1}-\vec{p}\,^{\prime}_{1}).$ (70)

States with a defined center of mass momentum $\vec{P}$,
$|\Psi_{\vec{P}}\rangle$, should ideally be normalized as

$\langle\Phi_{\vec{P}^{{}^{\prime}}}|\Phi_{\vec{P}}\rangle=2E\delta(\vec{P}-\vec{P}^{{}^{\prime}})$

For eq. 69 this imposes

$\int\frac{d^{3}p_{1}}{2E_{1}}\frac{d^{3}p_{2}}{2E_{2}}f_{P}^{*}(\vec{p}_{1},\vec{p}_{2})f_{P^{\prime}}(\vec{p}_{1},\vec{p}_{2})=2E\delta(\vec{P}-\vec{P}^{{}^{\prime}})$
(71)

We now consider the correlation function

$f_{i,\vec{P}}(\vec{p}_{1},\vec{p}_{2})=\delta(\vec{P}-\vec{p}_{1}-\vec{p}_{2})f_{i}(\frac{\vec{p}_{1}-\vec{p}_{2}}{2})$
(72)

where the index $i$ refers to possible internal states. Then, with
$\vec{q}=\frac{\vec{p}_{1}-\vec{p}_{2}}{2}$, the normalization reads:

$\displaystyle\langle\Phi_{i,P^{{}^{\prime\prime}}}|\Phi_{j,P^{\prime}}\rangle$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{2E_{1}}\frac{d^{3}p_{2}}{2E_{2}}\delta({\vec{P}}-{\vec{p}}_{1}-{\vec{p}}_{2})\delta({\vec{P}}^{{}^{\prime}}-{\vec{p}}_{1}-{\vec{p}}_{2})$
(73) $\displaystyle\times$ $\displaystyle
f_{i}^{*}(\frac{{\vec{p}}_{1}-{\vec{p}}_{2}}{2})f_{j}(\frac{{\vec{p}}_{1}-{\vec{p}}_{2}}{2})$
$\displaystyle=$
$\displaystyle\delta({\vec{P}^{\prime}}-{\vec{P}^{\prime\prime}})\int\frac{d^{3}q}{4E_{1}E_{2}}f^{*}_{j}(\vec{q})f_{i}(\vec{q})$
$\displaystyle=$
$\displaystyle\delta({\vec{P}^{\prime}}-{\vec{P}^{\prime\prime}})\int\frac{d^{3}q}{4(q^{2}+m^{2})}f^{*}_{j}(\vec{q})f_{i}(\vec{q}),$

where $m$ is the mass of the individual particles. In the last expression we
have assumed that the CM momentum is small compared to the mass of the heavy
quarks. Provided that the integral equals $\delta_{ij}E_{P}$ we obtain the
proper normalization. At small momentum $P$, one has $E_{P}\approx M_{i}$ –
the mass of the two-particle state –, so that it is always possible to impose
such constrain on the integral…but thus this ceases to be so for finite values
of $P$. To further discuss our approach we temporary neglect the transverse
degrees of freedom. This reduces our approach to one momentum space and one
energy dimension.

###  A.1 Variable Transformation

To advance towards fully relativistic definitions let us consider the
invariant measure

$\displaystyle I$ $\displaystyle=$
$\displaystyle\int\frac{dp_{1}}{2e_{1}}\frac{dp_{1}}{2e_{1}}=\int
dp_{1}de_{1}dp_{2}de_{2}\delta(p_{1}^{2}-m^{2})$ (74) $\displaystyle\times$
$\displaystyle\delta(p_{2}^{2}-m^{2})\theta(e_{1})\theta(e_{2}).$

We introduce the relative momentum
$\Delta=\frac{{p}_{1}-{p}_{2}}{2}=(q^{0},q)$ and the center of mass momentum
$\Sigma=p_{1}+p_{2}=(P^{0},P)$. With $d^{2}p_{1}d^{2}p_{2}=d^{2}\Sigma
d^{2}\Delta$, we obtain

$I=\int\frac{d^{2}\Sigma
d^{2}\Delta}{2}\delta(\frac{\Sigma^{2}}{4}+\Delta^{2}-m^{2})\delta(\Sigma\cdot\Delta)\theta(\frac{P^{0}}{2}-|q^{0}|).$
(75)

$\delta(\Sigma\cdot\Delta)$ implies

$|q^{0}|=\sqrt{\frac{(P^{0})^{2}-4(m^{2}+q^{2})}{(P^{0})^{2}-4q^{2}}}|q|\leq|q|$
(76)

for $|q|<P^{0}/2$. The $\delta$ functions takes therefore care that the
condition $|q^{0}|<P^{0}/2$ is always fulfilled. $\Sigma$ is a time-like
vector, therefore $\delta(\Sigma\cdot\Delta)$ implies that $\Delta$ has to be
a space-like vector. Introducing $Q^{2}=-\Delta^{2}$, we now perform a second
variable transformation

$\displaystyle P^{0}$ $\displaystyle=$ $\displaystyle\sqrt{s}\ \cosh Y\ ,\
P=\sqrt{s}\ \sinh Y$ $\displaystyle q^{0}$ $\displaystyle=$ $\displaystyle Q\
\sinh y\ ,\ q=\pm Q\ \cosh y$ (77)

with $d^{2}\Sigma d^{2}\Delta=\frac{ds}{2}QdQdYdy$, including both, the
$q=+Q\cosh y$ and the $q=-Q\cosh y$ sector. In these variables we find

$\delta(\Sigma\cdot\Delta)=\delta(\sqrt{s}Q\sinh(y\pm Y))=\frac{\delta(y\pm
Y)}{\sqrt{s}Q}$

where the $\pm$ stand for the $q=\pm Q\cosh y$ sectors, respectively. The 2nd
distribution $\delta(\frac{\Sigma^{2}}{4}+\Delta^{2}-m^{2})$ simply writes

$\delta\left(\frac{\Sigma^{2}}{4}+\Delta^{2}-m^{2}\right)=\delta\left(\frac{s}{4}-Q^{2}-m^{2}\right),$
(78)

implying that the condition $|q|<P^{0}/2$ is always satisfied. The integral I
(eq. 74) can thus be rewritten as

$\displaystyle I=I_{+}+I_{-}$ $\displaystyle=$ $\displaystyle\int
dY\frac{d|q|}{\sqrt{s}}(\Theta(q>0)+\Theta(q<0))$ (79) $\displaystyle=$
$\displaystyle\int dY\int_{-\infty}^{+\infty}\frac{dq}{\sqrt{s}},$

where the 2 sectors have been merged.

###  A.2 Relativistic Two-Particles States

These new variables allow for reformulating our relativistic state. For any
lab frame $S^{\prime}$, we assume that there exists a frame $S$ where the two-
particle state CM is nearly at rest and define $y_{\Phi}$ as the rapidity of
$S$ in $S^{\prime}$. In $S$, the state is defined as

$\displaystyle|\Phi_{i,P\approx 0}\rangle$ $\displaystyle=$
$\displaystyle\int\frac{dp_{1}dp_{2}}{2e_{1}2e_{2}}\delta({P}-p_{1}-p_{2})f_{i}(p_{1}-p_{2})|p_{1}\rangle|p_{2}\rangle$
(80) $\displaystyle=$
$\displaystyle\int\frac{dq}{\sqrt{s}}dY\delta({P}-p_{1}-p_{2})f_{i}(q)|p_{1}\rangle|p_{2}\rangle.$

Let $f_{i}$ be defined as a boost invariant wave function depending on the
relative momentum $q$ measured in the rest frame. Proceeding to a Lorentz
transform to the lab frame $S^{\prime}$, the two-particle CM state reads

$\displaystyle|\Phi_{i,y_{\Phi}}\rangle$ $\displaystyle=$
$\displaystyle\int\frac{dq}{\sqrt{s}}dY\delta(\sqrt{s}\sinh(y_{\Phi}-Y))f_{i}(q)|p_{1}^{\prime}\rangle|p_{2}^{\prime}\rangle$
(81) $\displaystyle=$
$\displaystyle\int\frac{dq}{s}dY\delta(y_{\Phi}-Y)f_{i}(q)|p_{1}^{\prime}\rangle|p_{2}^{\prime}\rangle$

where the momenta $p_{1}^{\prime}$ and $p_{2}^{\prime}$ are taken as
$(\sqrt{m^{2}+q^{2}},q)=m(\cosh(\hat{y}),\sinh(\hat{y}))$ in $S$ and then
boosted with a rapidity shift $+y_{\Phi}$, leading to

$\displaystyle p_{1}^{\prime}$ $\displaystyle=$ $\displaystyle
m(\cosh(y_{\Phi}+\hat{y}),\sinh(y_{\Phi}+\hat{y}))$ $\displaystyle
p_{2}^{\prime}$ $\displaystyle=$ $\displaystyle
m(\cosh(y_{\Phi}-\hat{y}),\sinh(y_{\Phi}-\hat{y})).$ (82)

It remains to be checked that the orthogonality conditions of the state (A) is
fulfilled. After some trivial calculation, one finds

$\displaystyle\langle\Phi_{i,y^{{}^{\prime}}_{\Phi}}|\Phi_{j,y_{\Phi}}\rangle$
$\displaystyle=$
$\displaystyle\int\int\frac{dp_{1}}{2E_{1}}\frac{dp_{2}}{2E_{2}}\times\frac{\delta(y^{{}^{\prime}}_{\Phi}-Y)}{\sqrt{s}}\frac{\delta(y_{\Phi}-Y)}{\sqrt{s}}$
(83) $\displaystyle\times$ $\displaystyle f_{i}^{\star}(q)f_{j}(q)$
$\displaystyle=$
$\displaystyle\int\frac{dqdY}{\sqrt{s}}\frac{\delta(y^{{}^{\prime}}_{\Phi}-Y)\delta(y_{\Phi}-Y)}{s}$
$\displaystyle\times$ $\displaystyle f^{\star}_{i}(q)f_{j}(q)$
$\displaystyle=$
$\displaystyle\delta(y^{{}^{\prime}}_{\Phi}-y_{\Phi})\int\frac{dq}{s^{3/2}}f^{\star}_{i}(q)f_{j}(q)$

and we can derive the invariant orthogonality relations provided we require

$\int\frac{dq}{s^{3/2}}f^{\star}_{i}(q)f_{j}(q)=\delta_{ij}.$

By introducing the non relativistic wave function

$\psi_{i}(q)=\frac{f_{i}(q)}{s^{\frac{3}{4}}(q)}$ (84)

we obtain the orthogonality relation for non relativistic wave functions $\int
dq\psi_{i}^{*}(q)\psi_{j}(q)=\delta_{ij}$.

Although we have been able to develop a prescription for the construction of
an orthogonal invariant base, there are still some comments to be made: first,
we do not strictly recover the anticipated normalization
$\langle\Phi_{y^{{}^{\prime}}_{\Phi},i}|\Phi_{y_{\Phi},i}\rangle=2E\delta(p-p^{{}^{\prime}})\delta_{ij}$.
This can be explained due to the fact that for a given state $i$ and a given
rapidity $y_{\Phi}$, the total energy depends on the relative momentum $q$.
This has consequences for the orthogonality relation A.2, which slightly
differs from 73. However, these states admit proper transformation laws under
Lorentz boosts and the relationship 83 allows to create a completely
orthogonal set of two-particle particle states

$|\Phi_{i}\rangle:=\int dy_{\Phi}g_{i}(y_{\Phi})|\Phi_{y_{\Phi},i}\rangle,$
(85)

where $g$ is an arbitrary function which satisfies the relation

$\langle\Phi_{i}|\Phi{{}^{\prime}}_{j}\rangle=\delta_{ij}\int
dy_{\Phi}g^{*}_{i}(y_{\Phi})g_{j}^{{}^{\prime}}(y_{\Phi}),$ (86)

with an easy connection to the rapidity spectrum:

$\frac{dN_{i}}{dy_{\Phi}}=|\langle\Phi_{y_{\Phi},i}|\Phi\rangle|^{2}=|g_{i}|^{2}.$
(87)

After having developed a method to build an orthogonal boost invariant basis
on which we can project our two-particle states, we turn to the building of
the Wigner function on this basis.

##  Appendix B Relativistic Wigner Function from a Given Basis in (1+1)D Case

The main difficulty for finding a relativistic Wigner function from an
orthogonal basis of the form (85), lies in the difficulty of defining a
conjugate variable of the rapidity. Only the option to take the relative
momentum $q$ in the center of mass as the conjugate variable allowed to arrive
at a satisfactory conclusion and at the same time to obtain results which have
a clear physical significance. Here we discuss this option $[Y,q]$.

To obtain a Wigner density in $[Y,q]$ for the basis of the form (85) it is
convenient to rewrite (87) in the form

$\frac{dN_{i}}{dy_{\Phi}}=Tr(\widehat{\rho}_{\Phi}\widehat{\rho}_{i,y_{\Phi}})$
(88)

where the trace is performed over the phase space variables ($[Y,q]$),
$y_{\Phi}$ being the rapidity of the two-particle system (quarkonium). The
density operators $\widehat{\rho}_{\Phi}$ and $\widehat{\rho}_{i,y_{\Phi}}$
have the form

$\begin{split}\widehat{\rho}_{\Phi}=|\Phi\rangle\times\langle\Phi|\quad;\quad\widehat{\rho}_{i,y_{\Phi}}=|\Phi_{i,y_{\Phi}}\rangle\times\langle\Phi_{i,y_{\Phi}}|\end{split}$
(89)

in which $|\Phi\rangle$ and $|\Phi_{i,y_{\Phi}}\rangle$ represents a generic
$(Y,q)$ two-particle wave function and the two-particle wave function for a
state i with the rapidity $y_{\Phi}$ respectively. Inserting the identity
operator

$\widehat{I}=\int dY\frac{d{q}}{\sqrt{s}}|1,2\rangle\langle 1,2|$ (90)

in the expression for the spectrum (88) we obtain

$\frac{dN_{i}}{dy_{\Phi}}=\int
dYdY^{{}^{\prime}}\frac{d{q}d{q^{\prime}}}{\sqrt{ss^{{}^{\prime}}}}\langle
Y,q|\widehat{\rho}|Y^{{}^{\prime}},{q^{\prime}}\rangle\langle
Y^{{}^{\prime}},{q^{\prime}}|\widehat{\rho}_{i,y_{\Phi}}|Y,q\rangle$ (91)

where $|Y,p_{r}\rangle$ refers to the $|1,2\rangle$ state with a total
rapidity Y and a relative momentum $p_{r}$ and $\sqrt{s}$ respectively
$\sqrt{s^{{}^{\prime}}}$ are the center of mass energy of the states.

Defining

$\overline{\rho}({Y},q;Y^{{}^{\prime}},{q^{\prime}})=\frac{\langle{Y}^{{}^{\prime}},{q^{\prime}}|\widehat{\rho}|{Y},q\rangle}{(ss^{{}^{\prime}})^{1/4}},$
(92)

we get

$\frac{dN_{i}}{dy_{\Phi}}=\int
dYdY^{{}^{\prime}}d{q}d{q^{\prime}}\overline{\rho}({Y^{\prime}},{q^{\prime}};{Y},q)\bar{\rho}_{i,y_{\Phi}}({Y},q;{Y}^{{}^{\prime}},{q^{\prime}})$
(93)

We introduce now, in preparation of the Wigner transformation, the auxiliary
variables $\overline{{Y}}=\frac{{Y}+{Y}^{{}^{\prime}}}{2}$,
$\overline{q}=\frac{q+{q^{\prime}}}{2}$, $\Delta{Y}={Y}-{Y}^{{}^{\prime}}$ and
$\Delta q=q-{q^{\prime}}$, which transform (93) to

$\displaystyle\frac{dN_{i}}{dy_{\Phi}}$ $\displaystyle=$ $\displaystyle\int
d\overline{Y}d\overline{q}d\Delta Yd\Delta
Y^{{}^{\prime}}d\Delta{q}d\Delta{q^{\prime}}$ (94) $\displaystyle\times$
$\displaystyle\overline{\rho}(\overline{{Y}}-\frac{\Delta{Y}}{2},\overline{{p}}_{r}-\frac{\Delta
q}{2};\overline{{Y}}+\frac{\Delta{Y}}{2},\overline{q}+\frac{\Delta q}{2})$
$\displaystyle\times$
$\displaystyle\overline{\rho}_{i,y_{\Phi}}(\overline{{Y}}+\frac{\Delta{Y}^{\prime}}{2},\overline{q}+\frac{\Delta
q^{\prime}}{2};\overline{Y}-\frac{\Delta{Y}^{\prime}}{2},\overline{q}-\frac{\Delta
q^{\prime}}{2})$ $\displaystyle\times$
$\displaystyle\delta(\Delta{Y}-\Delta{Y}^{{}^{\prime}})\delta(\Delta
q-\Delta{q^{\prime}}).$

The relationship between conjugate variables is given by

$\displaystyle\delta(\Delta q-\Delta{q^{\prime}})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\hbar}\int
dx_{r}e^{i\frac{{x}_{r}(\Delta{q}-\Delta{q^{\prime}})}{\hbar}}$
$\displaystyle\delta(\Delta{Y}-\Delta{Y}^{{}^{\prime}})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int
d{k_{3}}e^{i{k_{3}}(\Delta{Y}-\Delta{Y}^{{}^{\prime}})}$ (95)

where ${x}_{r}$ is the relative position measured in the $Q\bar{Q}$ center of
mass system and $k_{3}$ corresponds to the dimensionless eigenvalues of the
boost operator Durand and O’Raifeartaigh (1976). We recall as well the
definition of the Wigner function of a density operator
$\rho({r},{r}^{{}^{\prime}})$, for instance associated to some wave function
$\psi$ through
$\rho({r},{r}^{{}^{\prime}})=\psi({r})\psi^{\star}(r^{\prime})$:

$W(r,p)=\frac{1}{2\pi\hbar}\int
dye^{-i\frac{{py}}{\hbar}}\rho({r}+\frac{{y}}{2},{r}-\frac{{y}}{2}).$ (96)

Substituting the delta distributions relations (95) and comparing the
definition of the Wigner function with the factors in the equation (94) we
obtain the Wigner function for the two-particle density operator

$\displaystyle W_{i,y_{\Phi}}(\overline{{Y}},{k_{3}};q,{x}_{r})$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{2}\hbar}\int d\Delta
Y^{{}^{\prime}}d\Delta{q^{\prime}}$ (97)
$\displaystyle\times\overline{\rho}_{i,y_{\Phi}}(\overline{{Y}}+\frac{\Delta{Y}^{{}^{\prime}}}{2},\overline{{q}}+\frac{\Delta{q^{\prime}}}{2};\overline{{Y}}-\frac{\Delta{Y}^{{}^{\prime}}}{2},\overline{q}-\frac{\Delta{q^{\prime}}}{2})$
$\displaystyle\times e^{-i({k_{3}}\Delta{Y}^{{}^{\prime}}+{x}_{r}\frac{\Delta
q^{\prime}}{\hbar})}$

and

$\displaystyle W(\overline{{Y}},{k_{3}};{q},{x}_{r})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{2}\hbar}\int d\Delta Yd\Delta{q}$ (98)
$\displaystyle\times\overline{\rho}(\overline{Y}-\frac{\Delta
Y}{2},\overline{q}-\frac{\Delta q}{2};\overline{{Y}}+\frac{\Delta
Y}{2},\overline{{q}}+\frac{\Delta{q}}{2})$ $\displaystyle\times
e^{i({k_{3}}\Delta{Y}+{x}_{r}\frac{\Delta{q}}{\hbar})}.$

Substituting the equations for the Wigner functions, 98 and 97 in equation 94
we obtain

$\displaystyle\frac{dN_{i}}{dy_{\Phi}}$ $\displaystyle=$
$\displaystyle(2\pi)^{2}\hbar\int
d\overline{Y}d\overline{q}dk_{3}dx_{r}W(\overline{{Y}},{k_{3}};\overline{q},{x}_{r})$
(99) $\displaystyle\times$ $\displaystyle
W_{i,y_{\Phi}}(\overline{{Y}},{k_{3}};\overline{q},{x}_{r})$

so we can, as in the non relativistic case, consider
$\frac{dN_{i}}{dy_{\Phi}}$ as a convolution of the two Wigner densities. Next
we need to evaluate
$W_{i,y_{\Phi}}(\overline{{Y}},{k_{3}};\overline{q},{x}_{r})$. We start from
eq. 92 and obtain:

$\displaystyle\overline{\rho}_{i,y_{\Phi}}({Y},{q};{Y}^{{}^{\prime}},{q}^{{}^{\prime}})$
$\displaystyle=$
$\displaystyle\frac{f^{\star}_{i,y_{\Phi}}(q)f_{i,y_{\Phi}}(q^{\prime})}{(ss^{{}^{\prime}})^{\frac{1}{4}}}$
$\displaystyle=\frac{\delta(\sqrt{s}(Y-y_{\Phi}))\delta(\sqrt{s^{{}^{\prime}}}(Y^{{}^{\prime}}-y_{\Phi}))f^{{}^{\prime}\star}_{i}({q})f_{i}({q}^{{}^{\prime}})}{(ss^{{}^{\prime}})^{\frac{1}{4}}}$
$\displaystyle=\frac{\delta(Y-y_{\Phi})\delta(Y^{{}^{\prime}}-y_{\Phi})f^{{}^{\prime}\star}_{i}({q})f_{i}({q}^{{}^{\prime}})}{(ss^{{}^{\prime}})^{\frac{3}{4}}}.$

with $f_{i}(q)$, the ”wave function” in eq. 81. With this density matrix we
calculate now the Wigner density, eq. 97

$\displaystyle W_{i,y_{\Phi}}(\overline{{Y}},{k_{3}};\overline{q},{x}_{r})$
$\displaystyle=$ $\displaystyle\frac{1}{(2\pi)^{2}\hbar}\int\frac{d\Delta
Yd\Delta{q}}{(ss^{{}^{\prime}})^{3/4}}e^{-i{x}_{r}\frac{\Delta{q}}{\hbar}}\times$
(101) $\displaystyle\delta(\overline{{Y}}-{y}_{\Phi})\delta(\Delta
Y)f^{*}_{i}(\overline{{p}}_{r}+\frac{\Delta{q}}{2})f_{i}(\overline{{p}}_{r}-\frac{\Delta{q}}{2})$
$\displaystyle=\frac{\delta(\overline{{Y}}-{y}_{\Phi})}{(2\pi)^{2}\hbar}\int
d\Delta{q}\,\psi^{\star}_{i}(\overline{{p}}_{r}+\frac{\Delta{q}}{2})\psi_{i}(\overline{q}-\frac{\Delta{q}}{2})$
$\displaystyle\hskip 14.22636pt\times e^{-i{x}_{r}\frac{\Delta{q}}{\hbar}}$

where we have employed
$\delta(\overline{{Y}}+\frac{\Delta{Y}}{2}-{y}_{\Phi})\delta(\overline{{Y}}-\frac{\Delta{Y}}{2}-{y}_{\phi})=\delta(\overline{{Y}}-{y}_{\phi})\delta(\Delta{Y})$
and eq. 84. We realize that $\frac{1}{2\pi\hbar}\int
d\Delta{q}\psi^{\star}_{i}(\overline{{p}}_{r}+\frac{\Delta{q}}{2})\psi_{i}(\overline{{p}}_{r}-\frac{\Delta{q}}{2})e^{-i{x}_{r}\frac{\Delta{q}}{\hbar}}$
is nothing than the non relativistic Wigner function $W_{i,{\rm NR}}$
expressed in the coordinates of the center of mass. So we can finally
establish

$\displaystyle W_{i,y_{\Phi}}(\overline{{Y}},{k_{3}};\overline{q},{x}_{r})$
$\displaystyle=$
$\displaystyle\frac{\delta(\overline{Y}-{y}_{\Phi})}{2\pi}W_{i,{\rm
NR}}(\overline{q},{x}_{r})$ (102)

The equation for $W_{i,y_{\Phi}}$ does not depend on the boost operator
eigenvalue $k_{3}$, whose meaning in physical terms is that the states used in
the definition of $\bar{\rho}_{i,y_{\Phi}}$ are plane waves wrt the center of
mass motion and come with a single rapidity, leading to a $\delta(\Delta Y)$
in the Wigner calculation. $W_{i,y_{\Phi}}$ depends on the rapidity of the
center of mass of the quarkonia state in the computational frame. Equation
(102) describes how we can evaluate the Wigner function for a $Q\bar{Q}$ pair
in the bound state $i$ with $y_{\Phi}$. It is also important to realize that
the boosted Wigner function (102) inherits some properties of the non
relativistic Wigner function $W_{\rm NR}$, including the normalization.
Substituting equation (102) in equation (99), we obtain the rapidity
distribution of the $Q\bar{Q}$ pairs which are bound in a state $i$:

$\displaystyle\frac{dN_{i}}{dy_{\Phi}}$ $\displaystyle=$
$\displaystyle(2\pi)^{2}\hbar\int
d\overline{Y}d\overline{q}dk_{3}dx_{r}W(\overline{{Y}},k_{3};\overline{{q}},x_{r})W_{i,y_{\Phi}}(\overline{{Y}},\overline{{q}},{x}_{r})$
(103) $\displaystyle=$ $\displaystyle 2\pi\hbar\int
d\overline{Y}\delta(\overline{{Y}}-y_{\Phi})\int d\overline{q}dx_{r}W_{i,{\rm
NR}}(\overline{{q}},{x}_{r})$ $\displaystyle\times\underbrace{\int
dk_{3}W(\overline{{Y}},{k_{3}};\overline{{q}},{x}_{r})}_{\bar{W}(\overline{{Y}},\overline{{q}},{x}_{r})}$
$\displaystyle=$ $\displaystyle 2\pi\hbar\int d\overline{q}dx_{r}W_{i,{\rm
NR}}(\overline{{q}},{x}_{r})\overline{W}(y_{\Phi},\overline{{q}},{x}_{r})$

In equation (103) the Wigner function $W_{i,\rm NR}(\overline{{q}},{x}_{r})$
represents the probability density (non-relativistic) of formation of a
quarkonium state $i$. $\int
dk_{3}W(\overline{{Y}},{k_{3}};\overline{{q}},{x}_{r})$ represents the
probability density of finding a $Q\bar{Q}$ pair with relative momentum and
position $q$ and $x_{r}$ and with a center of mass rapidity $y_{\Phi}$. Since
$W$ references to free $Q$ and $\bar{Q}$ before they form a bound state, the
relationship with the two-particle probability density operator
$\rho_{Q\bar{Q}}$, can be traced back by using the Wigner function definition

$\displaystyle\overline{W}(\overline{Y},{q},{x}_{r})$ $\displaystyle=$
$\displaystyle\int dk_{3}W(\overline{Y},{k_{3}};\overline{{q}},{x}_{r})$
$\displaystyle=\frac{1}{(2\pi)^{2}\hbar}\int dk_{3}\int d\Delta
Yd\Delta{q}e^{-i{k_{3}}\Delta{Y}-i{x}_{r}\frac{\Delta{q}}{\hbar}}$
$\displaystyle\times\overline{\rho}(\overline{{Y}}+\frac{\Delta{Y}}{2},\overline{q}+\frac{\Delta{q}}{2};\overline{{Y}}-\frac{\Delta{Y}}{2},\overline{{q}}-\frac{\Delta{q}}{2})$
$\displaystyle=\frac{1}{2\pi\hbar}\int
d\Delta{q}\overline{\rho}(\overline{{Y}},\overline{{q}}+\frac{\Delta{q}}{2};\overline{{Y}},\overline{{q}}-\frac{\Delta{q}}{2})e^{-i{x}_{r}\frac{\Delta{q}}{\hbar}}.$

The integral of $\overline{W}$ over $x_{r}$ has the form :

$\int
dx_{r}\overline{W}(\overline{{Y}},\overline{q},{x}_{r})=\overline{\rho}(\overline{{Y}},\overline{q};\overline{{Y}},\overline{q})=\frac{\langle\overline{{Y}},\overline{q}|\overline{\rho}_{Q\bar{Q}}|\overline{{Y}},\overline{q}\rangle}{\sqrt{s}}$
(105)

which satisfies the following relationship:

$\displaystyle\int$ $\displaystyle d{q}d\overline{{Y}}\int
dx_{r}\overline{W}(\overline{{Y}},{q},{x}_{r})=\int\frac{d\overline{Y}d{q}}{\sqrt{s}}\rho_{Q\bar{Q}}(\overline{{Y}},{q})$
(106) $\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}d^{3}p_{2}}{2e_{1}e_{2}}\rho_{Q\bar{Q}}({p}_{1},{p}_{2})=N_{Q\bar{Q}}$

where $N_{Q\bar{Q}}=N_{Q}\times N_{\bar{Q}}$, is the total number of different
$Q\bar{Q}$ pairs that are present in the system. One can thus write the
following relationship for $\overline{W}(\overline{{Y}},\overline{q},{x}_{r})$

$\overline{W}(\overline{{Y}},\overline{q},{x}_{r})=\frac{d^{3}N_{Q\bar{Q}}}{d\overline{{Y}}d\overline{q}dx_{r}}$
(107)

which is simply the differential form of equation (106). The rapidity
distribution, eq. 103 can now be expressed as

$\frac{dN_{i}}{dy_{\Phi}}=2\pi\hbar\int
d\overline{q}dx_{r}\frac{d^{3}N_{Q\bar{Q}}}{dy_{\Phi}d\overline{q}dx_{r}}W_{i,\rm
NR}(\overline{q},{x}_{r})$ (108)

##  Appendix C Generalizing for the (3+1)D Case

The results derived in equation (108) can be generalized to the case (3 +1)D
if we start from a state generated according to equation (69)

$|\Phi_{i,\mathbf{u}_{\Phi}}\rangle=\int\frac{d^{3}u}{u^{0}}\frac{d^{3}q^{\rm
cm}}{\sqrt{s}}f_{i,\mathbf{u}_{\Phi}}(\mathbf{u},\mathbf{q}^{\rm
cm})|\mathbf{p}_{1},\mathbf{p}_{2}\rangle$ (109)

where $(u^{0},\mathbf{u})$ is the 4-velocity of the quarkonium state $\Phi$
and $\mathbf{q}^{\rm cm}$ is the relative momentum evaluated in the quarkonium
state CM (we use an explicit ”cm” superscript anticipating a similar
formulation in the lab-frame). It is also possible to express these quarkonium
states in the basis $(\mathbf{u_{T}},Y)$, if we separate $\mathbf{u}$ into its
transverse $\mathbf{u}_{T}$ and longitudinal $u_{z}$ components, where $Y={\rm
atanh}(\frac{u_{z}}{u^{0}})$, leading to the analogous of equation (81):

$\displaystyle|\Phi_{i,y_{\Phi},{\bf u}_{T,\phi}}\rangle$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}q^{\rm
cm}}{s}dYd^{2}{u}_{T}\delta(y_{\Phi}-Y)\delta({\bf u}_{T\phi}-{\bf u}_{T})$
(110) $\displaystyle f_{i}(q^{\rm cm})|p_{1}\rangle|p_{2}\rangle$

where ${\bf u}_{T,\phi}$ is the transverse component of the quarkonium
4-velocity. One can then generalize equation (108) by extending this formula
to the transverse component (3D case), we arrive at

$\displaystyle\frac{d^{3}N_{i}}{dy_{\Phi}d^{2}u_{T,\Phi}}=(2\pi)^{6}(\hbar)^{3}\int
dYd^{2}u_{T}d^{3}q^{\rm cm}d^{3}x^{\rm cm}_{r}$ $\displaystyle\times$
$\displaystyle\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm cm}d^{3}x^{\rm
cm}_{r}}W_{i,y_{\Phi},\mathbf{u}_{T,\Phi}}(Y,\mathbf{u}_{T},\mathbf{q}^{\rm
cm},\mathbf{x}^{\rm cm}_{r})$

where

$\displaystyle
W_{i,y_{\Phi},\mathbf{u}_{T,\Phi}}(Y,\mathbf{u}_{T},\mathbf{q}^{\rm
cm},\mathbf{x}^{\rm cm}_{r})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\delta(y-y_{\Phi})\delta^{(2)}(\mathbf{u}_{T,\Phi}-\mathbf{u}_{T})W_{\rm
NR}(\mathbf{q}^{\rm cm},\mathbf{x}^{\rm cm}_{r}).$

Equation (C) is the 3D generalization of equation (102). Inspection shows that
we have just to multiply by the factor
$\frac{\delta^{(2)}(\mathbf{u}_{T,\Phi}-\mathbf{u}_{T})}{(2\pi)^{2}}$. After
performing the integrals over the delta distributions we obtain:

$\displaystyle\frac{d^{3}N_{i}}{dy_{\Phi}d^{2}u_{T,\Phi}}$ $\displaystyle=$
$\displaystyle(2\pi\hbar)^{3}\int d^{3}q^{\rm cm}d^{3}x^{\rm
cm}_{r}\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm cm}d^{3}x^{\rm
cm}_{r}}$ (113) $\displaystyle\times$ $\displaystyle W_{i,\rm
NR}(\mathbf{q}^{\rm cm},\mathbf{x}^{\rm cm}_{r})$

in which the nine-fold distribution in the integral is taken at $Y=y_{\Phi}$
and ${\bf u}_{T}={\bf u}_{T,\Phi}$. This expression can be also expressed in
the coordinates of the computational frame

$\displaystyle\frac{d^{3}N_{i}}{dy_{\Phi}d^{2}u_{T,\Phi}}$ $\displaystyle=$
$\displaystyle h^{3}\int d^{3}q^{\rm lab}d^{3}x^{\rm
lab}_{r}\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm lab}d^{3}x^{\rm
lab}_{r}}$ (114) $\displaystyle\times$ $\displaystyle W_{i,\rm
NR}(\mathbf{q}^{\rm cm}(\mathbf{q}^{\rm lab}),\mathbf{x}^{\rm
cm}_{r}(\mathbf{x}^{\rm lab}_{r}))$

in which the coordinates $\mathbf{q}^{\rm cm}$ and $\mathbf{x}^{\rm cm}_{r}$
have to be expressed as a function of the $\mathbf{q}^{\rm lab}$ and
$\mathbf{x}^{\rm lab}_{r}$. From the previous equation, integrating by the
variables $u_{T,\Phi}$ and $y_{\Phi}$ we can obtain the absolute number of
states of quarkonium $i$-th

$\displaystyle N_{i}$ $\displaystyle=$ $\displaystyle h^{3}\int
dy_{\Phi}d^{2}u_{T}d^{3}q^{\rm lab}d^{3}x^{\rm
lab}_{r}\frac{d^{9}N_{Q\bar{Q}}}{dYd^{2}u_{T}d^{3}q^{\rm lab}d^{3}x^{\rm
lab}_{r}}$ (115) $\displaystyle\times$ $\displaystyle W_{i,\rm
NR}(\mathbf{q}^{\rm cm}(\mathbf{q}^{\rm lab}),\mathbf{x}^{\rm
cm}_{r}(\mathbf{x}^{\rm lab}_{r}))$

The nine-fold distribution has to be considered according to the physical
situation…In nucleus-nucleus collisions, we take all possible $(Q,\bar{Q})$
combinations into account in order to form $i$-th quarkonium state. In Monte
Carlo Implementation, equation 114 becomes

$\displaystyle\frac{dN_{i}}{dy_{\Phi}d^{2}u_{T,\Phi}}$ $\displaystyle=$
$\displaystyle h^{3}\sum^{N_{Q}\times
N_{\bar{Q}}}_{l=1}W_{i,y_{\Phi},\mathbf{u}_{T,\Phi}}(Y,\mathbf{u}_{T},\mathbf{q},\mathbf{x}_{r})$
$\displaystyle=h^{3}\sum^{N_{Q}\times
N_{\bar{Q}}}_{l=1}\delta(Y-y_{\Phi})\delta(\mathbf{u}_{T,\Phi}-\mathbf{u}_{T})W_{i,\rm
NR}(\mathbf{q}^{\rm cm},\mathbf{x}^{\rm cm}_{r})$

where the sum runs over all $N_{Q}\times N_{\bar{Q}}$ possible combinations
and where the $\\{\mathbf{q}^{\rm cm},\mathbf{x}^{\rm cm}_{r},\ldots\\}$ are
constructed for each pair. So the total number of states $\Phi_{i}$ is given
by the expression

$N_{i}=h^{3}\sum^{N_{Q}\times N_{\bar{Q}}}_{l=1}W_{i,\rm NR}(\mathbf{q}^{\rm
cm},\mathbf{x}^{\rm cm}_{r}).$ (117)

The Wigner distribution function $W_{i,\rm NR}$ has to be multiplied by the
factor $g=\frac{2s_{\Psi}+1}{8(2s_{Q}+1)^{2}}$ to take into account the spin
of the quarkonium state and of the quarks. The factor 1/8 is due to the fact
that in our approach only color neutral $Q\bar{Q}$ combinations can form a
quarkonium.

## References

  * Andronic et al. (2018) A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, Nature 561, 321 (2018), eprint 1710.09425.
  * Adam et al. (2017) J. Adam et al. (ALICE), Phys. Lett. B 766, 212 (2017), eprint 1606.08197.
  * Adare et al. (2011) A. Adare et al. (PHENIX), Phys. Rev. C 84, 054912 (2011), eprint 1103.6269.
  * Acharya et al. (2020a) S. Acharya et al. (ALICE), JHEP 10, 141 (2020a), eprint 2005.14518.
  * Matsui and Satz (1986) T. Matsui and H. Satz, Phys. Lett. B 178, 416 (1986).
  * Mocsy and Petreczky (2007) A. Mocsy and P. Petreczky, Phys. Rev. Lett. 99, 211602 (2007), eprint 0706.2183.
  * Digal et al. (2001) S. Digal, P. Petreczky, and H. Satz, Phys. Rev. D 64, 094015 (2001), eprint hep-ph/0106017.
  * Bala et al. (2021) D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf, and J. H. Weber (2021), eprint 2110.11659.
  * Rothkopf (2020) A. Rothkopf, Phys. Rept. 858, 1 (2020), eprint 1912.02253.
  * Xu et al. (2019) Y. Xu et al., Phys. Rev. C 99, 014902 (2019), eprint 1809.10734.
  * Cao et al. (2019) S. Cao et al., Phys. Rev. C 99, 054907 (2019), eprint 1809.07894.
  * Beraudo et al. (2018) A. Beraudo et al., Nucl. Phys. A 979, 21 (2018), eprint 1803.03824.
  * Cacciari et al. (2008) M. Cacciari, G. P. Salam, and G. Soyez, JHEP 04, 063 (2008), eprint 0802.1189.
  * Cacciari et al. (2012) M. Cacciari, G. P. Salam, and G. Soyez, Eur. Phys. J. C72, 1896 (2012), eprint 1111.6097.
  * Song et al. (2015) T. Song, H. Berrehrah, D. Cabrera, J. M. Torres-Rincon, L. Tolos, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 92, 014910 (2015), eprint 1503.03039.
  * Linnyk et al. (2008) O. Linnyk, E. L. Bratkovskaya, and W. Cassing, Int. J. Mod. Phys. E 17, 1367 (2008), eprint 0808.1504.
  * Andronic et al. (2021) A. Andronic, P. Braun-Munzinger, M. K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel, and V. Vislavicius, JHEP 07, 035 (2021), eprint 2104.12754.
  * Du and Rapp (2015) X. Du and R. Rapp, Nucl. Phys. A 943, 147 (2015), eprint 1504.00670.
  * Zhou et al. (2014) K. Zhou, N. Xu, Z. Xu, and P. Zhuang, Phys. Rev. C 89, 054911 (2014), eprint 1401.5845.
  * He et al. (2022) M. He, B. Wu, and R. Rapp, Phys. Rev. Lett. 128, 162301 (2022), eprint 2111.13528.
  * Akamatsu and Miura (2022) Y. Akamatsu and T. Miura, EPJ Web Conf. 258, 01006 (2022), eprint 2111.15402.
  * Delorme et al. (2022) S. Delorme, T. Gousset, R. Katz, and P.-B. Gossiaux, EPJ Web Conf. 258, 05009 (2022).
  * Blaizot and Escobedo (2018) J.-P. Blaizot and M. A. Escobedo, Phys. Rev. D 98, 074007 (2018), eprint 1803.07996.
  * Remler (1981) E. A. Remler, Annals Phys. 136, 293 (1981).
  * Gyulassy et al. (1983) M. Gyulassy, K. Frankel, and E. a. Remler, Nucl. Phys. A 402, 596 (1983).
  * Aichelin and Remler (1987) J. Aichelin and E. A. Remler, Phys. Rev. C 35, 1291 (1987).
  * Song et al. (2017) T. Song, J. Aichelin, and E. Bratkovskaya, Phys. Rev. C 96, 014907 (2017), eprint 1705.00046.
  * Drescher et al. (2001) H. J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog, and K. Werner, Phys. Rept. 350, 93 (2001), eprint hep-ph/0007198.
  * Werner et al. (2010) K. Werner, I. Karpenko, T. Pierog, M. Bleicher, and K. Mikhailov, Phys. Rev. C82, 044904 (2010), eprint 1004.0805.
  * Gossiaux and Aichelin (2008) P. B. Gossiaux and J. Aichelin, Phys. Rev. C 78, 014904 (2008), eprint 0802.2525.
  * Gossiaux et al. (2009) P. B. Gossiaux, R. Bierkandt, and J. Aichelin, Phys. Rev. C79, 044906 (2009), eprint 0901.0946.
  * Nahrgang et al. (2016) M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, Phys. Rev. C 93, 044909 (2016), eprint 1602.03544.
  * Gossiaux and Aichelin (2009) P. B. Gossiaux and J. Aichelin, J. Phys. G 36, 064028 (2009), eprint 0901.2462.
  * Nahrgang et al. (2015) M. Nahrgang, J. Aichelin, S. Bass, P. B. Gossiaux, and K. Werner, Phys. Rev. C 91, 014904 (2015), eprint 1410.5396.
  * Das et al. (2016) T. Das, D. K. Choudhury, and K. K. Pathak, Indian J. Phys. 90, 1307 (2016), eprint 1505.05634.
  * Lafferty and Rothkopf (2020) D. Lafferty and A. Rothkopf, Phys. Rev. D 101, 056010 (2020), eprint 1906.00035.
  * Katz (2015) R. Katz, Ph.D. thesis, SUBATECH, Nantes (2015).
  * Katz and Gossiaux (2014) R. Katz and P. B. Gossiaux, J. Phys. Conf. Ser. 509, 012095 (2014), eprint 1312.0881.
  * Wheeler and Feynman (1949) J. A. Wheeler and R. P. Feynman, Rev. Mod. Phys. 21, 425 (1949).
  * Sorge et al. (1989) H. Sorge, H. Stoecker, and W. Greiner, Annals Phys. 192, 266 (1989).