Thermal Charm and Charmonium Production in Quark Gluon Plasma

We study the effect of thermal charm production on charmonium regeneration in high energy nuclear collisions. By solving the kinetic equations for charm quark and charmonium distributions in Pb+Pb collisions, we calculate the global and differential nuclear modification factors $R_{AA}(N_{part})$ and $R{AA}(p_t)$ for $J/\Psi$s. Due to the thermal charm production in hot medium, the charmonium production source changes from the initially created charm quarks at SPS, RHIC and LHC to the thermally produced charm quarks at Future Circular Collider (FCC), and the $J/\Psi$ suppression ($R_{AA}<1$) observed so far will be replaced by a strong enhancement ($R_{AA}>1$) at FCC at low transverse momentum.

We study the effect of thermal charm production on charmonium regeneration in high energy nuclear collisions. By solving the kinetic equations for charm quark and charmonium distributions in Pb+Pb collisions, we calculate the global and differential nuclear modification factors R AA (N part ) and R AA (p t ) for J/ψs. Due to the thermal charm production in hot medium, the charmonium production source changes from the initially created charm quarks at SPS, RHIC and LHC to the thermally produced charm quarks at Future Circular Collider (FCC), and the J/ψ suppression (R AA < 1) observed so far will be replaced by a strong enhancement (R AA > 1) at FCC at low transverse momentum. Statistical Quantum Chromodynamics (QCD) predicts that, a strongly interacting matter will undergo a deconfinement phase transition from hadron matter to quark matter at finite temperature and density. It is expected that, this new state of matter, the so-called quark gluon plasma (QGP), can be created by liberating quarks and gluons from hadrons through high energy nuclear collisions. Since the QGP can only exist in the initial period and cannot be directly observed in the final state of the collisions, one needs sensitive probes to demonstrate the formation of this new state. J/ψ suppression has long been considered as such a probe since the original work of Matsui and Satz [1], and many progresses have been achieved both experimentally and theoretically, see for instance the recent review paper [2,3]. While the charmonium production mechanism changes from initial production at SPS energy [4][5][6] to initial production plus regeneration at RHIC and LHC energies [7][8][9][10][11][12][13][14], the charm quarks are all from the initial production.
Recently, the Future Circular Collider (FCC) at CERN is proposed to push the energy frontier beyond LHC, which includes the plan of Pb+Pb collision at √ s NN = 39 TeV [15].
What would we expect about the charmonium production at this new energy regime? Since a much more hot medium will emerge at FCC, gluons and light quarks inside the medium would be more energetic and denser. Therefore, the thermal production of charm quarks via gluon fusion and quark and anti-quark annihilation may have a sizeable effect on charmonium regeneration. For the in-medium charm quark production, there are already many studies, by considering leading order [16][17][18] and including next to leading order [19] QCD processes. Taking into account the quadratic dependence of the charmonium regeneration on charm quark density, we expect that, the extra increase of charm quark pairs via the thermal production in QGP will obviously enhance the charmonium yield at FCC. Since the very hot medium can eat up almost all the initially produced charmonia, the regeneration becomes the only source of the finally observed soft charmonia. This makes J/ψ more effective to probe the medium properties. In this paper, we focus on the effect of thermal charm production on charmonium production in heavy ion collisions at LHC and FCC energies. The full information of charm quarks in medium is contained in their distribution function f c (t, x, p) in phase space, its momentum integration is the number density n c (t, x) = d 3 p/(2π) 3 f c (t, x, p). When charm quarks approach kinetic equilibrium with the medium significantly fast in high energy nuclear collisions, only the evolution of the chemical abundance needs to be considered. By integrating out the charm quark momentum assuming thermal distribution, the Boltzmann equation for f c becomes the rate equation for n c , where n µ c = n c (1, v) is the charm current with medium velocity v, and the gain and loss terms r gain and r loss on the right hand side are respectively the charm quark production and annihilation rates inside QGP. The rates can be calculated through perturbative QCD.
It is convenient to use the Lorentz covariant variables η = 1/2 ln((t + z)/(t − z)) and τ = √ t 2 − z 2 to replace time t and longitudinal coordinate z. By using ∂ t = cosh η∂ τ − sinh η/τ∂ η and ∂ z = − sinh η∂ τ +τ cosh η∂ η , the rate equation is expressed as (2) with transverse medium velocity v T . From the experimentally observed large quench factor [20,21] and elliptic flow [22,23] for open charm mesons at RHIC and LHC energies, we assume that charm quarks are kinetically equilibrated with the medium during the whole evolution [24]. Therefore, the longitudinal motion of charm quarks will be consistent with the medium's Bjorken expansion [25] in mid rapidity region. To make this explicitly, we set n c = ρ c /τ with ρ c (τ, x T ) being the charm quark number density in transverse plane and controlled by the reduced rate equation at mid rapidity, Taking into account the nuclear geometry and nuclear shadowing effect on parton distributions in nuclei, the initial con-dition at time τ 0 for the rate equation in nuclear collisions at fixed impact parameter b can be written as, where dσ cc /dη is the rapidity distribution of charm quark production cross section in p+p collisions, and T A and T B are the thickness functions at transverse coordinate x T and x T − b for the two colliding nuclei. From the FONLL [26] simulation, we extract dσ cc /dη| η=0 = 0.7, 1.0 and 2.5 mb at colliding energy √ s NN = 2.76, 5.5 and 39 TeV. Considering gluon fusion as the dominant process of charm quark production in high energy p+p collisions, the cross section is multiplied by the shadowing modification factors R g (x 1 , x T ) for one gluon and R g (x 2 , x T − b) for the other to include the nuclear shadowing effect in A+B collisions, where x 1 and x 2 are the averaged gluon momentum fractions which can be es- quark transverse momentum [26]. The space dependence of the shadowing factor is taken as a linearized form [27], is the nuclear geometry factor, and R g (x) is the space independent shadowing factor which is taken from the EKS98 model [28] in the present study . Now we turn to the loss and gain terms for thermal charm production in medium. The general Lorentz-invariant form for a 2 → 2 process with initial particles 1 and 2 is where ν is the number of identical particles in the initial state, F 12 = (p 1 · p 2 ) 2 − m 2 1 m 2 2 is the kinetic flux factor, σ 12 is the corresponding cross section, and f 1,2 are the distribution functions of the initial particles 1 and 2. For the charm pair production, we take the NLO cross section derived by Mangano, Nason, and Ridolfi (MNR-NLO) [29,30] and choose the QCD running coupling constant α s at the renormalization scale µ = m c . For initial gluons and light quarks, we take their thermal masses [31] is the number of colors and N f = 3 the number of flavors. Most of the nonperturbative dynamics is contained in the temperature dependent coupling g(T ). We take it from Ref. [32] obtained by fitting the lattice QCD thermodynamics. Using fully thermal and chemical equilibrium distributions f eq g and f eq q for gluons and light quarks, one obtains the in-medium charm pair production rate r gain = r gg→cc(g) + r qq→cc(g) .
When charm quarks are dense enough, their annihilation starts to reduce the charm quark population. From the detailed balance between the production and annihilation processes, we can get the annihilation cross section. For the annihilation rate r loss = r cc(g)→gg + r cc(g)→qq , we need the charm and anti-charm quark distributions f c and fc. Taking into account the strong interactions between charm quarks and the medium, we assume charm quark thermalization and take f c(c) = n c(c) /n eq c(c) f eq c(c) , where f eq c(c) are the kinetically thermalized distribution functions for charm and anti-charm quarks and n c(c) /n eq c(c) = d 3 p/(2π) 3 3 f eq c(c) (p) are just the charm quark fugacity factors γ c(c) which measure how far the distributions are from the chemical equilibrium.
For the 3 → 2 annihilation processes which are at next to leading order, there is one more momentum integration for the initial gluon in the loss rate. We can effectively absorb this into the annihilation cross section and still take the same form (5). This can also be examined in terms of thermally averaged cross section [19].
The local temperature T (x) and fluid velocity u µ (x) appeared in the gluon, light quark and heavy quark distribution functions are controlled by the ideal hydrodynamics [33], where T µν is the energy-momentum tensor of the system, and we have neglected the baryon current at LHC and FCC energies. The initial time for the fluid evolution created in Pb+Pb collisions is chosen as τ 0 = 0.6 fm/c at colliding energy √ s NN = 2.76 and 5.5 TeV and 0.3 fm/c at 39 TeV.
The equation of state for the medium is needed to close the above hydrodynamical equations. In this work we use the Lattice QCD parametrization "s95p-v1" [34], which matches the recent Lattice QCD [35] simulation by HotQCD collaboration at high temperature (above T c ) to the Hadron Resnance Gas (HRG) at low temperature with a smooth crossover transition at temperature T c = 155 MeV.
The final charged multiplicity is used to determine the initial entropy density of the fluid. For Pb+Pb collisions at √ s NN = 2.76 TeV we have dN ch /dη = 1600 at mid rapidity from the ALICE collaboration [36]. At higher colliding energy, we parameterize dN ch /dη as (7) which leads to dN ch /dη = 2000 at √ s NN = 5.5 TeV and 3700 at 39 TeV. Assuming that the entropy is directly related to the final charged multiplicity, its initial configuration is usually estimated through the two-component model [37,38]. We take the mixing ratio α between the number density of participants n part and the number density of binary collisions n coll as 0.14, 0.16 and 0.2 for Pb+Pb collisions at √ s NN = 2.76,

and 39
TeV. Being different from usual treatment, we consider here the shadowing effect on the contribution from the hard processes, where the space averaged shadowing factors R A,B are estimated from the EKS98 model [28]. The p+p inelastic cross section σ NN is taken as 62, 72 and 100 mb at √ s NN = 2.76, 5.5 and 39 TeV.
For the freeze out of the fluid, we assume that the medium maintains chemical and thermal equilibrium until the energy density of the system dropping down to 60 MeV/fm 3 when the interaction among hadrons ceases and their momentum distributions are fixed. By solving the hydrodynamic equations (6) for the local temperature and medium velocity, and then substituting them into the rate equation (3), we can numerically obtain the time evolution of the charm quark yield in the medium.  duction inside QGP, the total charm quark number (solid lines) increases in the early period of the QGP phase and then decreases due to the charm and anti-charm quark annihilation in the later stage. At √ s NN = 2.76 TeV, the thermal production is not so strong and almost canceled by the inverse annihilation at the end of the QGP phase. Therefore, the consideration of thermal charm production in medium is negligible. However, at √ s NN = 5.5 TeV, the QGP medium becomes hotter and survives longer, and the thermal charm production becomes important and leads to a visible increase of the total charm quark number. At the end, the thermal production induced increase is still 15%. This sizeable enhancement agrees with the previous calculations like [19] using a fire-cylinder model for the medium evolution. At the FCC energy √ s NN = 39 TeV, the initial temperature of the medium in central collisions is T 0 = 840 MeV, and the densely populated gluons at so high temperature in QGP makes the thermal charm production more efficient. In this case, the total charm quark number increases exponentially in the very beginning of the QGP medium and keeps as almost a constant in the later evolution. At the end of the QGP, the enhancement reaches 50%. Due to the heavy mass, charmonium evolution in QGP can be controlled by classical transport approaches. At mid rapidity in heavy ion collisions, the Boltzmann equation can be used to describe the charmonium distribution function f Ψ (x T , p T , τ|b) in transverse phase space (x T , p T ) at time τ and fixed impact parameter b, where Ψ stands for J/ψ, χ c and ψ ′ to include the feeddown [39] contribution from excited states χ c and ψ ′ to ground state J/ψ, v Ψ = p T / p 2 T + m 2 Ψ is the Ψ transverse velocity. On the right hand side, the loss and gain terms α(x T , p T , τ|b) and β(x T , p T , τ|b) represent charmonium dissociation and regeneration in the hot medium. The elastic scattering is neglected here, since the Ψ masses are much larger than the typical medium temperature. Considering the gluon dissociation g + Ψ → c +c as the main dissociation process, the loss term is given by [6,10] where E g is the gluon energy and F gΨ = (pk) 2 − m 2 Ψ m 2 g = pk the flux factor. For the dissociation cross section σ gΨ (p, k, 0) in vacuum at T = 0, one can use the result from the operator production expansion (OPE) method with a perturbative Coulomb interaction [40][41][42][43]. When going to high temperature medium where the charmonium states become loosely bound, the OPE method is no longer valid. We estimate the temperature effect by taking into account the geometrical relation between the cross sections in medium and in vacuum [44,45], where r 2 Ψ (T ) is the averaged charmonium radius square which can be calculated by the potential model [46] with lattice simulated heavy quark potential [47] at finite temperature. The divergence of r 2 Ψ (T ) defines the charmonium dissociation temperature T d above which the charmonium state Ψ will melt due to the color screening.
Using detailed balance between the dissociation and regeneration processes [8], we can get the regeneration cross section. To obtain the regeneration rate β, we need the charm quark phase space distribution f c(c) which includes the contribution from the thermal charm production described above.
The initial condition for the charmonium transport equation is obtained [48] by taking the geometric superposition of free p+p collisions along with the modifications from the cold nuclear matter effects like Cronin effect [49] and nuclear shadowing effect [50]. Fig.2 shows the transverse momentum integrated (0 < p t < 30 GeV) J/ψ nuclear modification factor R AA as a function of the number of participant nucleons N part at mid rapidity in Pb+Pb collisions with colliding energy √ s NN = 2.76, 5.5 and 39 TeV. The thermal production of charm quarks does not change the initial charmonium production before the medium formation but enhances the charmonium regeneration inside the hot medium. At the current LHC energy √ s NN = 2.76 TeV, the weak thermal charm production slightly enhances the charmonium regeneration and in turn leads to a small J/ψ enhancement even in very central collisions, see the difference between the calculations with (solid line) and without (dashed line) thermal production shown in the upper panel of Fig.2. Since the degree of the nuclear shadowing is still an open question and there is not yet precise calculation of the shadowing factor, we show, as a comparison, also the calculation with thermal production but without shadowing in Fig.2 (dotted lines). Since the shadowing effect suppresses both the initial production and regeneration, the dotted line is always above the solid line. From the comparison with the ALICE data [51], the calculation without shadowing effect looks more close to the data. At √ s NN = 5.5 TeV, there exists a wide plateau of R AA in semi-central and central collisions, when the thermal charm production is excluded, similar to the structure at RHIC energy [44]. The thermal charm production which increases with centrality breaks the plateau structure, and the charmonium yield goes up sizeably with the number of participants. Considering the fact that the charmonium regeneration is proportional to the charm quark number square, a small change in the charm quark number by thermal production can lead to a remarkable charmonium enhancement. For instance, in very central collisions at √ s NN = 5.5 TeV, the charm quark enhancement is about 15%, but the charmonium enhancement becomes (0.37 − 0.26)/0.26 ∼ 40%. Very different from the case at colliding energies √ s NN = 2.76 and 5.5 TeV, the thermal charm production plays a dominant role in charmonium production at FCC energy √ s NN =39 TeV. Since the fireball created at FCC is extremely hot, the size is so large, and the life time is so long, the initially produced charmonia are almost all eaten up by the hot medium, and the charmonium production is controlled by the regener- ation, and therefore the contribution from the thermal charm production to the charmonium yield is largely amplified, see the bottom panel of Fig.2. In very central collisions, the thermal charm production makes the J/ψ R AA increase from 0.2 to 0.75, being enhanced by a factor of 3! When the shadowing is switched off, the R AA can even be larger than one. The other significant signal of the thermal charm production is the deep valley of R AA located at N pant ∼ 30, due to the competition between the initial production in peripheral collisions and strong thermal charm production in semi-central and central collisions.
The charmonium transverse momentum distribution is more sensitive to the production mechanism and the medium properties. The two charmonium production mechanisms, namely the initial production and the later regeneration, play roles in different p t regions. For the initially produced charmonia, the low p t part is almost all absorbed by the hot medium, while the high p t part can survive the medium due to the leakage effect [52,53]. Since the regeneration process happens in the later stage of the medium evolution, the regenerated charmonia carry low energy and mainly distribute in the low p t region. Therefore, the thermal charm production which contributes only to the regeneration will enhance the charmo-nium yield in low p t region. The p t dependence of the J/ψ R AA at mid rapidity is shown in Fig.3 for central Pb+Pb collisions (b=0) at different colliding energy. At the current LHC energy √ s NN = 2.76 TeV, the high p t region (p t > 5 GeV/c) is dominated by the initial production, and the regeneration sourced from those initially produced charm quarks controls the low p t region. The thermal charm production leads to a remarkable enhancement at very low p t , see the comparison between the dashed and solid lines in the top panel of Fig.3. Without considering the shadowing effect, the thermal production induced enhancement becomes more remarkable.
With increasing colliding energy, the initially produced charmonia are more suppressed by the hotter medium, which results in a smaller R AA at high p t . On the other hand, only those regenerated charmonia in the temperature region T c < T < T d can survive in the QGP phase, those regenerated charmonia above the dissociation temperature T d will be immediately dissociated by the hot medium. This is the reason why the R AA at low p t is very small (∼ 0.25) at FCC energy when the thermal charm production is excluded, see dashed lines in Fig.3. However, when the thermal charm production is turned on, the J/ψ yield at low p t goes up monotonously with increasing colliding energy. Especially, the J/ψ suppression (R AA < 1) at √ s NN = 2.76 and 5.5 TeV becomes enhancement (R AA > 1) at the FCC energy. When the shadowing is switched off, the R AA can even reach 1.6, the enhancement by the thermal charm production is tremendous, see the dotted line in the bottom panel.
In summary, the effect of thermal charm production in QGP phase on the charmonium production in high energy nuclear collisions at LHC and FCC energies is investigated. There might also be minor contributions from pre-thermal equilibrium stage for thermal charm production, here in this exploratory work we neglect it assuming very fast thermalization for the whole system. We calculated the global and differential nuclear modification factors R AA (N part ) and R AA (p t ) for J/ψs, by solving the kinetic equations for charm quarks and charmonia in the QGP phase. While the thermal charm production is still weak at the current LHC energy √ s NN = 2.76 TeV, it becomes sizeable at √ s NN = 5.5 TeV and significant at the FCC energy. As a consequence, the source for charmonium production changes from initially created charm quarks at SPS, RHIC and LHC to thermally produced charm quarks at FCC, and the charmonium production in Pb+Pb collisions at FCC is completely controlled by the thermal charm production. This is manifested in the following three aspects: 1) The J/ψ yield in central collisions is enhanced by a factor of 4, 2) There appears a deep valley of global R AA located at peripheral collisions, and 3) the J/ψ suppression (R AA < 1) at low transverse momentum at SPS, RHIC and LHC energies becomes J/ψ enhancement (R AA > 1) at FCC energy.