Scaling behaviors of heavy flavor meson suppression and flow in different nuclear collision systems at the LHC

We explore the system size dependence of heavy-quark-QGP interaction by studying the heavy flavor meson suppression and elliptic flow in Pb-Pb, Xe-Xe, Ar-Ar and O-O collisions at the LHC. The space-time evolution of the QGP is simulated using a (3+1)-dimensional viscous hydrodynamic model, while the heavy-quark-QGP interaction is described by an improved Langevin approach that includes both collisional and radiative energy loss inside a thermal medium. Within this framework, we provides a reasonable description of the $D$ meson suppression and flow coefficients in Pb-Pb collisions, as well as predictions for both $D$ and $B$ meson observables in other collision systems yet to be measured. We find a clear hierarchy for the heavy meson suppression with respect to the size of the colliding nuclei, while their elliptic flow coefficient relies on both the system size and the geometric anisotropy of the QGP. Sizable suppression and flow are predicted for both $D$ and $B$ mesons in O-O collisions, which serve as a crucial bridge of jet quenching between large and small collision systems. Scaling behaviors between different collision systems are shown for heavy meson suppression factor and the bulk-eccentricity-rescaled heavy meson elliptic flow as functions of the number of participant nucleons in heavy-ion collisions.


I. INTRODUCTION
Quantum Chromodynamics (QCD) predicts that at extremely high temperature and density, nuclear matter transits from the hadron state to a color-deconfined state, known as the quark-gluon plasma (QGP) [1]. This has been supported by a large number of experimental evidences from relativistic heavy-ion collisions at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) [2], among which anisotropic flow and jet quenching are considered as the two most important signatures of the formation of QGP. At low transverse momentum (p T ), hadrons emitted from the QGP exhibit strong anisotropy in their azimuthal angular distributions [3][4][5], which has been successfully described by relativistic hydrodynamic models [6][7][8][9][10]. The small values of the specific shear viscosity extracted from hydrodynamic calculations [11,12] clearly shows the strongly-coupled nature of this fluid-like QGP. At high p T , hadrons and reconstructed jets emanating from initial hard scatterings are significantly quenched after traversing the QGP medium [13][14][15][16][17][18]. The extracted large values of the jet quenching parameterq [19,20] indicate the quark-gluon degrees of freedom inside the dense nuclear matter.
Over the past few years, large anisotropic flows have been observed in small collision systems as well, such as deuteron-gold (d-Au) collisions at RHIC and proton-lead (p-Pb) collisions at the LHC [21][22][23][24]. Interestingly, these flow coefficients can also be successfully described by hydrodynamic calculations, implying the possible formation of mini-QGP in such small systems [25][26][27][28][29][30][31][32][33]. On the other hand, the other strong evidence of QGP, jet quenching, has not been observed so far. For instance, despite the large elliptic flow coefficient (v 2 ) of D mesons in p-Pb collisions [34], their nuclear modification factor (R pA ) is found to be consistent with unity [35]. This has triggered hot debates on whether the collectivity observed in small systems originates from final-state QGP effects or from initial-state gluon saturation effects [36][37][38][39][40][41][42][43]. One possible way of disentangling the initial-state and finalstate contributions to jet observables is to scan the jet quenching effect across various sizes of nuclear collision systems [44,45]. This would bridge the gap between large and small systems and may hopefully help identify the boundary across which QGP disappears. Along this direction, several theoretical efforts have been recently devoted to explore the nuclear modification effects on high p T hadrons in systems smaller than Pb-Pb collisions at the LHC energies [46][47][48][49], and how parton energy loss depends on the size of collision systems [50,51].
In this work, we aim at using heavy quarks to probe QGP with different sizes. Our state-of-the-art Langevinhydrodynamics framework is applied to calculate the nuclear modification factor (R AA ) and the elliptic flow coefficient (v 2 ) of heavy flavor mesons across Pb-Pb, Xe-Xe, Ar-Ar and O-O collisions. The (3+1)-dimensional viscous hydrodynamic model CLVisc [77][78][79][80] is adopted for simulating the realistic QGP profiles produced in these collision systems, while the improved Langevin approach [81,82] is used for describing both elastic and inelastic scatterings of heavy quarks through the QGP medium. Hadronization plays an important role in studying heavy flavor dynamics in heavy-ion collisions. In this work, heavy quarks exiting the QGP are converted to heavy flavor hadrons using our advanced coalescencefragmentation model [71], which has successfully predicted the heavy flavor hadron chemistry at RHIC and the LHC. These sophisticated models on heavy quark energy loss and hadronization are necessary for a more quantitative comparison to the experimental measurements. Within this framework, we predict both D and B meson R AA and v 2 for different collision systems and collision centralities, from which we investigate the hierarchy of heavy quark energy loss and its momentum anisotropy with respect to the medium size and geometric anisotropy. In particular, the scaling behaviors of heavy flavor meson observables between different collision systems are explored. We find that both R AA and the bulkeccentricity-rescaled elliptic flow (v 2 /ε 2 ) of heavy flavor mesons scale with the number of participant nucleons (N part ). These findings help to disentangle the effects of the overall intensity of medium modification and its geometric asymmetry on jet quenching observables, which can be tested by future measurements. This paper is organized as follows. In Sec. II, we present the CLVisc hydrodynamic model that we use to generate the QGP profiles produced in relativistic heavyion collisions. In Sec. III, we review our Langevin approach that describes both collisional and radiative energy loss of heavy quarks inside QGP. In Sec. IV, our numerical results on D and B meson R AA and v 2 are presented for different centrality regions across Pb-Pb, Xe-Xe, Ar-Ar and O-O collision systems, from which the hierarchy and scaling behaviors of these observables with respect to the system size, medium geometry and heavy quark mass will be investigated in detail. The conclusion of this study is presented in Sec. V.

II. HYDRODYNAMIC SIMULATION OF MEDIUM PROFILES
In this study, the dynamical evolution of the QGP medium is provided by the (3+1)-dimensional CLVisc hydrodynamic model [78,79]. The full initial entropy density distribution S(τ 0 , x, y, η s ) is constructed by folding the smooth entropy density s(x, y) in the transverse plane and the parametrized envelope function H(η s ) in the longitudinal direction at the initial proper time τ 0 , where K is a scale factor which can be adjusted from the final charged hadron spectra in the most central collisions [83,84]. The entropy density s(x, y) is generated by the Trento initial condition [85], in which the positions of nucleons within nucleus are first sampled using the Woods-Saxon distribution, in which ρ 0 denotes the nuclear density at the nucleus center, d is the surface thickness parameter, and R(θ) = R 0 (1 + β 2 Y 20 (θ) + β 4 Y 40 (θ)) is the nuclear radius with spherical harmonic functions Y nl (θ). Here, β 2 , β 4 and w parameters control the deviations from a spherical nucleus. The local entropy density s(x, y) can be then constructed from the generalized mean of the nuclear matter thickness function T A (x, y) and T B (x, y) as follows: where the thickness functions are obtained by summing over the Gaussian smearing functions (with width 0.5 fm) of the participant nucleons inside the two colliding nuclei (A and B). In this work, we choose p = 0 that corresponds to the IP-Plasma-model-like or EKRT-model-like entropy deposition. The envelope functions H(η s ) are chosen to describe the longitudinal profile [78], where we use η 0 = 2.23, σ ηs = 1.8 for Xe and η 0 = 1.7, σ ηs = 2.0 for other nuclei in this study. For each centrality interval of each collision system, we average over 5000 Trento events to get a smooth initial entropy distribution as our hydrodynamic input at the initial proper time τ 0 = 0.6 fm.
In the framework of the CLVisc hydrodynamic model [78], the equation of motion for the energymomentum tensor T µν and the dissipative equation for the shear stress tensor π µν are solved with the partial chemical equilibrium equation of state s95p-pce in the Milne coordinate using the Kurganov-Tadmor (KT) algorithm: where σ µν is the symmetric shear tensor and θ is the expansion rate. We set the specific shear viscosity η v /s = 0.16 and the relaxation time τ π = 3η v /(T s). After hydrodynamic evolution, the QGP is converted to hadrons via the Cooper-Frye formula with the switching temperature set as T sw = 137 MeV. With above setups, our hydrodynamic calculation provides reasonable descriptions of the soft hadron spectra in Pb-Pb collisions and Xe-Xe collisions. The QGP profiles of Ar-Ar and O-O collisions should be viewed as predictions at this moment.

III. HEAVY QUARK EVOLUTION INSIDE QGP
The time evolution of heavy quarks through the QGP medium is described using the modified Langevin equation [81] that simultaneously includes quasi-elastic scattering and medium-induced gluon bremsstrahlung processes of heavy quarks inside a thermal medium: In the above equation, the first two terms on the righthand side follow the classical Langevin equation, denoting the drag force and thermal random force experienced by a heavy quark while it frequently scatters with the constituents of a thermal medium. The thermal force ξ is assumed to be independent of the heavy quark momentum. Its strength is quantified by the correlation function of a white noise In addition to the drag and diffusion from the multiple scattering process, the effects of medium-induced gluon radiation is introduced into Eq. (7) as a recoil force f g = −d p g /dt exerted on heavy quarks while they emit gluons with momentum p g . The probability of gluon radiation during a time interval (t, t+∆t) is evaluated using the average number of emitted gluons during this time interval: In the calculation, we choose a sufficiently small ∆t to guarantee N g (t, ∆t) < 1, so that this average number can be utilized as a probability. In Eq. (8), the mediuminduced gluon spectrum is adopted from the higher-twist energy loss calculation [87][88][89]: in which x is the fractional energy taken by the emitted gluon from its parent heavy quark, k ⊥ is the transverse momentum of the gluon, α s is the strong coupling which runs with k 2 ⊥ at the leading order, P (x) is the Q → Qg splitting function, and τ f = 2Ex(1 − x)/(k 2 ⊥ + x 2 M 2 ) denotes the splitting time with E and M being the energy and mass of heavy quarks respectively. Hereq is the gluon transport coefficient which can be related to the quark diffusion coefficient viaq = 2κC A /C F , where C A and C F are color factors of gluon and quark respectively. Note that in our modified Langevin model, there is only one free parameter which we choose as the dimensionless quantity D s (2πT ). It is adjusted as D s (2πT ) = 4 [76] to provide a reasonable description of the heavy flavor meson observables in heavy-ion collisions at the LHC.
When we sample the energy-momentum of the medium-induced gluons according to Eq. (9), a lower cutoff is implemented for the gluon energy at ω 0 = x 0 E = πT , below which gluon is not allowed to form. Due to the lack of the gluon absorption process in the current implementation, this cut-off helps mimic the balance between gluon emission and absorption processes around the thermal energy scale. We have verified that an approximate, though not exact, thermal equilibrium of heavy quarks can be achieved after a sufficiently long time of evolution inside a thermal medium [81].
Using this Langevin framework, we can simulate the heavy quark evolution through the QGP. The realistic QGP medium is generated by the CLVisc hydrodynamic model as described in Sec. II. Meanwhile, the initial heavy quarks are sampled using the binary collision vertices from the Monte-Carlo Glauber model for their position space, and the fixed-order-next-to-leadinglog (FONLL) calculation [90][91][92] convoluted with the CT14NLO [93] parton distribution function for their momentum space. Then heavy quarks are placed into our Langevin model for their subsequent interaction with the QGP medium, which we assume to commence at the initial time (τ 0 = 0.6 fm) of the hydrodynamic evolution.
After heavy quarks travel outside the QGP boundary, i.e., the local temperature of the medium drops below T c = 160 MeV, they are converted to heavy flavor hadrons via an advanced hybrid fragmentationcoalescence model [71]. In this model, the coalescence probability between heavy quarks and thermal light quarks are calculated according to the wavefunction overlap between the free-quark state and hadronic bound state. Both s and p-wave hadronic states are included in our calculation, which naturally cover the majority of heavy flavor hadron states observed in the Particle Data Group [94]. Based on this probability, heavy quarks that do not hadronize through coalescence are fragmented into heavy flavor hadrons via Pythia [95] simulation. The heavy flavor hadrons produced from both coalescence and fragmentation processes are then utilized for analyzing the final-state observables.
As discussed earlier, a constant D s (2πT ) = 4 is used in this work since our main focus is on the system size dependence of heavy quark energy loss and heavy flavor suppression and flow at the LHC. One may refer to our previous study [76] for a detailed analysis of the systematic uncertainties introduced by various model ingredients, such as the initial heavy quark spectrum, the starting time of heavy-quark-medium interaction, the medium profile in the pre-equilibrium state, and the temperature dependence of the heavy quark diffusion coefficient, etc. Compared to the linear Boltzmann transport model used in our earlier work [51], the Langevin approach is expected to be applicable to quasi-particles with large masses inside a thermal medium. However, it is easier to include the non-perturbative interaction between low energy heavy quarks and the QGP in the Langevin approach than in the perturbative-based Boltzmann calculation. Note that both models use the same method to implement the radiative energy loss of heavy quarks [51,76].

IV. HEAVY FLAVOR MESON SUPPRESSION AND FLOW
In this section, we present numerical results on the nuclear modifications of D and B mesons, and compare them between different collision systems (Pb-Pb, Xe-Xe, Ar-Ar and O-O) at the LHC energies. The two most frequently quoted heavy flavor observables -nuclear modification factor (R AA ) and elliptic flow coefficient (v 2 )are utilized to quantify features of heavy quark energy loss inside the QGP. In the present study, they are extracted as follows from the final-state energy-momentum information of the heavy flavor mesons: where N coll denotes the average number of binary collisions in a given centrality bin of a given nucleus-nucleus collision system, and . . . represents the average over the final-state heavy flavor mesons generated in our simulations. Smooth hydrodynamic profiles are used in this work, in which the x-z axes define the event plane while the x-y axes define the transverse plane of nuclear collisions. Within this setup, the azimuthal angle φ in Eq. (11) is measured with respect to the +x direction. Note that each smooth hydrodynamic profile is generated from an initial entropy distribution that has been averaged over 5000 Trento events, whose participant planes have been individually rotated to the x-z plane of our computational frame. Therefore, such smooth hydrodynamic profile has captured key features of event-by-event fluctuations in the initial state and serve as a good approximation of direct event-by-event simulations of QGP for studying heavy flavor observables. Implementing full event-by-event calculations can lead to stronger energy loss of heavy quarks [96] and larger v 2 [97] than using the smooth profiles, though such difference is expected to be within 10%. In this work, we use the event plane method Eq. (11) to evaluate the heavy meson v 2 , following the ALICE [98] and STAR [99] collaborations. Note that the correlation method has also been applied in STAR [99] and CMS [100] measurements. As shown by Ref. [99], these two methods produce similar v 2 for heavy mesons.
With these setups, we first present in Fig. 1  Note that the impact parameter averaged nuclear shadowing parametrization is used in Fig. 1. The dependence of nuclear shadowing on the impact parameter could introduce additional dependence of nuclear modification on centrality. This will be included in our future study.
Since the EPPS16 parametrization does not cover all nucleus species that we investigate in the present work, we choose to exclude this cold nuclear matter effect for the rest of our calculation in order to conduct an unbiased comparison between different collision systems. With the heavy quark diffusion coefficient set as D(2πT ) = 4, our results on the D meson R AA and v 2 for Pb-Pb collisions are consistent with the data from ALICE and CMS collaborations [98,100,102]. This helps confirm the satisfactory path-length dependence of parton energy loss embedded in our transport model.
Comparing different centrality classes in each panel of Fig. 1, one can observe a clear hierarchy in the heavy meson R AA , i.e., larger heavy quark energy loss in more central collisions leads to a smaller nuclear modification factor. However, this hierarchy does not hold for the heavy meson v 2 which depends on the competing effects between the amount of energy loss and the geometric anisotropy of the medium. The former is stronger in more central collisions, while the latter is larger in more peripheral collisions. Therefore, one usually observes a maximum for elliptic flow v 2 in semi-central/peripheral collisions (e.g. 30-40%). Comparing upper panels and lower panels, one can observe that D mesons have smaller R AA and larger v 2 than B mesons because charm quarks have much smaller mass than bottom quarks.
Within the same framework, we investigate the nuclear modification of heavy flavor mesons in smaller systems at the LHC. Results are presented in Fig. 2  Comparing different collision systems (from Fig. 1 to Fig. 4), a general conclusion can be drawn, i.e., parton energy loss becomes weaker inside a smaller collision system, as suggested by the gradually larger heavy flavor meson R AA and smaller v 2 within the same centrality class as we move from Pb-Pb, Xe-Xe, Ar-Ar to O-O collisions. Such system size dependence of jet quenching effects provides a crucial bridge of jet-medium interac-   tion between large and small collision systems.
It is interesting to note that even in the relatively small system produced by O-O collisions, considerable amount of energy loss effects are found for both charm and bottom quarks in the most central collisions -the corresponding D and B meson R AA 's are significantly smaller than unity while their v 2 's have sizable values. As one moves from central to peripheral collisions, heavy flavor meson R AA increases and approaches unity at high p T in peripheral collisions. The rise-and-fall structure of R AA at low p T region is due to the coalescence mechanism in heavy hadron formation in the presence of QGP medium, which converts low p T heavy quarks into intermediate p T heavy flavor mesons. For the heavy flavor meson v 2 , it first increases and then decreases as a function of centrality class due to the competing effects between parton energy loss and geometric anisotropy of the collision zone.
To have a more quantitative understanding of how heavy quark energy loss depends on the system size of QGP, we present the participant number (N part ) dependence of the R AA and v 2 of D mesons in Fig. 5 and B mesons in Fig. 6. In each panel of these two figures, we present the p T -integrated observable as a function of N part for different collision systems. The upper panels are for 5 < p T < 8 GeV and the lower for 8 < p T < 16 GeV. In the left panels, we observe a stronger nuclear modification of heavy mesons with a larger N part . As previously discussed, one can find clear nuclear modification of both D and B mesons even in the small-size O-O collisions as long as N part is not small. In addition, a scaling behavior of the nuclear modification factor with respect to N part can be seen in the left panels: the heavy flavor meson R AA in different collision systems follow the similar N part dependence. In other words, heavy flavor mesons produced from different collision systems share a similar R AA as long as N part is fixed. The slight breaking of this N part scaling behavior shown in the figures could be due to different initial heavy quark spectra produced at different √ s NN for different collision systems.
Unlike R AA , the N part scaling behavior does not exist for v 2 , as shown in the right panels of Figs. 5 and 6. This is because v 2 is driven not only by the overall energy loss of heavy quarks that is determined by N part , but also by the geometric anisotropy of the medium. For the same centrality class, larger collision system (e.g. Pb-Pb) has higher N part than smaller system (e.g. O-O). In other words, for similar N part , larger system has stronger anisotropy. Therefore, one observes the hierarchy of Pb-Pb > Xe-Xe > Ar-Ar > O-O for the heavy meson v 2 in the right panels of these two figures.
To investigate how the heavy flavor observables rely on the geometry of the QGP, we present R AA and v 2 as functions of centrality in Figs. 7 and 8 for D and B mesons respectively. Similar to the N part dependence figures presented above, we show in each figure the p T -integrated observables within 5 < p T < 8 GeV in the upper panel and 8 < p T < 16 GeV in the lower panel. The left panels are for R AA and the right for v 2 . For a given collision system, we generally observe that the heavy flavor meson R AA increases from central co peripheral collisions due to smaller heavy quark energy loss in more peripheral collisions, while v 2 first increases and then decreases due to the competing effects between parton energy loss and medium geometry. The only exception here is the large D    Comparing different collision systems, one can clearly observe the hierarchies of both R AA and v 2 of heavy flavor mesons. As discussed earlier, for a given centrality class (or medium eccentricity), a larger collision system has a higher N part , resulting in stronger energy loss of heavy quarks through the medium. This yields Pb-Pb < Xe-Xe < Ar-Ar   v 2 than B mesons due to the mass dependence of charm and bottom quark energy loss through the QGP.
Since the elliptic flow coefficient v 2 of heavy flavor mesons depends on both the amount of heavy quark energy loss and the medium anisotropy, its scaling behavior between different collisions systems is hard to be displayed when plotting as a function of either N part or centrality. Interestingly, one may remove the medium anisotropy effect from the total contribution by rescaling the heavy meson v 2 with the bulk medium eccentricity ε 2 . As shown in Fig. 9, the rescaled v 2 /ε 2 is mainly determined by the amount of parton energy loss, thus scales with the system size or N part between different collision systems. This behavior is very similar to heavy meson R AA . Note that although the amount of heavy quark energy loss is the main source of the heavy meson v 2 after removing the bulk geometry effect, the coupling of heavy quark motion to the QGP flow and the hadronization process can also affect the final state heavy meson v 2 and break the scaling behavior of v 2 /ε 2 v.s. N part . Such breaking effect is more prominent for low energy heavy quarks and when the bulk radial flow effect is strong.

V. SUMMARY
Within our Langevin-hydrodynamics framework, we have performed a systematic study on the system size dependence of heavy quark energy loss in heavy-ion collisions at the LHC energies. The space-time evolution of the QGP produced in different collision systems is simulated using our (3+1)-dimensional CLVisc hydrodynamic model. The medium modification of the heavy quark energy-momentum is described by our modified Langevin equation that incorporates both elastic and inelastic scatterings of heavy quarks inside the QGP. By combining this Langevin model with the FONLL calculation for the initial heavy quark spectra and the fragmentationcoalescence model for hadronization, we have calculated the nuclear modification factor (R AA ) and elliptic flow coefficient (v 2 ) of D and B mesons in various centrality regions of Pb-Pb, Xe-Xe, Ar-Ar and O-O collisions at the LHC. The transverse momentum, participant number and centrality dependences of the heavy meson R AA and v 2 have been investigated in detail.
Our results show a clear system size dependence of the heavy meson R AA . For the same collision system, R AA increases from central to peripheral collisions. For the same centrality class, R AA increases (Pb-Pb < Xe-Xe < Ar-Ar < O-O) as the size of colliding nuclei decreases. We have demonstrated a clear scaling of the heavy meson R AA as a function of N part between different collision systems, which indicates a direct correlation between the amount of jet energy loss and the size of the QGP profiles. On the other hand, the heavy meson v 2 simultaneously depends on the size and anisotropy of the QGP. For the same collision system, v 2 first increases and then decreases from central to peripheral collisions. For the same centrality class, v 2 generally increases as the size of colliding nuclei increases, except for the relatively large v 2 in central O-O collisions due to the strong initial-state fluctuations of the small O nucleus. After eliminating the effects of different bulk medium anisotropy in different collision systems, the bulk-eccentricity-rescaled heavy meson elliptic flow (v 2 /ε 2 ) is found to scale with N part . This reveals a direct correlation between v 2 /ε 2 and the amount of heavy quark energy loss which depends on the overall size of QGP. Moreover, the comparison between D and B mesons demonstrates a clear mass dependence of parton energy loss that yields smaller R AA and larger v 2 of D mesons than B mesons for the same collisions system and the same centrality class.
The system size dependence of D and B meson observables discussed in this work provides a crucial bridge between large (Pb-Pb) and small (p-Pb) systems of relativistic nuclear collisions. Comparison between our numerical predictions here and future system-size-scan experiments on jet quenching is expected to help resolve several open questions in high-energy nuclear physics, such as the precise path-length dependence and mass dependence of parton energy loss, and the detailed correlation of collective flow coefficients between hard probes and soft hadrons. Interestingly, our calculation shows considerable amount of heavy quark energy loss even in the small O-O collisions, as suggested by the quenching effect on R AA as well as the finite v 2 of both D and B mesons in central O-O collisions. This further implies that R pA ∼ 1 in proton-nucleus collisions [35,41] is mainly due to the small size of the nuclear medium in these even smaller collision systems. We note that our earlier study [41] suggests that the strong elliptic flow [34] of D mesons in p-Pb collisions cannot be explained by the final state parton energy loss effects. In contrast, Refs. [40,103] show that the initial state gluon saturation effect can explain well the observed elliptic flow of heavy mesons in p-Pb collisions. Further investigations on both heavy and light flavor R AA and v 2 in large and small systems, and their scaling behaviors, may help to identify the boundary where QGP disappears.