Flavor hierarchy of parton energy loss in quark-gluon plasma from a Bayesian analysis

The quenching of light and heavy flavor hadrons in relativistic heavy-ion collisions probes the color and flavor dependences of parton energy loss through a color-deconfined quark-gluon plasma (QGP), and thus reveals the properties of QCD matter at extremely high density and temperature. By combining a next-to-leading order perturbative QCD calculation of parton production, a general ansatz of parton energy loss functions and parton fragmentation functions, we calculate the nuclear modification of various hadron species -- charged hadrons, $D$ mesons and $B$-decayed $J/\psi$ -- over a wide transverse momentum regime. Comparing our calculations to the experimental data using the Bayesian statistical analysis, we perform a first simultaneous extraction of the energy loss functions of gluons ($g$), light quarks ($q$), charm quarks ($c$) and bottom quarks ($b$) inside the QGP. We find that the average parton energy loss at high energies follows the expected hierarchy of $\langle \Delta E_g \rangle>\langle \Delta E_q \rangle \sim \langle \Delta E_c \rangle>\langle \Delta E_b \rangle$, while the parton energy loss distribution can further test the QCD calculations of parton interaction with the dense nuclear matter. We also find that the reduction of experimental uncertainties can significantly improve the precision of the extracted parton energy loss functions inside the QGP.

One of the central goals of studying jets in heavyion collisions is utilizing them to extract the properties of nuclear matter at different temperature scales.For instance, considerable efforts have been devoted to the extraction of the jet transport coefficient q in both hot and cold nuclear matter [46][47][48][49][50].This q parameter measures the transverse momentum broadening of an energetic parton due to its elastic scatterings with the medium [51,52].It is directly related to the density of medium constituents, and is a key quantity in estimating the intensity of the medium-induced gluon bremsstrahlung from the hard parton [53][54][55].Recently, it has been proposed that the enhancement of q near the critical endpoint may also be applied to probing the QCD phase diagram [56].Besides q, one may also directly extract the amount of parton energy loss (∆E) inside the QGP.For example, the energy loss distribution (or the quenching weight) was obtained in Ref. [57] by comparing a simplified model calculation with the experimental data on the nuclear modification factor (R AA ) of single inclusive hadrons.In this work, the form of the parton energy loss is taken from the BDMPS mediuminduced gluon spectrum [7,58,59], and both light flavor hadrons and prompt J/ψ's are assumed to be produced only from gluon fragmentation while D mesons are only from charm quark fragmentation.Later in Ref. [60], a more general ansatz of the jet energy loss distribution is assumed and convoluted with the medium-modified jet function, based on which the flavor-averaged jet energy loss function is extracted using the jet R AA data with the Bayesian statistic interference.This methodology is then extended in Ref. [61] for constraining the energy loss distribution functions of both gluons and charm quarks from the nuclear modification of J/ψ at high p T .However, a simultaneous data-driven constraint on the energy loss functions of all parton species -gluon (g), light quarks (q = u, d, s), charm quark (c) and bottom quark (b) -is still absent in literature.This is the focus of our present study, which is essential in unraveling the color, flavor and mass hierarchy of parton energy loss.
Considering the successful perturbative QCD (pQCD) description of different hadron species at high p T (≳ arXiv:2303.12485v2[hep-ph] 19 Apr 2024 10 GeV) in p+p collisions, we will convolute the well-established next-to-leading-order (NLO) perturbative calculations of parton production [62,63], a general ansatz of parton energy loss inside the QGP, and parton fragmentation [64][65][66].By comparing the nuclear modification of charged hadrons, D mesons and B-decayed J/ψ between our model calculation and the experimental data at the Large Hadron Collider (LHC) [14,67,68], we will perform a simultaneous extraction of the energy loss functions of g, q, c and b using the state-of-the-art Bayesian analysis method.The mean values of the parton energy loss at high energies exhibit a clear flavor hierarchy, ⟨∆E g ⟩ > ⟨∆E q ⟩ ∼ ⟨∆E c ⟩ > ⟨∆E b ⟩, and their distribution functions can provide a stricter test on QCD calculations on parton scatterings inside a hot nuclear matter.The uncertainties propagated from the experimental data to our result will also be explored.

II. HADRON PRODUCTION AND MEDIUM MODIFICATION
The differential cross section for high-p T hadron production in p+p collisions can be factorized as In the above equation, dσ pp→jX represents the cross section for inclusive parton (j) production and can be calculated by convoluting the parton distribution functions (PDFs) of the two colliding protons with the partonic scattering cross section; D j→h denotes the fragmentation function (FF) of parton j to a given hadron species h.In this work, the PDFs are taken from the CTEQ parameterizations [69], the partonic cross section is evaluated perturbatively at the NLO [62,63], and the FFs are taken from Ref. [64] for charged hadrons, Ref. [65] for D mesons, and Ref. [66] for B mesons.As shown in our earlier studies [18,70,71], this combination provides a good description of charged hadron, D meson and B meson spectra at p T ≳ 10 GeV in p+p collisions at the LHC energy.Note that within the NLO framework, one can consistently include both quark and gluon fragmentations to heavy and light flavor hadron productions.In particular, the gluon fragmentation dominates the charged hadron production at p T < 50 GeV, and contributes to almost 40% D meson and 50% B meson productions up to p T beyond 50 GeV.
In heavy-ion collisions, we assume the parton production from the initial hard scatterings is a superposition of ⟨N coll ⟩ p+p collisions, where ⟨N coll ⟩ represents the average number of (binary) nucleon-nucleon collisions in each A+A collision.The nuclear shadowing effect is taken into account by modifying the nucleon PDFs using the EPS09 parameterizations [72].The hard partons then interact with the QGP and lose transverse momentum ∆p T according to an energy loss function W AA (x) [60], where x = ∆p T / ⟨∆p T ⟩ is the ratio between the transverse momentum loss in a particular event and its mean value.Hadronization is assumed to take place outside the medium, and therefore the vacuum fragmentation function can be applied as in p+p collisions.With this setup, the binary-collision-number-rescaled cross section of hadron production reads: Here, dσ p ′ p ′ →jX represents the parton cross section after including the shadowing effect, and different species of hard partons lose different amount of ⟨∆p T ⟩ which is parameterized as where β g controls the overall magnitude of gluon energy loss, γ tunes its p T dependence, and C j represents the parton energy loss ratio relative to the gluon's -1 for gluon and C q , C c , C b for light, charm and bottom quarks respectively.The normalized distribution function of parton energy loss is given by with α being a model parameter, which in principle can be computed from jet energy loss theory.Generally, parameters γ and α above can also depend on the parton flavor.However, they were shown to be consistent between gluon and charm quark in an earlier study using the J/ψ data [61].Thus in the present work, they are set to be the same for g, q, c and b.Releasing this assumption will consume much longer computational time, and will be left for our future effort.Thus, we have in total 6 parameters -β g , C q , C c , C b , γ and α -in this analysis.In this work, we do not include a specific energy loss model, while uncertainties in our extracted parameters can arise from the nuclear PDFs and FFs we use.The energy loss function we obtain in the end also depends on the parameterization scheme designed here.However, different schemes should lead to comparable results if they capture the main feature of the parton energy loss and meanwhile provide good descriptions of the experimental data after the Bayesian calibrations.
The ratio of Eq. ( 2) to Eq. ( 1) gives the nuclear modification factor (R AA ) of a given type of hadron.By comparing our model calculation to the R AA data of charged hadrons, D mesons and B-decayed J/ψ, we will extract the above 6 parameters and obtain the energy loss functions of different parton species.Note that in order to exclude the effects of non-perturbative interactions of partons with the QGP at low p T [73][74][75], and the coalescence process during their hadronization [76,77] at low to intermediate p T , we restrict our current study to the high p T (> 10 GeV) regime where the perturbative calculation is reliable.Instead of directly using the B meson R AA to constrain the energy loss function of b quarks, we decay B mesons to J/ψ through Pythia [78], given that the experimental uncertainties of the latter is much smaller than the former [68,79].
The Bayesian inference method has been successfully employed in both soft and hard sectors of heavy-ion physics, such as the extractions of initial condition, bulk transport coefficients and the equation of state of the QGP medium [80][81][82][83][84][85][86], jet transport coefficient q [49], heavy quark diffusion coefficient D s [87], and the energy loss distribution of hard partons and jets [60,61].It has also been utilized to constrain the fluctuating structure of protons with experimental data on coherent and inherent diffractive J/ψ production in electron-proton collisions at HERA [88].In this work, we implement the Bayesian inference method to simultaneously constrain the energy loss functions of different partons (g, q, c and b) through the QGP medium.
Our model calibration on the parameter set, denoted by θ, is based on the following Bayes' theorem, P (θ|data) ∝ P (θ) P (data|θ) , where P (θ|data) on the left represents the posterior distribution of the parameter set given the knowledge of the experimental data, P (θ) on the right represents the prior distribution of θ without any knowledge of data, and P (data|θ) measures the likelihood of a given set of θ by comparing its model output with the data.In the present study, θ = (β g , C q , C c , C b , γ, α) is a 6-dimensional vector, and the likelihood function is assumed to follow the Gaussian form as where y i (θ) is the model output at a given data point i, and y exp i and σ i are respectively the mean value and the error of the experimental observation at the point i.The statistic error and the systematic error of the experimental data are combined in our analysis.
Using Eqs. ( 5) and ( 6), we perform a Markov-Chain Monte-Carlo (MCMC) random walk [89] in the parameter space according to the Metropolis-Hasting algorithm.Each step of this random walk updates the parameter set based on its current location in the parameter space.We start from a random location and execute 10 7 steps for the Markov-Chain to reach equilibrium, after which the locations given by further steps can constituent equilibrium distributions of the model parameters, which are the posterior distributions we seek in this work.We generate another 10 7 steps after reaching equilibrium, from which we draw these posterior distributions.To exclude correlation between adjacent steps, we draw one sample of θ from every 10 steps, and use 10 6 samples in the end to produce the probability distribution of the parameter space.
A uniform prior distribution P (θ) is assumed in the regions of β g ∈ (0, 10), C q ∈ (0, 1), C c ∈ (0, 1), C b ∈ (0, 1), γ ∈ (−0.15, 0.5) and α ∈ (0, 15), which are sufficiently wide for our model output to cover all the experimental data.In order to accelerate the calibration speed, a Gaussian Process (GP) emulator [90,91] is applied for the model-to-data comparison when we scan through the 6-dimensional parameter space.This emulator is trained with our model results on 600 design points uniformly distributed in the prior range of θ above, and then used as a fast surrogate of the real perturbative calculation during the Bayesian analysis on the hadron R AA .
To validate our setup of the GP emulator and the Bayesian analysis framework, we first conduct a closure test in Fig. 1 based on a set of pseudo-experimental data generated by the model calculation.We start with a particular design point -θ 0 = (2.35,0.55, 0.5, 0.2, 0.15, 7)as indicated by the black vertical lines in the figure, and use it to calculate the R AA of charged hadrons, D mesons and B-decayed J/ψ with the perturbative method developed in Sec.II.These model results are used to replace the mean values of the experimental data in 0-10% Pb+Pb collision at √ s NN = 5.02 TeV [14,67,68], while the error bars of these data points are taken from the original CMS measurements.Using these pseudo-data, we then examine whether the posterior distribution of θ extracted from the Bayesian analysis agrees with the "ground-truth" value θ 0 that we implant.The posterior distributions of the 6 model parameters are presented in Fig. 1  from the pseudo-experimental data, the third column from using the original uncertainties of the experimental data, and the fourth column from using the halved uncertainties.

IV. EXTRACTING PARTON ENERGY LOSS FUNCTIONS
Using the validated statistical analysis framework above, we now extract the distributions of parton energy loss from the real R AA data of charged hadrons, D mesons and B-decayed J/ψ in 0-10% Pb+Pb collisions at √ s NN = 5.02 TeV [14,67,68].
In Fig. 2, we present the posterior probability distributions of the 6 model parameters extracted from the Bayesian fit to the experimental data.The diagonal with   [14,67,68], middle column from using the original uncertainties of the experimental data, and right column from using the halved uncertainties.panels show the probability distributions of individual parameters while the off-diagonal ones show the twoparameter joint distributions that quantify their correlations.The lower triangle (blue) is from using the original error bars of the experimental data in our Bayesian analysis, while the upper triangle (red) from the halved error bars.The 90% C.R.'s of these parameters are summarized in Tab.II, including results from using both the full (middle column) and halved (right column) error bars.One can see that the model parameters can be well constrained by comparing our perturbative calculation to the experimental data.An inverse correlation between β g (or FIG. 3: (Color online) The nuclear modification factors of charged hadrons, D mesons and B-decayed J/ψ in 0-10% Pb+Pb collisions at √ sNN = 5.02 TeV, compared between the perturbative calculation using the posterior distributions of model parameters and the CMS data [14,67,68].The theoretical bands correspond to the ±1σ uncertainties around the mean values.
C q , C c , C b ) and γ is observed.This is consistent with earlier findings in Refs.[60,61], and can be understood with Eq. ( 3) where they both contribute positively to the parton energy loss.The posterior distribution of the parameter α here is also consistent with that obtained earlier in Ref. [61].In addition, reducing the uncertainties from the experimental data makes the posterior distributions of the model parameters narrower, thus helps better constrain the parton energy loss functions.Note that the ranges of C q obtained here are consistent with the ratio of the Casimir factors between quark and gluon.
Shown in Fig. 3 is our model result after the Bayesian calibration.We first calculate the R AA 's of different hadrons using parameters drawn from their posterior distributions.Then, we evaluate the mean values and standard deviations (σ) of R AA 's at each p T , and present their ±1σ bands in Fig. 3. Our calibrated model calculation provides a simultaneous description on the nuclear modification factors of charged hadrons, D mesons and B-decayed J/ψ measured by the CMS Collaboration.
Using the parameters constrained by the ±1σ band of the R AA results above, we generate the average fractional transverse momentum loss (⟨∆p T ⟩/p T ) according to Eq. (3), and compare results between g, q, c and b in Fig. 4. From the figure, one can clearly observe the flavor hierarchy of parton energy loss through the QGP: gluon loses about twice energy than quark does due to the larger color factor of the former.Light quarks and charm quark also suffer stronger energy loss than bottom quark does due to the larger mass of the bottom.Since we focus on the high p T regime where the perturbative calculation is reliable, the mass effect on the charm quark energy loss is small.As a result, the extracted values of ⟨∆p T ⟩/p T are similar between light and charm quarks.Such flavor hierarchy of parton energy loss from our data-driven analysis is consistent with expectation from perturbative QCD calculations on both elastic and inelastic parton energy loss inside the QGP [92,93].In the end, we present the distribution function of parton energy loss in Fig. 5.This distribution can in principle be calculated from jet quenching theory, thus can provide a stringent test on the QCD calculations of parton interactions with the QGP.

V. SUMMARY
We have conducted a systematic data-driven analysis on the parton energy loss inside a QGP medium.A perturbative framework has been developed for calculating the nuclear modification of high p T hadron produc-tion in relativistic heavy-ion collisions, which convolutes the cross section of parton production at NLO, a general ansatz of the parton energy loss function inside the QGP, and parton fragmentation functions into different hadrons.By comparing the nuclear modification factors of charged hadrons, D mesons and B-decayed J/ψ between our calculation and the LHC data, we have performed a first simultaneous extraction of energy loss of gluon, light quarks, charm quark and bottom quark using the Bayesian interference method.Our result shows a clear flavor hierarchy of parton energy loss at high energies, ⟨∆E g ⟩ > ⟨∆E q ⟩ ∼ ⟨∆E c ⟩ > ⟨∆E b ⟩ inside a hot nuclear matter, consistent with perturbative QCD expectation.We also find that a reduction of the data uncertainties can significantly improve the precision of the extracted parton energy loss functions.Our current study provides a model-independent method to quantitatively constrain the flavor dependence of parton energy loss from the experimental data.The extracted distribution function of parton energy loss also provides a stringent constraint on the perturbative QCD calculation of parton interactions with the QGP.

FIG. 1 :
FIG.1:(Color online) Closure test: posterior distributions of model parameters (diagonal panels) and their correlations (off-diagonal panels) extracted from pseudo-experimental data on the RAA's of charged hadrons, D mesons and Bdecayed J/ψ, with black vertical lines denoting the input ("ground truth") values of θ0 = (2.35,0.55, 0.5, 0.2, 0.15, 7) used to generate the pseudo-data.The lower triangle (blue) is obtained using the original error bars of the experimental data[14,67,68] while the upper triangle (red) is obtained using halved error bars.

FIG. 2 :
FIG.2:(Color online) Posterior distributions of model parameters (diagonal panels) and their correlations (off-diagonal panels) extracted from the CMS data on the RAA's of charged hadrons, D mesons and B-decayed J/ψ in 0-10% Pb+Pb collisions at √ sNN = 5.02 TeV[14,67,68].The lower triangle (blue) is obtained using the original error bars of the data while the upper triangle (red) is obtained using halved error bars.
, with diagonal panels showing the probability distributions of individual parameters and off-diagonal ones showing the correlations between different pairs of parameters.The lower triangle (blue) is obtained by using the original error bars of the CMS data in the Bayesian analysis, while the upper triangle (red) is obtained by using halved error bars.Meanwhile, we present the 90% credible regions (C.R.'s) of the extracted model parameters in Tab.I, which are obtained by discarding the lowest 5% and highest 5% areas of their distribution functions.
One can see that the probability distribution, or the 90% C.R., of each model parameter can be constrained to its pre-set "true" value within this analysis framework.Halving the error bars of our pseudo-data narrows the 90% C.R.'s of the posterior distributions, and thus improves the precision of the extracted parameters.

TABLE I :
The 90% C.R.'s of the model parameters extracted