text:
# Spin-Wave Description of Haldane-gap antiferromagnets

## Abstract

Modifying the conventional antiferromagnetic spin-wave theory which is plagued by the difficulty of the zero-field sublattice magnetizations diverging in one dimension, we describe magnetic properties of Haldane-gap antiferromagnets. The modified spin waves, constituting a grand canonical bosonic ensemble so as to recover the sublattice symmetry, not only depict well the ground-state correlations but also give useful information on the finite-temperature properties.

PACS numbers: 75.10.Jm, 05.30.Jp, 75.40.Mg

Haldane’s conjecture  that one-dimensional Heisenberg antiferromagnets should exhibit qualitatively different low-energy structures based on whether the constituent spins are integral or fractional sparked renewed interest in the field of quantum magnetism. An energy gap immediately above the ground state was indeed observed in a quasi-one-dimensional Heisenberg antiferromagnet Ni(C2H8N2)2NO2(ClO4)  and a rigorous example of such a massive phase was also found out . The energy gaps in magnetic excitation spectra, that is, spin gaps, are now one of the most attractive and important topics. In the context of theoretical progress, we may be reminded of quantized plateaux in the ground-state magnetization curves , a dramatic crossover from one- to two-dimensional quantum antiferromagnets , and an antiferromagnetic excitation gap accompanied by ferromagnetic background . From the experimental point of view, metal oxides such as spin-Peierls compounds Cu$`{}_{1x}{}^{}M_{x}^{}`$GeO3 ($`M=\text{Zn},\text{Mg}`$) , Haldane-gap antiferromagnets $`R_2`$BaNiO5 ($`R=\text{rare earth}`$)  and ladder materials Srn-1Cun+1O2n ($`n=3,5,7,\mathrm{}`$)  have significantly contributed to systematic investigations of the mechanism of gap formation.

In studying Haldane-gap antiferromagnets, numerical tools such as quantum Monte Carlo and density-matrix renormalization-group methods are indeed useful, but analytic approaches still play an important role. The nonlinear $`\sigma `$-model technique  is powerful enough to investigate low-energy structures with particular emphasis on the competition between massive and massless phases. It is the valence-bond-solid description  of integer-spin chains that enables us to readily understand the hidden order  inherent in the Haldane massive phase. However, these methods are not so useful for exploring thermal properties as they are in investigating the ground-state properties. Although we may consider applying the spin-wave scheme to Haldane-gap antiferromagnets in such circumstances, the difficulty of the zero-field sublattice magnetizations diverging in one dimension has been a problem in the conventional spin-wave theory  for years. Therefore, it was a major breakthrough that spin waves, being constrained to keep zero magnetization, succeeded in precisely describing the low-temperature thermodynamics of the one-dimensional Heisenberg ferromagnet . This new spin-wave scheme, which is now referred to as the modified spin-wave theory, was further applied to quantum antiferromagnets and ferrimagnets. The modified spin-wave scheme is highly successful for extensive ferrimagnets  and still applies well for two-dimensional antiferromagnets . As for one-dimensional antiferromagnets, there exists a pioneering argument , but it looks unsatisfactory for interpreting experimental and numerical observations.

In this article, we aim to do away with our vague but persistent impression that the spin-wave scheme hardly works for one-dimensional quantum antiferromagnets. Spin-$`2`$ Haldane-gap antiferromagnets , as well as those of spin $`1`$, have recently been synthesized and more explorations into novel quantum phenomena in one dimension are expected in the future. In such circumstances, we make our first attempt to construct a modified spin-wave theory for Haldane-gap antiferromagnets.

We consider integer-spin antiferromagnetic Heisenberg chains

$$=J\underset{j=1}{\overset{L}{}}𝑺_j𝑺_{j+1};𝑺_j^2=S(S+1).$$
(1)
We define bosonic operators for the spin deviation in each sublattice via

$$\begin{array}{ccc}S_{2n1}^+=\sqrt{2Sa_n^{}a_n}a_n,\hfill & S_{2n1}^z=Sa_n^{}a_n,\hfill & \\ S_{2n}^+=b_n^{}\sqrt{2Sb_n^{}b_n},\hfill & S_{2n}^z=S+b_n^{}b_n.\hfill & \end{array}$$
(2)
The Fourier-transformed operators are introduced as

$$\begin{array}{c}a_k=\frac{1}{\sqrt{N}}\underset{n=1}{\overset{N}{}}\mathrm{e}^{\mathrm{i}k(n1/4)}a_n,\hfill \\ b_k=\frac{1}{\sqrt{N}}\underset{n=1}{\overset{N}{}}\mathrm{e}^{\mathrm{i}k(n+1/4)}b_n,\hfill \end{array}$$
(3)
where twice the lattice constant $`2a`$ is set equal to unity and therefore $`k=2\pi n/N(n=0,1,\mathrm{},N1;N=L/2)`$. Through the Bogoliubov transformation

$$\begin{array}{ccc}a_k^{}\hfill & =\hfill & \alpha _k^{}\mathrm{cosh}\theta _k\beta _k\mathrm{sinh}\theta _k,\hfill \\ b_k^{}\hfill & =\hfill & \beta _k^{}\mathrm{cosh}\theta _k\alpha _k\mathrm{sinh}\theta _k,\hfill \end{array}$$
(4)
we obtain

$$=E_\mathrm{N}+E_1+E_0+_1+_0+O(S^1),$$
(5)
where

$`E_\mathrm{N}`$ $`=`$ $`2NJS^2,`$ (7)
$`E_1`$ $`=`$ $`4NJS(\mathrm{\Gamma }\mathrm{\Lambda }),`$ (8)
$`E_0`$ $`=`$ $`2NJ(\mathrm{\Gamma }\mathrm{\Lambda })^2,`$ (9)
$`_1`$ $`=`$ $`J{\displaystyle \underset{k}{}}[\omega _1(k)(\alpha _k^{}\alpha _k+\beta _k^{}\beta _k)`$ (12)
$`+\gamma _1(k)(\alpha _k\beta _k+\alpha _k^{}\beta _k^{})],`$
$`_0`$ $`=`$ $`J{\displaystyle \underset{k}{}}[\omega _0(k)(\alpha _k^{}\alpha _k+\beta _k^{}\beta _k)`$ (14)
$`+\gamma _0(k)(\alpha _k\beta _k+\alpha _k^{}\beta _k^{})],`$
with

$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{k}{}}\mathrm{cos}{\displaystyle \frac{k}{2}}\mathrm{sinh2}\theta _k,`$ (16)
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{k}{}}(\mathrm{cosh2}\theta _k1),`$ (17)
$`\omega _1(k)`$ $`=`$ $`2S\left(\mathrm{cosh2}\theta _k\mathrm{cos}{\displaystyle \frac{k}{2}}\mathrm{sinh2}\theta _k\right),`$ (19)
$`\omega _0(k)`$ $`=`$ $`2(\mathrm{\Gamma }\mathrm{\Lambda })\left(\mathrm{cosh2}\theta _k\mathrm{cos}{\displaystyle \frac{k}{2}}\mathrm{sinh2}\theta _k\right),`$ (20)
$`\gamma _1(k)`$ $`=`$ $`2S\left(\mathrm{cos}{\displaystyle \frac{k}{2}}\mathrm{cosh2}\theta _k\mathrm{sinh2}\theta _k\right),`$ (22)
$`\gamma _0(k)`$ $`=`$ $`2(\mathrm{\Gamma }\mathrm{\Lambda })\left(\mathrm{cos}{\displaystyle \frac{k}{2}}\mathrm{cosh2}\theta _k\mathrm{sinh2}\theta _k\right).`$ (23)
The naivest diagonalization of the Hamiltonian (5), whether up to $`O(S^1)`$ or up to $`O(S^0)`$, results in diverging sublattice magnetizations even at zero temperature. Although eq. (2) assumes that $`a_n^{}a_n2S`$ and $`b_n^{}b_n2S`$, the conventional spin-wave theory cannot reasonably control the boson numbers. Then we consider introducing a grand canonical constraint to the noncompact Hamiltonian (5).

Isotropic magnets should lie in the state of zero magnetization $`_jS_j^z=0`$ and the minimization of the free energy under such a condition indeed yields an excellent description of the low-temperature thermodynamics of ferromagnets . However, in the cases with antiferromagnetic exchange interactions, the zero-magnetization constraint, claiming that $`_n(a_n^{}a_nb_n^{}b_n)=0`$, still fails to overcome the divergence of the numbers of the sublattice bosons. Hence we minimize the free energy constraining the sublattice magnetizations to be zero:

$$\underset{n}{}a_n^{}a_n=\underset{n}{}b_n^{}b_n=SN.$$
(24)
Within the conventional spin-wave theory, spins on one sublattice point predominantly up, while those on the other predominantly down. This condition (24) restores the sublattice symmetry. In order to enforce the constraint (24), we first introduce a Lagrange multiplier and diagonalize

$$\stackrel{~}{}=+2J\lambda \underset{k}{}(a_k^{}a_k+b_k^{}b_k).$$
(25)
Then the ground-state energy and the dispersion relation are obtained as

$`E_\mathrm{g}=E_\mathrm{N}+\stackrel{~}{E}_1;\stackrel{~}{E}_1=E_1+4NJ\lambda \mathrm{\Lambda },`$ (26)
$`\omega (k)=\stackrel{~}{\omega }_1(k);\stackrel{~}{\omega }_1(k)=\omega _1(k)+2\lambda \mathrm{cosh2}\theta _k,`$ (27)
within the linear modified spin-wave scheme and as

$`E_\mathrm{g}=E_\mathrm{N}+\stackrel{~}{E}_1+E_0,`$ (28)
$`\omega (k)=\stackrel{~}{\omega }_1(k)+\omega _0(k),`$ (29)
within the up-to-$`O(S^0)`$ interacting modified spin-wave scheme. For eqs. (26) and (27), $`\theta _k`$ is given by $`\gamma _1(k)2\lambda \mathrm{sinh2}\theta _k\stackrel{~}{\gamma }_1(k)=0`$, whereas for eqs. (28) and (29), $`\theta _k`$ may be determined in two ways. One idea is the perturbational treatment of $`_0`$, which is referred to as the perturbational interacting modified spin-wave scheme, where $`\theta _k`$ is still given by $`\stackrel{~}{\gamma }_1(k)=0`$ and the $`O(S^0)`$ quantum correction is the $`O(S^1)`$-eigenstate average of $`_0`$. The other is the full diagonalization of $`_1+_0`$, which is referred to as the full-diagonalization interacting modified spin-wave scheme, where $`\theta _k`$ is given by $`\stackrel{~}{\gamma }_1(k)+\gamma _0(k)=0`$. Once $`\theta _k`$ is given, we calculate the free energy and obtain the optimum thermal distribution functions as

$$\alpha _k^{}\alpha _k=\beta _k^{}\beta _k\overline{n}_k=\frac{1}{\mathrm{e}^{J\omega _k/k_\mathrm{B}T}1},$$
(30)
where $`\lambda `$ is self-consistently determined by the condition

$$\underset{k}{}\left(2\overline{n}_k+1\right)\mathrm{cosh2}\theta _k=(2S+1)N.$$
(31)
First, let us evaluate the ground-state energy. We compare the modified spin-wave calculations with the highly accurate quantum Monte Carlo estimates  in Table I. The modified spin-wave findings are generally in good agreement with the quantum Monte Carlo results. The interacting modified spin waves describe the ground-state correlations much better than the linear ones and their description becomes increasingly refined with increasing $`S`$. The interacting modified spin-wave findings miss the correct value by 0.5 percent for $`S=1`$ and by only 0.008 percent for $`S=3`$.

Secondly, we consider the Haldane gap $`\omega (ak=\pi )\mathrm{\Delta }(T)`$. Table II shows that the modified spin-wave scheme, in contrast with the conventional spin-wave theory, succeeds in generating the gap but considerably underestimates it. Unavailability of the absolute energy level structure is an inevitable consequence of our employing the effective Hamiltonian (25). Then we examine the present scheme by scaling $`\mathrm{\Delta }(T)`$ to its zero-temperature value $`\mathrm{\Delta }(T=0)\mathrm{\Delta }_0`$. Such an argument is quite usual with field-theoretical investigations, which do not lead to an estimate of the normalization factor but derive finite-temperature expressions involving only ratios such as $`k_\mathrm{B}T/\mathrm{\Delta }_0`$. In Fig. 1, we compare the modified spin-wave calculations with numerical , field-theoretical  and experimental  findings. Now the superiority of the modified spin-wave scheme is clear at a glance. Of all the theoretical tools, the full-diagonalization interacting modified spin-wave approach is the most successful to reproduce the observed upward behavior of $`\mathrm{\Delta }(T)`$ with increasing temperature . The nonlinear $`\sigma `$-model treatment  is justified well in the low-temperature region $`k_\mathrm{B}T\mathrm{\Delta }_0`$, while the maximum-entropy technique  works less with increasing temperature. Nickel compounds such as Y2BaNiO5  and Ni(C2H8N2)2NO2(ClO4)  are good candidates for spin-$`1`$ Haldane-gap antiferromagnets, but magnetic anisotropy and interchain interaction are not negligible there. They split the lowest excitation gap of the isotropic chain (1) into several levels and smear the intrinsic behavior of the ideal integer-spin chains. Considering such practical factors, the present theory satisfactorily interprets the observations.

Lastly, we show the modified spin-wave calculations of the magnetic susceptibility in Fig. 2. Considering the total breakdown of the conventional spin-wave theory in one-dimensional thermodynamic calculations, the modified spin-wave achievement is highly successful. All the calculations are guaranteed to reproduce the paramagnetic susceptibility $`\chi /Lg^2\mu _\mathrm{B}^2=S(S+1)/3k_\mathrm{B}T`$ at high temperatures. Since the interacting modified spin-wave scheme gives better estimates of the gap than the linear modified spin-wave one (Table I), it is somewhat better at describing the low-temperature behavior. With increasing $`S`$, the activation-type initial behavior is suppressed and the antiferromagnetic peak is broadened. The $`S=1`$ modified spin-wave calculations are in fine agreement with the quantum Monte Carlo findings over a wide temperature range, while those for $`S=2`$ look somewhat poorer at intermediate temperatures. It may be closely related to the fact that the modified spin-wave estimates of the $`S=2`$ gap are worse than those of the $`S=1`$ gap (Table II). However, Table II suggests that the validity of the modified spin-wave scheme for excitations significantly improves with increasing $`S`$, possibly in a staggered way. Thermodynamic calculations, whether by quantum Monte Carlo or density-matrix renormalization group, for systems with large degrees of freedom are less feasible numerically, in particular, at low temperatures. The present scheme has the advantage of saving time and computational resources.

We have demonstrated the applicability of the new spin-wave scheme to Haldane-gap antiferromagnets. This is the first comprehensive attempt to describe one-dimensional spin-gapped antiferromagnets in terms of spin waves. The modified spin waves are free from their thermal as well as quantum divergence and can therefore microscopically interpret various magnetic properties. Besides the magnetic susceptibility, the spin correlation function and the nuclear spin-lattice relaxation time can be revealed. Unfortunately, we have less information on the bare energy spectrum, because we get rid of the quantum divergence at the cost of the original Hamiltonian. However, with the zero-temperature spin gap, which is readily and precisely available through numerical calculations within a canonical ensemble , we can still argue the energy structure quantitatively.

Haldane-gap antiferromagnets in a magnetic field provide further interesting issues. With the increase of an applied field, the gap is reduced and the ground state is mixed increasingly with the first excited state. Indeed the nuclear spin-lattice relaxation is accelerated in the vicinity of the critical field, but its overall behavior as a function of an applied field and temperature  is far from understandable at a glance. We may expect the present new scheme to open the way for the total understanding of low-dimensional spin-gapped antiferromagnets.


format: text

text:
# On the 𝑥_𝐹 distribution of 𝐽/𝜓’s produced in heavy ion collisions

## I Introduction

Over the past fifteen years charmonium suppression has been considered as one of the best signatures of quark gluon plasma (QGP) formation. Recently this belief was questioned by some works. $`J/\psi `$ suppression in relativistic heavy ion collisions is based on the following simple argument, presented in the original work by Matsui and Satz . Lattice simulations show that the heavy quark-antiquark potential becomes screened at high enough temperatures. As a consequence, when the range of the potential becomes smaller than the $`J/\psi `$ Bohr radius the bound state can no longer be formed. On the other hand, detailed simulations  of a population of $`c\overline{c}`$ pairs traversing the plasma suggested that, given the large number of such pairs, the recombination effect of the pairs into charmonium Coulomb bound states is non-negligible and might even lead to an enhancement of $`J/\psi `$ production. This conclusion received support from the calculations of , where a two component model for $`J/\psi `$ production was proposed.

Taking the existing calculations seriously, it is no longer clear that an overall suppression of the number of $`J/\psi `$’s will be a signature of QGP. A more detailed analysis is required and a more complex pattern can emerge. In particular, we might have suppression in some regions of the phase space and enhancement in others. Indeed, this is the result of the analysis presented in . In that work the authors study the fate of $`c\overline{c}`$ pairs produced in the early stage of heavy ion collisions, comparing the case where they have to traverse QGP with the case where they have to traverse ordinary nuclear matter. One of the merits of that paper is to emphasize that very interesting information can be extracted from the scaled momentum ($`x_F`$) distribution of the produced $`J/\psi `$’s.

Motivated by these observations, in this work we address the $`J/\psi `$ Feynman momentum distributions. In some aspects we follow refs.  and  with one important difference: we include $`J/\psi `$ production within the plasma. This is usually neglected because the number of $`c\overline{c}`$ pairs produced in the plasma is believed to be small. However, as it will be discussed in section III, a closer look into the existing estimates shows discrepancies of two orders of magnitude. In , for example, it was estimated to be of the order of $`1\%`$ of the total number of charm quark pairs. In , with the inclusion of thermal parton masses this fraction was estimated to reach $`2030\%`$. Finally, in  this number could, in some cases, be equal to the number of directly produced pairs. No systematic effort was made to reconcile these different estimates mainly because the screening mechanism was believed to destroy all bound states. Therefore in-plasma charm production has been an issue of open charm physics whereas the main focus has always been on hidden charm production and suppression.

Although well established in several lattice calculations , the color screening mechanism may be not so effective in real collisions. The existing lattice results are valid for an infinite mass, static, charge-anticharge pair and we know that this is not very close to the real situation, where charges are not so heavy and are moving. Recently, a calculation of the screening effect in moving charges  showed that, in this case the screening is not so strong. With the improvements of lattice gauge determinations  of the temperature-dependent quark-antiquark potential new calculations of the quarkonium spectrum in QGP were performed. In  it was found that at up to temperatures of the order of $`1.10T_c`$ the $`J/\psi `$ can survive as a bound state in the plasma. At higher temperatures the potential is too shallow to hold a bound charmonium state.

We shall assume that screening is not so strong and, in a narrow temperature window just above the phase transition, will allow for charmonium (Coulomb) bound states which “survive” from QGP. More precisely, we will allow up to one in hundred bound states to survive in the region where they should be screened and destroyed. As it will be seen, these “survivors” tend to escape with $`x_F0`$, thus giving us a new kind of QGP signal. Moreover, this peak at low $`x_F`$ will fill the dip predicted in .

We will consider $`J/\psi `$ production in nucleus-nucleus collisions by three different mechanisms: i) direct (primordial QCD parton fusion) ii) thermal (statistical coalescence at hadronization) and iii) QGP (in-plasma parton fusion). As already discussed in  i) and ii) are just two extreme cases. In realistic simulations, for a fixed number of $`c\overline{c}`$ pairs and of $`J/\psi `$’s, these last are gradually destroyed, giving origin to the “regenerated” component ii). In the particular calculation presented in , for central collisions and for high enough energies a near to complete replacement of the initial by the final thermal $`J/\psi `$’s was found. Component iii) was considered long time ago , before the work of Matsui and Satz, and then it was left in oblivion. Once produced the charmonium state will suffer interactions with the partons in the plasma, or with hadrons in hot and dense medium and in a later stage with comoving hadrons .

Before performing detailed calculations one could try to guess the shape of the outcoming $`J/\psi `$’s $`x_F`$ distribution in the three mechanisms mentioned above. QCD production has been studied in the past and the naive extrapolation of the spectra measured in $`pp`$ and $`pA`$ collisions would lead us to conclude that the $`J/\psi `$’s produced in i) should follow a distribution of the type :

$$(1x_F)^c$$
(1)
where $`c56`$. As it will be discussed in the next section, perturbative QCD (pQCD) calculations within the framework of the color evaporation model (CEM) confirmed this behavior in leading order and next-to-leading order . No explicit calculation of this $`x_F`$ distribution was performed with the color octet model (COM), but a recent study of asymmetries in open charm production  suggests that a similar smooth behavior at $`x_F0`$ might be found for the $`J/\psi `$ differential cross section. Of course, this has to be verified. Thermal production follows a Bose-Einstein distribution , which integrated in $`p_T^2`$ gives:

$$exp\left[\frac{1}{T}\left(\sqrt{x_F^2s/4+m_{J/\psi }^2}\right)\right]\left[1+\frac{1}{T}\sqrt{x_F^2s/4+m_{J/\psi }^2}\right]$$
(2)
where $`x_F`$ and $`\sqrt{s}`$ are the Feynman momentum and the c.m.s. nucleon-nucleon invariant energy. Comparing (1) with (2), we can see that the latter flattens out at low values of $`x_F`$, because of the mass $`m_{J/\psi }`$, forming a “plateau”. A spectrum from mechanism iii) is a little more difficult to predict. It involves the convolution of two distributions of the type (2) (for the two colliding massless partons) and therefore we expect it to be more steeply falling with $`x_F`$. As it will be seen, within certain approximations, it is given by:

$$\frac{1}{x_F}$$
(3)
This behavior was found with the CEM. It remains to be shown that it is not an artifact of the model and that it would also be found in the thermal version of the COM. Comparing (1), (2) and (3) we see that, if QGP is formed and if the multiplicity of “in-plasma-formed” $`J/\psi `$’s is large enough, we may expect to see an enhancement (or even a peak) at the origin of the $`x_F`$ distribution. This signal is interesting because it was not predicted by any other production mechanism. However it relies on a very accurate knowledge of the pQCD, thermal and QGP abundancies, which is very difficult, as it will be dicussed in detail later.

At this point one might argue that after being produced by the usual pQCD mechanism, the $`c\overline{c}`$ bound state suffers interactions either with the plasma or with comovers in a hadron gas. These interactions (called here Final State Interactions) will distort the initial $`x_F`$ distribution of the bound state. In Ref.  it was shown that FSI (either with a plasma or a hadron gas) will suppress the charmonium yield at $`x_F=0`$, giving origin to a dip in the central $`x_F`$ region.

We are suggesting that a peak will arise at $`x_F=0`$ as a consequence of in-plasma production. This is a signal which will emerge from the background composed by initial pQCD production. If we ignore the final state interactions we will calculate a background which already contains an enhancement of $`J/\psi `$ production in the low $`x_F`$ region, making it harder to detect the signal. In other words, if we treat correctly the FSI this will lead to a background with a dip around $`x_F=0`$ instead of an enhancement. Therefore, at least for the sake of our argument, neglecting FSI is being conservative. This is the correct procedure in order to isolate the new effect, which we are looking for. Nevertheless, at the end of this work we will include FSI quantitatively and investigate their effect on the $`J/\psi `$ momentum distribution.

The text is organized in the following way: in section II we describe direct production; in sections III and IV we address production in the plasma and final thermal production respectively; in section V we discuss final state interactions; in section VI we present numerical results and in the last section we make some concluding remarks.

## II Initial $`J/\psi `$ production in nucleus-nucleus collisions

Nowadays perturbative QCD calculations of $`c\overline{c}`$ production can be found in textbooks. Nevertheless, we shall, in what follows, give some formulas to introduce the notation, to stress the role played by shadowing and to obtain expressions in terms of $`x_F`$, the variable of interest. One popular approach to charm production is the color evaporation model . We are aware of the discrepancies between different implementations of perturbative QCD in charm production . In particular, there is another model, namely the color octet model , which gives a more sophisticated description of the hadronization process in terms of non-perturbative matrix elements. While the COM seems to give a very good description of high $`p_T`$ data ($`p_T>2`$ GeV), especially those from the Fermilab Tevatron, it is not so reliable at smaller $`p_T`$  ($`p_T<2`$ GeV), where the bulk of charmonium production takes place and where we wish to make a comparison between QCD and plasma production. For our purpose of studying the $`x_F`$ distributions, as can be inferred from the calculations presented in , the COM prediction is compatible with a smooth behavior of the cross section $`d\sigma /dx_F`$ at the origin $`x_F0`$, much like the CEM. Therefore, for simplicity, we shall use this latter model in what follows. Our results will certainly depend on the choice of approach. Inspite of the uncertainties, we will assume that, since we are interested in the low $`x_F`$ region, pQCD is enough and nonperturbative effects are small, in sharp contrast to what happens at large $`x_F`$ .

In the CEM, charmonium is defined kinematically as a $`c\overline{c}`$ state with mass below the $`D\overline{D}`$ threshold. In leading order (LO) the cross section is computed with the use of perturbative QCD for the diagrams of the elementary processes $`q\overline{q}c\overline{c}`$ and $`ggc\overline{c}`$ convoluted with the parton densities in the projectile and in the target. Calling $`x_F`$ the fractional momentum of the produced pair (with respect to the momentum of a projectile nucleon in cm frame) and $`\sqrt{s}`$ the cm energy of a nucleon-nucleon collision, the cross section for production of a $`c\overline{c}`$ pair with mass $`m`$ is just given by:

$`{\displaystyle \frac{d\sigma ^{ppc\overline{c}}}{dx_Fdm^2}}`$ $`=`$ $`{\displaystyle _0^1}𝑑x_1𝑑x_2\delta (x_1x_2sm^2)\delta (x_Fx_1+x_2)H(x_1,x_2;m^2)`$ (4)
$`=`$ $`{\displaystyle \frac{1}{s\sqrt{x_F^2+4m^2/s}}}H(x_{01},x_{02};m^2);x_{01,02}={\displaystyle \frac{1}{2}}(\pm x_F+\sqrt{x_F^2+4m^2/s})`$ (5)
where $`x_1`$ and $`x_2`$ are the nucleon momentum fractions carried respectively by partons in the projetile and target. The function $`H(x_1,x_2;m^2)`$, which represents the convolution of the elementary cross sections and parton densities is given by:

$`H(x_1,x_2;m^2)`$ $`=`$ $`f_g(x_1,m^2)f_g(x_2,m^2)\widehat{\sigma }_{gg}(m^2)`$ (6)
$`+`$ $`{\displaystyle \underset{q=u,d,s}{}}[f_q(x_1,m^2)f_{\overline{q}}(x_2,m^2)+f_{\overline{q}}(x_1,m^2)f_q(x_2,m^2)]\widehat{\sigma }_{q\overline{q}}(m^2)`$ (7)
with the parton densities $`f_i(x,m^2)`$ in the nucleon computed at the scale $`m^2=x_1x_2s`$.

The LO elementary cross sections in terms of the pair invariant mass ($`m`$) are given by :

$$\widehat{\sigma }_{gg}(m^2)=\frac{\pi \alpha _s^2(m^2)}{3m^2}\left\{\left(1+\frac{4m_c^2}{m^2}+\frac{m_c^4}{m^4}\right)ln\left[\frac{1+\lambda }{1\lambda }\right]\frac{1}{4}\left(7+\frac{31m_c^2}{m^2}\right)\lambda \right\}$$
(8)
$$\widehat{\sigma }_{q\overline{q}}(m^2)=\frac{8\pi \alpha _s^2(m^2)}{27m^2}(1+\frac{2m_c^2}{m^2})\lambda ;\lambda =[1\frac{4m_c^2}{m^2}]^{1/2}$$
(9)
where $`m_c`$ is the mass of the $`c`$ quark.

The production cross section of the charmed state $`i`$ ($`=J/\psi ,\psi ^{}`$ or $`\chi _{cJ}`$) $`\sigma _i`$, is then finally obtained by integrating the free pair cross section $`c\overline{c}`$ over the invariant mass $`m`$ starting from the production threshold $`2m_c`$ up to open charm production threshold $`2m_D`$:

$$\frac{d\sigma ^{ppJ/\psi }}{dx_F}=F_{J/\psi }_{4m_c^2}^{4m_D^2}𝑑m^2\frac{d\sigma ^{ppc\overline{c}}}{dx_Fdm^2}$$
(10)
where $`F_{J/\psi }`$ is the fraction of $`\sigma ^{c\overline{c}}`$ which contains the corresponding $`c\overline{c}`$ resonance.

This model describes well the experimentally measured $`x_F`$ distribution of hidden charm both with LO and NLO cross sections, provided that $`F_i^{LO}`$ is defined as $`F_i^{NLO}`$ multiplied by a theoretical factor $`\kappa `$, which is equal to the ratio of the NLO and LO cross sections ($`F_{J/\psi }^{NLO}2,54\%`$) .

In what follows we shall use the CEM to study perturbative $`J/\psi `$ production in nucleus-nucleus collisions at RHIC. As it is well known, in nuclear collisions and in processes involving small values of $`x`$, shadowing plays an important role (see, for example, ). In this case expression (7) is rewritten as:

$`H_{AB}(x_1,x_2;m^2)`$ $`=`$ $`AB\{f_g(x_1,m^2)R_g^A(x_1,m^2)f_g(x_2,m^2)R_g^B(x_2,m^2)\widehat{\sigma }_{gg}(m^2)`$ (11)
$`+`$ $`{\displaystyle \underset{q=u,d,s}{}}[f_q(x_1,m^2)R_q^A(x_1,m^2)f_{\overline{q}}(x_2,m^2)R_{\overline{q}}^B(x_2,m^2)`$ (12)
$`+`$ $`f_{\overline{q}}(x_1,m^2)R_{\overline{q}}^A(x_1,m^2)f_q(x_2,m^2)R_q^B(x_2,m^2)]\widehat{\sigma }_{q\overline{q}}(m^2)\}`$ (13)
where

$$R_i^A(x,m^2)=\frac{f_i^A(x,m^2)}{f_i(x,m^2)}$$
(14)
with $`f_i^A(x,m^2)`$ being the $`i`$ parton momentum distribution in a nucleon inside the nucleus $`A`$. Replacing $`H`$ by $`H_{AB}`$ in (5) we obtain the cross section for $`c\overline{c}`$ production in $`AB`$ collisions:

$`{\displaystyle \frac{d\sigma ^{ABc\overline{c}}}{dx_F}}`$ $`=`$ $`\kappa {\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle \frac{1}{s\sqrt{x_F^2+4m^2/s}}}H_{AB}(x_{01},x_{02};m^2)`$ (15)
and the $`J/\psi `$ production cross section as:

$`{\displaystyle \frac{d\sigma ^{ABJ/\psi }}{dx_F}}`$ $`=`$ $`\kappa F_{J/\psi }^{NLO}{\displaystyle _{4m_c^2}^{4m_D^2}}𝑑m^2{\displaystyle \frac{1}{s\sqrt{x_F^2+4m^2/s}}}H_{AB}(x_{01},x_{02};m^2)`$ (16)
In a central $`A+A`$ collision the number of $`J/\psi `$’s produced directly is related to the cross section by :

$$\frac{dN_{direct}^{AAJ/\psi }}{dx_F}\frac{1}{\pi R_A^2}\frac{d\sigma ^{AAJ/\psi }}{dx_F}$$
(17)
where $`R_A`$ is the radius of nucleus $`A`$.

The expression above would give the final distribution if there were no absorption due to further interactions with partons in a QGP or hadrons in a comoving fireball. The effect of these final state interactions will be discussed later.

Since we are going to compare the total number of initially produced $`c\overline{c}`$ pairs and $`J/\psi `$’s with those produced in the plasma it is useful to introduce the compact notation

$$N_{QCD}^{c\overline{c}}=\frac{1}{\pi R_A^2}_0^1𝑑x_F\frac{d\sigma ^{AAc\overline{c}}}{dx_F}$$
(18)
for the number of $`c\overline{c}`$ pairs and

$$N_{QCD}^{J/\psi }=_0^1𝑑x_F\frac{dN_{direct}^{AAJ/\psi }}{dx_F}$$
(19)
for the number of $`J/\psi `$’s.

Expressions (18) and (19) are calculated with the nuclear parton distributions which take shadowing into account. While this makes calculations more realistic, it also introduces some model dependence in the results. In order to have a baseline for comparisons it is useful to introduce the equivalent definitions of $`N_{QCD}^{c\overline{c}}`$ and $`N_{QCD}^{J/\psi }`$ without shadowing. For central $`AA`$ collisions these are:

$$N_{QCD}^{c\overline{c}}=T_{AA}(b=0)_0^1𝑑x_F\frac{d\sigma ^{ppc\overline{c}}}{dx_F}\text{ (no shadowing)}$$
(20)
and

$$N_{QCD}^{J/\psi }=F_{J/\psi }^{NLO}T_{AA}(b=0)_0^1𝑑x_F\frac{d\sigma ^{ppc\overline{c}}}{dx_F}\text{ (no shadowing)}$$
(21)
where $`T_{AA}`$ is the usual nuclear overlap function. In Table I we present some quantitative results for (18), (19), (20) and (21). As it can be seen, these numbers change with different choices for the charm quark mass and for the parton distribution. Inclusion of shadowing reduces the number of charm quark pairs and $`J/\psi `$’s in about 10 %, choosing different parton densities may change these numbers in 40 % but what changes the most our results is the adopted value for the charm quark mass. This choice may alter the numbers by a factor 4. Our number of $`c\overline{c}`$ pairs is smaller but still compatible with other estimates presented, for example, in . On the other hand, the number of $`J/\psi `$’s can be estimated taking into account the recent PHENIX data  on $`pp`$ collisions. These estimates  indicate that $`N_{QCD}^{J/\psi }0.1`$, which is close to the number of $`J/\psi `$’s in Table I.

## III $`J/\psi `$ production in the plasma

### A The number of $`c\overline{c}`$ pairs

As in the previous section we reproduce below (to a great extent) textbook material , with the purpose of defining our notation. The derivation presented here is also useful because, relaxing the screening hypothesis, we calculate the production of charm bound states in the plasma and we integrate the rate on all variables except $`x_F`$. All this requires some straightforward but not very often shown manipulations.

The computation of the in-plasma $`c\overline{c}`$ pair production rate goes back to the late eighties , was discussed in short papers (for example ), included in comprehensive review articles  during the nineties and experienced improvements due to advances in thermal field theory .

Assuming that QGP is formed, we then have a gas of quarks and gluons with momenta obeying respectively Fermi-Dirac and Bose-Einstein distributions that can collide producing $`c\overline{c}`$ pairs. This is the same mechanism, which in the strange sector causes the strangeness enhancement. In the charm case, there will be always a “charm enhancement” but, because of screening, the charmed quarks will mostly form open charm. In the present calculation we will let part of the thermally produced charm Coulomb bound states escape. Since the temperature is high, there will be a significant number of parton-parton collisions at a cms energy high enough to produce charm quark pairs, which may form charmonium states. We will now estimate their production rate using the CEM in a thermal environment.

The charmonium production rate in the reaction $`ggc\overline{c}`$, at temperature $`T`$ is given by :

$$\frac{dN^{ggc\overline{c}}}{dtd^3x}=\frac{1}{2}\frac{1}{(2\pi )^6}g_g^2d^3p_1d^3p_2f_g(E_1)f_g(E_2)\widehat{\sigma }_{gg}^{LO}(m^2)v_{12}$$
(22)
where $`g_g`$ is the gluon statistical factor (number of colors $`\times `$ number polarization states) $`v_{12}`$ is the relative velocity between colliding gluons with energies $`E_1`$ and $`E_2`$ and three momenta $`\stackrel{}{p_1}`$ and $`\stackrel{}{p_2}`$, $`\widehat{\sigma }_{gg}^{LO}`$ is the elementary gluon-gluon cross section (8) and $`f_g(E_i)`$ the usual thermal distribution function:

$$f_g(E_i)=\frac{1}{e^{E_i/T}1}$$
(23)
We now introduce the charm pair four momentum $`p`$ with help of the delta function $`\delta ^4[p(p_1+p_2)]`$:

$$\frac{dN^{ggc\overline{c}}}{d^4p}=\frac{1}{2}\frac{1}{(2\pi )^6}g_g^2𝑑td^3xd^3p_1d^3p_2f_g(E_1)f_g(E_2)\widehat{\sigma }_{gg}^{LO}(m^2)v_{12}\delta ^4[p(p_1+p_2)]$$
(24)
In the expression above we have:

$$p_1(E_1,\stackrel{}{p}_{T_1},p_{z_1});p_2(E_2,\stackrel{}{p}_{T_2},p_{z_2});p(E,\stackrel{}{p}_T,p_z)$$
(25)
$$m^2=(p_1+p_2)^2;v_{12}=\frac{1}{2}\frac{m^2}{E_1E_2}$$
(26)
where $`m`$ is the invariant mass of the pair. We next decompose the delta function into temporal, transverse and longitudinal components:

$`{\displaystyle \frac{dN^{ggc\overline{c}}}{dm^2d\stackrel{}{p}_Tdp_z}}`$ $`=`$ $`{\displaystyle \frac{1}{4(2\pi )^6}}g_g^2{\displaystyle 𝑑td^3x\frac{1}{E_1}d^3p_1\frac{1}{E_2}d^3p_2f_g(E_1)f_g(E_2)\widehat{\sigma }_{gg}^{LO}(m^2)m^2}`$ (27)
$`\times `$ $`\delta [m^2(p_1+p_2)^2]\delta ^2[\stackrel{}{p}_T(\stackrel{}{p}_{T_1}+\stackrel{}{p}_{T_2})]\delta [p_z(p_{z_1}+p_{z_2})]`$ (28)
Making use of the identity

$$\frac{1}{2E_i}d^3p_i=d^4p_i\delta (p_i^2)\theta (E_i)=𝑑E_i𝑑\stackrel{}{p}_{T_i}𝑑p_{z_i}\delta [E_i^2p_{t_i}^2p_{z_i}^2]\theta (E_i)$$
(29)
we arrive at:

$`{\displaystyle \frac{dN^{ggc\overline{c}}}{dm^2d\stackrel{}{p}_Tdp_z}}`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^6}}g_g^2{\displaystyle 𝑑td^3x_0^{\mathrm{}}𝑑E_1_{\mathrm{}}^{\mathrm{}}𝑑\stackrel{}{p}_{T_1}_{\mathrm{}}^{\mathrm{}}𝑑p_{z_1}\delta [E_1^2p_{t_1}^2p_{z_1}^2]\theta (E_1)}`$ (30)
$`\times `$ $`{\displaystyle _0^{\mathrm{}}}𝑑E_2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\stackrel{}{p}_{T_2}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p_{z_2}\delta [E_2^2p_{t_2}^2p_{z_2}^2]\theta (E_2)f_g(E_1)f_g(E_2)`$ (31)
$`\times `$ $`\widehat{\sigma }_{gg}^{LO}(m^2)m^2\delta [m^2(p_1+p_2)^2]\delta ^2[\stackrel{}{p}_T(\stackrel{}{p}_{T_1}+\stackrel{}{p}_{T_2})]\delta [p_z(p_{z_1}+p_{z_2})]`$ (32)
Integrating in $`p_{z_2}`$, $`p_{T_1}`$, $`p_{T_2}`$ and $`E_2`$ and defining $`\beta `$ as the angle between $`\stackrel{}{p}_T`$ and $`\stackrel{}{p}_{T_1}`$, we finally obtain:

$`{\displaystyle \frac{dN^{ggc\overline{c}}}{dp_z}}`$ $`=`$ $`{\displaystyle \frac{\pi }{4(2\pi )^6}}g_g^2{\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle _0^{\mathrm{}}}𝑑p_T^2{\displaystyle _0^{2\pi }}𝑑\beta {\displaystyle 𝑑td^3x_0^{\mathrm{}}𝑑E_1_{\mathrm{}}^{\mathrm{}}𝑑p_{z_1}f_g(E_1)f_g(EE_1)}`$ (33)
$`\times `$ $`\widehat{\sigma }_{gg}^{LO}(m^2){\displaystyle \frac{m^2}{E}}\delta [m^22EE_1+2p_T(E_1^2p_{z_1}^2)^{1/2}\mathrm{cos}(\beta )+2p_zp_{z_1}]`$ (34)
$`\times `$ $`\theta (E_1)\theta (EE_1)\theta (E_1^2p_{z_1}^2)`$ (35)
$`=`$ $`{\displaystyle \frac{\pi }{4(2\pi )^6}}g_g^2{\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle _0^{\mathrm{}}}𝑑p_T^2{\displaystyle _0^{2\pi }}𝑑\beta {\displaystyle 𝑑td^3x_0^{\mathrm{}}𝑑E_1f_g(E_1)f_g(EE_1)}`$ (36)
$`\times `$ $`\widehat{\sigma }_{gg}^{LO}(m^2){\displaystyle \frac{m^2}{E}}\left[{\displaystyle \frac{1}{|H(h_1)|}}\theta (E_1^2h_1^2)+{\displaystyle \frac{1}{|H(h_2)|}}\theta (E_1^2h_2^2)\right]\theta (E_1)\theta (EE_1)`$ (37)
where $`E`$, $`H(h_{1,2})`$ and $`h_{1,2}`$ are given by:

$`E=[m^2+p_T^2+p_z^2]^{1/2};H(h_{1,2})=2p_z{\displaystyle \frac{2p_T\mathrm{cos}(\beta )}{(E_1^2h_{1,2}^2)^{1/2}}}h_{1,2}`$ (38)
$`h_{1,2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{p_z^2+p_T^2\mathrm{cos}^2(\beta )}}\{p_z(2EE_1m^2)\{p_T^2\mathrm{cos}^2(\beta )`$ (39)
$`\times `$ $`[4E_1^2(p_T^2\mathrm{cos}^2(\beta )+p_z^2)(m^22EE_1)^2]\}^{1/2}\}`$ (40)
In the first line of (37) we have $`\delta [m^22EE_1+2p_T(E_1^2p_{z_1}^2)^{1/2}\mathrm{cos}(\beta )+2p_zp_{z_1}]`$. It is easy to see that taking either $`E_1p_{z_1}`$ or $`E_1>>p_{z_1}`$ and using the $`\delta `$ function to perform the integral in $`p_{z_1}`$ we obtain a factor $`1/p_z`$. This factor survives and, at the end, gives the $`1/x_F`$ structure (3) mentioned in the introduction. The actual numerical calculation involves no approximation and has a similar behavior. Notice that, because of the kinematical cuts in the invariant mass integral, we are now restricting the result to the rate to charmonium states. A similar expression can be derived for $`c\overline{c}`$ production originating from quark-antiquark annihilation:

$`{\displaystyle \frac{dN^{q\overline{q}c\overline{c}}}{dp_z}}`$ $`=`$ $`{\displaystyle \frac{\pi }{2(2\pi )^6}}g_{q(\overline{q})}^2{\displaystyle \underset{q=u,d,s}{}}{\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle _0^{\mathrm{}}}𝑑p_T^2{\displaystyle _0^{2\pi }}𝑑\beta {\displaystyle 𝑑td^3x_0^{\mathrm{}}𝑑E_1f_q(E_1)f_{\overline{q}}(EE_1)}`$ (41)
$`\times `$ $`\widehat{\sigma }_{q\overline{q}}^{LO}(m^2){\displaystyle \frac{m^2}{E}}\left[{\displaystyle \frac{1}{|H(h_1)|}}\theta (E_1^2h_1^2)+{\displaystyle \frac{1}{|H(h_2)|}}\theta (E_1^2h_2^2)\right]\theta (E_1)\theta (EE_1)`$ (42)
where $`g_{q(\overline{q})}`$ is the statistical factor for quarks (antiquarks) (number of colors $`\times `$ number of spin states), $`\widehat{\sigma }_{q\overline{q}}^{LO}(m^2)`$ is given by (9) and $`f_{q,\overline{q}}(E_i)`$ is the usual quark (antiquark) thermal distribution function:

$$f_{q,\overline{q}}(E_i)=\frac{1}{e^{E_i/T}+1}$$
(43)
In order to account for expansion effects we shall assume that the system cools down following Bjorken hydrodynamics, in which temperature and proper time are related by :

$$\frac{T(\tau )}{T_0}=\left(\frac{\tau _0}{\tau }\right)^{1/3}$$
(44)
where $`T_0`$ is the initial temperature and $`\tau _0`$ is the thermalization time, which marks the beginning of the hydrodynamical expansion. With the help of this relation we can perform the following change of variables:

$$dtd^3x=dx_{}\tau dyd\tau =3dV\tau _0T_0^3T^4dT𝑑td^3x=3_{T_f}^{T_0}\tau _0T_0^3T^4V(T)𝑑T$$
(45)
Finally we introduce the energy scale $`\sqrt{s}`$ (the invariant energy of a single nucleon-nucleon collision):

$$p_z=x_F\frac{\sqrt{s}}{2}$$
(46)
Expressions (37) and (42) can thus be rewritten as:

$`{\displaystyle \frac{dN^{ggc\overline{c}}}{dx_F}}`$ $`=`$ $`{\displaystyle \frac{3\pi }{8(2\pi )^6}}g_g^2\sqrt{s}\kappa {\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle _0^{\mathrm{}}}𝑑p_T^2{\displaystyle _0^{2\pi }}𝑑\beta {\displaystyle _{T_f}^{T_0}}\tau _0T_0^3T^4V(T)𝑑T{\displaystyle _0^{\mathrm{}}}𝑑E_1`$ (47)
$`\times `$ $`f_g(E_1)f_g(EE_1)\widehat{\sigma }_{gg}^{LO}(m^2){\displaystyle \frac{m^2}{E}}\left[{\displaystyle \frac{1}{|H(h_1)|}}\theta (E_1^2h_1^2)+{\displaystyle \frac{1}{|H(h_2)|}}\theta (E_1^2h_2^2)\right]`$ (48)
$`\times `$ $`\theta (E_1)\theta (EE_1)`$ (49)
$`{\displaystyle \frac{dN^{q\overline{q}c\overline{c}}}{dx_F}}`$ $`=`$ $`{\displaystyle \frac{3\pi }{4(2\pi )^6}}g_{q(\overline{q})}^2\sqrt{s}\kappa {\displaystyle \underset{q=u,d,s}{}}{\displaystyle _{4m_c^2}^{\mathrm{}}}𝑑m^2{\displaystyle _0^{\mathrm{}}}𝑑p_T^2{\displaystyle _0^{2\pi }}𝑑\beta {\displaystyle _{T_f}^{T_0}}\tau _0T_0^3T^4V(T)𝑑T{\displaystyle _0^{\mathrm{}}}𝑑E_1`$ (50)
$`\times `$ $`f_q(E_1)f_{\overline{q}}(EE_1)\widehat{\sigma }_{q\overline{q}}^{LO}(m^2){\displaystyle \frac{m^2}{E}}\left[{\displaystyle \frac{1}{|H(h_1)|}}\theta (E_1^2h_1^2)+{\displaystyle \frac{1}{|H(h_2)|}}\theta (E_1^2h_2^2)\right]`$ (51)
$`\times `$ $`\theta (E_1)\theta (EE_1)`$ (52)
The volume of the system evolves in time according to:

$$V(\tau )=V_0\frac{\tau }{\tau _0}V(T)=V_0\left(\frac{T_0}{T}\right)^3$$
(53)
where $`V_0=\pi R_A^2\tau _0`$, $`R_A=(1.18A^{1/3}0.45)`$ fm and the $`\kappa `$ factor was introduced explicitly.

Summing Eqs. (49) and (52) and integrating over $`x_F`$ we obtain the total number of in-plasma produced $`c\overline{c}`$ pairs, which we call $`N_{QGP}^{c\overline{c}}`$.

### B Comparison with other works

We shall now compare the number of $`c\overline{c}`$ obtained in this work with the number of charmed pairs obtained in other works. There are some known papers on the subject . From the reading of these papers, we conclude that there are large discrepancies in the numbers and in the way to obtain them. The sources of these discrepancies are:

a) initial temperature of the plasma, $`T_0`$ (ranging from $`300`$ MeV to $`550`$ MeV)
b) degree of parton equilibration (described by the fugacity factors)
c) initial volume and/or thermalization time, $`V_0`$ and $`\tau _0`$
d) total energy contained in the fireball
e) use or not of a $`\kappa `$ factor (=2) in computing thermal rates
f) use or not of temperature dependent $`\alpha _s`$
g) mass of the charm quark $`m_c`$ (going from $`1.2`$ to $`1.5`$ GeV)
h) use or not of thermal masses for gluons and quarks in the reactions $`g+gc+\overline{c}`$ and $`q+\overline{q}c+\overline{c}`$.

Depending on the choices that one has to make in dealing with a) $``$ h) the final value of in-plasma produced $`c\overline{c}`$ pairs can change by orders of magnitude, going roughly from $`N_{QGP}^{c\overline{c}}=0.02`$ in  to up $`15`$ in . In a comprehensive analysis this same variation was found in . In , the authors arrive to the conclusion that in-plasma production was only a factor two smaller than the initial production. However in that paper the factor $`\kappa =2`$ was not included in calculation. If it were, then $`N_{QGP}^{c\overline{c}}N_{QCD}^{c\overline{c}}`$.

The list of uncertainties given above is not unique and could be enlarged. Each of the items implies taking a decision or making an assumption. Since the list is already large it seems that a high precision calculation is hopeless. In fact the situation is not so bad since our knowledge will increase both because of new experimental results and because of lattice results. RHIC data will impose severe constraints on items from a) to d). For example, taking together data from all the four collaborations, we may expect to know with sufficient accuracy the rapidity distribution of charged particles, from which we can have a good knowledge of the total energy contained in the fireball (item d)). The global analysis of data on rapidity, $`p_T`$ spectra, abundancies, eliptic flow and HBT interferometry, will eventually rule out several initial conditions used in hydrodynamical models and we will have much less uncertainty in the initial temperature of the plasma. On the theoretical side, the perturbative QCD analysis of other processes like, for example, $`J/\psi `$ photoproduction or $`J/\psi `$ production in $`e^{}e^+`$ collisions may significantly reduce the uncertainty in the charm quark mass (item g)). Lattice calculations at finite temperature will hopefully reduce the freedom in the choices of items f) and h). A review of what we may learn in the field in a near future can be found in the summary talk presented by Pisarski in the last Quark Matter meeting .

We have explicitly investigated the effect of changing the mass of the quark $`c`$ and the differences which arise when we use the coupling constant running with $`m^2`$:

$$\alpha _s(m^2)=\frac{12\pi }{(332N_F)\mathrm{ln}\frac{m^2}{\mathrm{\Lambda }_{QCD}^2}}$$
(54)
or with $`T`$ :

$$\alpha _s(T)=\frac{6\pi }{(332N_F)\mathrm{ln}\frac{19T}{\mathrm{\Lambda }_{\overline{MS}}}}$$
(55)
where $`N_F=3`$, $`\mathrm{\Lambda }_{QCD}=230`$ MeV and $`\mathrm{\Lambda }_{\overline{MS}}=80`$ MeV.

In Table II we present our results for the number of “in-plasma born” $`c\overline{c}`$ pairs for different values of couplings and charm quark masses. All the calculations were done with an initial plasma temperature of $`T_0=550`$ MeV, $`\alpha _s(M^2)`$ and $`\alpha _s(T)`$ are given by (54) and (55) respectively. The numbers inside parenthesis correspond to the choice $`\kappa =1`$. Otherwise the numbers were obtained with $`\kappa =2`$.

The fourth line of Table II shows that our results may change by a factor $`30`$ depending on the inputs used. Since there is no strongly preferred value for the $`c`$ quark mass, neither for $`\kappa `$ or for the functional form of $`\alpha _s`$, we are not able to a priori discard any of these choices and our final results will reflect these uncertainties. In some cases a direct comparison with other works is possible. For example, in the second column, comparing the numbers in parenthesis, we notice that we obtain nearly five times more pairs than in . A similar excess is observed comparing the numbers in parenthesis in the fourth column. Comparing our work with Ref.  we can see that the pair production mechanism is quite similar but the treatment given to the plasma expansion is different. Whereas we have used the standard Bjorken hydodynamics, in Ref.  a new hydodynamical model was introduced. Comparing the details of both approaches we concluded that in Ref.  the expansion and cooling of the system is much faster than in Bjorken hydrodynamics. Consequently the system stops much earlier to create $`c\overline{c}`$ pairs and the final yield will be smaller. We think that the expansion of the plasma must be included in any serious calculation, but hydrodynamics is in itself a very complex subject. A prudent strategy, which we adopt here, is to study the desired effect (charm production) first with the most standard and simple hydrodynamical model and then, as a second step, plugg the charm production formalism into a state-of the-art hydrodynamical code. We will leave this last step for the future. The comparison of the other results in Table II with other works shows that for similar inputs we obtain numbers which are compatible with those presented in Refs.  and .

### C The number of $`J/\psi `$’s

The number of $`J/\psi `$’s produced in the plasma can be obtained from Eqs. (49) and (52). We must change the upper limit of integration in $`m^2`$ introducing a kinematical cut-off, i.e., making the replacement $`\mathrm{}4m_D^2`$. In doing so, we rename the superscripts in Eqs. (49) and (52) to $`ggJ/\psi `$ and $`q\overline{q}J/\psi `$. Moreover we introduce the CEM multiplicative factor $`F_{J/\psi }`$ and arrive at:

$$\frac{dN_{QGP}^{AAJ/\psi }}{dx_F}=f_sF_{J/\psi }[\frac{dN^{ggJ/\psi }}{dx_F}+\frac{dN^{q\overline{q}J/\psi }}{dx_F}]$$
(56)
and

$$N_{QGP}^{J/\psi }=_0^1𝑑x_F\frac{dN_{QGP}^{AAJ/\psi }}{dx_F}$$
(57)
The color evaporation model has an ”intrinsic efficiency” given by the fractional factor $`F_{J/\psi }`$ above, which was fixed  in the analysis of $`pp`$ reactions to be $`F_{J/\psi }0.02=2\%`$. It is by no means obvious that the same value should hold for $`AA`$ collisions. For simplicity, we shall assume that $`F_{J/\psi }`$ is universal and holds even for in-plasma production. The “screening factor” $`f_s`$ is thus the only free number introduced here. It gives the probability that a $`J/\psi `$ formed inside the plasma survives the passage through the medium. In other words, $`f_s`$ accounts for dynamical screening, in which gluons in the plasma destroy the charmonium bound state. It varies from $`0`$ to $`1`$. A literal interpretation of  or  would imply $`f_s=0`$. We will however tolerate a small value of $`f_s`$ and examine the consequences. In fact, it is impossible to say what a ”realistic” value of $`f_s`$ would be. Our choices, based on some numerical estimates, are in the following interval:

$$10^3f_s\mathrm{\hspace{0.17em}10}^2$$
(58)
If $`f_s<10^3`$ we do not observe any visible effect of the in-plasma production. On the other hand, $`f_s=10^2`$, which can be interpreted as meaning that 1 % of the in-plasma born $`J/\psi `$’s can survive as bound states, can be taken as an upper limit, beyond which, rather than taking small fluctuations into account we would be really challenging the well established concept of screening.

## IV  Final thermal $`J/\psi `$ production

We will also, for completeness, consider the case where the $`J/\psi `$’s produced in the early stage of the collision traverse the plasma being destroyed and then “regenerated” . We call them “thermal”. As has been discussed in the literature , regenerated $`J/\psi `$’s follow a thermal distribution with the temperature of the quark-hadron phase transition $`T_c170`$ MeV:

$`{\displaystyle \frac{dN_{thermal}^{AAJ/\psi }}{dx_F}}`$ $`=`$ $`{\displaystyle \frac{\pi \sqrt{s}V_0A}{2}}{\displaystyle _0^s}𝑑p_T^2exp\left[{\displaystyle \frac{1}{T_c}}(x_F^2s/4+p_T^2+m_{J/\psi }^2)^{1/2}\right]`$ (59)
$`=`$ $`constexp\left[{\displaystyle \frac{1}{T_c}}\left(\sqrt{x_F^2s/4+m_{J/\psi }^2}\right)\right]\left[1+{\displaystyle \frac{1}{T_c}}\sqrt{x_F^2s/4+m_{J/\psi }^2}\right]`$ (60)
where the constant in front of the integral will be fixed later through the normalization condition:

$$N_{th}^{J/\psi }=_0^1𝑑x_F\frac{dN_{thermal}^{AAJ/\psi }}{dx_F}$$
(61)
## V Final state interactions

In sections II, III and IV we discussed $`J/\psi `$ production. After being produced these states will suffer interaction with the free partons of the plasma and later with the hadronic comovers. In , the final state interactions (FSI) were incorporated (to the rate analogous to our (17)) by the introduction of the multiplicative suppression factors $`S_{FSI}^{COM}`$ and $`S_{FSI}^{QGP}`$, for interactions with comovers and QGP partons respectively. We shall adopt here the same procedure. In fact, in the case of the initial QCD production, since we are using the same formalism, for the sake of comparison we shall later borrow the final corrected expression from . In the case of (57), the $`J/\psi `$ dissociation by partons within the plasma is already taken into account by the “screening factor” $`f_s`$. If the charmonium state survives the plasma phase, it will interact with hadronic comovers and (57) will be multiplied by the suppression factor :

$$S_{FSI}^{COM}=exp\left[\sigma _{co}n_{co}\mathrm{ln}\left(\frac{n_{co}}{n_{fo}}\right)\right]$$
(62)
where $`n_{co}`$ is the comover density, $`n_{fo}`$ is the freeze-out density ($`=1.15fm^2`$) and $`\sigma _{co}14`$ mb is the charmonium-hadron cross section . The above expression depends on the impact parameter, on the collision energy and on time at which the charmonium interacts with a given comover. Moreover, since $`n_{co}`$ is a function of the rapidity, it will depend also on $`x_F`$. However, restricting the analysis to central collisions and to a fixed ($`\sqrt{s}`$) energy, the suppression factor tends to be constant, as concluded in . For simplicity we shall take here the average of the values quoted in these works:

$$S_{FSI}^{COM}=0.25$$
(63)
Final state interactions with comovers will be included in the thermal production differential rate (60) in the same way and with the use of (63).

## VI Feynman momentum distributions

In Fig. 1 we show, for $`AuAu(\sqrt{s}=200`$ GeV) (1a) and $`PbPb(\sqrt{s}=17`$ GeV) (1b) collisions, the results obtained with (17) (dash-dotted line), (56) (solid line) and (60) (dashed line). For the QCD distribution, we have used the nuclear parton distribution functions (PDF) parametrized as in . In doing so, the proton PDF’s were taken from  (GRV98 LO). We also used $`m_c=1.2`$ GeV, $`m_D=1.87`$ GeV and $`\kappa =2`$. With these choices the total number of produced $`J/\psi `$’s is determined. In computing (56) we have used $`\tau _0=0.7`$ fm and $`T_0=550`$ MeV, for $`Au+Au`$ and $`\tau _0=1`$ fm and $`T_0=300`$ MeV for $`Pb+Pb`$ collisions. $`T_c`$ was taken to be $`170`$ MeV.

In each figure all curves are normalized to the same number of produced $`J/\psi `$’s, which is given by (19). In Fig. 1 we have used the parton distributions of , obtaining $`N_{J/\psi }0.09`$ at RHIC and $`N_{J/\psi }0.001`$ at SPS.

Although the normalization is still artificial (it overestimates the number of $`J/\psi `$’s produced in QGP) it anticipates our main claim: QGP production will create a peak at $`x_F0`$. Even if suppressed, it will leave a signal. Moreover, this may happen at RHIC and at SPS as well.

In order to further investigate the differences between the three $`J/\psi `$ production mechanisms we will study the relation between the measured number of $`J/\psi `$’s and the naive expectation based on what we know from $`pp`$ collisions. Experimentally this ratio is easily constructed from the measured $`x_F`$ spectra in $`A+A`$ and $`p+p`$ collsions. Theoretically, the numerator can be written in terms of our theoretical prejudices.

As a starting point, we assume that all $`J/\psi `$’s are produced directly and we show in Fig. 2 the following ratio:

$$R_I(x_F)=\frac{\frac{1}{\pi R_A^2}\frac{d\sigma _{QCD}^{AAJ/\psi }}{dx_F}}{T_{AA}(b=0)\frac{d\sigma ^{ppJ/\psi }}{dx_F}}=\frac{\frac{dN_{direct}^{ABJ/\psi }}{dx_F}}{T_{AA}(b=0)\frac{d\sigma ^{ppJ/\psi }}{dx_F}}$$
(64)
where the denominator is given by the product of (10) with the nuclear overlap function $`T_{AA}(b=0)`$ $`=29.9`$ mb<sup>-1</sup> at RHIC and $`=30.4`$ mb<sup>-1</sup> at SPS. The numerators were obtained with (16), (17) and the nuclear parton distribution functions (PDF) parametrized as in  (solid line) and as in  (dashed line). In the first case the proton PDF were taken from  (GRV98 LO) whereas in the second case they were taken from  (MRST LO). With this last set of parton distributions $`N_{J/\psi }0.06`$ at RHIC and $`N_{J/\psi }0.001`$ at SPS.

$`R_I`$ shows essentially the effect of shadowing in the low $`x`$ region. We can see that if no effect is present, other than shadowing, $`J/\psi `$ spectrum produced in $`AA`$ collisions is just a constant times the corresponding $`pp`$ spectrum. As it can be seen this constant may change a lot with the PDF. This constancy in the region $`x_F<0.1`$ will be important for the subsequent discussion. This figure shows also that the value of the constant decreases when we go from SPS to RHIC collisions. This same behavior was found in .

In Figure 3 we present the $`x_F`$ distribution of the in-plasma produced $`J/\psi `$’s. The screening factor will be fixed to $`f_s=10^2`$. The curves are obtained with expression (56) normalized to the largest (solid line) and to the smallest (dashed line) number of $`J/\psi `$’s. These normalization constants correspond to the choices 2.7 and 0.16, in the last line of Table II. Fig. 3a and 3b correspond to RHIC and SPS collisions respectively.

If we believe that all production is mainly direct plus some small plasma component, then (64) aquires the form:

$$R_{II}(x_F)=\frac{\frac{dN_{direct}^{AAJ/\psi }}{dx_F}+\frac{dN_{QGP}^{AAJ/\psi }}{dx_F}}{T_{AA}(b=0)\frac{d\sigma ^{ppJ/\psi }}{dx_F}}$$
(65)
where, for the moment, we neglect the final state interactions discussed in the previous section. If the screening in the plasma is so strong that $`f_s=0`$, then, $`R_{II}=R_Iconst`$. If, however, there is a small $`J/\psi `$ survival probability, for example $`f_s=10^2`$, there will be a noticeable change in $`R_{II}`$. This is illustrated in Fig. 4, where we see a pronounced deviation from a flat behavior in the low $`x_F`$ region if there is a contribution from the plasma. This interesting feature is due to nature of plasma production, which is peaked at $`x_F0`$. In Fig. 4 the two upper lines correspond to the largest plasma $`J/\psi `$ production ($`N_{QGP}^{J/\psi }=2.7`$) and the lower two lines to the smallest plasma contribution ($`N_{QGP}^{J/\psi }=0.16`$). The scrrening factor is kept fixed to $`f_s=10^2`$. In all cases solid and dashed lines mean the same PDF’s as in Fig. 2.

Figs. 4a and 4b refer to $`Au+Au`$ at RHIC and to $`Pb+Pb`$ at SPS respectively. We can clearly see that deviations from unity happen more strongly for RHIC collisions.

The above situation has to be regarded as an extreme case, in which $`J/\psi `$’s are produced via parton fusion (with shadowing taken into account) in the initial state of the nucleus-nucleus collision and then just escape from the fireball, without further nuclear or comover suppression.

We now consider the opposite extreme case, where all $`J/\psi `$’s are formed in the final state. Final state here means that the original $`J/\psi `$’s are totally destroyed by plasma screening and recriated at the final stage of the fireball life (at $`T=T_c170`$ MeV) by coalescence and then just “emerge ready” during the fase transition, obeying the thermal distribution (60). This scenario is supported by the calculations performed in . Assuming this production mechanism the ratio defined above is modified to:

$$R_{III}(x_F)=\frac{S_{FSI}^{COM}\frac{dN_{Thermal}^{AAJ/\psi }}{dx_F}+S_{FSI}^{COM}\frac{dN_{QGP}^{AAJ/\psi }}{dx_F}}{T_{AA}(b=0)\frac{d\sigma ^{ppJ/\psi }}{dx_F}}$$
(66)
where $`S_{FSI}^{COM}=0.25`$. When discussing thermal production a crucial aspect is how to normalize the distribution (60). The most straightforward procedure would be to normalize it in such a way that the total number of $`c`$ and $`\overline{c}`$ quarks in the final state would match the number of charm quarks pairs produced initially and computed with the help of perturbative QCD. However, as we can see from Table I and also from the recent and comprehensive analysis performed in Ref. , already at the $`p+pc\overline{c}+X`$ level there is some uncertainty in the cross section coming from choices of the quark masses, renormalization scales and parton density parametrizations. A further source of uncertainty arises when we go from $`pp`$ to $`AA`$ collisions and introduce shadowing effects. Moreover, in a very recent calculation of thermal charm production , already taking into account the first PHENIX data , it was shown that, to explain the overall magnitude of the data, one needs to increase the NLO QCD $`c\overline{c}`$ yield by a factor $`2.8`$. In view of this lack of precise knowledge, in our study we will normalize the integral of Eq. (60), $`N_{J/\psi }^{th}`$ in two different ways:

$$I)N_{th}^{J/\psi }=3N_{QCD}^{J/\psi };II)N^{J/\psi }_{th}=0.5N^{J/\psi }_{QCD}$$
(67)
where $`N_{QCD}^{J/\psi }`$ is given by Eq. (19), which, in turn, will be evaluated with two different nuclear parton distributions

In Fig. 5 we plot the ratio $`R_{III}`$. The two lower (upper) pannels refer to $`Pb+Pb`$ ($`Au+Au`$) data measured at SPS (RHIC). In the two left (right) pannels we use the Eskola et al.  (Hirai et al. ) nuclear parton distributions. Inside each of the pannels, the two upper curves are obtained with I) $`N_{th}^{J/\psi }=3N_{QCD}^{J/\psi }`$ and the two lower curves with II) $`N_{th}^{J/\psi }=0.5N_{QCD}^{J/\psi }`$. Finally, solid and dashed lines represent the maximal and the minimal contribution from the QGP component Eq. (56). Although $`f_s=10^2`$ it is easy to see that choosing a smaller value for the screening factor would move the solid lines to the dashed ones, since they contain already a very small number of “in-plasma” born $`J/\psi `$’s.

In Fig. 5 we see that there is already a rise of $`R_{III}`$ at decreasing values of $`x_F`$, but this rise has the shape of a plateau (as anticipated in the introduction), starting from $`x_F0.01`$. If the screening is less effective, $`f_s=10^2`$, then we find a deviation from the plateau at RHIC energies, as in Figs. 5a and 5b, but do not observe any visible effect at SPS energies, as in Figs. 5c and 5d. From these last figures we conclude that in-plasma production can hardly be detected at SPS since the $`x_F`$ distribution is always dominated by the thermal (plateau-like) contribution. On the other hand, it should be remarked that if a plateau is found experimentally, this is interesting in itself, being a strong evidence in favor of the statistical hadronization model, which implies plasma formation.

Now we come back to the prediction made in , namely that the initially produced $`J/\psi `$’s after the interaction either with hadronic comovers or with the plasma will disappear around $`x_F0`$ giving place to a dip in this region. We can split Eq. (65) into two pieces, the first one containing initial production and the second containing the plasma contribution. In order to include the suppression resulting from the final state interactions predicted in , we replace the first piece by a parametrization of $`R_{AB}(x_F,b=0)`$, Eq. (3) of , keeping our plasma contribution the same and construct a new version of (65):

$$R_{IV}(x_F)=R_{AB}(x_F)+\frac{S_{FSI}^{COM}\frac{dN_{QGP}^{AAJ/\psi }}{dx_F}}{T_{AA}(b=0)\frac{d\sigma ^{ppJ/\psi }}{dx_F}}$$
(68)
This quantity, plotted in Fig. 6 is a more realistic version of $`R_{II}`$ (plotted in Fig. 4). As in previous figures solid and dashed lines correspond to the nuclear parton densities of Refs.  and  respectively. Figs. 6a and 6b correspond to RHIC and SPS collisions respectively. Following the same convention employed in Fig. 4, in each pannel the two upper curves represent the largest QGP contribution, whereas the two lower curves represent the smallest QGP contribution. Comparing Figs. 4 and 6 we can see that the inclusion of FSI changes the results, decreasing the strength of the low $`x_F`$ enhancement. Nevertheless, the effect of in-plasma production remains quite visible. Notice that the lowest line crosses the origin of the horizontal axis ($`x_F0`$) at $`0.16`$. This is still far above the point predicted in , of $`0.055`$.

We discuss now two more conservative scenarios, i.e., with less QGP production. First, we consider the lower part of Fig. 6a, amplifying the scale and allowing for smaller values of the screening factor $`f_s`$ ($`f_s=10^3`$). This amplification is shown in Fig. 7a). Inspite of the uncertainties in the calculations, our results in this figure, being still larger than the one found in  (the dash-dotted line), suggest that an enhancement in $`R_{IV}`$ in the low $`x_F`$ region might be seen if QGP would be formed. Finally, we compute again the QGP production rate, using the same inputs except for the initial temperature of the plasma, which we take to be $`T_0=300`$ MeV. At these smaller temperatures the screening factor must be larger and we choose it to be $`f_s=0.1`$. This is shown in Fig. 7b) with the same conventions used in the previous figures. As it can be seen, our lowest line is still above the result of .

## VII Conclusions

In this paper we have considered $`J/\psi `$ production in nucleus-nucleus collisions by three different mechanisms: direct (primordial parton fusion) ii) thermal (statistical coalescence at hadronization) and iii) QGP (in-plasma parton fusion). As already discussed in  i) and ii) are just two extreme cases. In realistic simulations, for a fixed number of $`c\overline{c}`$ pairs and of $`J/\psi `$’s, these last are gradually destroyed, giving origin to the “regenerated” component ii). Component iii) is reconsidered here after being forgotten for long time. Our main point was that, even being small, contribution iii) is very strongly peaked around $`x_F0`$ and can thus become visible if enough plasma is formed.

With this in mind we defined the ratio (64) and computed it assuming mechanism i) plus a small component of iii) (ratio $`R_{II}`$) and assuming mechanism ii) plus a small component of iii) (ratio $`R_{III}`$). In the first case it is completely flat without plasma contribution. A QGP contribution creates a big bump at $`x_F<0.01`$. In the second case, without plasma contribution, we observe a step structure, with a plateau at $`x_F<0.01`$. Switching on the plasma contribution creates a steeply falling curve. We believe that any of these features can be measured and will be very interesting for charm physics at RHIC.

In  it was suggested that the $`x_F`$ distribution to be measued at RHIC is flat and has a dip in the region $`x_F0`$. Our calculations indicate that the inclusion of the “in-plasma born” $`J/\psi `$’s will fill this dip.

Of course, there are several improvements that one could do in the in-plasma production, as, for example, the use of a more sophisticated hydrodynamical expansion. Moreover, final state interactions of this component should be considered. Work along these lines is already in progress. Our findings in this work encourage us to pursue this program.

Acknowledgements: This work was partially supported by FAPESP under contract 00/04422-7. We are grateful to R. Vogt, P. Levai, X.N. Wang and W. Zajc for useful discussions.