R_{AA} of J/\psi near mid-rapidity in heavy ion collisions at \sqrt{s_{NN}}=200 GeV

of near mid-rapidity in heavy ion collisions at GeV


We build up a model to reproduce the experimentally measured of near midrapidty in Au+Au collision at GeV. The model takes into account the suppression from the quark-gluon plasma and hadron gas as well as the nuclear absorption of primordial charmonia and the regeneration effects at the hadronization stage, and hence is a generalization of the two component model introduced by Grandchamp and Rapp. The improvements in this work are twofold; the addition of the initial local temperature profile and a consistent use of QCD NLO formula for both the dissociation cross section in the hadron gas and the thermal decay widths in the quark-gluon plasma for the charmonium states. The initial local temperature profile is determined from the assumption that the local entropy density is proportional to a formula involving the number densities of the number of participants and of the binary collisions that reproduces the multiplicities of charged particles at chemical freeze-out. The initial local temperature profile brings about a kink in the curve due to the initial melting of . The initially formed fireball, composed of weakly interacting quarks and gluons with thermal masses that are extracted from lattice QCD, follows an isentropic expansion with cylindrical symmetry. The fit reproduces well the Au+Au as well as the Cu+Cu data. The same method is applied to predict the expected from the Pb+Pb collision at LHC energy.


I Introduction

Ever since the seminal work by Matsui and Satz (1), has been investigated for a long time as a diagnostic tool to probe the properties of hot nuclear matter created in the early stages of a relativistic heavy ion collision. The original claim stated that cannot be formed in the quark-gluon plasma(QGP) due to color screening between the charm and anti-charm quark and that such mechanism will lead to the suppression of production in a heavy ion collision if the QGP is formed. However, the situation was found to be more involved as recent lattice QCD results showed that might survive past the critical temperature ( and dissolve only at a higher temperature ()(2); (3). If so, it is important to know the detailed temperature dependencies of the properties of the above , as large thermal width for example might still lead to a very small survival rate of the (4); (5); (6). Moreover, if the quark-gluon plasma is formed and the number of charm quarks are large, there is an additional production mechanism that has to be estimated and comes from the formation of through the regeneration of charm and anti-charm that are well described by statistical approaches(7); (8); (9); (10). A successful phenomenological model that is able to explain the signals coming from p+p to heavy ion collisions at high energy should consistently include all the above mentioned ingredients.

The commonly measured and calculated nuclear modification factor shows whether the is suppressed or enhanced in the heavy ion collision (11). For instance, if , the number of produced in A+A collision is equal to the number of produced in p+p collision at the same energy, times the number of binary collisions of nucleons from the two colliding nuclei. If , the number of dissolved in hot nuclear matter is larger (smaller) than that regenerated from the quark-gluon plasma at the hadronization point. Presently of is below unity up to the maximum energy of RHIC. (12).

Phenomenological models have been developed to describe of from SPS to RHIC, and to predict the upcoming LHC heavy ion data. One of the successful models to explain the heavy ion data is the thermal model (13). Here the production follows the statistical description at the chemical freeze out point, and thus can be thought of as being regenerated from the quark-gluon plasma at the hadronization point (). In the so called two-component model by Grandchamp and Rapp(GR) (14), one additionally takes into account the fact that the survives past and dissociates only at a higher temperature so that some of produced at initial stage survives the high temperature phase and reaches the final stage; these together with those regenerated from the quark-gluon plasma at the hadronization point constitutes the two components of the observed (15). In a third model the continuous dissociation and regeneration of is calculated using the Boltzman transport equation from to evolved using the hydrodynamics equation (16).

Crucial quantities in the above models are the ’s for the charmonia; these are the melting temperatures for the Deby screened Coulomb type of charmonium bound states above . While, initial calculations seemed to suggest that for and lower for the excited states(2); (3); (17), recent studies claim that these could be lower due to the large thermal width the states acquire (4); (18). The current experimental error for of as a function of the number of participants are still too large to draw any definite conclusions about these different temperatures. However, one finds interesting characteristics if one takes the central values of the experimental data for the curve of Au+Au collision at GeV; there seems to be two sudden drops (12). One occurs initially at a very small number of participants and the other somewhere between 170 and 210. Assuming that the initial temperatures of the fireball is correlated to the number of participants, the two drops might be related to two distinct initial temperatures. Indeed, as we will show in our two component model, the first sudden drop at small number of participants can be related to the dissociation of excited charmonia such as and ; the initial dissociation of these states will suppress the expected 30 to 40% feedback production of from its excited states. The second sudden drop is found to be related to the initial dissociation of . In geometrical terms, the first drop occurs when there appears regions where the initial temperature is larger than the of the excited charmonia, and the second drop when it is larger than of . Moreover, if the thermal decay width of is not so large, these effects become prominent.

In this work, we generalize the two-component model to take into account the initial hot region mentioned above, and calculate the of at midrapidity in Au+Au collision at GeV. We estimate the initial local temperature with the assumption that local entropy density is simply the linear combination of number density of the number of participants and of binary collision, a formula that well reproduce the multiplicities of charged particles at chemical freeze-out stage in the scenario of isentropically expanding fireball. The average temperature of the fireball is obtained from equating the entropy density to that of the quark-gluon plasma calculated using a non interacting gas of quarks and gluons with effective thermal masses extracted from lattice QCD calculation for the energy density and pressure (19). The baryon chemical potential is obtained from the ratio of the proton and the antiproton, and other chemical potentials from the conservation of respective flavors.

In the two-component model, the primordial charmonia, which have not evolved to the final charmonia yet, first undergo nuclear absorption. Then after they form into the initially produced charmonia, they undergo through thermal decay in the quark-gluon plasma, and then hadronic decay in the hadron gas until the chemical freeze-out stage. On the other hand, the regenerated charmonia only go through the hadronic decay in the hadron gas. Consequently, another important quantity in the model is the thermal widths of charmonia in the quark-gluon plasma and in hadron gas. We improve previous calculations by using the results obtained using perturbative QCD up to the next to leading order (NLO) in the coupling constant(30). In our work, this coupling constant is not a free parameter, because it is related to the screening mass of quark-gluon plasma and sequential melting temperatures of charmonia. The details is mentioned in chapter VII. The relaxation factor of charm quarks in quark-gluon plasma is calculated up to leading order in perturbative QCD; this factor determines the fraction of thermalized charm quarks and used to estimate the regeneration of charm quarks.

The paper is organized as follows. In section II, we discuss nuclear absorption and introduce the necessary concepts. In section III, we introduce the initial local temperature profile. In section IV, we show how the thermodynamic parameters are determined. In section V and VI we respectively discuss the thermal decay and regeneration of charmonia. In the last two sections, we give the results and conclusions. In the appendix, we summarize the perturbative NLO formula for the dissociation cross section of charmonium.

Ii nuclear absorption

Within the Glauber model, heavy ion collision can be described with collisions of the nucleons. The model has two important scales, the number of participants and the number of binary collisions. Participants mean the nucleons in both colliding nuclei that go through inelastic scattering at least once. The rest are called spectators. Usually the amount of bulk matter created from heavy ion collision is proportional to the number of participants. The number of binary collisions counts only primary collisions of two nucleons. Usually hard particle such as the charmonium is produced through primary collision, because it requires large energy transfer; hence, its production number is proportional to the number of binary collisions. The two numbers are calculated in Glauber model as follows (20):


where is the impact parameter of the two colliding heavy nuclei, and A, B are their mass numbers. is the inelastic cross section of two nucleons, which is about at GeV (11). is called the thickness function defined as


and is called the overlap function. Here, is the distribution function of one nucleon in nucleus A(B), for which we use the following Woods-Saxon model:


For the Au nucleus, (21). The normalization factor is determined by the requirement .

The charm and anticharm quark pair produced through primary collision is called the primordial charmonia before they form into on shell states. Such primordial charmonia can still be absorbed by nucleons passing through it. The survival rate at transverse position from the nuclear absorption is


Here, is the absorption cross section of primordial charmonia by a nucleon, and is about at GeV (15). and are the positions where the primordial charmonium is created in beam direction from the center of nucleus A and from the center of nucleus B respectively. The second line in Eq. (4) is the absorption factor by nucleus A, and the third by nucleus B. Mass number A and B in the exponents of the second and third line are both subtracted by 1, because one nucleon of nucleus A and one nucleon of nucleus B take part in producing a charm pair, and cannot take part in absorbing the pair.

Iii initial local temperature

To build up a model for the suppression, we will closely follow the space time picture of nucleus nucleus collision generally accepted from a hydrodynamic simulations(22). After the heavy ions pass through each other, hot nuclear matter is created with small net baryon density between the two receding nuclei. The hot matter is expected to thermalize early; while the exact thermalization time could be controversial, we will just take the typical time accepted in hydrodynamic simulation given as fm/c(23); (24).

The total entropy produced in heavy ion collision is estimated from the multiplicity of produced particles. The multiplicity of charged particles near midrapidity was found to be well described under the assumption that it is proportional to a linear combination of the number of participants and the number of binary collisions. Because both the multiplicity of particles and the entropy are extensive thermal quantities, it is assumed that the entropy per pseudorapidity is also proportional to the same linear combination of number of participants and binary collisions as the particle multiplicity.


In the linear combination for multiplicities of charged particles, is 0.09 at GeV and 0.11 at GeV; the same values are used for entropy (25); (26). The overall factor of 30.3 is determined so that multiplicities of charged particles are well reproduced under the condition that the ratio between entropy and multiplicity is determined from the statistical model at the chemical freeze out point; , where are the total entropy and multiplicity (density) respectively. The upper figure of Fig. 4 shows that Eq. (5), which is the solid line, reproduces the data well, when isentropic expansion of hot nuclear matter is assumed.

Dividing Eq. (5) by the volume, entropy density per pseudorapidity has the form


where , are volume densities of the number of participants and binary collisions respectively. These densities at thermalization time , which is taken to be 0.6 in this work, are given respectively as,


where is the position vector on the transverse plane, and homogeneous densities in beam direction is assumed.

Near midrapidity, the rapidity is approximated to the longitudinal velocity,


and is further assumed to be equal to the pseudorapidity. Supposing that the rapidity interval of interest is , the total entropy within the interval is . Ignoring longitudinal acceleration, longitudinal size of nuclear matter within is at thermalization time . Therefore, entropy density at thermalization time is


where A is the transverse area.

The equation of state is required to extract the temperature from the entropy density. In the quasi-particle picture, the strongly interacting massless particles are replaced by the noninteracting massive particles, whose properties are determined from the condition that they reproduce the thermal quantities such as energy density and pressure of the strongly interacting massless particles. The thermal quantities of strongly interacting quarks and gluons are simulated by lattice QCD. In this work, we will use the parametrization for the effective thermal masses of quark and gluon extracted by Levai and Heinz from lattice QCD data for pure gauge, for and for (19). They fix the degrees of freedom of quarks and gluon to those in the vacuum, and use the forms of LO perturbative QCD for their masses, parameterizing the the nonperturbative effects into a temperature dependent coupling :




, 140, 170 MeV and , 1.03, 1.05 for pure gauge, and cases respectively. Here we choose the case of for thermal masses, because its critical temperature is more reasonable than that of .

With these effective thermal masses and the equation of state obtained by assuming a gas of weakly interacting quarks and gluons, the entropy density is related to temperature as given below,


The entropy density of quark-gluon plasma has a term proportional to the derivative of the squared thermal mass with respect to temperature, because thermal masses depend on temperature. But the term is to be canceled by the term induced from bag pressure to maintain thermodynamic consistency (27). The summation includes gluon and quark of three flavors, although the lattice data is extracted for , we use the realistic number of flavors to get a more realistic magnitude for the entropy. The signs in the denominators are minus for gluon and plus for quarks. Dashed line and dotted line in Fig. 1 show the entropy density as a function of the temperature in quark-gluon plasma respectively for parameters extracted from lattice QCD for and .

Figure 1: Entropy density vs. temperature in the hadron gas (solid line) and in the quark-gluon plasma with parameters extracted from lattice QCD with two flavors (dotted line) and with four flavors (dashed line)

From the entropy density of Eq. (6) and the relation between entropy density and temperature given in Eq. (13), the initial local temperature distribution at is obtained as given in Fig. 2 and 3. y-axis is the direction of impact parameter and x-axis is perpendicular to beam axis and y-axis. Fig. 2 shows temperature profiles along the y-axis at various impact parameters. As we can see, the maximum temperature becomes lower as impact parameter increases. Fig. 3 shows isothermal lines on transverse plane. In head-on collision case, the transverse radius of quark-gluon plasma is about 7.0 fm at thermalization time. This size decreases as the collision becomes peripheral. More importantly, the equithermal line of MeV exists at fm, but does not at fm. If the melting temperature of is 380 MeV, the region where cannot be formed exists at and at fm, but does not at fm, which, as will be seen later, leads to a sudden drop to curve of as a function of the number of participants.

Figure 2: Profiles of temperature along y-axis at various impact parameters

Figure 3: Isothermal lines on transverse plane at with impact parameters,

Iv thermal quantities in expanding hot nuclear matter

The hot nuclear matter created through a heavy ion collision expands in time; assuming isentropic expansion, the entropy density and the temperature decreases along with it. Once the critical temperature is reached, the phase of the matter changes from quark-gluon plasma to the hadron gas. It is known that the phase transition is a crossover at small baryon chemical potential and that there is a critical point where the phase transition changes from a crossover to the first order as baryon chemical potential increases. While the exact location of the point is not certain yet, the phase transition is expected to be still at the rapid crossover region at RHIC and at LHC energies.

The statistical model is very successful in reproducing the observed particle ratios with two free parameters: the temperature and the baryon chemical potential at the chemical freeze-out stage. The other chemical potentials are obtained from the condition that each flavor must be conserved, which can be written as follows:


Here is the volume of the hot nuclear matter of interest, and the number of protons and neutrons in that volume respectively, and the constituents of the matter: the constituents are the quarks and gluons in the quark-gluon plasma, while they are the mesons and baryons in the hadron gas. For hadron gas, all mesons and baryons whose masses are less than 1.5 GeV and 2.0 GeV respectively are included in the sum. These upper limits for the hadron masses are taken such that inclusion of hadrons with masses higher than the upper bounds will not change the result; for example, if the upper limit of meson mass is set at 2.0 GeV, much more mesons are to be included, but their contribution is small. is the number density of particle in the grandcanonical ensemble. , and are baryon number, the third component of isospin and strangeness of particle respectively.

In our work, will represent the volume in the midrapidity region, which is defined by . To determine all the chemical potentials from Eq. (14), we have to know the two numbers and . The two numbers are determined from two constraints: the first is that the ratio should be the same as in the original combined nuclei, and the second is that the observed ratio of proton to antiproton within the rapidity is well reproduced. The second constraint needs some more explanation. Instead of solving Eq. (14) for each , we will try to find the best fit for the experimentally observed ratio as a function of by assuming that appearing in the right hand side of the first equation in Eq. (14) is equal to the total number of participants in the collision divided by a constant number; the constant number is found to be 20.5.

The relation between and the ratio of proton and antiproton is as following. If is very large, the nuclear matter with volume will have much more baryons than antibaryons and thus baryon chemical potential will be large and positive; in this case, the proton to antiproton ratio will be much larger than 1. On the other hand, if is very small, the ratio will be close to 1. As one can see in the lower figure of Fig. 4, the best fit shown as the solid curve well reproduces the experimental data.

Figure 4: multiplicities of charged particles per pseudorapidity divided by a half number of participants (upper) and ratios of proton and antiproton (lower) near midrapidity at GeV

For the volume , the following cylindrical form is used for simplicity:


The factor 2 is multiplied to take into account the forward and the backward expansion in the beam direction. is the longitudinal velocity of nuclear matter, which is approximately equal to the rapidity near midrapidity. is thermalization time, and is the initial radius of quark-gluon plasma on the transverse plane. As can be seen in Fig. 3, the shape of quark-gluon plasma is not a circle but almond-like except for the head-on collision case.

In this work, the almond shape is transformed to a circle with the same area. is the transverse acceleration of quark-gluon plasma, which is set to fm (15). The acceleration is turned off, when the transverse velocity reaches 0.6 c. The initial total entropy of the quark-gluon plasma is calculated from a modified version of Eq. (5), where the and are not the total number of participants and total number of binary collisions but are respectively the number of participants and the number of binary collisions within quark-gluon plasma. Then, the entropy density is obtained by dividing the total entropy of the quark-gluon plasma by the initial volume as given in Eq. (15) at . Finally, the entropy density at a later invariant time is obtained by assuming the expansion to be isentropic and that the volume expands as given in Eq. (15).

With the chemical potentials obtained from Eq. (14), the entropy density determines temperature. The evolution of temperature of the quark-gluon plasma is obtained by equating the expression for the entropy density given in Eq. (13) to the previously determined entropy density at a later time. In contrast to the case for the quark-gluon plasma, it is assumed that the masses of all hadrons do not change at finite temperature. As we discussed before, the order of the phase transition that we probe at RHIC seems to be a strong cross over. However, the equation of states for the quark-gluon plasma and the hadron gas that we use gives a 1st order transition. Nevertheless, we will still use these equation of states because they are simple, physically intuitive and easy to manipulate. Moreover, although at the physical quark masses, the order of phase transition is a rapid cross over, the change occurs at a very small temperature range so that effectively we can approximate the transition with a simple 1st order transition.

The solid line in Fig. 1 shows the correspondence between the entropy density to the temperature in the hadron gas. As can be seen in Fig. 1, the entropy density of the hadron gas is higher than that of quark-gluon plasma at high temperature. These curves cross at around MeV below which the entropy density of hadron gas is lower than that of quark-gluon plasma. But the curve of hadron gas and curves of quark-gluon plasma are very close to each other, which means that the phase transition from quark-gluon plasma to hadron gas is fast. As hot nuclear matter expands, its entropy density and temperature decreases along the curve of quark-gluon plasma in Fig. 1. At critical temperature, the nuclear matter transfers from the curve of quark-gluon plasma to that of hadron gas. If two curves meet at critical temperature and their derivatives with respect to temperature are the same, the phase transition is a crossover. If two curves meet but their derivatives are different, the phase transition is second order. But two curves in Fig. 1 do not meet at critical temperature so the phase transition is first order. However, the phase transition takes short time, because two curves are close at critical temperature. This time interval is the period of mixed phase where two phases coexist. The upper figure of Fig. 5 shows time dependence of the average temperatures of the hot nuclear matter in head-on collision.

The lifetime of quark-gluon plasma phase in head-on Au+Au collision at GeV is 5.63 fm/c. Mixed phase lasts for 1.44 fm/c. Temperature for chemical freeze-out is set at 161 MeV, as motivated by thermal models. After chemical freeze-out, the number of particles of each species is assumed to remain constant. Multiplicities of charged particles and proton to antiproton ratio are all calculated at this temperature. Hadron gas phase lasts for 1.52 fm/c to chemical freeze-out.

Figure 5: The time evolution of average temperature(upper figure) and chemical potentials(lower figure) of hot nuclear matter at fm. The solid line, dotted line, and dashed line in the lower figure are baryon, strangeness and isospin chemical potentials respectively.

Charmonia produced outside hot region where initial temperature is higher than dissociation temperature of charmonia go through thermal decay caused by interactions with surrounding particles, which are the quarks and the gluons in the quark-gluon plasma, and the hadrons in the hadron gas or both in the mixed phase. Once thermal quantities of hot nuclear matter are given as functions of time, one needs to know the properties of charmonia in hot nuclear matter in order to calculate how many charmonia survive the thermal dissociation.

V Thermal decay of charmonia

The thermal decay width of charmonia at temperature can be calculated using the following factorization formula (28):


where repents the quark or the gluon in the quark-gluon plasma, or the baryons and the meson in the hadron gas. is the degeneracy factor of particle , the number density of particle in the grand canonical ensemble, the relative velocity between the charmonium and the particle , and the dissociation cross section of charmonium by particle . It is assumed that thermal width in the mixed phase is a linear combination of contributions from the quark-gluon plasma and from the hadron gas as given in the following form:


where is the fraction of the hadron gas in the mixed phase.

v.1 dissociation cross section

The crucial quantity in Eq. (16) is the dissociation cross section. As emphasized before, one of the improvements in our work over the previous two component model calculations is the consistent application of the NLO perturbative formula to calculate the dissociation of charmonia both in the quark-gluon plasma (29) and in the hadron gas (30). In fact, it is difficult to describe the dissociation of charmonia both in the quark-gluon plasma and in the hadron gas with the same approach. As an example, GR in ref.(14) uses a meson exchange model for estimating the dissociation of charmonia in the hadron gas and a quasifree particle approximation for charm quark inside the in the quark-gluon plasma. Binding energies of charmonia become very small at high temperature so that the charm and anti-charm quarks inside the charmonia indeed can be approximated by a quasi free particles. With this idea, GR approximate the dissociation cross section of charmonia with the elastic cross section of charm or anti-charm quark, once the energy transfer is larger than the binding energy of charmonium.

However, as we show in the appendix, such contribution constitutes only the leading monopole contribution, whose contribution from the charm and anticharm quarks cancel as they have opposite charge. A consistent calculation shows that the dissociation cross section is of the dipole type and only sensitive to the size of the wave function as is expected of a system composed of a quark and an antiquark system with opposite charge.

v.2 Thermal width in the hadron gas

The dissociation cross section of charmonium by hadron can be calculated by the following factorization formula:


where is the distribution function of parton in hadron . is the longitudinal momentum fraction of parton in hadron , which is a value between 0 and 1. is the renormalization scale of the parton distribution function. Suppose that a gluon is emitted from a light quark in a hadron. If transverse momentum of the gluon is smaller than , parton distribution function absorbs this splitting process. But if the transverse momentum is larger than , this splitting process has to be calculated in of Eq. (18). In other word, dissociation cross section of charmonia by hadron , the left hand side of Eq. (18), is composed of nonperturbative part, which is the parton distribution function, and perturbative part, which is the dissociation cross section of charmonia by parton . That is why Eq. (18) is called factorization formula.

The factorization formula also shows how collinear divergence is removed. If a massless gluon is emitted in the same or opposite direction from the original light quark or a massless gluon, then the propagator will have a pole and induce a collinear divergence. In this case, transverse momentum of the emitted gluon is smaller than and thus the divergence can be absorbed by the parton distribution function. The elementary cross section as well as the parton distribution function depends on the scale but the final result, does not depend on the scale in principle. There are several parton distribution functions available. However the minimum scales above which they are defined, are all larger than the binding energies of charmonia, which is the renormalization scale for . For the case of heavy quark scattering, the separation scale can be taken to be the transferred energy or the mass of heavy quark. But in the case of heavy quarkonium, the separation scale is the binding energy of heavy quarkonium.

We use the Bethe-Salpeter amplitude to describe the bound state of quarkonium in which the generated potential becomes Coulombic in the heavy quark limit; in principle, a consistent counting and renormalization is possible in such limit. Although the phenomenological heavy quark potential for strong force is composed of the Coulombic potential part and the linearly rising string potential energy (31), in the heavy quark limit, the distance between quark and antiquark is very small and the effect of the linearly rising string part of the potential is not felt by the quarks. In such a limit, the binding energy of , which is 1S state, can be estimated using the mass difference to its first excited 2S state, which will also be Coulombic. Although the charm quark mass is not so large, we can still approximate the states to be Coulombic and obtained the binding energy of to be 780 (32); (33). Unfortunately, no parton distribution functions are defined in such a low .

For our purpose, we just take the MRSSI parton distribution function (34) of the pion at GeV and use it for all hadrons as the pion is the most abundant particle in the hadron gas stage. The minimum scale of MRSSI is still much larger than the binding energy of . In fact, due to this reason, it was found that combining the parton distribution function at larger value than the scale of the scattering cross section in Eq. 18 we found inconsistent results for the charmonium states where the cross section becomes negative at some incoming energies(30). Therefore, we will first apply the approach to calculate the dissociation cross section for the bottom system, perturbative QCD approach always gives positive cross section(30) and then extrapolate the result to the case by the following dipole approximation,


Here, are the radii of and respectively.

In Coulombic potential, the radius , where is mass of constituent heavy quark of quarkonium, GeV for charm, GeV for bottom, and is the binding energy of quarkonium. The difference in the renormalization scales of the couplings appearing in the and cross section appearing in Eq. (19) is neglected. Two different energy scales are assumed to have the relation , where are masses of constituent charm and bottom quarks.

The upper figure of Fig. 6 shows the dissociation cross section of by the obtained using the method above. The present cross section is still smaller than that estimated using the meson-exchange model (35); (36) or the quark-exchange model (37).

The cross section obtained above is now put to Eq. (16) to get the thermal width of by . For more simplicity, the thermal width in hadron gas is simplified as following:

Figure 6: Upper figure: the dissociation cross section of by . Lower figure: the thermal widths of (solid line), (dotted line) and (dashed line) in the hadron gas and in the quark-gluon plasma obtained in perturbative QCD approach.

where is relative velocity between charmonium and pion in rest frame of charmonium, and is dissociation cross section of charmonium by the pion. A is the ratio of number density integrated in momentum space of all hadrons, whose masses are less than 1.5 GeV for mesons and 2.0 GeV for baryons, and that of pion. The integrated number density of pion is about at temperature of chemical freeze-out, and that of all hadrons is about .

Figure 7: Wavefunctions of (upper), (middle), and (lower) in momentum space, at screening masses, (solid line), (dotted line), and (dashed line), and (dot-dashed line).

In order to calculate the dissociation cross section of charmonia in perturbative QCD approach, we have to know their wavefunctions (30); (29). For in hadron gas, wavefunction of 1S Coulombic bound state is used (30). For excited charmonia in hadron gas, it is assumed that the cross section is proportional to squared radius. According to the screened Cornell potential (31), vacuum screening masses for , and are , and GeV, from which the radii are found to be , and fm respectively. As a result, the dissociation cross section of () is (7.1) times larger than that of in the hadron gas.

v.3 Thermal width in the QGP

The thermal width depends on the charmonium wave function. Here we use the Cornell potential with temperature dependent screening mass. The form of screened Cornell potential is given as


where is and . The screening mass in quark-gluon plasma depends on temperature. In the limit that goes to zero, the screened potential becomes the normal Cornell potential. Fig. 7 shows the wavefunctions of charmonia in momentum space at various screening mass. We can see that the wavefunctions of and are sensitive to the change of screening mass, while that of is not and is more localized at the origin in momentum space.

For the temperature dependence of the screening mass, we simply use the leading order pQCD result, . The coupling constant is in the case of quark-gluon plasma, which is obtained from the best fit for data in our formalism. The details for in quark-gluon plasma will be explained later. In the case of quark-gluon plasma, the dissociation cross section does not have divergence, because the effective thermal masses of partons play the role of cutoff in momentum space; in a sense the quasi particles are the asymptotic states of the QGP that can perturbatively scatter with the charmonium states. As long as the temperature scale is smaller than the separation scale, such description can be thought of as an effective factorization formula, where all the soft physics is included in the quasi particles.

The lower figure of Fig. 6 shows the thermal widths of charmonia both in the hadron gas and in the quark-gluon plasma, calculated using our NLO formula. We see that the thermal widths of charmonia suddenly increase by about an order of magnitude at the critical temperature. There are several reasons for such a sudden increase of the thermal widths. One of them is the difference in kinetic energies of partons. The kinetic energy of a parton has to be larger than the binding energy of the charmonium states for the break up to take place. Such energy are readily available for massive quarks and gluons at finite temperature. However, for pions, the temperatures have to be higher. Moreover, there is a large increase in the degree of freedom as phase transition occurs.

Another important factor for survival rate of is the feed-down contribution from its excited states such as the , . Some is produced directly but some of them is produced through the decay of its excited states. It is not clear yet what fractions of are produced through feed-down. Here it is assumed that 25% and 8% of are fed down from and respectively, which is the world average (38). Then the survival rate of from thermal decay is



The superscript in the lower equation is species of charmonia, and , the upper limit of the integration, is proper time for the chemical freeze-out. is the hot region where initial temperature is over the dissociation temperature of each charmonia.

Vi the regeneration of charmonia

As discussed before, results from hydrodynamics suggest that hot nuclear matter created in relativistic heavy ion collision is thermalized very fast. Moreover particle yields and their ratios from relativistic heavy ion collision are well reproduced using the grand canonical ensemble. If the number of charm and anti-charm quarks are also thermalized, the total number of charm pairs produced from the collision of nucleus A and B would satisfy the following equation (8):


Here, and are the number density of open and hidden charms respectively, and is the volume of hot nuclear matter; all are a function of the proper time .

In general, if creation or annihilation of charm pairs is considered during the expansion of hot nuclear matter, the number of charm quark pairs will be a function of proper time . Moreover, if the lifetime of hot nuclear matter is sufficiently longer than the thermalization time, the number of charm pairs will become thermalized and satisfy Eq. (23); thermalization will proceed through the processes such as . However, estimates based on perturbative QCD suggests that the lifetime of hot nuclear matter is not long enough for the number of charm pairs to be thermlized, because the cross section for creation or annihilation of charm pairs is small with respect to the life time of the fireball.

Therefore, taking into account the off-equilibrium effects for the heavy quarks, the number of charm pairs given in Eq. (23) will be modified as follows,


where is called the fugacity (8). If the number of charms is more (less) then what it would be if the charm quarks are thermalized, fugacity is larger (smaller) than unity. In the central collision of RHIC the initial temperature of created hot nuclear matter is very high and the fugacity at is smaller than 1. As the system cools down, the fugacity increases. This fact can in fact be simply understood qualitatively as follows; assuming isentropic expansion, where is the density of light partons, substituting this into Eq. (24) gives , since , where is the mass of charm quark, and assuming the number of charm quarks do not change during the expansion, we can conclude that the fugacity should increase as the temperature decreases during the expansion. The final value of fugacity at as a function of is shown in the lower figure of Fig. 8.

In Eq. (24) the number density of open charm is multiplied by and hidden charm by , because open charm has one charm or anti-charm quark while hidden charm has both charm and anti-charm quarks. It may be assumed that the number of charm pairs does not change during the expansion of hot nuclear matter, because the creation and annihilation cross section of the pair is very small. Then the left hand side of Eq. (24) dose not change and its initial value is calculable with Glauber model as following:


where is the cross section for pair production in p+p collision, which is per rapidity near midrapidity at GeV in perturbative QCD (39); (13).

Another correction is about system size. So far grandcanonical ensemble is used for number densities of particles. Canonical ensemble, however, is to be used technically, because quantum numbers of the system are to be conserved. Quantum numbers of the system are fixed in canonical ensemble, while they are not in grandcanonical ensemble. The function of probability density with respect to quantum number is a delta function in canonical ensemble, while it is a gaussian type distribution in the grandcanonical ensemble. As the system size is large, the width of function of probability density in grandcanonical ensemble is narrow. In the thermodynamical limit, where the system size is infinite, the function becomes like a delta function. Therefore, grandcanonical ensemble and canonical ensemble are equal in thermodynamical limit. The equality seems to be true in central collisions of heavy ions, because the number of produced charm or anticharm are not small. However it is questionable for peripheral collisions. The quantum number in grandcanonical ensemble can be converted to that in canonical ensemble by simply multiplying modification factors as follows;


where and are modified Bessel functions and their ratio that depends on the number of charm quarks is called the canonical ensemble correction (9).

The upper figure of Fig. 8 is the canonical ensemble corrections at hadronization. In central collisions, the number of produced charm quarks are not small so that the correction is almost unity. However the correction rapidly decreases in peripheral collisions because the number of produced charm quarks becomes small. In other words, the difference between grandcanonical ensemble and canonical ensemble becomes serious.

Fugacity is calculated from Eq. (25) and (26). The lower figure of Fig. 8 shows the fugacity at hadronization. The Fugacity is about at central collision and remains constant until decreases to about 100, and then rapidly increases. The rapid decrease of canonical ensemble correction causes fugacity to rapidly increase at small .

Figure 8: Canonical suppression factor(upper) and fugacity(lower) with respect to number of participants

So far it is assumed that charm quarks are thermalized kinetically. In other words, it is assumed that the kinetic energy of the charm quarks follow thermal distribution. This thermalization is attained through elastic scattering of charm quarks with the surrounding partons. However the elastic, as well as the inelastic, cross sections of the heavy quarks are smaller than those of the light quarks. The formation of charmonium states through recombination are more probable when the charm and anti-charm quarks are randomly distributed in momentum space. Therefore, to estimate the regeneration effect, it is important to know how much charm quarks and anti-charm quarks are thermalized kinetically at the hadronization point.

To estimate this effect, one defines the relaxation factor as follows, (14):


where is the proper time the quark-gluon plasma phase terminates. Relaxation time, , is defined as follows:


where is the number density of parton in the quark-gluon plasma, and is the elastic cross section of charm or anti-charm quark by parton in the charm or anti-charm rest frame. The cross section multiplied by number density of partons is the mean free path of charm or anti-charm quark. Therefore is the average elastic collision time of charm or anti-charm quarks with the light partons. The elastic cross section is calculated to leading order in perturbative QCD. The effective thermal masses of partons extracted from lattice by Levi and Heinz are used (40); (19). The relaxation factor is about at central collision and decreases as the collision becomes peripheral.

The number of regenerated can now be written as


where and are branching ratios from to and from to respectively. Canonical ensemble correction is not multiplied, because charmonia are hidden charms. , and are respectively the survival rates of , and from thermal decay in hadron gas defined as follows:


Here, is the proper time the hadronization is completed.

Vii results

is a massive particle produced through hard collision of two nucleons from each colliding nucleus. If there is no nuclear modification, the number of produced in a heavy ion collision will be proportional to the number of binary collisions as follows,


Here, is the cross section for production in p+p collision, which is recently measured to be per rapidity near midrapidity at GeV (41); (13). is the thickness function defined in Eq. (1).

The nuclear modification factor is defined as follows,


The numerator in the right hand side of the above equation is the actual number of produced in A+A collision, while the denominator is the expected number of production in the same collision. However, nuclear absorption and thermal decay of charmonia suppress the actual number of produced so that these effects supress . On the contrary, regeneration of charmonia enhances the number so that it raises . Considering these effects altogether, the nuclear modification factor is from Eq. (4), (22) and (32)


The first term of the right hand side represents nuclear absorption and thermal decay of , which includes the convolution in the transverse plane that takes into account the initial suppression if the local temperature is larger than the charmonium dissociation temperature, and the second term regeneration of .

Figure 9: The upper figure compares from the experiment (12) to that calculated from Eq. (33). The dashed line is the nuclear absorption and thermal decay, the dotted line the regeneration effect and the solid line their sum. As for the lower figure, the experimental data and the solid line is the same as in the upper figure, but the dashed line shows the result when the initial suppression of in the hot region is not taken into account.

vii.1 results for Au+Au

The upper figure of Fig. 9 shows that Eq. (33), with a best fit coupling , reproduces the experimental data of (12) well. The dashed line is , the dotted line and the solid line their sum. The coupling constant is an important parameter in the result. It determines the thermal widths of charmonia, their dissociation temperatures and the relaxation factor . When the coupling constant increases, the thermal widths of charmonia become wider, the dissociation temperatures of charmonia become lower and the relaxation factor larger, and vice versa. While perturbative QCD formulas are used, it is known that the quark-gluon plasma from up to about is nonperturbative. In this sense, the coupling constant is kind of a free parameter fitted to reproduce the experimental data.

In our approach the value of the coupling also determines the position of a kink on curve. At present, the experimental uncertainties in is still too large to extract any detailed structure as a function in . Nevertheless, taking the central experimental data of seriously, there seems to be two drops. The first one appears at small number of participants, and the second one between and .

There are two important thermal effect on charmonia production: one is the initial suppression which occurs when the local initial temperature of the hot region is higher than the dissociation temperature, the other is thermal decay of charmonia through interaction with surrounding particles. Two effects occur at different temperatures, which are well distinguished in the lower figure of Fig. 6. From the critical temperature to the dissociation temperature, which are respectively identified by the sudden increase and the divergence of the thermal width, the thermal width of charmonium increases monotonically; in the lower figure of Fig. 6, the divergence is represented by the termination point. For the initial suppression effect, the survival rate of charmonium is 0, while for the effect coming from thermal decay, it depends on the temperature and lifetime of the QGP relevant for thermal decay.

It should be noted that the survival rate will change appreciably only when the thermal width becomes large, because the lifetime of quark-gluon plasma in heavy ion collision is not so long. In fact, as can be seen from Fig. 5, the temperature region near is more important, which last almost 2 to 3 fm/c. This means that the thermal width has to be of the order 100 MeV to have an appreciable suppression effect. Therefore, there will be a critical change in the survival rate depending on whether or not the initial temperature is larger than the dissociation temperature; because that condition corresponds to having either an infinite or finite thermal width. This critical change in survival rate will appear as a sudden change of along the axis of number of participants, because the initial temperature rises as the number of participants increases. Therefore, if at some future point, the error bars in the experimental data are reduced and one indeed finds the sudden kink, it would strongly imply that the dissociation takes place at the corresponding point. Moreover, an initial kink could appear when the initial temperature is higher than the dissociation temperatures of excited charmonia, such as and . It should be kept in mind however that when the thermal decay width is in the order of few hundred MeV just above the deconfinement, the sudden kink could also have occurred due to the thermal decay.

Again, taking the central values of seriously, the second kink can be translated to a dissociation temperature of ; this can be obtained from Fig. 2 and the Glauber model given in Eq. (1) that links to the impact parameter. Now, solving Schrodinger equation with the screened Cornell potential, one finds that the dissociation occurs when the screening mass is about . Since the screening mass is simply in our calculation, the coupling constant has to be fixed at to reproduce the screening mass at the dissociation temperature . This coupling relates to which is not so different from the lattice calculation for the asymptotic value of the coupling reached at higher temperature(42) and to the coupling given in Eq. (11).

With the same value of , one finds the dissociation temperatures of and to be and respectively. The initial drop of the solid line in both the upper and lower figures of Fig. 9 is coming from the initial suppression of both and . The dashed line of the lower figure of Fig. 9 corresponds to the case where the initial suppression of is not taken into account. We can see how much is improved by considering the effects of the initial hot region.

vii.2 understanding the results

To gain some more insights into what causes the characteristic features of the solid line in Fig. 9. Let us artificially change some parameters and investigate the changes in the prediction. The results are shown in Fig. 10. The solid line is the best fit value using the same method as described above but obtained with and with a constant thermal width of 10 MeV. The dashed and dot-dashed curves are when the thermal width is changed to 100 MeV and zero respectively. One sees that the slope of the curve critically depends on the thermal width. While the increase in the width suppresses , one should also note that an increase in the width also implies a larger coupling to the medium and hence an increase in the relaxation factor. This will enhance the regeneration effects and thus compensate the suppression.

In the dotted curve, the mass of is artificially changed by -100 MeV. This is to see if the sudden mass shift of recently advocated by Morita and one of us(5) will have any effect on . As one can see, due to the enhancement of the statistical factor in regeneration effect, the curve is pushed up by a nontrivial amount. Therefore, the effects of a sudden decrease in the mass and a sudden increase in the width at the phase transition point, as advocated in Ref.  (6); (43), might just compensate each other in , making it difficult to identify such effects from a measurement of alone.

Figure 10: when the thermal width and mass are artificially changed. The different curves are explained in the text.

Figure 11: The upper figure is of in Cu+Cu collision at GeV. The solid line for fm and the dotted line for fm are overlapped. The lower figure is of in Pb+Pb collision at TeV. The dashed line is after nuclear absorption and thermal decay, the dot-dashed line is from regeneration effect and the solid line is their sum.

vii.3 results for Cu+Cu

The same method is applied to Cu+Cu collision at GeV. The only difference from the previous analysis is the geometry of the colliding nuclei. and for Cu are set at fm and fm respectively in Eq. (3) (21). The upper figure of Fig. 11 shows of in Cu+Cu collision at GeV. Transverse acceleration of hot nuclear matter in Cu+Cu collision will be less than that in Au+Au collision, so that two different transverse accelerations, fm and fm, are used in the calculation for . However there is almost no difference in so that the two curves almost overlap in Fig. 11. The calculation for the Cu+Cu collision slightly overestimates the experimental data at central collisions while it underestimates at peripheral collisions. In Cu+Cu collision at GeV, the maximum temperature is lower than the dissociation temperature of even at the most central collision, hence no kink appears in .

vii.4 results for LHC

We can attempt to predict of in Pb+Pb collision at LHC as well. Because the collision energy is much higher than at RHIC, parameters will have to be changed from those used before. We change only some of them. and for Au are set at fm and fm respectively (21). and are use for the cross sections for production and for pair production per rapidity in p+p collision (13). The constant for produced entropy in Eq. (5) is replaced by (13). The variable in Eq. (5) will be larger than , because increases as collision energy increases. However we use the same value for lack of information. All chemical potentials are set at zero, because midrapidity region will almost be baryon-free at LHC energy. All other parameters including transverse acceleration and its terminal velocity of hot nuclear matter are set at the same values as at RHIC for simplicity and for lack of knowledge.

The lower figure of Fig. 11 shows the prediction for at LHC energy. In contrast to the case at RHIC energy, the survival rate from thermal decay is very low and the regeneration effect is dominant except in the peripheral collision. This is so because the lifetime of quark-gluon plasma is much longer than at RHIC. According to our rough estimation, the lifetime at the most central collision is twice that of the RHIC case. The long lifetime of quark-gluon plasma suppresses the survival rate of charmonia through thermal decay and enhances the relaxation factor of charm and anti-charm quarks. Abundance of charm quarks is another reason for the dominance of the regeneration effect.

The same method seems to be applicable to lower collision energies such as SPS. In this case, however, the baryon chemical potential is too high to make use of effective thermal masses by Levai and Heniz, which are extracted from the lattice QCD data with zero baryon chemical potential.

Viii conclusion

We introduce a generalized two-component model to calculate the of near midrapidity in Au+Au collision at GeV. In the model, suppression is caused by the nuclear absorption of primordial charmonia and by the thermal decay of charmonia in the quark-gluon plasma and in the hadron gas. Enhancement is caused by the regeneration from the QGP. One of our emphases is the use of a consistent perturbative QCD approach, up to next to the leading order, to estimate the thermal decay widths of charmonia both in the quark-gluon plasma and in the hadron gas. While the coupling is considered to be a free parameter to be determined by experiment, it is also related to the screening mass and subsequently dissociation temperatures of charmonia in screened Cornell potential model.

Another emphasis of our approach is the initial temperature profile which leads to initial suppression of charmonia if the local temperature is higher than the dissociation temperature. When the initial temperature profiles, and thus the initial hot regions are considered, the calculated seems to show two sudden drops; one at very small number of participants and the other at around . The first sudden drop of at very small number of participants is caused when the maximum temperature of hot nuclear matter reaches above the dissociation temperatures of excited charmonia and the second sudden drop when it reaches the dissociation temperature of . On the other hand, if the thermal width of is large compared to the inverse time scale of the QGP, the second sudden drop of disappears. Moreover, some effects can cancel other effects. For example, a increased thermal width of will reduce the , while a decreased mass of will enhance it.

The present experimental data from RHIC are still large to discriminate all the different and competing mechanisms discussed here. However, the central values seem to be consistent with the two drops. In this sense, the maximum temperature of SPS seems to be lower than the melting temperature of , because the second drop is not seen in there (44). Similar conclusion could also be made from the data for the Cu+Cu collision at GeV. However refined data for the Au+Au collision at GeV are required to make a definite conclusion. Moreover fitting a single data can not disentangle all the different suppression mechanisms separately. In that sense, the of at LHC energy should provide a very interesting information. As shown in our rough estimate, at LHC is very different from that at RHIC, in that it is dominated by the regenerated from the QGP. Moreover, within the thermal decay width given in the paper, the life time of the plasma seems to be much longer than the inverse of the thermal decay with. Hence, initially formed will be mostly suppressed through thermal decay. Therefore, the upcoming LHC experiment will provide a totally new information about in relativistic heavy ion collision and thus the nature of QGP and its formation.


This work was supported by KOSEF(R01-2006-000-10684-0) and by KRF-2006-C00011. We thank K.Morita for useful discussions.


Figure 12: one diagram of quark-induced dissociation at next to the leading order in perturbative QCD, and its decomposition

The upper figure of Fig. 12 is a diagram for the dissociation by a light quark at next to the leading order in perturbative QCD. The double line represents external line of and the adjacent small circle the Bethe-Salpeter amplitude. The Bethe-Salpeter vertex represents bound state of charm quark and anti-charm quark and has the following form in the heavy quark limit and in rest frame of :


where is the spin index of and its binding energy. is the relative three momentum of and , , and the wavefunction of . is number of color. The charm quark propagator between Bethe-Salpeter amplitude and gluon vertex is reexpressed as


If the binding energy of is small, the momentum of charm quark propagator will be slightly off-shell. By using Eq. (35), the invariant amplitude for Fig. 12 is


where is the invariant amplitude for with spin of incoming charm quark in the approximation of quasifree particle (14). The lower figure of Fig. 12 shows the diagrammatic decomposition of Eq. (36). In heavy quark limit, the condition of energy conservation,


gives the following order counting


because the binding energy is of order (30). Under this order counting, the denominator of charm quark propagator can be approximated as follows,