Heavy quark production at RHIC and LHC within a partonic transport model

Heavy quark production at RHIC and LHC within a partonic transport model

Jan Uphoff uphoff@th.physik.uni-frankfurt.de Institut für Theoretische Physik, Johann Wolfgang Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany    Oliver Fochler Institut für Theoretische Physik, Johann Wolfgang Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany    Zhe Xu Frankfurt Institute for Advanced Studies, Ruth-Moufang-Str. 1, D-60438 Frankfurt am Main, Germany Institut für Theoretische Physik, Johann Wolfgang Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany    Carsten Greiner Institut für Theoretische Physik, Johann Wolfgang Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt am Main, Germany
July 13, 2019
Abstract

The production and space-time evolution of charm and bottom quarks in nucleus-nucleus collisions at RHIC and LHC are investigated with the partonic transport model BAMPS (Boltzmann Approach of MultiParton Scatterings). Heavy quarks, produced in primary hard parton scatterings during nucleon-nucleon collisions, are sampled using the Monte Carlo event generator PYTHIA or the leading order mini-jet model in conjunction with the Glauber model, revealing a strong sensitivity on the parton distribution functions, scales, and heavy quark mass. In a comprehensive study exploring different charm masses, K factors, and possible initial gluon conditions, secondary production and the evolution of heavy quarks are examined within a fully dynamic BAMPS simulation for central heavy ion collisions at RHIC and LHC. Although charm production in the quark-gluon plasma can be neglected at RHIC, it is significant at LHC but very sensitive to the initial conditions and the charm mass. Bottom production in the quark-gluon plasma, however, is negligible both at RHIC and LHC.

pacs:
25.75.-q, 25.75.Bh, 25.75.Cj, 12.38.Mh, 24.10.Lx

I Introduction

Experimental observations from relativistic heavy ion collisions at the BNL Relativistic Heavy Ion Collider (RHIC) indicate the production of a hot and dense partonic medium Adams:2005dq (); Adcox:2004mh (); Arsene:2004fa (); Back:2004je (), commonly referred to as the quark-gluon plasma (QGP). A good agreement between ideal hydrodynamic simulations Kolb:2000sd (); Heinz:2001xi (); Huovinen:2001cy () and measurements Adams:2003am (); Adler:2003kt () of the elliptic flow hints at fast thermalization of this system, which behaves like a nearly perfect fluid, that is, possesses a very small shear viscosity to entropy density ratio Romatschke:2007mq (); Xu:2007jv ().

Heavy quarks, that is, charm and bottom, are a unique probe for this medium. Due to their large mass (, ), a large amount of energy is needed to produce heavy quarks. Such high energy densities are primarily found at the early stage of heavy ion collisions: in hard scatterings of partons in the nucleons of the heavy ions or during the early phase of the QGP. In that energy domain, the running coupling of the strong interaction is small and nearly constant Bethke:2006ac (). Therefore, heavy quark production should be describable within the framework of perturbative QCD (pQCD) Nason:1987xz (); Nason:1989zy (), even for small transverse momenta.

Another theoretically proposed implication of the large mass is the “dead cone effect” Dokshitzer:2001zm (); Zhang:2003wk (), which implies a smaller energy loss via gluon radiation compared to light quarks and delays the thermalization of heavy quarks by a factor of Rapp:2008qc (), resulting in a thermalization time of about the lifetime of the QGP. In contrast, the experimentally observed energy loss Abelev:2006db (); Adare:2006nq () and elliptic flow Adare:2006nq () of heavy flavor is comparable to that of light quarks, the reason for this puzzle being under investigation Armesto:2005mz (); vanHees:2005wb (); Moore:2004tg (); Wicks:2005gt (); Peigne:2008nd (); Gossiaux:2008jv (). Due to the – in principle – unique identification of heavy quarks because of their flavor and roughly known production time during the early stage of the collision, their distributions can reveal information about the interaction history, rendering them as an ideal probe of the medium.

In principle, heavy quarks can be produced at three stages of the collision: during initial hard parton scatterings, in the QGP, or during the hadronic phase. As we will show in this article, most of the heavy quarks are created in hard parton scatterings during initial nucleon-nucleon collisions in the heavy ion collision at RHIC. However, at the CERN Large Hadron Collider (LHC) secondary production during the QGP phase becomes important. Here, one often distinguishes between prethermal and thermal production depending on whether the partonic medium is already thermalized or not. In the present article, we focus on the first two possibilities of heavy quark production and neglect the hadronic phase, which hardly contributes to the heavy quark yield.

An interesting signature in heavy ion collisions is the suppression Matsui:1986dk () or enhancement of heavy quarkonia like (hidden charm) Andronic:2006ky (). In this article, however, we examine only open heavy quark production and postpone the investigation of hidden charm or bottom to future investigations. For the simulation of the heavy quark production in the QGP, we use the partonic transport model called the Boltzmann Approach of MultiParton Scatterings (BAMPS) and study the impact of various initial gluon distributions such as from PYTHIA Sjostrand:2006za (), the mini-jet model Kajantie:1987pd (); Eskola:1988yh (), and the color glass condensate Iancu:2003xm (); Drescher:2006ca (). The initial heavy flavor yield in hard parton scattering is obtained with PYTHIA and compared with that from leading order (LO) pQCD.

Most of the time we will only refer to charm quarks instead of heavy quarks in general. However, the majority of the concepts and findings apply qualitatively also to bottom quarks if one takes the larger mass into account. For LHC calculations, we will explicitly mention bottom production and make predictions on the bottom yield.

This article is organized as follows. First we introduce our model BAMPS and the partonic cross sections for heavy quark production. In Section III we use this model to estimate the chemical equilibration time of a static system filled with gluons and charm quarks and compare it to the analytic solution. The next two sections address charm and bottom production in heavy ion collisions. After investigating the initial heavy quark distributions from primary nucleon-nucleon scatterings in Section IV, we study the heavy quark production during the QGP phase in Section V. Here, the charm quark fugacity is also calculated and related to the estimated chemical equilibration time from Section III. Finally, we conclude with a short summary in Section VI.

Ii Parton cascade BAMPS

For the simulation of the QGP phase we use the partonic transport model BAMPS Xu:2004mz (); Xu:2007aa (), which stands for Boltzmann Approach of MultiParton Scatterings. BAMPS simulates the fully space-time evolution of the QGP at heavy ion collisions by solving the Boltzmann equation,

((1))

dynamically for on-shell partons with a stochastic transport algorithm and pQCD interactions. are the relevant collision integrals, and the one-particle distribution function of species , since light quarks are not included yet. Bottom quarks can also be investigated if one changes and . In addition to the binary collisions , also scatterings for the gluons are possible. That is, the following processes are implemented in BAMPS:

((2))

Within the stochastic method, the probability for a collision of two particles during a time step in a volume element can be obtained from the collision term of the Boltzmann equation Xu:2004mz ():

((3))

where Cleymans:1992je ()

((4))

stands for the relative velocity, and for the binary cross sections. The cross section for as well as the treatment of the interactions are given in Xu:2004mz ().

The most dominant charm production process is gluon fusion with a differential cross section of

((5))

is the averaged matrix element, that is, the averaged (summed) over color and spin of the incoming (outgoing) particles. It can be expressed in terms of the Mandelstam variables , and Levai:1994dx (); Combridge:1978kx ():

((6))

denotes the mass of the charm (or bottom) quarks. If not otherwise specified, we use for charm and for bottom quarks, which are widely adopted in the literature Muller:1992xn (); Levai:1994dx (); Gavai:1994gb (); Levai:1997bi (); Cacciari:2005rk (); Andronic:2006ky (); Rapp:2008qc (); Nason:1999ta (); Frixione:1997ma (); Vogt:2001nh (). For the coupling of the strong interaction we choose to take a constant value of .

After substituting with the relation , the total cross section is obtained by integrating Equation (5), Combridge:1978kx (); Gluck:1977zm (); Babcock:1977fi (); Barger:1981rx (); Matsui:1985eu ()

((7))

with the abbreviation

((8))

The angular distribution of the charm and anti-charm quarks after the collision is sampled using Equation (5).

The back reaction is negligible in heavy ion collisions at RHIC and LHC Andronic:2006ky () due to the small number of produced charm quarks compared to the number of gluons. However, we want to implement the back reaction in a box model in Section III in order to estimate the chemical equilibration time scale. The cross section of this process can be obtained through detailed balance,

((9))

The factor comes into play due to gluons being identical particles. takes the color and spin averaging and summation into account whereas is a kinematical factor.

The process for charm production by light quark and anti-quark annihilation is not yet included in the cascade. However, we use this process to estimate the initial charm quark yield in heavy ion collisions (see Section IV.3). The cross section for that process is given, for instance, in Combridge:1978kx (); Gluck:1977zm (); Babcock:1977fi (); Barger:1981rx (); Matsui:1985eu ().

Iii Charm production in a static medium

In this section, we test our implementation with BAMPS by comparing our numerical results to the analytic solution of a rate equation. To keep the problem simple and find an analytic solution, we consider a box of gluons and charm quarks, in which only charm-anti-charm production through gluon fusion and the back reaction

((10))

are allowed, whereas all other interactions between gluons and charm quarks are forbidden. However, we checked that adding other possible processes like , , , and does not have an impact on our findings. The additional processes ensure the kinetic equilibration and maintains the chemical equilibration of the gluons, but they do not affect the chemical equilibration time scale of the charm quarks.

For the system with only processes from (III), we write down a rate equation for the evolution of the charm density Matsui:1985eu (); Biro:1993qt (); Levai:1994dx (); Zhang:2008zzc ()

((11))

where and denote the rates and ( being the Lorentz factor) the four-velocity of the considered volume element with charm density , which is identical to the charm pair density since charm and anti-charm quarks are always produced in pairs.

The rates are given by Zhang:2008zzc (); Biro:1993qt (); Levai:1994dx ()

((12))
((13))

In Equation (12), stands for the gluon density and the factor is needed to take into account that gluons are identical particles. and are the mean cross sections weighted with the relative velocity defined in Equation (4). In all calculations we use a constant coupling of for the cross sections, which are given in equations ((7)) and ((9)).

iii.1 Analytic solution of the rate equation

For a static box the four-velocity of each volume element is and the rate equation ((11)) simplifies to

((14))

In chemical equilibrium the rates are equal, . From that relation one can obtain the charm density in chemical equilibrium,

((15))

where denotes the constant total particle density.

Initially the box may only contain gluons which are chemically and thermally equilibrated. That is, the initial gluon density for an initial temperature of reads

((16))

(gluons are treated as Boltzmann particles, is the degeneracy factor for gluons), whereas for charm quarks .

Solving the rate equation ((11)) provides the time evolution of the charm density, from which one can read off the chemical equilibration time scale. Taking into account that the total particle number is constant, the charm quark density as a function of time is given by

((17))

where the abbreviations

((18))
((19))

have been introduced.

The solution is implicitly dependent on the temperature via and and is only valid assuming that the temperature stays constant over the whole period of time. As we will see in the next section, due to the mass creation this is, however, not exactly but in a good approximation the case.

iii.2 Comparison between analytic and numerical solution

The numerical solution is obtained with the parton cascade BAMPS (see Section II). In order to compare it to the analytic solution from the previous section only the two processes from III are allowed. Again, the initial charm density is and the gluons are sampled thermally employing according to Equation (16) with temperatures of 400 MeV and 800 MeV, which correspond to the expected temperatures of the quark-gluon plasma at RHIC Thews:2000rj (); Soff:2000eh (); Fries:2002kt (); Turbide:2005fk (); Rapp:2000pe (); Cooper:2002td (); Wang:1996yf () and LHC Fries:2002kt (); Cooper:2002td (); Wang:1996yf (), respectively (cf. also Figures 19 and 21).

iii.2.1 Charm production at RHIC temperature

The initial temperature of the gluons is taken as Thews:2000rj (); Soff:2000eh (); Fries:2002kt (); Turbide:2005fk (); Rapp:2000pe (); Cooper:2002td (); Wang:1996yf (). Being precisely, during the evolution of the box the temperature is ill-defined because of the gluons being not in thermal equilibrium (due to the neglect of elastic scatterings among the gluons). However, one can extend the temperature definition to this non-equilibrium domain, introducing an “effective gluon temperature”

((20))

Here () stands for the total gluon energy (number) in the box, and and for the respective densities.

The time evolution of the gluon temperature is illustrated in Figure 1.

Figure 1: (Color online) Time evolution of the gluon temperature in a static system with an initial temperature of . The temperature decreases due to the creation of massive charm quarks.

The temperature declines to about since part of the kinetic energy is used for the production of massive charm quarks. Consequently, this implies a problem in comparing the numerical results with the analytic solution, which is only valid for a constant temperature. Therefore, we will compare the analytic solutions for initial and final temperatures against the numerical results at the beginning and at the end, respectively.

Figure 2 shows the time evolution of the charm pair density .

Figure 2: (Color online) Time evolution of the charm pair density in a static medium with an initial temperature of . For comparison, the analytic solutions from ((17)) for initial and final temperature are also shown.

In addition, the analytic solutions of the rate equation for the charm density from Equation (17) for initial and final temperature are also plotted. In the beginning, the numerical results are in very good agreement with the analytic solution for . Thereafter, the temperature of the gluon plasma decreases (cf. Figure 1), and the numerical results must be compared with the analytic solution for the stationary final temperature of , both being in good agreement.

The rates for and the back reaction are shown in Figure 3.

Figure 3: (Color online) Time evolution of the rates in a static system with an initial temperature of . Again, the analytic solutions for initial and final temperatures are also plotted.

As one would expect from the principle of detailed balance, the rates converge to a mutual constant equilibrium value with time. The chemical equilibration time scale may be defined as . By fitting the curve in Figure 2 with

((21))

we obtained an approximate value of for . Of course, this result of the chemical equilibration time was estimated within a simple static model without expansion. However, the order of this result should be roughly the same as in realistic heavy ion collisions, which indicates that not very many charm quarks are produced during the QGP phase at RHIC considering this large time scale of chemical equilibration compared to the lifetime of the QGP, as we will explicitly show in Section V.1.

An obvious expression for the temperature dependence of can be obtained by differentiating Equation (21) by , equating with Equation (14), and evaluating at :

((22))

It is related to from Equation (19), which does not have a direct physical meaning, by

((23))

The temperature dependence of can be obtained analytically with the assumption of a constant temperature throughout the evolution of the system. In order to distinguish it from the obtained numerically with BAMPS, in which the temperature drop due to the mass creation is considered, we label the analytic solution as . Figure 4 depicts the temperature dependence of .

Figure 4: Temperature dependence of the chemical equilibration time for charm production in a static medium of a constant temperature.

Compared to the results from BAMPS, these values for a constant temperature are slightly larger. For high temperatures, which can be present at LHC at a very early stage of the QGP phase (cf. Figure 21), the chemical equilibration time lies in the same range as the lifetime of the QGP. Consequently, a substantial charm production in the QGP at LHC can be expected, if the initial temperature of the medium is large. We will address this in more detail in Section V.2.

The fugacity, which is defined by

((24))

for particle species , is an interesting variable for investigating the chemical equilibration. Figure 5 shows the fugacities of gluons and charm quarks in our simulation as a function of time.

Figure 5: (Color online) Evolution of gluon and charm quark fugacities in a static system. The final value is above 1 due to the fixed number of particles and the temperature drop.

Due to detailed balance, they adopt the same value in equilibrium, although that differs from 1, a consequence of the total particle number being constant. As a note, explicitly allowing inelastic processes () as in the full transport simulation (see Section V) ensures chemical equilibration and leads to a final fugacity of 1 both for gluons and charm quarks.

iii.2.2 Charm production at LHC temperature

For the initial gluon temperature, a value of Fries:2002kt (); Cooper:2002td (); Wang:1996yf () (cf. also Figure 21) is chosen, which drops to about as shown in Figure 6 due to the mass creation.

Figure 6: (Color online) Evolution of the gluon temperature in a static medium with an initial temperature of .

Figure 7 depicts the numerical results of the charm pair density evolution, which is in excellent agreement with the analytic solutions taking the temperature drop into account.

Figure 7: (Color online) Charm quark density as a function of time in a static medium with an initial temperature of . The analytic solutions for initial and final temperature are also shown.

From that figure, using Equation (21), the time scale of chemical equilibration can be estimated to about , which is much smaller than for RHIC temperature, but still sizable.

The rates for both considered processes can be found in Figure 8.

Figure 8: (Color online) Evolution of the rates in a static medium with .

Iv Initial parton distribution in heavy ion collisions

The initial distributions of the partons play a crucial role for the dynamics of the heavy ion collision. In this section, we want to outline a prescription to describe the heavy ion collision – according to the Glauber model – as a superposition of nucleon-nucleon collisions, which are sampled with the event generator PYTHIA Sjostrand:2006za (). In that framework, we also study charm production and compare the yields with experimental data from nucleon-nucleon collisions. To be able to examine the impact of the initial conditions on charm production during the QGP phase in Section V, we also discuss other models for the initial parton distributions such as the mini-jet model and the color glass condensate.

The prescription for the position sampling of the partons according to a geometric model is described in great detail in Xu:2004mz ().

iv.1 Parton and momenta sampling with PYTHIA

For our simulation we use PYTHIA 6.4 Sjostrand:2006za (), allow hard and soft QCD interactions, and turn off the hadronization. PYTHIA distinguishes between soft and hard events. Therefore, we define particles stemming from a hard (soft) event as being hard (soft), regardless of their momenta. The only exception to this rule is that beam remnants such as diquarks are always considered as soft.

For nucleon-nucleon collisions at RHIC with a center of mass energy of , PYTHIA yields the following results. On average, 53 % of all processes are hard and 47 % soft. All particles produced in hard processes are partons. These hard partons account for about 58 % of all particles created in nucleon-nucleon collisions. On the contrary, in soft processes, most formed particles are non-partonic, for instance, diquarks or excited nucleons. Their fraction is about 28 % of all particles, whereas partons which are produced in soft events only account for 14 % of all particles. The total number of produced particles in one nucleon-nucleon collision averages to 7.5, whereof 5.4 are partons. The total energy is, of course, , but all partons together possess only . Consequently, non-partonic particles account for 85 % of the total energy. The reason for that is their large energy per particle ratio, which results from the non-partonic particles either being beam remnants or stemming from elastic or diffractive hadronic processes and, therefore, carrying a huge amount of energy.

The reason for the separation in soft and hard particles lies in the different scaling behavior from nucleon-nucleon to heavy ion collisions. According to the Glauber model, hard processes scale with the number of binary collisions, wong ()

((25))

where is the impact parameter, the nuclear overlap function for collisions of two nuclei and , and for RHIC and for LHC the inelastic p+p cross sections. For central collisions, we approximate Sarcevic:1994ma (); Muller:1992xn (), which leads to for Au+Au collisions at RHIC being in excellent agreement with a numerical calculation using misko_overlap (). Therefore,

((26))

is the number of binary collisions at RHIC for and . Taking shadowing into account reduces this number to about Adams:2004cb (), which we use as a scaling factor for hard partons, . To obtain the scaling factor for soft particles, one can make use of energy conservation,

((27))

where stands for the total energy available in the heavy ion collision (at RHIC ) and for the total energy of all hard (soft) particles in one nucleon-nucleon collision. Solving this equation for leads to , which is in the same order as the number of nucleons, the usually used scaling factor for soft particles.

Employing the scaling prescription obtained from Equation (27) results in the following particle yields for Au+Au collisions at RHIC with : Partons from hard processes account on average for about 93 % of the produced particles. 4.5 % of all particles are non-partonic particles from soft processes, whereas 2.5 % are soft partons. The total number of produced particles is about 4600, of which 4400 are partons. However, for the energy distribution the picture looks a bit different. Hard partons carry just 53 % of the total energy and non-partonic particles with 44 % nearly the same amount. Soft partons only account for 3 % of the total energy. Consequently, the whole energy deposited by partons, which is available for the parton cascade, is about 56 % of the total energy.

Figure 9 depicts the rapidity distribution of the particle number and of their transverse energy in a Au+Au collisions at RHIC.

(a) (b)
Figure 9: (Color online) Rapidity distribution of particle number (a) and transverse energy (b) in central Au+Au collisions at RHIC. As explained in the text, the distributions are obtained by simulating nucleon-nucleon collisions with PYTHIA and scaling them to heavy ion collisions. CTEQ6l is used for the parton distribution functions. “others” denotes non-partonic particles such as diquarks or protons from soft processes. The distribution of soft gluons is almost zero and is, therefore, not shown here or in the following figures.

Gluons from hard processes dominate the spectrum at mid-rapidity in both distributions. In addition, hard quarks take a considerable fraction of the transverse energy. Partons from soft processes have lost much of their significance compared to unscaled nucleon-nucleon collisions due to the smaller scaling factor of soft particles.

In Figure 10, the transverse momentum spectra are shown.

Figure 10: (Color online) Transverse momentum spectra of particles at rapidity in central Au+Au collisions at RHIC.

Of course, hard partons also dominate these spectra. The transverse momentum of soft quarks is always smaller than 2 GeV because of a momentum cut-off in PYTHIA. For partons from hard events, we did not introduce a cut-off. However, PYTHIA has an internal cut-off for so called semi-hard scatterings, which we also consider as hard events. The cut-off is dependent on and lies for RHIC energy at 141 MeV causing a jump at this value in the spectra of hard partons.

For LHC, the initial parton distribution is also sampled with PYTHIA and scaled by using the same prescription. For central Pb+Pb collisions at , the overlap function is and proton-proton cross section , which determines the number of binary collisions to be about 2000. Actually, we reduce this value to to take also shadowing into account Emel'yanov:1999bn (); Accardi:2004be (). Figure 11 depicts the rapidity distributions of the particle number and their transverse energy for heavy ion collisions at LHC.

(a) (b)
Figure 11: (Color online) As in Figure 9, but for central Pb+Pb collisions at LHC ().

Currently, BAMPS does not include light quarks. Therefore, to conserve energy and particle number, light quarks from PYTHIA are converted to gluons before we evolve the system with the parton cascade. This conversion does not account for effects such as the faster thermalization of gluons compared to light quarks and the larger cross section of than of . However, since gluons dominate the early phase of the QGP, where most of the heavy quark production takes place, we do not expect an impact on the production of heavy quarks. Nevertheless, this will be checked in the near future after light quarks are included in BAMPS.

iv.2 Other models for the initial parton distribution

Two other models for the initial parton distribution were already used in the parton cascade BAMPS: the mini-jet model Xu:2004mz (); Xu:2007aa () and a color glass condensate (CGC) inspired model El:2007vg (); Xu:2008zi ().

In the mini-jet model Kajantie:1987pd (); Eskola:1988yh (), the initial parton distribution is given by several independent 2-jet events, which are sampled according to Wang:1991hta ()

((28))

where denotes the transverse momentum, the rapidity and the Bjorken variable. The cross section is calculated in LO pQCD and a factor is introduced to account for higher orders.

To avoid problems at low momenta where pQCD is not valid anymore, we employ a momentum cut-off, which is set to for RHIC in order to get the final transverse energy distribution of the gluons in agreement with data Xu:2005wv (). For LHC we choose to get a gluon yield, which is comparable to the PYTHIA result.

The number of produced partons is given by

((29))

To generate the CGC initial conditions we use the model from Drescher:2006pi (); Drescher:2006ca (). In that model gluons are sampled using the factorization ansatz Gribov:1984tu ()

((30))

where stands for the number of colors, for the light cone momentum fractions and for the center of mass energy.

According to the KLN (Kharzeev-Levin-Nardi) ansatz Kharzeev:2000ph (); Kharzeev:2002ei (), the unintegrated gluon distribution function is related to the saturation scale via

((31))

The normalization of Equation (31) and therefore also ((30)) is determined by the experimentally measured particle multiplicity at mid-rapidity in central collisions at RHIC. For LHC it is normalized in such a way that it is in agreement with the yield from PYTHIA.

Figure 12 compares the initial particle and energy distributions of all partons for the three used models as a function of rapidity.

(a) (b)
Figure 12: (Color online) Rapidity distribution of parton number (a) and their transverse energy (b) for PYTHIA, mini-jet and CGC initial conditions in a central Au+Au collision at RHIC.

Although the particle number at mid-rapidity in all three models is nearly the same, the energy distributions differ considerably. The small value from PYTHIA is a consequence of omitting beam remnants like diquarks, since BAMPS works on a partonic level. A future project will be to include string fragmentation of diquarks and partons in order to transport part of their energy into the partonic phase.

Figure 13 shows the same distributions for LHC.

(a) (b)
Figure 13: (Color online) The same as in Figure 12 for central Pb+Pb collisions at LHC.

The model for the CGC only samples gluons with rapidity smaller than 6, which is, of course, not a physical effect. However, since we are interested in the mid-rapidity region that is not a problem. PYTHIA yields also for LHC the smallest initial energy distribution. The particle number distributions are comparable with various other models Armesto:2008fj (). However, there are huge uncertainties due to shadowing of  % Emel'yanov:1999bn (); Accardi:2004be (), which also affects secondary charm production in the QGP due to the unclear initial properties of the medium. On the other hand primary charm production in initial hard scatterings is influenced by shadowing as well (cf. next section).

iv.3 Charm production in hard processes

The most general expression for the double differential cross section for charm pair production at a collision of the hadrons and is given by the partonic cross section for the process convoluted with the parton distribution functions (PDF) in the hadrons Gavai:1994gb ():

((32))

The indices and denote the partons from and . and are the factorization and renormalization scale, respectively, which are mostly chosen to be equal, , where Gavai:1994gb (); Smith:1996sb (); Eskola:2003fk () or Cacciari:2005rk (); Vogt:2007aw (); Eskola:2003fk () are common choices.

iv.3.1 Leading order charm production within the mini-jet model

In leading order charm quarks are produced in the two processes

((33))

Experimentally, due to the confinement one cannot measure single charm quarks or but only mesons with charm content , or . The invariant cross section for an open charm production process with two charmed mesons in the final state reads Gavai:1994gb ()

((34))

with

((35))

Here is the center of mass energy of the partons, which is related to the center of mass energy of the hadrons via . The partonic cross sections are given in Section II.

The fragmentation function describes the hadronization of the charm quarks, with . Since fragmentation only affects the momentum distribution but not the total cross section Gavai:1994gb (), in which we are interested, the momentum of is assumed to be equal to the charm quark momentum, that is, and . Thus, taking energy conservation into account, Equation (34) can be written as

((36))

being the rapidity. In the center-of-mass frame, the charm and anti-charm quarks hold the same transverse momentum .

Integration of the previous equation leads to the total cross section

((37))

The factor results from the method of charm number counting, which is widely adopted in the literature. The charm quark number corresponding to the cross section should not be the total number of particles (charm plus anti-charm quarks) but the number of charm pairs. The number of produced charm pairs and the cross section is simply related via Vogt:2001nh (); Andronic:2006ky ()

((38))

The limits of the integrals over the rapidities in Equation (37) can be obtained from the definition of the rapidity and some kinematic considerations,

((39))

with .

From equations ((34)) and ((35)) it is obvious that the charm production cross section is fundamentally dependent on the PDFs. Therefore, using different PDFs provided by the Les Houches Accord Parton Density Function (LHAPDF) group Whalley:2005nh () we explore this dependency in detail, as listed in Table 1.

PDF Scale
CTEQ6m 1.2 160 38
1.5 72 19
1.2 140 36
1.5 79 20
PYTHIA 540 130
CTEQ6l 1.2 230 57
1.5 90 25
1.2 280 68
1.5 120 31
PYTHIA 370 91
MRST2007lomod 1.2 290 65
1.5 120 30
1.2 320 66
1.5 150 35
PYTHIA 370 89
GRV98lo 1.2 190 38
1.5 78 17
1.2 220 43
1.5 97 20
PYTHIA 120 30
PHENIX
STAR
Table 1: LO pQCD cross sections for charm pair production in nucleon-nucleon collisions at with a running , , Bethke:2006ac () and a factor of for various parton distribution functions, renormalization () and factorization () scales and charm masses . Results from PYTHIA and experimental data :2008asa_PHENIX (); Adare:2006hc_PHENIX_dsigmadY (); Adams:2004fc_STAR_dcsdY_cstot () are also listed. In PYTHIA is used.

For that study we also vary the charm mass and the factorization and renormalization scale. The calculations are done with a running coupling . In contrast, a fixed coupling of yields much smaller results. Additionally, results from PYTHIA (see next section) and experimental data :2008asa_PHENIX (); Adare:2006hc_PHENIX_dsigmadY (); Adams:2004fc_STAR_dcsdY_cstot () are shown in the table. The reason for the strong deviation between both experiments is currently under intensive investigation.

Our findings are in good agreement with the literature. In Levai:1994dx (), the charm production cross section in LO is calculated to using the HIJING model and taking shadowing into account. Vogt:2001nh () finds a value between 133 and depending on the choice of the parton distribution functions. In Pop:2009sd (), the cross section lies between 300 and , according to whether shadowing and/or strong color electric fields are incorporated.

Our results in Table 1 as well as those from the previously mentioned LO models are strongly dependent on various parameters. Mainly, the dependence on the scales shows that LO calculations are not sufficient to describe these processes. Indeed, compared to experimental data, the LO results are much too small. In contrast, next-to-leading order (NLO) calculations are closer to the data Cacciari:2005rk (). However, errors due to uncertainties in PDFs, scales and charm mass are still significantly high.

Table 2 lists our predictions for the charm production cross section of nucleon-nucleon collisions at LHC energy. Again, there is a considerable dependency on the aforementioned parameters.

PDF Scale
CTEQ6m 1.2 1600 170
1.5 1100 130
1.2 770 83
1.5 690 78
PYTHIA 2600 300
CTEQ6l 1.2 5300 640
1.5 3200 400
1.2 3500 420
1.5 2400 310
PYTHIA 2500 310
MRST2007lomod 1.2 6600 730
1.5 3900 460
1.2 4700 480
1.5 3100 350
PYTHIA 2700 320
GRV98lo 1.2 7600 890
1.5 4100 500
1.2 - -
1.5 - -
PYTHIA - -
Table 2: As in Table 1, but at LHC energy of . The results for GRV98lo with a factorization scale of and for PYTHIA cannot be calculated, because these PDFs are not designed for such large scales.

Levai:1994dx () calculates the cross section for LHC to and Pop:2009sd () to . In Vogt:2001nh () it varies between 2000 and , depending on the PDFs chosen.

iv.3.2 Charm production with PYTHIA

We also used PYTHIA to estimate the number of produced charm quarks in a nucleon-nucleon collision. The corresponding cross sections are listed in Tables 1 and 2. Applying the scaling prescription introduced in Section IV.1, the number of initial charm quarks in a heavy ion collision can be computed from this value. Since they are produced in hard processes, their number scales with the number of binary collisions.111Actually, the scaling factor is a bit smaller if one takes shadowing into account as is done in this calculation (cf. Section IV.1). Table 3 lists the number of initially produced charm quarks in Au+Au collisions at RHIC according to PYTHIA for various parton distribution functions.

PDF Reference Number of charm pairs
CTEQ5l (LO) Lai:1999wy ()
CTEQ6l (LO) Pumplin:2002vw_CTEQ6 ()
CTEQ6m () Pumplin:2002vw_CTEQ6 ()
MRST2001LO Martin:2002dr ()
MRST2007LOmod Sherstnev:2007nd ()
HERAPDF01 Nagano:2008ip ()
GJR08 (FF LO) Gluck:2007ck (); Gluck:2008gs ()
GRV98 (LO) Gluck:1998xa ()
Table 3: Number of charm pairs produced in primary hard scatterings in central Au+Au collisions at RHIC for some parton distribution functions by sampling nucleon-nucleon collisions with PYTHIA and scaling to Au+Au collisions (cf. Section IV.1).

Again, the uncertainties of using different PDFs are reflected. The mean value lies at about 9 pairs.

According to Duraes:2004zt ()  charm quarks are produced in initial hard collisions at RHIC in LO, dependent on the PDFs, charm mass and shadowing. Muller:1992xn () estimates this number to 2 using the mini-jet model and a phenomenological factor of . Gavin:1996bx () extrapolates from p+p collisions in NLO calculations 8.7 produced pairs. The authors also calculate the number of charm quarks at mid-rapidity to about 3 pairs, which is in good agreement with our results of about 2 for CTEQ6l and 3 for CTEQ6m. In another NLO calculation Vogt:2002vx () the total number of charm quarks varies between 8 and 13. Thus, recent results from the literature are in good agreement with our results from PYTHIA.

Experimental data for the differential charm production cross section at mid-rapidity in a nucleon-nucleon collision at is available both from STAR and PHENIX. The former one measured a value of Adams:2004fc_STAR_dcsdY_cstot () and the latter of Adare:2006hc_PHENIX_dsigmadY (), again showing big deviations. The rapidity distribution of the charm production cross section simulated with PYTHIA is depicted in Figure 14 together with the experimental data points and the pQCD calculations from Section IV.3.1.

Figure 14: (Color online) Charm production cross section as a function of rapidity in a nucleon-nucleon collision at RHIC energy simulated with PYTHIA and pQCD, respectively, for the PDFs CTEQ6l and CTEQ6m together with experimental data Adams:2004fc_STAR_dcsdY_cstot (); Adare:2006hc_PHENIX_dsigmadY (). The pQCD calculation is done in LO with , , , Bethke:2006ac () and .

The charm distribution from the LO pQCD calculation lies far below the experimental data. Hence, higher order corrections should be taken into account or the introduction of a phenomenological factor of is necessary. The results from PYTHIA agree well with the PHENIX data point, which is a bit peculiar, since PYTHIA is based on LO pQCD cross sections. However, PYTHIA is tuned with a running coupling and factors in order to describe experimental data well.

Although CTEQ6m reproduces the data better than CTEQ6l, we will use in the following the latter PDF set, because it is designed for LO event generators such as PYTHIA Pumplin:2002vw_CTEQ6 (). Therefore, we can use PYTHIA also for a reliable sampling of the initial gluon distribution.

Heavy quark pairs
PDF Reference Charm Bottom
CTEQ6l (LO) Pumplin:2002vw_CTEQ6 () 62 7.2
CTEQ6m () Pumplin:2002vw_CTEQ6 () 66 6.9
MRST2007LOmod Sherstnev:2007nd () 67 8.9
Table 4: As in Table 3, but for central Pb+Pb collisions at LHC.

Table 4 lists the number of charm and bottom quarks produced during initial nucleon scatterings in Pb+Pb collisions at LHC according to PYTHIA. The results from different PDFs deviate not as much as for collisions at RHIC. Nevertheless, there are still big systematic errors due to uncertainties in mass, shadowing, factorization and renormalization scale, although these are not reflected in the table.

Muller:1992xn () predicts 34 produced charm pairs in Pb+Pb collisions at LHC within the mini-jet model. The NLO value from Gavin:1996bx () was estimated to 450, but reduced to after taking newer PDFs and shadowing into account Vogt:2002vx (). Vogt:2002vx () also makes a prediction for bottom quarks, of which about 5 pairs should be produced at LHC.

The charm production cross sections as a function of rapidity for nucleon-nucleon collision at LHC energy simulated with PYTHIA and LO pQCD, respectively, with different PDFs are shown in Figure 15.

Figure 15: (Color online) As in Figure 14, but for LHC.

It is surprising that the curves from PYTHIA and LO pQCD do not differ very much for CTEQ6L. For CTEQ6M, however, the curve from PYTHIA is nearly identical to that from CTEQ6L, but it differs by more than a factor of 3 from the LO pQCD result for the same PDF. From this surprising result, which can also be observed in Table 4, one sees again the big uncertainties due to different parameters, models, PDFs and shadowing. The NLO calculation including shadowing from Vogt:2001nh () lies about a factor of 2 above the PYTHIA yield with CTEQ6L.

V Heavy quark production in the QGP

v.1 Charm production at RHIC

In this section, charm production during the quark gluon plasma phase in central Au+Au collisions at RHIC will be investigated within the framework of the parton cascade BAMPS (cf. Section II). In contrast to Section III, where we limited ourselves to a static medium with just two processes, we consider now the full BAMPS simulation of the expanding fireball with all interactions from II. In that framework, because of the processes, rapid thermalization Xu:2007aa (), for instance, or the build-up of the elliptic flow Xu:2008av () of the gluonic medium can be explained.

The initial distribution of the charm quarks is simulated using PYTHIA with the parton distribution function CTEQ6l and scaled to heavy ion collisions as described in Section IV.1. Since light quarks are not implemented in BAMPS yet, they are treated as massless gluons to take conservation of energy and particle number into account. In this gluonic medium charm quarks are produced by gluon fusion and interact with gluons in elastic scatterings within our model. For both processes, LO cross sections (cf. Section II) and a constant coupling of are used. First, we study charm production without multiplying the cross section by any factor, although we will occasionally employ later. Charm annihilation is very unlikely due to the small number of charm quarks which are produced and can be neglected Andronic:2006ky ().

The primary sources of produced charm quarks at RHIC are initial hard parton scatterings during nucleon-nucleon collisions. With PYTHIA and CTEQ6l the initial yield amounts to 9.2 charm quark pairs as was shown in the previous section. The time evolution of the number of charm quarks during the quark gluon plasma phase within our BAMPS simulation is depicted in Figure 16.

Figure 16: (Color online) Number of charm quark pairs produced in a central Au+Au collision at RHIC according to BAMPS. The initial parton distribution is obtained with PYTHIA. Note that only a small interval of the -axis is shown. Therefore, we also depict the full range in the inset. In addition, the experimental value for the number of final charm pairs is plotted Adler:2004ta ().

The number of charm pairs increases slightly in the first by 0.3 pairs on average due to the high initial temperature of the plasma. Subsequently, the number saturates at a value of 9.5 charm pairs after , which indicates a marginal charm production of only 3 % of the total final charm quarks in the QGP. To illustrate this, the full range of the -axis is shown in the plot in the inset. Additionally, the experimental value of charm pairs, which is obtained from analyzing meson decays, is also shown. However, we calculated and its errors from and since these are the only values given in Adler:2004ta (). The number of charm pairs estimated with BAMPS and PYTHIA is below the experimental value, but still within the errors. This deviation results from the initial charm number obtained with PYTHIA and CTEQ6l, which is also smaller than the experimental value in proton-proton collisions (cf. Figure 14).

In addition, experimental data show that the number of produced charm quarks at RHIC scales with the number of binary collisions Adler:2004ta (); Adams:2004fc_STAR_dcsdY_cstot (); :2008hja (). That is, charm quarks are mainly created in initial hard parton scatterings and not in the QGP, which is in good agreement with our findings.

The small charm yield in the QGP at RHIC confirms our results from Section III.2.1, where we estimated a huge chemical equilibration time scale for charm quarks in a gluonic box at RHIC temperatures.

The number of charm quarks produced in the QGP is very sensitive to the initial gluon distribution. Therefore, we employed in addition to PYTHIA the color glass condensate (CGC) and the mini-jet model for the initial gluon distribution (cf. Section IV.2). Both models lead to an increase of a factor of 2.5 in the number of charm pairs produced in the QGP compared to PYTHIA, as is shown in Figure 17.

Figure 17: (Color online) As in Figure 16, but with PYTHIA, CGC, and mini-jet initial conditions for gluons. For better comparison, the initial charm distribution from PYTHIA is used for all models.

This larger yield is a result of the larger initial gluon energy density of both models compared to PYTHIA, as discussed in Section IV.2.

Often, the charm production cross section is multiplied by a factor of 2, which leads (naturally) to twice as many charm quarks. Nearly the same effect is attained by altering the charm mass from 1.5 GeV to 1.3 GeV, as depicted in Figure 18.

Figure 18: (Color online) As in Figure 16, but varying the charm mass and the factor for the cross section of .

Either changing the initial gluon distribution, choosing a smaller charm mass, or employing a factor can enhance the charm production by about a factor of 2, but the ratio of charm quarks produced in the QGP and total charm yield is still just 6 %. A combination of all these parameter changes (minijet initial conditions with high initial energy density, small charm mass of and factor of 2) yields 3.4 charm quarks which are produced during the QGP phase. This corresponds to 27 % of the final charm quarks, a value being the upper limit of charm production in our model.

In Table 5 our results are summarized.

IC
PYTHIA 1.5 1
Mini-jets 1.5 1
CGC 1.5 1
PYTHIA 1.5 2
PYTHIA 1.3 1
Mini-jets 1.3 2
Table 5: Comparison of the total charm pair yield and the number of charm pairs at mid-rapidity () after the QGP phase at central Au+Au collisions at RHIC simulated with BAMPS. Masses are in units of GeV. stands for factor, and IC for initial conditions.

In addition to the total charm numbers, we list the numbers at mid-rapidity. Table 6 gives an overview about charm production at RHIC in other models.

Reference IC Comments
Present work 1.5 1 PYTHIA Gluonic medium,
Present work 1.5 1 Mini-jets Gluonic medium,
Present work 1.5 1 CGC Gluonic medium,
Present work 1.5 2 PYTHIA Gluonic medium,
Present work 1.3 1 PYTHIA Gluonic medium,
Present work 1.3 2 Mini-jets Gluonic medium,
Levai:1997bi () 1.5 1 HIJING Gluonic medium
Levai:1997bi () 1.5 1 HIJING Quarks and gluons with thermal masses
Levai:1997bi () 1.2 1 HIJING Gluonic medium
Levai:1997bi () 1.2 1 HIJING Quarks and gluons with thermal masses
Levai:1997bi () 1.5 1 Minijets Gluonic medium
Gavin:1996bx () 1.2 1 Hydro
Duraes:2004zt () 1.5 1 Therm. running coupling
Duraes:2004zt () 1.2 2 Therm. running
Duraes:2004zt () 1.5 1 Therm. constant
Duraes:2004zt () 1.2 2 Therm. constant
Muller:1992xn () 1.5 2 HIJING
Levai:1994dx () 1.5 2 HIJING
Table 6: Comparison of charm pair yield in the QGP at RHIC from various models. Masses are given in GeV. stands for factor and IC for initial conditions.

The big deviations of the various models result from the influence of many parameters, which are often chosen differently in different models. For instance, there is disagreement regarding the temperature of the QGP, thermalization time scale, constant or running , thermal masses of quarks and gluons, volume, energy density, factor or charm mass. Within our results we also saw the sensitivity of the production rates concerning some of these parameters.

Another interesting topic we want to address is the question of whether charm quarks are chemically equilibrated. For that, we investigate the fugacity defined in Equation (24) for charm quarks which are located in central tubes with radius in transverse direction and longitudinal boundaries at space-time rapidity . The equilibrium value for the charm quark number is computed with the effective temperature of the gluonic medium. Both variables are shown in Figure 19 for all three initial models.

(a) (b)
Figure 19: (Color online) Evolution of the effective gluon temperature (a) and the charm quark fugacity (b) in the center of the collision ( and ) at RHIC for different initial models. The insets show the evolution of both quantities for PYTHIA initial conditions until when the temperature drops below the phase transition temperature of about Cheng:2006qk (); Aoki:2006br (); Aoki:2009sc ().

The temperature is quite large for RHIC energies, because we consider only a gluonic medium. If we took also light quarks into account, the temperature of the medium would be smaller.

The number of charm quarks is below the equilibrium value and increases with time. As we saw in Section III.2.1, however, their equilibration time scale is by far too large to see a significant rise before hadronization. In addition, the production rate decreases with time due to the accompanied temperature decline Zhang:2008zzc (). If the system evolves for a longer time, as is shown in the insets in Figure 19, the fugacity increases dramatically because of the decreasing temperature. At the temperature of the fugacity is about 4, which is of the same order as the result from Andronic:2006ky (); Andronic:2007ff () for the same temperature. However, in Andronic:2006ky (); Andronic:2007ff () the fugacity is computed on a hadronic level from the equilibrium number of charmed hadrons, whereas in our approach the fugacity is calculated on a quark level from the equilibrium number of charm quarks. Taking this and the difference of the considered volume into account, it is understandable that the values of the fugacity are not exactly the same.

At this low effective temperature of , the energy density in the considered volume is about (in a totally equilibrated gluonic medium it would be ). As a note, if the energy density in a cell in BAMPS drops below the critical energy density of the particles stream freely and do not interact anymore. In the present simulation this energy density in the considered volume corresponds to (cf. in an equilibrated medium).

v.2 Charm production at LHC

Due to the higher center of mass energy at LHC than at RHIC, more charm quarks are produced during the QGP phase and the yield from initial nucleon-nucleon scatterings is also much higher at LHC. However, since the charm production in the QGP increases almost exponentially with the QGP temperature but the initial charm yield only logarithmically with the collision energy Zhang:2008zzc (), the importance of the former rises significantly.

For the initial gluon distribution we use again PYTHIA, the CGC, and the mini-jet model. In Section IV.3.2 we estimated with PYTHIA and CTEQ6l about 62 charm pairs for the initial yield at LHC, which evolve in the QGP as a function of time as shown in Figure 20.

Figure 20: (Color online) Number of charm pairs in a central Pb+Pb collision at LHC simulated with BAMPS with PYTHIA, CGC, and mini-jet initial conditions for gluons.

With initial conditions from PYTHIA about 11 charm pairs are produced, which corresponds to about 15 % of the total final charm quarks. In contrast to RHIC, where nearly all charm quarks are produced in initial hard parton scatterings, at LHC a sizable fraction is created in the QGP. For initial gluon distributions from the CGC or the mini-jet model this fraction is even higher. In the latter framework the charm yield in the QGP is actually comparable to the initial yield (47 % of the total charm comes from the QGP). The introduction of a factor of 2 or the lowering of the charm mass to increases these values again by a factor of 2.

Table 7 summarizes our results and Table 8 compares them to other models.

IC
PYTHIA 1.5 1
Mini-jets 1.5 1
CGC 1.5 1
PYTHIA 1.5 2
PYTHIA 1.3 1
Table 7: As in Table 5, but for central Pb+Pb collisions at LHC.
Reference IC Comments
Present work 11 1.5 1 PYTHIA Gluonic medium,
Present work 55 1.5 1 Mini-jets Gluonic medium,
Present work 26 1.5 1 CGC Gluonic medium,
Present work 23 1.5 2 PYTHIA Gluonic medium,
Present work 20 1.3 1 PYTHIA Gluonic medium,
Present work 38 1.5 2 PYTHIA Gluonic medium,
Levai:1997bi () 43 1.5 1 HIJING Gluonic medium
Levai:1997bi () 94 1.5 1 HIJING Quarks and gluons with thermal masses
Levai:1997bi () 102 1.2 1 HIJING Gluonic medium
Levai:1997bi () 245 1.2 1 HIJING Quarks and gluons with thermal masses
Levai:1997bi () 21 1.5 1 Minijets Gluonic medium
Gavin:1996bx () 23 1.2 1 Hydro.
BraunMunzinger:2000dv () 5 1.5 1 SSPC & HIJING , , ,
Zhang:2008zzc () (7) 1.3 1 Therm. , NLO, at mid-rapidity
Zhang:2008zzc () (33) 1.3 1 Therm. , NLO, at mid-rapidity
Table 8: As in Table 6, but for central Pb+Pb collisions at LHC.

As for RHIC (cf. previous section), most of these predictions from other models are only for thermally produced charm quarks, whereas in our calculations the prethermal charm production within the first