Charmonium production in a thermalizing heat bath

TL;DR

This paper extends the Remler formalism to study charmonium production in thermalizing environments, incorporating temperature-dependent effects and quark-antiquark separation dynamics to align with statistical model predictions.

Contribution

It introduces two additional rates to the Remler formalism, accounting for temperature effects and quark-antiquark separation, improving its applicability to heavy-ion collision scenarios.

Findings

The modified rates align the formalism with statistical model results.

Temperature dependence of Wigner density affects charmonium formation.

Quark-antiquark separation dynamics are crucial for accurate modeling.

Abstract

Using the Remler formalism for the creation of composed particles, we study charmonium production both in thermalized and thermalizing boxes, which contain charm and anticharm quarks. The thermalizing box studies include the lowering of the box temperature, the spatial diffusion of charm and anticharm quarks, which are initially confined in the central region, as well as the combination of both, what imitates heavy-ion collisions. Comparing numerical and analytical results we demonstrate that the rate of the original Remler formalism has to be supplemented by two rates to obtain, for $t\to \infty$, results, which are consistent with the statistical model predictions: i) a rate, which takes into account the temperature dependence of the Wigner density of the quarkonium during the expansion and, in the case that a heavy quark potential is not implemented in the Monte Carlo approach, ii) a rate which comes from the change of the relative distance between the heavy quark and antiquark. These results provide the basis for future applications of the Remler formalism to heavy-ion collisions.

Equations

$$P_{\Phi }(t\to \infty )=\mathrm{Tr}[{\rho }_{\Phi }{\rho }_N(t\to \infty )],$$

$$\underset{t\to 0}{\lim }{P}_{\Phi }=0,$$

$$\displaystyle\Gamma=\frac{dP_{\Phi}}{dt}={\rm Tr}\bigg{(}\frac{d\rho_{\Phi}}{dt}\rho_{N}\bigg{)}+{\rm Tr}\bigg{(}\rho_{\Phi}\frac{d\rho_{N}}{dt}\bigg{)}$$

$$\displaystyle={\rm Tr}\bigg{(}\frac{d\rho_{\Phi}}{dt}\rho_{N}\bigg{)}-i{\rm Tr}\bigg{(}\rho_{\Phi}[H,\rho_{N}]\bigg{)}\equiv\Gamma_{\rm local}+\Gamma_{\rm coll}.$$

$$H$$

$$H_{1,2}$$

$$H_{N-2}$$

$$\displaystyle\Gamma_{\rm coll}=-i{\rm Tr}\bigg{(}\rho_{\Phi}[H,\rho_{N}]\bigg{)}=i{\rm Tr}\bigg{(}\rho_{N}[H,\rho_{\Phi}]\bigg{)},$$

$$[H_{N-2},{\rho }_{\Phi }]=0,[H_{1,2},{\rho }_{\Phi }]=0,$$

$$\displaystyle\Gamma_{\rm coll}=-i\sum_{i\geq 3}{\rm Tr}\bigg{(}\rho_{\Phi}[V_{1i}+V_{2i},\rho_{N}]\bigg{)}.$$

$$\displaystyle\rho_{\Phi}\rightarrow W_{S}(r,p)=8\exp\bigg{[}-\frac{r^{2}}{\sigma^{2}}-\sigma^{2}p^{2}\bigg{]},$$

$${\rho }_N\approx \prod _{i=1}^N{h}^{3N}\delta (r_i-r_i^{\ast }(t))\delta (p_i-p_i^{\ast }(t))$$

$$\frac{d{\rho }_N}{dt}\approx \sum _i{v}_i\cdot {\nabla }_r{\rho }_N$$

$$+\sum _{i>j}\sum _{\nu }\delta (t-t_{ij}(\nu ))\{{\rho }_N(t+\epsilon )-{\rho }_N(t-\epsilon )\}$$

$${\Gamma }_{\mathrm{coll}}(t)$$

$$P_{\varphi }(t^{\prime })-P_{\varphi }(0)\approx {\int }_0^{t^{\prime }}dt\{{\Gamma }_{\mathrm{local}}(t)+{\Gamma }_{\mathrm{coll}}(t)\}.$$

$$N_{\Phi }=\frac{d}{d_1{d}_2}\int \frac{d^3{p}_1{d}^3{p}_2{d}^3{r}_1{d}^3{r}_2}{(2\pi {)}^6}$$

$$\times f_Q(r_1,p_1)f_{\bar{Q}}(r_2,p_2)W_S(r,p)$$

$$=\frac{d}{d_1{d}_2}\int \frac{d^{3}P{d}^{3}p{d}^{3}R{d}^{3}r}{(2\pi {)}^6}$$

$$\times f_Q(r_1,p_1)f_{\bar{Q}}(r_2,p_2)W_S(r,p),$$

$$(2\pi {)}^3\frac{d{N}_S}{V{d}^{3}P}=\frac{d}{{\pi }^3{d}_1{d}_2}\int d^{3}p{d}^{3}r{f}_Q(p_1)f_{\bar{Q}}(p_2)e^{-\frac{r^2}{{\sigma }^2}-{\sigma }^2{p}^2}$$

$$\displaystyle=\bigg{(}\frac{\sigma^{2}}{\pi}\bigg{)}^{3/2}\frac{d}{d_{1}d_{2}}\int d^{3}pf_{Q}(p_{1})f_{\bar{Q}}(p_{2})e^{-\sigma^{2}p^{2}},~{}~{}~{}$$

$$\int d^{3}r{e}^{-r^2/{\sigma }^2}=2\pi {\sigma }^3\Gamma (3/2)=(\pi {\sigma }^2{)}^{3/2}.$$

$$\frac{1}{d_1{d}_2}{f}_Q(p_1)f_{\bar{Q}}(p_2)e^{-{\sigma }^2{p}^2}\approx e^{-E_1/T-E_2/T-{\sigma }^2{p}^2}$$

$$\displaystyle\approx\exp\bigg{[}-\bigg{(}2m_{Q}+\frac{p_{1}^{2}}{2m_{Q}}+\frac{p_{2}^{2}}{2m_{Q}}\bigg{)}\bigg{/}T-\sigma^{2}p^{2}\bigg{]}$$

$$\displaystyle=\exp\bigg{[}-\bigg{(}2m_{Q}+\frac{P^{2}/2+2p^{2}}{2m_{Q}}\bigg{)}\bigg{/}T-\sigma^{2}p^{2}\bigg{]}.$$

$$\displaystyle(2\pi)^{3}\frac{dN_{S}}{Vd^{3}P}=d\exp\bigg{[}-\bigg{(}M+\frac{P^{2}}{2M}\bigg{)}\bigg{/}T\bigg{]}$$

$$\displaystyle\times\bigg{(}\frac{\sigma^{2}}{\pi}\bigg{)}^{3/2}\int d^{3}pe^{-\{\sigma^{2}+1/(m_{Q}T)\}p^{2}}$$

$$\displaystyle=d~{}e^{-E/T}\bigg{(}\frac{\sigma^{2}}{\sigma^{2}+1/(m_{Q}T)}\bigg{)}^{3/2}$$

$$(2\pi {)}^3\frac{d{N}_S}{V{d}^{3}P}\approx d{e}^{-E/T}.$$

$$m_Q\gg \frac{3}{8⟨r^2⟩T}.$$

$$\displaystyle P_{\Phi}={\rm Tr}[\rho_{\Phi}\rho_{N}]={\rm Tr}\bigg{(}|\Phi\rangle\langle\Phi|$$

$$\displaystyle\times|Q_{1}\bar{Q}_{2}...Q_{N-1}\bar{Q}_{N}\rangle_{i}\rho_{ij}\langle Q_{1}\bar{Q}_{2}...Q_{N-1}\bar{Q}_{N}|_{j}\bigg{)}$$

$$\displaystyle=\bigg{|}\langle\Phi|Q_{1}\bar{Q}_{2}...Q_{N-1}\bar{Q}_{N}\rangle_{i}\bigg{|}^{2}\rho_{ii},$$

$$\mathcal{W}(r_1^{\ast },r_2^{\ast },p_1^{\ast },p_2^{\ast };t)\equiv \sum _{i=1,2}\sum _{j\ge 3}\int d^3{r}_1{d}^3{p}_1...d^3{r}_N{d}^3{p}_N$$

$$\times (2\pi {)}^{3N}{\rho }_{\Phi }(r_1,p_1;r_2,p_2){\rho }_N(t),$$

$$\mathcal{W}(r_1^{\ast },r_2^{\ast },p_1^{\ast },p_2^{\ast };t)={\mathcal{W}}_p(p_1^{\ast },p_2^{\ast };t){\mathcal{W}}_r.$$

Full text

I Introduction

Quarkonium is a bound state of a heavy quark and its antiquark. Since the object is flavorless, it is often called hidden heavy flavor meson. Several decades ago Matsui and Satz suggested the suppression of quarkonium production in heavy-ion collisions as a signature for the formation of a quark-gluon plasma (QGP), because the color screening, which exists only in a deconfined phase, prevents the heavy quark pair from forming a bound state Matsui and Satz (1986).

Quarkonium suppression was indeed later observed at the Super Proton Synchrotron (SPS) and the Relativistic heavy-ion Collider (RHIC) Alessandro et al. (2005); Adare et al. (2007). As the collision energy increases, however, more and more heavy quarks are produced and the possibility of regeneration of quarkonia emerges also because the c and $\bar{c}$ from different primary vertices may form a quarkonium. In fact, the nuclear modification factor of $J/\psi$, which measures the normalized ratio of $J/\psi$ produced in heavy-ion reactions as compared to pp collisions, is larger at the Large Hadron Collider (LHC) than at RHIC and larger in mid-rapidity than in forward and backward rapidities, though the temperature is higher at LHC and in the mid-rapidity region  Manceau (2013). Therefore the study of quarkonium production in heavy-ion collisions is not simply focused on the dissociation in a hot dense nuclear matter but also on regeneration in a QGP Grandchamp and Rapp (2002); Yan et al. (2006); Linnyk et al. (2008); Song et al. (2011) or when the QGP hadronizes.

In the seventies Remler devised a formalism to study the production of composite particles in heavy-ion collisions by using the Wigner representation of density operators Remler and Sathe (1975); Remler (1975, 1981), which was successfully applied to deuteron production in heavy-ion collisions Remler and Sathe (1975); Gyulassy et al. (1983).

Recently an attempt has been made to use the Remler formalism to study quarkonium production in heavy-ion collisions Villar et al. (2022); Villar (2021) as well as in p+p collisions Song et al. (2017). A distinguished feature of this approach is that the temperature and therefore in a QGP time dependence of the Wigner density of the quarkonium is taken into account. Such a dependence is predicted by lattice gauge calculations Lafferty and Rothkopf (2020), which show that below the dissociation temperature (above which a $J/\psi$ is not stable anymore) the root-mean-square (rms) radius of a $J/\psi$ depends on the temperature.

In this study we apply the Remler formalism for quarkonium production in a thermal box by using Monte-Carlo methods. Contrary to simulations of heavy-ion collisions, box simulations have the advantages that everything is controllable and analytic solutions are available. The equilibration of the quarkonia in the box is achieved by scattering with virtual particles at a given temperature. The asymptotic distribution can then be compared with statistical model predicitions. The statistical model, which has successfully described the particle production in heavy-ion collisions, provides reliable solutions for the quarkonium production Andronic et al. (2006, 2007).

We first briefly review the Remler formalism in Sec. II and present the solutions for quarkonium production in a thermalized box in Sec. III. In the next section we carry out box simulations for four different initial conditions and present our results. In Sec. V the necessity of an additional term, which is responsible for spatial diffusion, is explained. Finally, a summary is given in Sec. VI.

II Remler formalism

In the Remler formalism the multiplicity of a quarkonium state $\Phi$, $P_{\Phi}$, in a $N$-body system, composed of heavy (anti)quarks, is asymptotically (for $t\to\infty$) obtained by

[TABLE]

where $\rho_{\Phi}$ (which is assumed to be time independent) and $\rho_{N}$ are, respectively, the density operator of the quarkonium and of the $N$ heavy (anti)quarks.

In practice we cannot solve the time evolution of the quantal N-body density matrix without approximations. In the past it turned out that very satisfying results are obtained, if one does not study the time evolution of the density matrix itself but that of the Wigner density distribution, the Fourier transform of the density matrix, and approximates, guided by the mathematical work on the solution of the quantal Vlasov equation, the quantal Wigner density distribution by an ensemble of classical phase space distributions of point like particles. Averaging over this ensemble, one can calculate observables, which agree remarkable well with the experimental results. This procedure is the background of the Boltzmann-Uehlung -Uhlenbeck (BUU) (cf. Aichelin and Bertsch (1985); Cassing (2021)) and the Vlasov-Uehling-Uhlenbeck (VUU) Kruse et al. (1985) approach, which are widely used to describe the results of heavy-ion collisions with center of mass energies between few GeV and several TeV (cf. Petersen et al. (2019); Linnyk et al. (2008); Bleicher and Bratkovskaya (2022)).

In these approaches the particles scatter and move on curved trajectories due to a mean potential, created by the fellow particles. If this potential is absent or neglected one talks about a cascade approach. If heavy (anti)quarks are described in the cascade mode, where only scattering and free streaming but no potential interaction is present, they will have diverging trajectories. Therefore, even if the heavy (anti)quarks are initially confined to a space region, asymptotically we will find

[TABLE]

because two and more body correlations are lost in this approach. To overcome this problem one introduces a rate

[TABLE]

The Hamiltonian is decomposed into

[TABLE]

where $K_{i}$ and $V_{ij}$ are, respectively, the kinetic and interaction terms and

[TABLE]

Using the cyclic property of traces,

[TABLE]

and supposing for simplicity that $\Phi$ contains the particles 1 and 2 we obtain

[TABLE]

because $H_{N-2}$ does not affect $\rho_{\Phi}$ and $\rho_{\Phi}$ is a eigenfunction of $H_{1,2}$. The collision term in Eq. (3) is then simplified to Gyulassy et al. (1983); Villar et al. (2022)

[TABLE]

For a S-state the density operator of the quarkonium, expressed in Wigner representation, is approximated by

[TABLE]

where $r$ and $p$ are, respectively, relative distance and relative momentum in the center-of-mass frame. The width $\sigma$ is given by the rms radius of the quarkonium.

The classical phase space distribution of point like particles can be expressed as

[TABLE]

where $r_{i}^{*}(t)$ and $p_{i}^{*}(t)$ is the phase space trajectory of particle $i$. The time derivative of the density matrix is then given by

[TABLE]

where $t_{ij}(\nu)$ is the time for the $\nu$ th collision between the particles $i$ and $j$. The first term implies free streaming between instant scatterings, which are described by the second term. The second term in Eq. (11) is equivalent to Eq. (8). Therefore,

[TABLE]

and we obtain for the multiplicity at time t’

[TABLE]

We note that $\Gamma_{\rm local}(t)$ contributes only if $\sigma$ in Eq. (9) changes with time Villar et al. (2022).

III quarkonium in thermal box

The multiplicity of a quarkonium $\Phi$ can be obtained in the coalescence approach by projecting the phase space distribution of the charmed quarks onto the Wigner density of the $\Phi$ state Song and Berrehrah (2016); Song and Coci (2022)

[TABLE]

where $f_{Q}(r_{1},p_{1})$ and $f_{\bar{Q}}(r_{2},p_{2})$ are, respectively, the heavy quark (Q) and heavy antiquark ($\bar{Q}$) distribution functions, $r=r_{1}-r_{2}$, $R=(r_{1}+r_{2})/2$, $p=(p_{1}-p_{2})/2$ and $P=p_{1}+p_{2}$, and $d_{1}$, $d_{2}$ and $d$ is, respectively, the spin-color degeneracy of $Q$, $\bar{Q}$ and quarkonium. $r$ and $p$ in the Wigner functions are the distance and the relative momentum of the heavy quark pair in their center-of-mass frame. Since we are studying heavy quarks with a mass much larger than the temperature, we can safely use the Galilean transformation instead of the Lorentz transformation.

We assume an uniform distribution of $Q$ and $\bar{Q}$ in space,

[TABLE]

because

[TABLE]

In the nonrelativistic (or heavy quark) limit one can take a Boltzmann distribution for $f_{Q}(p_{1})$ and $f_{\bar{Q}}(p_{2})$

[TABLE]

Substituting Eq. (17) into Eq. (15) we obtain

[TABLE]

where $M=2m_{Q}$. Assuming $m_{Q}\gg 1/(\sigma^{2}T)$ we find,

[TABLE]

Since $\sigma^{2}=8/3\langle r^{2}\rangle$ with $r$ being the quarkonium radius Song and Berrehrah (2016); Song and Coci (2022), the assumption is justified if

[TABLE]

It is hence valid even close to the critical temperature $T_{c}$, supposing $\sqrt{\langle r^{2}\rangle}\sim$ 0.5 fm. We note that from Eq. (18) that the charmonium abundance at $T=200$ MeV is about 77 % of the statistical model abundance for $m_{Q}=1.5$ GeV.

Therefore we expect that in a heat bath the statistical model and the coalescence approach yield a similar multiplicity. Strictly speaking, there must be an attractive force between the heavy quark and the heavy antiquark to form a quarkonium and the mass of the quarkonium is therefore less than twice the heavy quark mass. In this study, however, we assume that the binding energy is small and therefore the $J/\psi$ is weakly bound in a QGP Lee et al. (2014); Gubler et al. (2020); Song et al. (2021). The inclusion of heavy quark potential in the Remler formalism Villar et al. (2022) will be discussed in a future publication.

The coalescence model of Eq. (14) is closely related to Eq. (1) Song and Berrehrah (2016), since

[TABLE]

where $\rho_{ij}$ is the density matrix element of the $N-$body density matrix $\rho_{N}$ and the eigenstates $|\Phi\rangle$ are taken for a basis of the $c\bar{c}$ states.

The Remler formalism starts from the same equation (1). Therefore we can test this formula in a thermal box, which is completely controllable and for which explicit and analytical solutions of Eq. (19) can be obtained. That the numerical realization of the Remler formalism gives the correct result in a box is a prerequisite for its application in numerical simulation of heavy-ion collisions.

IV box simulations

In this section the Remler formalism is tested in a box in which the heavy (anti)quarks are in thermal equilibrium as well as for three scenarios in which their initial momentum distribution and/or their initial spatial distribution does not correspond to the equilibrium distribution.

IV.1 static thermal box

We prepare a box of size $100^{3}~{}{\rm fm^{3}}$ in which we place charm quarks and charm antiquarks in thermal equilibrium at $T=$ 200 MeV. Charm (anti)quarks scatter off artificial partons which have a thermal momentum of $T=$ 200 MeV and are not affected by the scattering. The interaction rate is fixed to 1.0 c/fm. To remove boundary effects, we extended the box by 3 fm in each direction but for the analysis we exclude the extended volume.

Figure 1 compares the charmonium multiplicity as a function of time from statistical model calculations, assuming that the quarkonium mass is twice the heavy quark mass, as in Eq. (19), from Eq. (14) and from Eq. (13), which are named in the figure ”statistical model”, ”Wigner density” and ”time integration of $\Gamma$”, respectively. The charm quark mass is assumed to be 1.5 GeV and and we take a charmonium radius of 0.5 fm. The results are an ensemble average of 200 events. The colored band indicates the statistical error. One can see a slight difference between the statistical model and the Wigner projection, which originates from the approximations in Eqs. (17) and (19).

One can see that the time integration of $\Gamma_{\rm coll}$ does not deviate from the dashed and solid blue lines within the statistical error. The reason we will explain now.

We define the Wigner projection at $t$,

[TABLE]

where $r_{1}^{*},r_{2}^{*},p_{1}^{*},p_{2}^{*}$ characterize trajectories of (anti)cham quarks in Eq. (10) projected to quarkonium state, see Eq. (10).

Since the projection is carried out in a homogeneous box, one can separate the spatial dependence, which will be a constant in time, such that

[TABLE]

Then Eq. (13) is expressed as

[TABLE]

where $t=0$ is the initial projection time where only the production term appears, which corresponds to $P_{\phi}(0)$ in Eq. (II). We number the scatterings of heavy quarks or antiquarks by ”$i$”. $t_{i}$ is the time of the $i$’th scattering in which either the charm quark or the anticharm quark is involved and both, production and destruction term, appear. We note that $\mathcal{W}_{p}(p_{1}^{*},p_{2}^{*};t)$ in the above equation means $\mathcal{W}_{p}(p_{1}^{*}(t),p_{2}^{*}(t))$.

Since there is no scattering of charm quarks with a light quark between $t=0$ and $t=t_{1}-\varepsilon$,

[TABLE]

The same applies to each following time interval between collisions:

[TABLE]

Thus, most of the terms cancel and only the Wigner projection at the end of the last collision remains:

[TABLE]

where $p_{1}^{*}$ and $p_{2}^{*}$ are still the thermal momenta of the charm quark and anticharm quark. The time integration of $\Gamma_{\rm coll}$ fluctuates around the Wigner projection for a system in thermal equilibrium at T=200 MeV.

One can see that the statistical error in Fig. 1 increases with time. The reason is as following: Whenever scattering happens, a new Wigner projection is added and the old Wigner projection is subtracted as in Eq. (24). Since the box is already in thermal equilibrium both the addition and the subtraction are random fluctuations. As a result, like for all random walks, some events deviate far away from the thermal average as time passes. That is why the statistical error increases with time, while the average value stays near the thermal equilibrium.

IV.2 cooling down of (anti)charm quarks

Now we apply the Remler formalism to charm and anticharm quarks in a box of the same size in which the initial temperature of (anti)charm quarks is 400 MeV, but their number is same as in the previous subsection. In other words, only their thermal momentum changes. Initially it is given by a thermal distribution at $T=$ 400 MeV and then cools down to 200 MeV through the scattering with the artificial partons, which are assumed to have a thermal distribution with $T=$ 200 MeV.

In this case Eq. (26) is still valid. The only difference is that the initial momentum distributions of $p_{1}^{*}$ and $p_{2}^{*}$ correspond to the thermal distribution of $T=$ 400 MeV. Then they gradually change with time, through scattering with the artificial partons, to the distribution corresponding to $T=$ 200 MeV. Simply expressed, though it is not quite true, Eq. (24) will be like

[TABLE]

However, if the radius or $\sigma$ in the Wigner function of Eq. (9) depends on the temperature, Eqs. (25) and (26) are not valid any more:

[TABLE]

because the temperature at $t=t_{i}+\varepsilon$ is different from that at $t=t_{i+1}-\varepsilon$. Therefore, it is necessary to add the rate $\Gamma_{\rm local}$ as in Eqs. (3) and (13), which is expressed by Villar et al. (2022); Villar (2021)

[TABLE]

Then Eq. (24) changes to

[TABLE]

Since nothing changes between $t=0$ and $t=t_{1}-\varepsilon$ except of the temperature dependent $\sigma$

[TABLE]

ignoring the temperature change between $t=t_{1}-\varepsilon$ and $t=t_{1}$. Therefore, one arrives the same result as in Eq. (27):

[TABLE]

The upper panel of Fig. 2 shows the time evolution of the effective temperature, which is defined by

[TABLE]

where $\langle E_{kin}\rangle$ is the mean value of the (anti)charm kinetic energy. We note that due to the approximation of Eq. (34) the initial and the final effective temperatures are a bit higher than the real initial and final temperatures, which are, respectively, 400 MeV and 200 MeV. One can see that the temperature reaches its final value before $t=$ 100 fm/c. The radius of quarkonium is simply modeled as

[TABLE]

such that $\sqrt{\langle r^{2}\rangle}=0.5$ fm at $T=$ 200 MeV, a reasonable approximation to the lattice results Lafferty and Rothkopf (2020). The lower panel of Fig. 2 corresponds to Fig. 1 . We have added the magenta line, which is the result if we apply only the collisional rate, $\Gamma_{\rm coll}$, whereas the orange line is obtained if we take the sum of both rates, $\Gamma_{\rm local}$ and $\Gamma_{\rm coll}$. The multiplicity of charmonia starts from a lower value than in equilibrium at $T=$ 200 MeV because a larger thermal momentum lowers the Wigner projection. The orange line recovers the equilibrium multiplicity at about the same time when the temperature of the box reaches its final value. One can clearly see that the inclusion of $\Gamma_{\rm local}$ is necessary to obtain a results consistent with the assumed equilibrium.

IV.3 expanding (anti)charm quarks

Now we confine the initial (anti)charm quarks within a smaller box of size $50^{3}~{}{\rm fm^{3}}$ in the center of the large box, assuming a momentum distribution of the (anti)charm quarks, corresponding to $T=$ 200 MeV. The density of the charm and anticharm quarks in this smaller box is therefore 8 times higher as compared to the above discussed configuration.

As time passes the charm density decreases and the number of charmonia converges to that expected for a system in equilibrium in the full volume. One can see in Fig. 3 that the multiplicity, which is given by the solid blue line, indeed converges to the number expected for a system in equilibrium. The time integration of $\Gamma_{\rm coll}$, however, does not catch up with the decrease and remains higher. The reason is that $\mathcal{W}_{r}$ in Eq. (24) now depends on time and we obtain

[TABLE]

Eqs. (25) and (26) are still valid, but the two terms do not cancel any longer:

[TABLE]

because $\mathcal{W}_{r}(t_{i}+\varepsilon)$ is larger than $\mathcal{W}_{r}(t_{i+1}-\varepsilon)$ due to the spatial diffusion of (anti)charm quarks. This is why the equilibrium distribution of Eq. (27) cannot be achieved in this scenario.

IV.4 cooling and expanding (anti)charm quarks

Now we imitate heavy-ion collisions by combining the two previous scenarios, in other words, heavy (anti)quarks are initially at a high temperature and densely populated in a small volume. As before the initial temperature is given by 400 MeV and the initial volume by $50^{3}~{}{\rm fm^{3}}$. Then they cool down and expand in space.

Therefore the projection probability first increases with time due to the momentum loss, as in Fig. 2, and then decreases due to the spatial diffusion, as in Fig. 3. However, one can see in Fig. 4 the multiplicity, obtained by the time integration of $\Gamma_{\rm coll}$, differs from that obtained by applying the Wigner projection directly and from that, which is given by the statistical model.

V spatial diffusion term

The discrepancies between the statistical model and the time integration of $\Gamma$ in Figs. 3 and 4 originate from the calculations of $d\rho_{N}/dt$. Since $\rho_{N}$ is the density operator of $N$ (anti)charm quarks, it includes both spatial and momentum information. But only momentum space information has been taken from the comparison between Eq. (8) and Eq. (11). In other words, only the interaction terms are taken into account and the kinetic terms (the free streaming) are neglected in Eq. (11). The spatial diffusion is attributed to the kinetic terms in $K_{1}$ and $K_{2}$ of Eq. (5). In principle, it must not contribute to $d\rho_{N}/dt$ because heavy quark and heavy antiquak are bound by the potential $V_{12}$ and move together as shown in Eq. (7). However, in standard cascade simulations we have only free propagation and instant scatterings Gyulassy et al. (1983) and it is still challenging to properly implement microscopic potentials in numerical simulations Villar et al. (2022); Krenz et al. (2018).

We therefore propose that the following term should be added to Eq. (3) if a standard cascade approach is employed

[TABLE]

and for an S-state it will be (see Eq. (9)) of the form

[TABLE]

where $\vec{v}=d\vec{r}/dt$ with $\vec{r}$ being $\vec{r}_{Q}-\vec{r}_{\bar{Q}}$ in their center-of-mass frame. Then the inequality of Eq. (37) can be removed through

[TABLE]

In fact, if the bound state is perfectly described in the simulations, even $\Gamma_{\rm local}$ is not needed, because the equality of Eq. (29) will dynamically be restored:

[TABLE]

As temperature decreases, the binding of quarkonium will be stronger so that the relative distance $r$ decreases and the relative momentum $p$ increases, which compensates the change of $\sigma$ with temperature.

Figure 5 shows the time evolution of the multiplicity of $\Phi$ including $\Gamma_{\rm diff}$ in comparison the cases already studied in figures 3 and 4. One finds that including $\Gamma_{diff}$ the time integration of $\Gamma$ is now in good agreement with the statistical model predictions. Considering the results of our study of the $\Phi$ multiplicity in a box, the best and simplest method to obtain the asymptotically correct values is to add to $\Gamma$ a diffusion rate $\Gamma_{\rm diff}$

[TABLE]

where $\nu$ means $\nu$ th scattering of $i=1$ or of $i=2$.

Consequently, to be consistent with statistical model predictions, one has to supplement the old projection probability of particles $i=1,2$ (the second term in the bracket) by a new one (the first term in the bracket) whenever a scattering happens to $i=1$ or to $i=2$. Then the change of the temperature, reflected in the change of $\sigma$, and of the spatial separation of (anti)charm quarks between $t=t_{ij}(\nu-1)$ and $t=t_{ij}(\nu)$ will completely be canceled. This approach is also more natural than Eq. (II), which assumes an instant interaction between $t-\varepsilon$ and $t+\varepsilon$, because it increases $\varepsilon$ to the time between two scatterings. We note that the combination of Eq. (32) and Eq. (40) yields

[TABLE]

VI summary

The Remler formalism has been advanced to study the production of $J/\psi$ in a thermalized expanding system. In this study we have tested the Remler approach for $J/\psi$ production in thermalized and thermalizing boxes, composed of c and $\bar{c}$ quarks. The goal was to verify whether the numerical, Monte Carlo based, realization of the Remler algorithm gives the right asymptotic solution for $t\to\infty$ , which can be calculated analytically.

We started out with calculations in a completely thermalized box, where we find that the Remler formula produces the results which are consistent with that of a statistical model calculation. Then three different types of thermalizing boxes have been investigated.

In the first scenario the initial temperature of the charm quarks is high. They cool down with time through collisions with background particles and finally reach thermal equilibrium at a lower temperature. The multiplicity of the charm quarks is kept constant. We have found that the temperature derivative of the Wigner function is required to assure that the Remler approach agrees for large times with statistical model calculations.

In the second scenario charm and anticharm quarks are initially concentrated in a smaller box and then diffuse in space. The last scenario, which we studied, is the combination of the first and second one: The initial temperature of the charm quark is higher than that of the background particles and (anti-)charm quarks are concentrated initially in the central region of the box. Then they cool down and diffuse in space. This is a simple model for the expansion of the midrapidity QGP, which is created in heavy-ion collisions. We have found that in the second and third scenarios the discrepancy between the multiplicities, calculated in the Remler approach and in a statistical model, does not disappear even for $t\to\infty$.

We identified the origin of this discrepancy: it is caused by the fact that in the numerical realization of the Remler algorithm, as presented in Gyulassy et al. (1983) for deuterons, the expansion of the system between two subsequent collisions is not taken properly into account. Neglecting the potential, it treats the c and $\bar{c}$ quarks as freely moving particles between two collisions whereas in reality they are bound when they form a quarkonium. Introducing a diffusion rate, which adds to the local rate and the collision rate, this discrepancy disappears. We note that recently also an approach has been advanced which includes the $c\bar{c}$ potential in an approximate way Villar et al. (2022). It would be interesting to verify whether there equilibrium is obtained without a spacial diffusion rate. In conclusion, we have found that also for an expanding system, which cools down, the Monte Carlo realization of the Remler formalism describes correctly the approach to equilibrium .

Acknowledgements

The authors acknowledge valuable discussions with P.-B. Gossiaux. We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the grant CRC-TR 211 ’Strong-interaction matter under extreme conditions’ - Project number 315477589 - TRR 211. This work is supported by the European Union’s Horizon 2020 research and innovation program under grant agreement No 824093 (STRONG-2020). The computational resources have been provided by the LOEWE-Center for Scientific Computing and the ”Green Cube” at GSI, Darmstadt and by the Center for Scientific Computing (CSC) of the Goethe University.