Introduction

In nuclear–nucleus collisions at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC), a system of color deconfinement, known as Quark-Gluon Plasma (QGP), is expected to be produced [1, 2]. QGP is an almost perfect fluid in the sense that its shear viscosity to entropy density ratio is close to the minimum value [3, 4]. Because the volume and lifetime of the QGP may be much larger and longer than those given by the confinement scale ($\approx 1/{\Lambda }_{\text{Q}CD}$), it is expected that the collective flow of partons should arise and have observable consequences that provide information on the dynamics in the QGP [5]. Studies on the anisotropic flow of particle production have contributed significantly to the characterization of the system created at RHIC and LHC [5]. Using a general Fourier series decomposition of the azimuthal distribution of the emitted particles, $\frac{\text{d}N}{\text{d}\phi }\propto 1+2{\displaystyle \sum _{n=1}^{\infty }}{v}_n\cos [n(\phi -{\Psi }_n)]\mathrm{,}$ (1) anisotropic flow is quantified using coefficients vn and the corresponding symmetry planes Ψn [6], where φ represents the azimuthal angle of the emitted particles in the momentum space. Systematic measurements of the correlations and fluctuations of vn coefficients and Ψn have been performed in Refs. [7-13]. Comprehensive comparisons with theoretical model calculations provide critical information on the event average initial density distribution in the nuclear overlap region, as well as its connection to QGP dynamics [14-19].

Charm and bottom quarks are important probes for the QGP. They were mainly produced in the initial state of the hard-scattering process prior to QGP formation and survived throughout the evolution stage of the QGP. Quarkonia are bound states of heavy flavor quark–antiquark pairs, which are believed to be suppressed because of the color screening effect in the deconfined interior of the interaction region [20]. They offer a complementary way to study the interaction of the heavy flavor with the medium and, thus, independently shed light on the properties of QGP [21-23]. Dozens of measurements in this direction have been carried out, and rich physical information has been delivered [24-31], with a clear future direction [21].

To study the characteristics of the QGP and its evolution toward chemical and thermal equilibrium, it is necessary to derive transport equations for quarks and gluons when they are in equilibrium and pre-equilibrium [32]. There have been several studies on deriving relativistic abelian and non-abelian transport and constraint equations for the relativistic gauge covariant Wigner operator for color particles interacting via SU(N) gauge fields. This indicates that viscous hydrodynamic models can describe the transverse momentum spectra of light-flavor hadrons and distribution in the azimuthal angle for low p_T [33, 16]. On the other hand, transport models have been applied to study the evolution of QGP medium [34]. Owing to its relatively small cross section, a large event sample is required to study charm or even bottom quark transport in the QGP within hydrodynamics calculations or transport model simulations. In this study, we implement an additional heavy quark–antiquark pair production trigger in a multiphase transport (AMPT) model to study the dynamics of charm and bottom quarks in a dense light quark soup. This significantly improves the charm or bottom quark simulation efficiency, and it will be closer to the data at LHC energies because of a large charm quark production cross section [21].

Transport Models of Heavy Flavor

We employed the AMPT model with the string melting version (v2.265) to simulate the collective motion of heavy quarkonium in Pb–Pb collisions at LHC energies. The AMPT model is a hybrid model consisting of fluctuating initial conditions, parton elastic scattering, quark hadronization, and hadronic interaction [34]. It has been used extensively to describe the dynamics of heavy-ion collisions at RHIC and LHC energies in some recent publications [35-42]. For a detailed description of the AMPT model, please refer to Ref. [34] and Ref. [43] for a summary of recent developments.

In the AMPT model, the initial production of heavy quarks is handled using the HIJING two-component model [44]. It includes pair productions ($q+\overline{q}\to Q+\overline{Q}$, $g+g\to Q+\overline{Q}$) and gluon splitting ($g\to Q+\overline{Q}$) with a cross section in perturbative quantum chromodynamics (QCD) at leading order: $\frac{\text{d}{\sigma }^{Q\overline{Q}}}{\text{d}{p}_T^2\text{d}{y}_1\text{d}{y}_2}\text{}=\text{}K\text{}{\displaystyle \sum _{a\mathrm{,}b}}{x}_1{f}_a\left(x_1\mathrm{,}\text{}{\mu }_{\text{F}}^2\right)x_2{f}_b\left(x_2\mathrm{,}\text{}{\mu }_{\text{F}}^2\right)\times \frac{\text{d}{\sigma }^{ab\to Q\overline{Q}}}{\text{d}\widehat{t}}\mathrm{.}$ (2) Here, y₁ and y₂ are the rapidities of the two produced partons, the K factor aims to account for higher-order corrections of heavy quark production, a and b refer to the types of interacting partons in the initial state, x denotes the nucleon momentum fraction taken by the interaction parton, μ_F represents the factorization scale, fa is the parton distribution function of parton type a in a nucleon, and ${\sigma }^{ab\to Q\overline{Q}}$ is the cross section for parton types a and b to produce the heavy quark pair [44]. Because the cross section of heavy quark production, especially for the bottom quark, is several orders of magnitude smaller than that of light quarks in heavy-ion collisions at LHC energies, it calls for a trigger for heavy quark production to increase the simulation efficiency. The algorithm was originally developed by Wang and Gyulassy [44] and was extended to study the bottom quark dynamics in the QGP in this study. We refer to the details in Ref. [44] and Ref. [45] and do not repeat them here.

The final state parton interactions can be described by the equations of motion for the quark and gluon Wigner operators [46, 47], which are closely connected with the classical distribution function. For example, the relativistic gauge covariant Wigner operator for spin-$\frac{1}{2}$ particles in coordinate and momentum space [46]: $\widehat{W}\left(x\mathrm{,}p\right)={\displaystyle \int }\frac{{\text{d}}^{4}y}{{\left(2\pi \right)}^4}{e}^{-ip\cdot y}\overline{\psi }\left(x\right)e^{\frac{1}{2}y\cdot D_x^{†}}\otimes e^{-\frac{1}{2}y\cdot D_x}\psi \left(x\right)\mathrm{,}$ (3) where $D^{\mu }\equiv {\partial }^{\mu }+ig{A}^{\mu }$ is the covariant derivative, and $\widehat{W}$ is a matrix in spinor (4×4 components) as well as in the color space (N×N components).

By using the properties of the link operators, one can derive the quantum transport equation for the QCD Wigner operator [46, 47]. Under semiclassical and Abelian dominance approximations, in the weak-field limit, the equations of motion become a set of equations that resemble classical transport equations [48]. This motivated the study of the final state interaction by solving the Boltzmann equations for quarks and gluons. If only two-body interactions are considered, the Boltzmann equations can be reduced to the following: $p^{\mu }{\partial }_{\mu }f(x\mathrm{,}p\mathrm{,}t)\propto {\displaystyle \int }\sigma f(x_1\mathrm{,}{p}_1\mathrm{,}t)f(x_2\mathrm{,}{p}_2\mathrm{,}t)\mathrm{,}$ (4) where f(x,p,t) represents the distribution function of the parton at time t in the phase space, and $\sigma$ is the cross section for the two-body scattering of partons. Zhang's parton cascade (ZPC) model was developed to study the final-state parton interactions, which solved the Boltzmann equations for partons using the cascade method [48]. In this model, two partons scatter when their closest distance is less than $\sqrt{\sigma /\pi }$. The scattering angle is determined by the differential parton–parton scattering cross section. Because the cascade prescription for solving the semiclassical transport equation does not assume any particular form of the phase space distribution, it can be used to study systems far from local thermal equilibrium, and it automatically takes into account the finite mean free path effect [49]. The current study focuses on heavy quark dynamics at LHC energies, where the bottom quark may be far from thermalization and the charm quark may be close to thermalization [27, 28, 31].

Results and Discussion

To follow the parton scattering history in the AMPT model, we define N_scat. as the number of two-body elastic scatterings suffered by a parton. The N_scat. distribution of the active partons (N_scat. > 0) is shown in Fig. 1. The distributions from the AMPT model with an enhanced heavy quark trigger are flatter than those in the normal AMPT calculation, which suggests a larger scattering possibility for partons with enhanced heavy quark triggers in the event. The effect is rather small for the charm quark scenario and more significant for the bottom quark scenario.

Fig. 1Full-screen View

AMPT simulations of: (a) normalized parton scattering distributions of active quarks; (b) mean transverse velocity of freeze-out quarks in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV and 5%–60% centrality, where normal, w/ c-trig., and w/ b-trig. represent normal AMPT, normal AMPT with triggering an additional charm quark–antiquark pair production, and normal AMPT with triggering an additional bottom quark–antiquark pair production, respectively.

Figure 1 also shows the approximate average transverse radial velocity $\langle {\beta }_{\perp }\rangle \equiv \langle {\widehat{r}}_{\perp }\cdot {\widehat{p}}_{\perp }\rangle$ of freeze-out partons with different flavor in the AMPT model as a function of N_scat., where ${\widehat{r}}_{\perp }$ and ${\widehat{p}}_{\perp }$ are the transverse radial position and momentum unit vectors. The ${\beta }_{\perp }$ is not exactly zero at N_scat.=0 because partons can form only after a finite formation time over which the displacement of a parton depends on its momentum [34]. One sees that ${\beta }_{\perp }$ of heavy quarks are larger in magnitude and almost independent of N_scat., whereas it is smaller in magnitude and drops for light quarks as N_scat. increases. This is understood because heavy quarks are formed earlier, and they are more likely to inherit the direction of mother partons, which will give them larger $\langle {\beta }_{\perp }\rangle$. Compared with light quarks, it is harder to push the heavy quarks away from their original trajectory, their $\langle {\beta }_{\perp }\rangle$ are independent of N_scat.. As the mid-rapidity hadron density of a Pb–Pb collision at $\sqrt{s_{\text{NN}}}=5.02$ TeV reaches as high as a few thousand, depending on the collision centrality [50], triggering on an additional heavy quark–antiquark pair production has a negligible effect on the hadron density, but may affect the evolution dynamics significantly. We use these trigger events to study the properties of QGP matter created in Pb–Pb collisions at LHC energies.

Figure 2 presents v₂, v₃ of a charged hadron in mid-central Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV from AMPT model. To do a fair comparison with experimental data [51], the results are calculated using the two-particle correlation method [52]. The result from the reaction plane method, as described in Eq. (1), has been checked and is found to be consistent with Fig. 2. Overall, the AMPT calculations with the normal, enhanced heavy quark trigger scenario reproduce the experimental data reasonably well. Larger elliptic flow is seen in the enhanced heavy quark trigger scenario. The enhancement is approximately 5% and p_T independent in AMPT with bottom quark trigger (cf. the bottom panel of Fig. 2). Including the discussion in Fig. 1 where the AMPT model with enhanced bottom quark–antiquark pair production presents a larger scattering possibility for active quarks, which tends to generate larger elliptic flow and is seemingly closer to v₂ data as pT increases.

Fig. 2Full-screen View

Charged particles elliptic flow (v₂) and triangular flow (v₃) in mid-central Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV from the AMPT model with normal, enhanced heavy quark production trigger settings. The filled points represent experimental data for comparison [51].

Figure 3 highlights the v₂(p_T) of quarkonia in middle and forward rapidity of Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV. The uncertainties in calculations are determined by statistics in different centrality and different rapidity bins. Our results are consistent with the experimental data considering current uncertainties. Overall, AMPT simulations on v₂ of $\Upsilon$ are smaller than those of J/Ψ. One should note that in the current version of the AMPT calculation, regeneration of J/Ψ or dissociation of $\Upsilon$ are missing, only continued two-body elastic scattering is involved, which seems to generate the finite v₂ of quarkonia. The magnitude difference between $\Upsilon$ and J/Ψ may be due to the mass of different quarks associated with each calculation. We have learned that it is also one of the key parameters in hydrodynamics calculations [33, 16].

Fig. 3Full-screen View

Heavy flavor quarkonia v₂(p_T) of mid-rapidity and forward rapidity in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV from the AMPT model with enhanced heavy quark–antiquark pair trigger settings. The centralities were chosen to be the same as those in the experimental measurement [53, 31].

To study this further, we examined the differential rapidity (pseudorapidity) dependence on v₂, as shown in Fig. 4. Because v₂(η) data for a charged hadron in Pb–Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV has not yet been published, we chose the data [54] at $\sqrt{s_{\text{NN}}}=2.76$ TeV for comparison. Thus, the AMPT calculation for the charged hadron is performed on $\sqrt{s_{\text{NN}}}=2.76$ TeV accordingly. Our calculations with normal or enhanced heavy quark–antiquark pair production scenarios describe the data fairly well. Although the uncertainties presented in the measurement do not allow us to distinguish between different scenarios, the difference in the calculations is of a few per cent. The value of $\Upsilon$ v₂ is positive and at a maximum at mid-rapidity, and close to zero at higher rapidity. The calculation appears to agree with the data at forward rapidity [31]. The value of J/Ψ v₂ is larger than $\Upsilon$ and exhibits similar rapidity dependence. However, it underpredicts the data at forward rapidity. This deviation may originate from the missing dissociation and regeneration mechanisms for heavy quarkonium in the AMPT. Further studies are required to implement these effects in the AMPT model.

Fig. 4Full-screen View

AMPT model simulations on hadron v₂(η) in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76$ and 5.02 TeV. The centralities and p_T bins were chosen to be the same as the experimental data [31, 54]. Different model configurations were applied in the charge hadron calculation, whereas only the heavy quark–antiquark pair trigger settings were used for studying J/Ψ and $\Upsilon$.

The enhanced heavy quark–antiquark pair production event can also be applied to study the direct flow (v₁) of quarkonium. This represents a new probe for studying the strong electromagnetic field effect of hot QCD matter via v₁ measurement, particularly on the difference in v₁(η) between the particle and its antiparticle [56, 57, 58]. Measurements of the open charm meson v₁(η) show a difference in the slopes of D⁰ and ${\overline{D}}^0$, but with a different sign from RHIC to LHC [59, 55, 60]. The origin of this physics is still debated. Because the electromagnetic field effect has not been implemented in the AMPT model, we consider the v₁(η) of J/Ψ and $\Upsilon$ instead of the v₁ difference between the particle and antiparticle in our current framework. In this aspect, we used the coordinate system of the model with the x-axis along the impact parameter direction ($\overrightarrow{b}$) of each event in the transverse plane, with the center of the projective and target nuclei at the transverse coordinates (b/2, 0) and (-b/2, 0). The direction of the projectile momentum specifies the z-axis; consequently, the total angular momentum of the system is along the -y direction. Figure 5 shows our results for J/Ψ and $\Upsilon$ together with the charged hadron in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV. It can be seen that quarkonia present an order of magnitude larger v₁(η) than charged hadrons, presumably because of the early formation and larger mass. In comparison with the data [55], the v₁(η) of the charged hadron from the AMPT model is on the correct sign but an order of magnitude lower. Because the dynamics that generates v₁ is quite complicated, it may depend on the equation of state, mean-field potentials, particle scattering, shadowing from spectator nucleons, and tilt of the hot dense QCD matter in the x–z plane [58]. Our results rely on two-body parton elastic scattering, and naive quark hadronization can be understood. Future studies with more realistic dynamics are of interest.

Fig. 5Full-screen View

AMPT model calculations on hadron v₁(η) in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV. The charged hadron result is scaled by a factor of 10 for comparison with the experimental data [55].

Summary

In summary, we have implemented an additional heavy flavor quark–antiquark pair production trigger in the AMPT model to study their dynamics in hot and dense media created in high-energy heavy-ion collisions. In this study, we have focused on the heavy flavor quarkonium anisotropy flows in Pb–Pb collisions at LHC energies. The new implementation presents a systematic increase in the flow coefficients as a function of p_T and η in comparison with the result from the normal AMPT model and provides an efficient way to study quarkonium dynamics with relatively small event samples. It describes the charged particle, J/Ψ, and $\Upsilon$ v₂ reasonably well. Future efforts are required to include more realistic dynamics in the AMPT model, such as the dissociation and regeneration of heavy quarkonium production at LHC energies.