Production of Charm Quarks in a Parton Cascade Model for Relativistic Heavy Ion Collisions at $\sqrt{s_{\textrm {NN}}}$= 200 GeV

TL;DR

This paper investigates charm quark production in relativistic heavy ion collisions using a parton cascade model, highlighting the roles of rescattering and gluon radiation in different collision scenarios.

Contribution

It introduces a detailed parton cascade model for charm production, analyzing the effects of rescattering and gluon radiation in heavy ion collisions at 200 GeV.

Findings

Multiple collisions significantly enhance charm production.

Gluon radiation contributes notably to charm quark dynamics.

Rescattering effects are crucial in understanding charm quark evolution.

Abstract

We study the production and dynamics of heavy quarks in the parton cascade model for relativistic heavy ion collisions. The model is motivated by the QCD parton picture and describes the dynamics of an ultra-relativistic heavy-ion collision in terms of cascading partons which undergo scattering and multiplication while propagating. We focus on the dynamics of charm quark production and evolution in p+p and Au+Au collisions for several different interaction scenarios, i.e. collisions only between primary partons without radiation of gluons, multiple collisions without radiation of gluons and multiple collisions with radiation of gluons, allowing us to isolate the contributions of parton rescattering and radiation to charm production. We also discuss results of an eikonal approximation of the collision which provides a valuable comparison with mini-jet calculations and clearly brings out the importance of multiple collisions.

Tables

Table 1. Table 1:

The number of scatterings and time like branchings involving different sub-processes in $pp$ collisions at $\sqrt{s_{NN}}$ = 200 GeV.
process	Only Primary-Primary	Eikonal	Multiple	Multiple collisions
	collisions	approximation	collisions	and radiation
q+q $\to$ q+q	0.204	0.255	0.269	0.364
q+$\overline{q}$ $\to$ q+$\overline{q}$	0.002	0.002	0.002	0.012
q+$\overline{q}$ $\to$ g+g	0.006	0.008	0.008	0.032
q+$\overline{q}$ $\to$ g+$\gamma$	0.0	0.0	0.0	0.0
q+$\overline{q}$ $\to$ 2 $\gamma$	0.0	0.0	0.0	0.0
q+g $\to$ q+g	1.062	1.447	1.617	2.874
q+g $\to$ q+$\gamma$	0.0	0.001	0.001	0.004
g+g $\to$ q+$\overline{q}$	0.008	0.011	0.014	0.052
g+g $\to$ g+g	0.985	1.533	1.852	4.814
total hard a+b $\to$ c+d	2.267	3.257	3.764	8.152
q $\to$ q+g				0.710
g $\to$ g+g				2.989
g $\to$ q+$\overline{q}$				0.255
q $\to$ q+$\gamma$				0.0
Total a $\to$ b+c				3.954
$N_c$ (=$\overline{N_c}$)/event	4.98 $\times$ 10-4	7.62 $\times$ 10-4	9.08 $\times$ 10-4	34.95 $\times$ 10-4

Table 2. Table 2:

The number of scatterings and time like branchings involving different sub-processes in $Au+Au$ collisions at 200A GeV.
process	Only Primary-Primary	Eikonal	Multiple	Multiple collisions
	collisions	approximation	collisions	and radiation
q+q $\to$ q+q	111.2	185.3	227.8	591.5
q+$\overline{q}$ $\to$ q+$\overline{q}$	1.1	1.4	2.2	34.0
q+$\overline{q}$ $\to$ g+g	3.7	5.4	7.6	86.7
q+$\overline{q}$ $\to$ g+$\gamma$	0.0	0.1	0.0	0.8
q+$\overline{q}$ $\to$ 2 $\gamma$	0.0	0.0	0.0	0.0
q+g $\to$ q+g	527.7	1056.9	1538.3	6370.0
q+g $\to$ q+$\gamma$	0.3	0.4	0.7	11.8
g+g $\to$ q+$\overline{q}$	3.7	8.0	16.3	164.8
g+g $\to$ g+g	432.7	1109.8	2020.6	13304.0
total hard a+b $\to$ c+d	1080.6	2367.3	3813.4	20563.7
q $\to$ q+g				787.0
g $\to$ g+g				3868.2
g $\to$ q+$\overline{q}$				345.5
q $\to$ q+$\gamma$				0.0
Total a $\to$ b+c				5000.8
$N_c$ (=$\overline{N_c}$)/event	0.33	0.82	1.33	5.21

Equations

$$q_i{q}_j$$

$$q_i{\bar{q}}_i$$

$$q{\bar{q}}_i$$

$$q_{i}g$$

$$gg$$

$$g^{\ast }$$

$$\frac{d\hat{\sigma }}{d\hat{t}}=\frac{1}{16\pi {\hat{s}}^2}\sum |\mathcal{M}{|}^2,$$

$$\sum |\mathcal{M}{|}^2={\pi }^2{\alpha }_s^2(Q^2)\left[a_1+a_2+a_3+a_4+a_5+a_6\right]$$

$$a_1$$

$$a_2$$

$$a_3$$

$$a_4$$

$$a_5$$

$$a_6$$

$$\hat{\sigma }(\hat{s})=\frac{1}{16\pi {\hat{s}}^2}{\int }_{M^2-\hat{s}(1+\chi )/2}^{M^2-\hat{s}(1-\chi )/2}d\hat{t}\sum {\left|\mathcal{M}\right|}^2,$$

$$\chi =\sqrt{1-\frac{4M^2}{\hat{s}}}.$$

$${\hat{\sigma }}_{\mathrm{gg}\to \mathrm{Q}\bar{\mathrm{Q}}}(\hat{s})$$

$$\sum |\mathcal{M}{|}^2$$

$${\hat{\sigma }}_{\mathrm{q}\bar{\mathrm{q}}\to \mathrm{Q}\bar{\mathrm{Q}}}(\hat{s})=\frac{8\pi {\alpha }_s^2(Q^2)}{27{\hat{s}}^2}(\hat{s}+2M^2)\chi .$$

$$Qq(\bar{q})\to Qq(\bar{q})\&Qg\to Qg.$$

$$\frac{d\hat{\sigma }}{d\hat{t}}=\frac{1}{16\pi (\hat{s}-M^2{)}^2}\sum |\mathcal{M}{|}^2.$$

$$\sum |\mathcal{M}{|}^2$$

$$b_1$$

$$b_2$$

$$b_3$$

$$b_4$$

$$b_5$$

$$b_6$$

$${\hat{\sigma }}_{tot}=\frac{1}{16\pi (\hat{s}-M^2{)}^2}{\int }_{-(\hat{s}-M^2{)}^2/\hat{s}}^{0}d\hat{t}\sum |\mathcal{M}{|}^2.$$

$${\mathcal{P}}_a=\sum _{b,c}\frac{{\alpha }_{abc}}{2\pi }{P}_{a\to bc}\frac{d{Q}^2}{Q^2}dz$$

$$\frac{d\sigma }{d{p}_T^{2}d{y}_{1}d{y}_2}$$

$$R_{\text{AA}}=\frac{d{N}_{\text{AA}}/d{p}_{T}dy}{N_{\text{coll}}\times d{N}_{\text{pp}}/d{p}_{T}dy}$$

Full text

I Introduction

The quark-gluon plasma (QGP), a deconfined strongly interacting matter, which filled the nascent universe a few micro-seconds after the Big Bang, is now being routinely produced and studied in relativistic heavy ion collisions at RHIC and LHC. Its presence has been confirmed through a host of previously predicted phenomena Gyulassy:2004zy ; Muller:2006ee ; Muller:2012zq ; Wang:2016opj , as well as through novel and unexpected discoveries, such as strong elliptic flow Kolb:2003dz ; Song:2010mg , initial state fluctuations leading to higher order flow Alver:2007qw ; Alver:2008zza ; Alver:2010gr of hadrons as well as phenomena related to parton recombination Fries:2003vb ; Greco:2003xt .

Heavy quarks serve as excellent probes of the QGP fireball as they are primarily produced by early-state hard scatterings and thus have the potential to probe the whole space-time history of the transient matter. It was generally expected that heavy quarks do not lose as much energy as light partons while traversing the quark gluon plasma due to their large mass. However, early estimates of drag and diffusion coefficients Svetitsky:1987gq ; GolamMustafa:1997id as well as an early estimate of heavy-quark energy loss due to radiation of gluons Mustafa:1997pm suggested that heavy quarks could possibly lose as much energy as light quarks and gluons during their passage through the QGP. These predictions have been subsequently confirmed by multiple other calculations GolamMustafa:1997id ; Gyulassy:2000fs ; Djordjevic:2003zk ; Moore:2004tg ; Peigne:2008nd ; Das:2010tj ; Gossiaux:2010yx ; Abir:2012pu ; Das:2015ana , with tools ranging from simple phenomenological models used to estimate the medium modification of charm production for heavy ion collisions Younus:2010sx , to more detailed comparisons with models of energy loss embedded in a hydrodynamic evolution of the plasma Cao:2013ita . Studies including the later hadronic stages of the evolution have also indicated that charm-mesons may lose energy due to (resonant) scatterings with hadrons vanHees:2007me ; Rapp:2009my . The experimental observation of large elliptic flow for D and B mesons as well as the small value of $R_{AA}$, the nuclear modification factor, have added experimental evidence to these findings Andronic:2015wma ; Xie:2016iwq . Several of these aspects have been studied in great detail using the Parton Hadron String Dynamics Model at energies reached at RHIC and LHC Cassing:2000vx ; Bratkovskaya:2003ux ; Bratkovskaya:2004ec ; Linnyk:2008uf ; Linnyk:2008hp ; Song:2015sfa ; Song:2015ykw .

In the high momentum domain the production of heavy quarks via the interaction between (mini)-jets following the initial scattering as well as due to interaction of thermal partons at high incident energies Levai:1994dx ; Lin:1994xma ; Younus:2010sx has been studied. An estimate at next to leading order using the formalism developed by Mangano et al. Mangano:1991jk ; Frixione:1997ma was also used to suggest Younus:2011mn that angular correlations of charm-anticharm quarks in p+p collisions would be drastically different for various production processes, e.g., gluon fusion or quark-antiquark annihilation, gluon splitting and when a gluon is radiated off one of the legs of the scattering diagram. This angular correlation is expected to get even more interesting for nucleus-nucleus collisions because of the interaction of the heavy quarks with the medium.

State of the art approaches for the investigation of heavy quarks in a hot and dense QCD medium need to incorporate a dynamical treatment for both, the heavy quarks as well as the QCD medium. This can be either accomplished via hybrid approaches that utilize a hydrodynamic evolution for the medium in concert with a dynamical interaction model for the heavy quarks Moore:2004tg ; He:2011qa ; Gossiaux:2011ea ; Cao:2013ita or in the context of a purely microscopic transport for both, the heavy quarks and the medium Uphoff:2012gb ; Scardina:2014lda ; Uphoff:2014hza . Parton Cascade Models Geiger:1991nj ; Geiger:1992si ; Geiger:1992ac ; Geiger:1993py ; Geiger:1993ix ; Geiger:1994he ; Bass:2002fh as used in this study, fall into the latter category. Over the past decade, PCMs have been used for multiple studies of heavy quark dynamics in a QCD medium, ranging from the study of heavy quark energy loss in infinite matter Younus:2013rja , to studies of equilibration and elliptic flow built-up at RHIC and LHC energies Uphoff:2012gb ; Uphoff:2014hza . However, most of these studies focused on the low to intermediate momentum regime (mostly due to the availability of data in that regime).

As more heavy quark data at high transverse momenta become available, it is now the right time to ask whether approaches for the dynamics of heavy quarks based purely on perturbative QCD are able to describe the experimental observations in that domain. In addition, even if pQCD approaches cannot fully account for the observed features, it remains an open question to what extent the observed features of final state observables are already imprinted on the heavy quarks prior to the formation of the QGP, i.e. via initial state effects as well as contributions from early non-equilibrium evolution. To address these questions we use a Monte Carlo implementation of the PCM, VNI/BMS Bass:2002fh , that is based on a Boltzmann transport description with pQCD matrix elements for parton-parton interactions. This particular PCM implementation also contains parton shower emission. In the present work we discuss the implementation of heavy quarks into the PCM and calculate charm production in relativistic collisions of gold nuclei at $\sqrt{s_{\textrm{N}N}}$ = 200 GeV. We consider several different interaction scenarios for pp and nucleus-nucleus systems, namely, an implementation of the eikonal approximation for a direct comparison with minijet calculations and to clearly demonstrate the dynamics which emerges as the partons change their momenta after collisions (which is neglected in the eikonal approximations), a scenario which involves only collisions between primary partons to clearly bring out the consequences of multiple scattering, and finally a scenario which additionally includes gluonic multiplication by radiation of gluons following scattering. The last scenario considerably increases the number of scatterings and partons (including the heavy quarks) are capable of losing energy by radiation of gluons.

II Formulation

The details of the Monte Carlo implementation of the parton cascade model have been discussed in  Bass:2002fh . However, given how our understanding of the hot and dense QCD matter created in relativistic heavy-ion collisions has changed since the inception of the PCM, a brief discussion of its past and future uses and relevancy is in order: at the time when the first Parton Cascade Models were developed, the notion of a Quark-Gluon-Plasma was still in its infancy and it was thought of as a weakly interacting gas of quark and gluons that could be described using perturbative QCD. Having a space-time evolution based on the Boltzmann equation with quark and gluon degrees of freedom and pQCD cross sections was thought to provide a full picture of the dynamics of the deconfined system up to hadronization. With the discovery of the near perfect fluidity of the QGP at RHIC and LHC, this picture has been found to be inadequate and most PCM implementations have failed to generate the observed amount of collective flow, while utilizing pQCD based interactions for quark and gluon degrees of freedom 111the one exception is the BAMPS implementation .

However, the production of hard probes – jets, photons and heavy quarks – can be perfectly well understood in terms of perturbative QCD production cross sections and their dynamical evolution in the early reactions stages prior to the formation of the QGP at $\tau_{0}\approx 1$ fm/c is dominated by interactions with fairly large momentum transfers. Parton Cascade Models are thus ideally suited to describe the early out-of-equilibrium evolution of hard probes prior to the formation of a thermalized Quark-Gluon-Plasma. Their ability to describe the dynamical evolution of hard probes during the early times of the collision evolution, including multiple interactions (rescattering and gluon splitting), provides significant added insight beyond what can be gained from initial state production Monte-Carlo codes such as PYTHIA or Herwig 222historically, VNI/BMS is based on PYTHIA 6 with the addition of a Boltzmann equation solver for the space-time evolution of all partons.

Thus, we shall utilize the PCM for the study of heavy quarks during early times of Au+Au collisions at RHIC. In the following we shall first describe the implementation of heavy quark (Q) production and interactions before moving on to our analysis.

II.1 Parton cascade model involving light quarks and gluons

The $2\rightarrow 2$ scatterings included in the version of the parton cascade model utilized for our work, VNI/BMS, are:

[TABLE]

The $2\rightarrow 3$ reactions are included via time-like branchings of the final-state partons:

[TABLE]

For details we refer the reader to Ref Bass:2002fh . We add that the gluon splitting in the above, includes the splitting $g^{*}\rightarrow Q\bar{Q}$ where $Q$ stands for heavy quarks (charm or bottom) if the virtuality is large enough to admit this process. We shall see later that this process plays an important role, which is likely to increase with the increasing centre of mass energy of the collision. Due to the specifics of the implementation of the collision term in VNI/BMS the reverse process, parton fusion, is not included, since it would lead to the propagation of partons off-shell, which is not included in the VNI/BMS formulation. Technically, parton fusion with on-shell partons can be formulated via a $3\rightarrow 2$ process that can be implemented using collision rates in the collision term. This has been successfully done in the BAMPS model Xu:2004mz . The main drawback of the lack of parton fusion in VNI/BMS is that it violates detailed balance in the collision term and thus leads to an incorrect equilibrium limit for the bulk in the infinite time and size limit. While this is a valid concern if we were to study the equilibration of bulk QCD matter in our approach, our focus on high-momentum heavy quarks that are not thought to thermalize fully Cao:2011et reduces the systematic uncertainty induced by the omission of these processes.

II.2 Heavy quark production by gluon fusion and annihilation of light quarks

At lowest order pQCD, heavy quarks are produced from fusion of gluons ($gg\rightarrow Q\bar{Q}$) and annihilation of light quarks ($q\bar{q}\rightarrow Q\bar{Q}$).

The partial differential cross-section can be written as

[TABLE]

where the summed spin and coloured averaged squared matrix element for the process $gg\rightarrow Q\bar{Q}$ is given by

[TABLE]

where

[TABLE]

The total cross-section $\hat{\sigma}(\hat{s})$ is obtained from the above by integrating over $\hat{t}$ :

[TABLE]

where

[TABLE]

We note that due to the mass of the heavy quark $M$ this cross-section remains finite and the total cross-section reduces to

[TABLE]

Similarly, the summed spin and coloured averaged squared matrix element for the process $q\bar{q}\rightarrow Q\bar{Q}$ is given by:

[TABLE]

so that the total cross-section becomes.

[TABLE]

Once again, these cross-sections remain finite due to the mass of the heavy quarks.

II.3 Heavy quark production due to flavour excitation

The implementation of heavy quark production due to flavour excitation has been subject of significant debate in the community. In principle, the (off-shell) heavy sea-quarks in the target (projectile) can be brought on-shell by scattering with other off-shell partons from the projectile (target) (or by scattering with on-shell partons that are produced during an earlier partonic interaction). This contribution was referred to as flavour excitation in the original paper of Combridge Combridge:1978kx . Its inclusion in the parton cascade model, where all parton scatter are treated at lowest order, has been questioned Lin:1994xma . It is also seen to give rise to a substantial production of charm and bottom quarks Geiger:1993py , as was also pointed out in the original work of Combridge. It has been argued that this contribution should be strongly suppressed if higher order coherent interference terms are accounted for (see Ref. Collins:1985gm and references there in).

In our calculation, we have excluded this process for the production of charm or bottom quarks. This is done in two steps. Firstly, we initialize the quark and gluon distribution functions at a rather low $Q_{0}^{2}$ of 0.589 GeV2 (see Bass:2002fh ) which effectively eliminates the presence of heavy sea-quarks. In addition, we disallow any partonic interaction if one of the partons is a heavy (charm or bottom) quark from the sea.

The flavour excitation process as it is understood now is included in NLO calculations, e.g., Ref. Mangano:1991jk ; Frixione:1997ma , and proceeds via a splitting of a gluon into a $Q\bar{Q}$ pair followed by scattering of one of the heavy quarks with another parton in the system. The contribution from this process to heavy quark production is not large. We re-emphasize that as mentioned earlier, the splitting of the final state gluon into a charm- anti charm pair (or for that matter bottom-anti bottom pair as well) is included in our calculations. We also add that Soft Collinear Effective Theory has been used recently to estimate the gluon splitting into heavy quarks Kang:2016ofv , which suggests that this contribution can be large at LHC energies. Once produced, the heavy quarks are allowed to scatter with quarks and gluons. We shall see later that this plays an important role in their dynamics.

II.4 Scattering of light quarks and gluons with heavy quarks

We consider elastic scattering of heavy quarks with light quarks and gluons, i.e,

[TABLE]

The differential scattering cross section is defined as:

[TABLE]

The summed spin and coloured averaged squared matrix element for the process $Qq(\bar{q})\rightarrow Qq(\bar{q})$ have been derived by Combridge Combridge:1978kx :

[TABLE]

The corresponding matrix element for the process $Qg\rightarrow Qg$ is given by:

[TABLE]

where

[TABLE]

The total cross section is formally defined as:

[TABLE]

We immediately note that the total cross-sections for both the processes diverge because of the pole at $\hat{t}$ equal to zero. We reduce the upper limit of integration to $\hat{t}=-Q_{0}^{2}$, where $Q_{0}^{2}=p_{0}^{2}$ and $p_{0}$ is the cut-off used in the parton cascade model to regularize the $2\rightarrow 2$ pQCD matrix element, as in PYTHIA.

II.5 Radiation of gluons by heavy quarks

The parton cascade model includes final state radiation of gluons or photons following a hard ($2\rightarrow 2$) scattering using the well established approach of parton showers Sjostrand:1985xi as implemented in PYTHIA, in a leading logarithmic approximation. It incorporates a sequence of nearly collinear splittings, $a\rightarrow bc$, where the initial parton $a$ is called mother parton and $b$ and $c$ are called daughters, which can split further and populate a tree-like structure. Strong and electromagnetic interactions allow for several different possibilities of splittings, which are all included in the PCM. The differential probability for a splitting to occur is given by:

[TABLE]

where the sum runs over all the allowed splittings, $\alpha_{abc}$ is equal to $\alpha_{em}$ for emission of photons and for QCD splittings it is $\alpha_{s}$. Further, $Q^{2}$ is the momentum scale of the splitting, $z$ gives fraction of energy carried by the daughter $b$ and $1-z$ is the remaining fraction, carried by the daughter $c$. The kernels describing the splitting $P_{a\rightarrow bc}$ are taken from Altarelli and Parisi Altarelli:1977zs .

The parton cascade model uses a cut-off $\mu_{0}$, taken as $\geq$ 1 GeV, to regulate the collinear singularities, by terminating the splitting once the virtuality of the time-like parton drops to $\mu_{0}$. The soft gluon interference is accounted for by implementing an angular ordering of the radiated gluons. The mass of the partons is accounted for by modifying the cut-off $\mu_{0}$ such that $\mu_{0}^{2}$= 1 GeV2 for gluons and $\mu_{0}^{2}=m_{q}^{2}+1$ (GeV2) for quarks Geiger:1991nj and thus radiation of gluons from heavy quarks is at the same footing as for that from light quarks. This value of $\mu_{0}$ is same as that used in PYTHIA for regularizing the final (and initial) state radiations. Increasing this will lead to a hardening of the transverse momentum spectra as well as a reduction in the number of collisions due to a reduction in number of gluons, while taking a smaller value would have an opposite effect. We shall see later that we get a reasonable description of the spectra for charm quarks both for $pp$ and $AA$ collisions and thus have kept it fixed, though exploratory calculations were done to verify the comments made above.

II.6 Production of heavy quarks using the minijet formalism

In order to validate the Monte-Carlo implementation of the PCM that we are using, we perform a separate calculation of heavy quark production using the minijet formalism Eichten:1984eu . The cross-section for the production of a heavy-quark pair ($Q\bar{Q}$) is written in terms of the rapidities of the two quarks $y_{1}$ and $y_{2}$ and their transverse momentum $p_{T}$ as

[TABLE]

In the above $f_{i}^{(a)}$ denotes structure function of the parton $i$ for the nucleon $a$, etc., and $x_{a}$ is the fraction of linear-momentum of the nucleon $a$ carried by the parton, $Q^{2}$ is the momentum scale, in the standard notation. However, in order to be able to directly compare to the parton cascade model, we fix the fragmentation scale at the momentum cut-off $p_{0}^{2}$ used in the calculation for the momentum cut-off and the coupling constant is kept at a fixed value of $\alpha_{s}$.

III Results

We shall first discuss our results for $pp$ collisions at $\sqrt{s_{NN}}$ = 200 GeV for several different parton interaction scenarios. Next we give our findings for $Au+Au$ collisions at the same centre-of mass energy for central collisions. Finally we compare the proton-proton results to the Au+Au results to get an idea about the emerging dynamics of the propagation of charm quarks, subject to semi-hard scatterings and radiation of gluons in the QCD medium produced by relativistic heavy ion collisions.

III.1 $pp$ collisions at $\sqrt{s_{NN}}=$ 200 GeV

In order to test the accuracy of our implementation of heavy quark production in VNI/BMS we first study collisions of protons at $\sqrt{s_{NN}}=$ 200 GeV.

We calculate the production of charm quarks due to $gg\rightarrow Q\bar{Q}$ and $q\bar{q}\rightarrow Q\bar{Q}$ processes using a mini-jet calculation with a GRV-HO parametrization of the parton distribution function. In order to compare it with our implementation of the PCM we fix $\alpha_{s}$ as 0.3 and keep $Q_{0}^{2}$ fixed at 0.589 GeV2 for the renormalization and factorization scales. Note that for a pure theory to theory comparison the specifics of the parton distribution function are not of particular importance, which is why we choose the well-established (even though somewhat dated) GRV-HO parametrization.

The minijet calculation is compared to a VNI/BMS calculation using an eikonal approximation, where the final momenta and flavours of the particles produced in a hard scattering are replaced by the momenta and flavours of the corresponding initial state partons. This procedure thus mimics the eikonal approximation. A comparison of these two calculations can be found in figures 1 and 2: the rapidity integrated transverse momentum spectrum and the $p_{T}$ integrated rapidity spectrum of the produced charm (or anti-charm) are given in Fig 1 and the $p_{T}$ spectra at several rapidities are given in Fig. 2. We observe satisfactory agreement between the VNI/BMS calculation and the minijet reference. The minor differences at very low and very high $p_{T}$ arise due to the $p_{T}$-cut-off in the parton cascade model, which is also reflected in $x_{\textrm{min}}$ which enters in the parton cascade model. This difference is again reflected in the rapidity distribution as well.

Having established reasonable quantitative agreement between the two independent calculations and thus confirming that heavy quark interactions are correctly implemented in the Monte Carlo scheme of VNI/BMS, we proceed to study the dynamics of the cascade evolution. As a first step we consider a calculation in which only primary-primary parton interactions are permitted, i.e. each parton is allowed to interact only once. A comparison of these calculations with the results of the eikonal approximation then immediately highlights the effect of multiple parton interactions in the eikonal approximation, already in $pp$ collisions. Again we show the rapidity integrated transverse momentum spectrum and the $p_{T}$ integrated rapidity spectrum of the produced charm (or anti-charm) in Fig 3. The spectra are quite similar in shape, but differ by the factor which accounts for the increased number of collisions in the eikonal approximation. Details of the evolution of the cascade are quantified in Table 1. We see from Table 1 that the leading process which contributes to the production of charm quarks, i.e., $gg\rightarrow c\overline{c}$, is about 50% more likely in the eikonal approximation compared to the case when only primary-primary collisions are allowed.

Now we proceed to results with the proper accounting of multiple scatterings along with the changes in the flavour and momenta of the produced partons, before they scatter again. These studies, we believe, will help us identify the effects of multiple scatterings already at the level of p+p collisions.

Looking at Table 1, we find there is an increase in the number of collisions by about 20% compared to the eikonal approximation. Compared to the primary-primary collision calculation, the increase is about 60%. While the increase in the number of collisions for the case of multiple scatterings over and above that of the case of primary-primary collisions is expected, the increase in the number of collisions above that for the eikonal approximation comes as a surprise. This may be due to interactions that impart a large momentum to partons that previously would have been below the momentum cutoff for an interaction now acquiring sufficient transverse momentum for subsequent interactions.

The rapidity integrated $p_{T}$ and the $p_{T}$ integrated rapidity spectra are shown in Fig. 4. We see that a proper accounting of the changing momenta and the flavours of the partons following a hard scattering enhances the production of charm quarks at moderate transverse momenta. This effect will undoubtedly be even more prominent in $Au+Au$ collisions, where the probability of multiple scatterings increases considerably. However there are already indications of non-uniform variation of the factor by which the two calculations differ as a function of $p_{T}$ and $y$, presumably due to the larger density of interacting partons at midrapidity.

Finally, we include time-like branchings of partons in our calculation, which leads to a rapid increase in the number of partons at midrapidity and subsequent multiple scatterings. As expected, we observe a considerable increase in the number of hard collisions, fueled by the time-like branchings which leads (mainly) to gluon multiplication. Once again we notice that the increase in $gg\rightarrow q\bar{q}$ is reflected in the increase of the production of charm quarks. This close correspondence may be due to fact that the medium generated in $pp$ collisions is not yet dense and long lived enough to have several scatterings by partons which can reduce their momenta such that further collisions do not produce heavy quarks. The momenta of charm quarks at medium and large $p_{T}$ are also greatly affected (see Fig. 5) while a comparison of the rapidity distribution suggests increased interactions at central rapidities (Fig. 5).

Even though not substantial, we also observe perceptible change in the momentum distribution obtained in eikonal approximation compared to the case for proper multiple collisions, as these actually evolve. This would imply a slight correction to the standard calculations for the proton-proton baseline for the medium modification of particle production if estimated within the eikonal approximation. However, one would require very high precision data to quantitatively investigate this effect.

Charm production in $pp$ collisions has been measured by STAR experiment Ye:2014ska by measuring $D^{0}$ production and other $D$ mesons. We can compare our results for the distribution of charm quarks by accounting for the fragmentation ratio ($\approx 0.565$) for the $c\rightarrow D^{0}$ fragmentation as also in the above reference for a comparison results from with PYTHIA (tuned) and FONLL calculations. This amounts to assuming that the fragmentation function varies as $\delta(1-z)$ (see later). Our results for the calculations incorporating multiple scattering and radiations are shown in Fig. 6. In view of the $p_{T}^{\text{cut-off}}$ and $\mu_{0}$ used to regularize the matrix elements and fragmentation functions used in our work we show the results for our calculations only for $p_{T}\geq$ 2 GeV. A fair agreement is seen without any attempt to adjust any parameters.

III.2 $Au+Au$ collisions at $\sqrt{s_{NN}}=$ 200 GeV

When and where are most of the charm quarks produced in a nucleus-nucleus collision? How often do they undergo collisions? How often do they radiate gluons? How does the charm production depend on gluon multiplication? These and many other interesting details of charm production are easily discernible in the space-time description of the cascading partons in our calculations. We proceed to address these questions in the following sections.

III.2.1 Dynamics of charm production in cascading partons

We start by calculating the spatial distribution of the production vertices of the charm quarks along the beam axis (Fig. 7) for the different interaction scenarios introduced in the previous section. We see that for all cases studied (primary interactions only, secondary rescattering and rescattering plus radiation) the production of charm quarks is symmetrically concentrated around the region of complete overlap ($z$=0). We further note that charm-quark production is mostly limited to the zone around $\pm$ 0.5 fm along the beam axis, due to the Lorentz-contraction of the colliding nuclei.

Focusing on the case where multiple collisions along with the fragmentation of partons is permitted, we note that the number of produced charm quarks is significantly larger than in the other cases. However, the fraction of charm quarks produced at later separations is now much smaller than that for the case when only multiple scatterings are permitted (see Fig. 7). Multiple collisions occur predominantly in the region of complete overlap. The colliding partons radiate gluons and rapidly lose energy, while the parton clouds continue to propagate through each other. However by the time the gluons undergo third or fourth interaction their energy has dropped significantly and subsequent interactions may not have enough energy to produce heavy quarks. On the other hand, if the partons were not allowed to radiate gluons as in the case of multiple scatterings without radiations, even the later collisions have enough energy to produce heavy quarks. This scenario is further reinforced by the observations of Fig. 8, which shows the distribution of the production times of charm quarks.

Charm production generally commences by the time the nuclei start touching other, increases rapidly as they interpenetrate and peaks when they are on top of each other and falls off as the nuclei disengage. A similar trend is seen in the behavior of the $z$-distribution of the production vertices. When the multiple scatterings are permitted (without radiations), the production of charm quarks continues (albeit at a very small rate) till almost 3 fm/$c$ after the complete over-lap of the nuclei, though most of them are produced by about 0.5 – 1.0 fm/$c$. For the last case, when the radiation of gluons is permitted, a rapid multiplication of partons and a considerable increase in multiple scatterings is observed. As a consequence, there is a large increase in production of charm quarks as well. However, this additional production is now mostly limited to early times, say about 0.3–0.5 fm/$c$. We attribute this restriction to early times to the energy lost by partons due to radiation of gluons following multiple scatterings: collisions at later times involve energies which are below the threshold for production of charm. We emphasize however that even though the production of charm quarks is essentially over by about 0.5 fm/$c$ the collisions (and radiations) will continue till about 3–4 fm/$c$ and thus charm quarks (and other high energy quarks or gluons) will continue to lose energy due to semihard processesBass:2002fh .

In Fig 9 we show the actual distribution of $t$ and $z$ positions of the production vertices (for about 2700 charm quarks) for the three cases discussed above. We notice that the $t$ and $z$ co-ordinates of the production vertices for the primary-primary collision case are located close to the world lines describing the motion of the centers of the two nuclei and are limited to the time of overlap. For the case of multiple collisions only, the production continues to late times and fills up larger distances from the point of complete overlap. For the case when both multiple scatterings and radiation (of gluons) are permitted, we have a concentration of production vertices around the time of overlap, as at later times the partons do not have enough energy to cross the threshold for the production of charm-anticharm pairs. The corresponding $x$ and $y$ co-ordinates for the three cases are shown in Fig. 10. We have verified that these distributions closely follow the distribution of binary collisions $n_{\textrm{BC}}(x,y)$ obtained using a Glauber model Kolb:2003dz .

How often do the charm quarks scatter after production and how often do they radiate gluons? The distribution of the number of scatterings suffered by charm quarks for the two cases that allow rescattering after production is given in Fig. 11. The produced charm quarks, being secondary particles, will not scatter when only primary-primary parton interactions are considered. We also show the distribution of the number of gluon splittings by the charm quarks (Fig. 12) for the case when we consider both multiple scattering and multiplication of gluons due to their radiation off scattered partons.

When only multiple scatterings are permitted, charm quarks will rarely scatter more than 3–4 times. We attribute this limitation to the number of partons remaining limited and thus the possibility of an interaction decreasing rapidly as the parton clouds disengage. However, when the fragmentation of final state partons produced in a hard scattering is allowed, the number of gluons increases rapidly and leads to a significantly larger number of scatterings and a significant production of gluons due to radiation. Now we find that charm quarks may scatter up to 20 times, with 8–10 scatterings not being too uncommon. Despite the restriction of the PCM to interactions above a $p_{T}$ cut-off, this large number of semi-hard scatterings has the potential to start propelling charm quarks towards thermalization!

Finally in Fig. 12 we provide results for the occurrences of the splitting of gluons off charm quarks following a semi-hard scattering (note that we only count the number of initial splittings, not the total number of gluons that may be emitted in each of these processes). While one splitting is most frequent, we do observe some occurrences of multiple splittings.

Table 2 provides a quantitative summary of the above analysis. For example, we find that total number of semi-hard scatterings for all partons increases by a factor of 2.2 when we perform calculations in the eikonal approximation, a factor of about 3.5 when we perform the calculation using multiple scatterings and by a factor of about 19 when we perform calculations allowing for multiple scatterings and multiplication of partons due to time-like branchings or fragmentations – in comparison to the calculation involving only primary-primary interactions. The number of collisions for the leading process for charm production included in the counting of, $gg\rightarrow q\bar{q}$, rises by a factors of 2.2, 4.3 and 44 respectively.

The production of charm quarks shows increases by factors of 2.5, 4, and 15.2 respectively. We note that while the rise in charm production in going from primary-primary scatterings to the eikonal approximation to multiple scatterings closely follows the trend for the rise in the $gg\rightarrow q\bar{q}$ process, the rise in charm production when we consider multiple scatterings with parton multiplications is much less (only a factor of about 15) compared to the rise in such processes (about a factor of 44) for light parton interactions. As discussed previously, we attribute this reduced increase to gluon multiplication leading to an ensemble of less energetic partons so that despite the increase in the number of multiple collisions that partons suffer, they may not be energetic enough for the production of charm in later interactions.

III.2.2 Spectra of charm quarks in $Au+Au$ collision at $\sqrt{s_{\textrm{NN}}}$ = 200 GeV

In this section we proceed to discuss our results for spectra of charm quarks produced in central collisions for $Au+Au$ at $\sqrt{s_{\textrm{NN}}}$ = 200 GeV, again using the same set of interaction models as in the previous section.

As a first step we compare our results for the implementations of the eikonal approximation with those for primary-primary collisions (Fig. 13). We find that $p_{T}$ spectra are similar but not identical in shape and differ by factors of 2–5 varying with $p_{T}$ and also $y$.

The results for a comparison of the multiple scattering to that of the eikonal approximation are shown in (Fig. 14). The $p_{T}$ spectra for both calculations are quite similar for large $p_{T}$. However, we observe considerable differences at intermediate transverse momenta: while the results for the eikonal approximation mimic a power-law, those for the case of multiple scattering exhibit more of an exponential shape. This is a consequence of multiple scatterings altering the parton momenta, which are properly accounted for in the parton cascade model. This observation is also supported by estimates of the production of charm quarks due to multiple scatterings among partons produced in semi-hard (mini)-jets by several authors Younus:2010sx ; Levai:1994dx ; Lin:1994xma . We note that the production of charm is proportionately higher at central rapidities as compared to larger rapidities in the case of multiple scattering.

Now we study the effect of gluon splitting on the charm spectra: for this purpose we compare results for the production of charm quarks in the multiple scattering mode to that in which multiple scatterings as well as fragmentation of partons is enabled (Fig. 15). We observe a substantially larger production of charm at lower $p_{T}$ and a suppressed production of charm at larger $p_{T}$. The increase in the production of charm at lower $p_{T}$ is a consequence of the large increase in the population of low momentum partons due to gluon multiplication. In addition, high momentum charm quarks loose energy via gluon splitting and multiple scatterings with lower momentum partons.

The rapidity distribution of charm shows an enhanced production at central rapidities: while the overall increase in production of charm is by a factor of about four, at central rapidity the rise is by a factor of six.

Comparing data on final state hadrons, such as D-mesons, to the PCM is inherently challenging, due to the limitation of the PCM to parton scattering above a momentum cut-off $p_{T}^{\rm min}$ as well as due to the systematic uncertainties inherent in the modeling of hadronization. A realistic hadronization model, accounting for bulk hadronization, parton recombination and hadronization via fragmentation is currently out of reach for this particular PCM implementation. Nevertheless, we shall attempt a series of schematic comparisons, by assuming a simple hadronization scheme for D-mesons, based solely on the fragmentation of c-quarks. Starting with a fragmentation function for the production of $D$ mesons from charm quarks as $\delta(1-z)$ where $z$ is the fraction of momentum of the charm quark carried by the $D$ meson, we can assume that the $p_{T}$ of the D-mesons will be very similar to the $p_{T}$ of the charm quarks themselves. This ansatz has widely been used in the literature (see e.g., Ref. Vogt:1993kn ; Brodsky:1991dj ). Please note, however, that more accurate fragmentation functions would show a narrow peak around the ratio of the masses of the charm quark and $D$ mesons $M_{c}/M_{D}$ (Ref. Peterson:1982ak ). With this in mind we compare the results of our calculations for the $p_{T}$-spectrum of charm quarks at central rapidity, with the results for $D^{0}$ mesons obtained by the STAR Collaboration Adamczyk:2014uip for the most central collisions (Fig. 16) after adjusting it for the fragmentation ratio for the $c\rightarrow D^{0}$ fragmentation. We show our results only for $p_{T}>2$ GeV, since our calculations do not include hydrodynamic flow, which is known to modify the spectra at lower transverse momenta and we have used a $p_{T}^{\text{cut-off}}$ and $\mu_{0}$ to regularize the pQCD matrices and $P_{a\rightarrow c}$ fragmentations. We see that the shape for the calculated spectrum is quite similar to the experimental result above a $p_{T}$ of 2-3 GeV, where fragmentation should start to significantly contribute to charm quark hadronization. We do add that while our calculations are for zero impact parameter, the experimental data are for 0–10% central collisions which have about a 10% smaller number of collisions.

III.2.3 Nuclear Modification Factor for Charm Production

The nuclear modification factor is defined by:

[TABLE]

where $N_{\text{coll}}$ is the number of binary nucleon-nucleon collisions for a given centrality. $R_{\textrm{AA}}$ is a measure of modification of the production processes due to the presence of the medium produced in the wake of the nucleus-nucleus collisions and has been used extensively to study the phenomenon of jet-quenching.

As a first step we compare calculations for the rapidity integrated $R_{\textrm{AA}}$ for different interaction assumptions (Fig. 17). Unsurprisingly, we find that the results for the eikonal approximation are consistent with one, as one would expect from these calculations. The results for multiple scattering suggests a considerable increase in the production of charm quarks at intermediate $p_{T}$ beyond that which one would have expected on the basis of a mere superposition of $pp$ collisions and is consistent with our findings discussed in previous sections.

The results for the full calculation, including multiple scatterings as well as parton fragmentations (mostly radiation of gluons) provide some interesting findings: in the low $p_{T}$ domain, we observe an enhancement of charm production of more than 50% beyond what one would expect from mere superposition of $pp$ collisions. Noting that the $p_{T}$ of charm quarks can get strongly affected by hydrodynamic flow Cao:2014fna , these result are expected to change if a more realistic bulk medium evolution is included in the calculation. At large $p_{T}$ we see a substantial suppression of charm production. Thus we find that while the inclusion of mere multiple scattering can lead to increase in production of charm at intermediate $p_{T}$, the inclusion of gluon emission rapidly depletes the momenta of the partons. Please note that this significant energy-loss is due to vacuum radiation matrix elements, since the PCM does not contain any medium modified matrix elements. Please also note that the presence of hydrodynamic flow would significantly alter the shape of $R_{AA}$ at low transverse momenta.

In Fig 18 we show a comparison of the nuclear modification factor $R_{\text{AA}}$ for central collisions at central rapidity. Interestingly we find a significant suppression, even down to intermediate $p_{T}$, indicating that charm interactions in the early non-equilibrium phase of the reaction play a significant role. Clearly, the PCM calculation of $R_{\text{AA}}$ suppression is less than observed in the data, i.e. it cannot account for the full suppression observed. Also, the radial flow bump at low $p_{T}$ is missing in the PCM calculation. Yet our results indicate that a comprehensive calculation of heavy flavour energy loss would need to account for both, the early time dynamics as described by the PCM as well as the later time dynamics as described by some form of hydro-kinetic transport Rapp:2009my ; Gossiaux:2010yx ; Cao:2014fna ; Das:2015ana ; Song:2015ykw . Please note that the experimental data are for 0–10% centrality while our calculations are for impact parameter equal to zero.

IV Summary

We have studied the production and scattering of heavy quarks in parton cascade model and presented results for production of charm quarks in $pp$ and $AuAu$ collisions at $\sqrt{s_{\textrm{NN}}}$ = 200 GeV. We have verified the PCM implementation by comparing an eikonal approximation implementation against an independent calculation using the formalism of minijets. The effect of charm-medium interactions is studied by successively opening up reaction channels, starting with only primary-primary collisions at first, then allowing for multiple scatterings and finally switching on multiple scatterings along with parton splittings. Interestingly, while the last named calculations lead to a very large increase in multiple collisions, the charm production does not increase by the same factor, mainly due to their production threshold as the initial light partons lose energy while radiating gluons. The space-time distribution of the production vertices of charm quarks provide valuable insight into the production process. The transverse momentum spectrum obtained using the complete calculation is found to have a shape similar to the experimental results, confirming the importance of radiative energy loss suffered by heavy quarks at intermediate and large $p_{T}$. Results for different systems, energies and cetralities would provide a very valuable confirmation of our procedure and findings. The considerable suppression of production of charm quarks having large transverse momenta due to the multiple collisions and radiations seen in the present work suggests that a significant contribution to jet-quenching could also arise from the same mechanism. This is under investigation.

Acknowledgments

DKS gratefully acknowledges the support by the Department of Atomic Energy, India. This research was supported in part by the ExtreMe Matter Institute EMMI at the GSI Helmholtzzentrum für Schwerionenforschung, Darmstadt, Germany. We thankfully acknowledge the High Performance Computing Facility of Variable Energy Cyclotron Centre Kolkata for all the help during these calculations performed over several months. SAB acknowledges support by US Department of Energy grant DE-FG02-05ER41367.