Constraining the equation of state with heavy quarks in the quasi-particle model of QCD matter

In a quasi-particle model of QCD matter at finite temperature with thermal masses for quarks and gluons from hard thermal loops, the equation of state (EOS) can be described by an effective temperature dependence of the strong coupling $g(T)$. Assuming the same effective coupling between the exchanged gluon and thermal partons, the EOS can also be related to parton energy loss.} Based on the quasi-particle linear Boltzmann transport (QLBT) model coupled to a (3+1)-dimensional viscous hydrodynamic model of the quark-gluon plasma (QGP) evolution and a hybrid fragmentation-coalescence model for heavy quark hadronization, we perform a Bayesian analysis of the experimental data on $D$ meson suppression $R_{\rm AA}$ and anisotropy $v_2$ at RHIC and the LHC. We achieve a simultaneous constraint on the QGP EOS and the heavy quark transport coefficient, both consistent with the lattice QCD results.

Introduction.Since the discovery of the quark-gluon plasma (QGP) [1,2] at the beginning of this century, the precise extraction of its properties has been one of the key goals of the relativistic heavy-ion collision programs at RHIC and the LHC [3].It is now generally accepted that the QGP in heavy-ion collisions is a strongly coupled system that behaves like a nearly perfect fluid.This allows a hydrodynamic description of the QGP [4,5], based on which one may obtain crucial properties of the system from experimental data on the low transverse momentum (soft) hadrons emitted from the QGP with the help of Bayesian interference.These include the shear and bulk viscosities of the QGP [6][7][8][9], and its equation of state (EoS) [10]; the latter agrees well with the lattice QCD results.
The strength of jet-medium interactions can be quan-tified by a set of transport coefficients, such as the jet quenching parameter q [94,95] that characterizes the transverse momentum broadening of jet partons inside the QGP.Over the past decade, considerable efforts have been dedicated to developing theoretical and statistical tools for extracting q from the jet quenching data, ranging from a traditional χ 2 -fit method that determines its values at RHIC and LHC separately [96], to a Bayesian statistical method that provides a continuous function of q with respect to the medium temperature and the parton energy [97], and a recent information field based global Bayesian approach that starts with no specific ansatz of its functional form [98]. Similarly, the diffusion coefficient of heavy quarks D s has also been successfully obtained using various approaches and shown in agreement with the lattice QCD results [38,[99][100][101][102].
In the above studies, EOS parameterized from lattice QCD results, transport coefficients and initial conditions determined from the soft hadron spectra are used in the viscous hydrodynamics for the QGP evolution.Jet transport coefficients are then constrained separately by jet quenching data with the QGP evolution given by the calibrated hydrodynamic models.A closer investigation on the medium properties is performed in Ref. [103] where the density of scattering centers inside the QGP has been extracted from the jet quenching data and shown in qualitative consistency with the lattice EoS.In this study, we present a direct Bayesian extraction of the QGP EOS using heavy flavor observables within the quasi-particle linear Boltzmann transport (QLBT) model [104] which couples to the (3+1)-D CLVisc hydrodynamic model [105,106] of the QGP evoluton and a hybrid fragmentation-coalescence approach for heavy flavor hadronization [107].In QLBT model, the EOS of the QCD matter can be described by an effective tem-arXiv:2304.08787v2[hep-ph] 11 Dec 2023 perature dependence of the strong coupling g(T ) or the thermal mass of quasi-particles (quarks and gluons).Assuming the same effective coupling g(T ) between the exchanged gluon and thermal partons, the EOS can also be related to parton energy loss and therefore can be constrained by the jet quenching observables.We can achieve a simultaneous constraint on QGP EOS and the heavy quark diffusion coefficient, both consistent with the lattice QCD results.The shear viscosity from this constrained quasi-particle model is also consistent with the values used in CLVisc and from other Bayesian inferences.
Heavy quark interaction in QGP.In the QLBT model [104], QGP medium is viewed as a gas of massive quarks and gluons (quasi particles) [23,24,108], and the elastic and inelastic scatterings between heavy quarks and the medium are described by the linear Boltzmann transport model [63,109].The interaction among medium partons is encoded in their thermal masses [108], for quarks and gluons respectively, where N c = 3 is the number of colors and N f = 3 is the number of flavors.Since the interactions between the medium partons are at the thermal scale, we assume the coupling strength runs with the medium temperature as: in which T c denotes the transition temperature between the QGP and the hadronic matter, a, b, c and d are four parameters to be determined.Although the above equations are inspired by perturbative QCD (pQCD) calculations, it has been shown [104] that they are able to capture key features of the QCD EOS, which is nonperturbative in nature.However, unlike Ref. [104] where g 2 (T ) is directly fitted from the lattice EOS, the inverse question is addressed in the present study -whether one can constrain the QCD EOS with the heavy flavor observables in heavy-ion collisions.
Given thermal masses of these quasi-particles, interactions between them and heavy quarks (denoted by "a") can be simulated via the Boltzmann equation, where p i (i = a, b, c, d) is the four-momentum of parton i (b denotes the thermal parton with a spin-color degeneracy factor γ b , c and d are the final states of a and b respectively), and f i is the phase-space distribution which is taken as f i = 1/(e pi•u/T ± 1) for thermal partons with the local temperature T and flow velocity u provided by the hydrodynamic simulation of the QGP evolution.The summation over b, c and d takes into account all possible elastic scattering processes whose scattering amplitudes |M ab→cd | 2 are calculated at the leading order of pQCD [110].A factor . As discussed in Ref. [104], the thermal masses of medium partons enter in three locations of Eq. ( 4): the thermal distribution, the Mandelstam variables (ŝ, t, û), and the kinematic constraint (S 2 × δ 4 -function).The inelastic scattering process is related to the medium-induced gluon emission from heavy quarks, where the gluon spectrum is taken from the higher-twist energy loss calculation [30,111].There are two types of coupling vertices in the heavy-quark-QGP interactions.In order to simplify the Monte-Carlo simulation, we assume the coupling between the exchanged gluon and thermal partons as the same as that between thermal partons in Eq. ( 3) that leads to the quasi-particle masses in Eqs. ( 1) and ( 2); For vertices directly connecting with heavy quarks, we assume the interactions rely on the heavy quark energy scale and therefore parameterize the corresponding coupling strength as: with A and B as two additional parameters.Therefore, we have a total of six parameters in QLBT.The energy of heavy quark here is evaluated in the local rest frame of the QGP medium.Note that Eq. ( 3) can be reduced to Eq. ( 5) if the energy scale of thermal partons is significantly larger than T c .We use the (3+1)-dimensional CLVisc hydrodynamic model [105,112] initialized with the AMPT model [113] to simulate the spacetime evolution of the QGP fireballs in heavy-ion collisions, from which quasi-particle distributions are extracted using the local temperature and flow velocity.The starting time of the hydrodynamic expansion (τ 0 = 0.6 fm/c) and the specific shear viscosity η/s = 0.08 have been adjusted to provide reasonable descriptions of the soft hadron spectra including the anisotropic flows.Heavy quarks are initialized using the Monte-Carlo Glauber model for their initial position and the LO pQCD calculations of the pair production and flavor excitation processes for their initial momentum distributions.Their subsequent interactions with the QGP are simulated using Eq. ( 4).On the chemical freezeout hypersurface (T c ), they are converted to hadrons using a hybrid coalescence-fragmentation hadronization model [107] that has been constrained by the heavy flavor hadron chemistry measurements.Interactions before (τ < τ 0 ) and after (T < T c ) the QGP stage are neglected in this work.
Bayesian extraction of the model parameters.To simultaneously constrain the six model parameters with the heavy flavor observables, we employ the Bayesian statistical analysis method, in which the posterior distribution of parameters (θ = {a, b, c, d, A, B}) given the experimental observation (data = {y exp i }) is proportional to the product of the likelihood of a given set of θ and its distribution: The likelihood function can be estimated with a Gaussian form that quantifies the similarity between our model output ({y i }) with a given set of θ and the data: where the index i runs over all data points included in our analysis, and σ i is the standard deviation at each data point that combines the experimental error and the interpolation error from our Gaussian process emulator (GP) [119] which is trained with our QLBT calculations using over 800 sets of θ and then applied as a fast surrogate of QLBT when we scan across the parameter space to constrain its posterior distribution.
With the above setup, the prior distribution of θ is assumed uniform in the regions summarized in Tab.I.Note that in principle, the transition temperature T c could also be a model parameter to be constrained from data.However, in order to compare our extracted EoS to two different lattice QCD calculations -Wuppertal-Budapest (WB) [120] and HotQCD (HQ) [121], here we use two separate values of T c from these two work -150 MeV for WB and 154 MeV for HQ.The ranges in Tab.I are chosen such that they can produce a sufficiently wide prior range of the EoS with various shapes, a reasonable coverage of the experimental data from our model calculation, and a well constrained posterior distribution of the model parameters via the Bayesian calibration, as will be shown later.
Starting with these prior distributions of θ, we calculate the nuclear modification factor R AA and the elliptic flow coefficient v 2 of D mesons using the QLBT model in Au-Au collisions at 200 AGeV and Pb-Pb collisions at 5.02 ATeV with different centralities.The GP is trained with 821 parameter sets for the T c = 150 MeV setup and 815 sets for T c = 154 MeV.These design points are sampled according to the Latin-Hypercube algorithm [122] and iterated for multiple rounds of trials to ensure both efficiency and stability in obtaining our results.With the GP, we employ the maximum a posterior (MAP) method in the PyMC library [123] to perform the Markov Chain Monte Carlo (MCMC) estimation of the parameters with Metropolis-Hastings random walk in the parameter space according to the Bayesian statistics.We first implement 320,000 steps of random walks to approach the equilibrium of the posterior distribution, and then run another 320,000 steps from which we take 32,000 sets of parameters, each drawn from every 10 steps in order to exclude correlations from adjacent steps.These 32,000 sets represent the posterior distribution of θ and are used to evaluate the QGP EoS in the end.We use the integrated autocorrelation time method in emcee [124] to quantify the effects of sampling error on our results.The integrated autocorrelation time τ f is summarized in Tab.II.We can see that the chains contain 32,000 sets of parameters, longer than 100 τ f for all parameters.This means that the chains are sufficiently converged.Shown in Fig. 1 is the model calibration against the D meson data at RHIC and LHC, left for R AA and right for v 2 .The upper panels are from our model calculations using the prior distributions of the parameters while the lower panels using their posterior distributions.One observes that after the Bayesian calibration, our QLBT model is able to provide a reasonable description of the D meson observables in heavy-ion collisions.The remaining discrepancies between model and data may result from the limits existed in our ansatz of the coupling strength based on its perturbative form, and the energy loss and hadronization models.Here we present results  [114,115] and LHC [116,117], upper panel for calculations with the prior distributions of parameters and lower panel with the posterior distributions after the calibration.The inserted figures are RAA and v2 at the LHC in the same pT range as RHIC results (i.e., pT = 0 ∼ 12 GeV).The ALICE data [118] are also presented for comparison, though they have not been included in our model calibration.with T c = 150 MeV.Similar results can be obtained for T c = 154 MeV.We expect that heavy flavor observables at high p T are roughly determined by the heavy quark energy loss, which is related to the overall magnitude of the couplings g(T ) and g(E).However, at low to intermediate p T , these observables also depend on the medium flow, whose evolution profile is related to both the magnitude and the shape of the EoS.Therefore, comparing to high p T region, heavy flavor observables at low to intermediate p T region are more sensitive to the EoS.The extracted model parameters are presented in Fig. 2, in which the diagonal panels show their posterior distributions and off-diagonal panels show their correla-tions.Reasonable constraints on these model parameters have been obtained, and certain sensitivity of the extracted parameters on the value of T c can also be observed in the figure (the upper triangle for T c = 154 MeV and the lower triangle for T c = 150 MeV).This will affect our constraint on the QGP EoS later.
QGP EoS and transport coefficients.Using the six parameters constrained by data on the D meson R AA and v 2 , we can evaluate the EoS of the quark-gluon plasma and the transport coefficients of heavy quarks and shear viscosity of the quasi-particle system.
With the temperature-dependent thermal masses of quarks and gluons that are determined by parameters (a, b, c, d) via g 2 (T ), one may calculate the pressure of the relativistic gas system as where we sum over all parton species with E i (p, T ) = p 2 + m 2 i (T ) being their kinetic energy, f i (p, T ) being the Bose/Fermi distribution for gluon/quark and γ i being their spin-color degeneracy factors.Similarly, the energy density is obtained as Although a temperature-dependent bag term B(T ) is introduced in the above two equations to generate a phase transition between the QGP and the hadronic gas, it has no contribution to the entropy density: s(T ) = [ϵ(T ) + P (T )]/T .Since this bag term cannot be determined from the coupling strength g(T ), it has not been included in the Boltzmann equation in this work.Shown in Fig. 3 is the T 3 -rescaled entropy density as a function of the temperature in the QGP phase.The  left and middle panels show the variations of EOS from prior distributions of the parameters (a, b, c, d) in the range given by Tab.I, while the right panel is obtained from their posterior distributions.One observes that the heavy flavor observables do provide constraints on the QGP EoS.Its 95% Credible Region (C.R.) in the right panel collapses into two bands depending on the value of T c we use.For a larger T c that terminates the heavy-quark-QGP interaction earlier, a higher entropy density of the medium is required in order to produce the same amount of nuclear modification on heavy flavor hadrons observed by experiments.The EoS we extract with T c = 150 MeV agrees well with the WB lattice data that shares the same T c , though some deviation can be observed from HQ data.We note that deviation exists between WB and HQ results on s(T )/T 3 as a function of either T or T /T c .One can also extract the interaction strength between heavy quarks and the QGP, which can be quantified by the transverse transport coefficient q of heavy quarks.It is defined as the transverse momentum broadening square of a heavy quark per unit time, d⟨k 2 ⊥ ⟩/dt, due to elastic scatterings and can be derived from Eq. ( 4) as This depends on both the temperature through thermal masses and the coupling between quasi-particles g(T ) and the heavy quark energy through the coupling between heavy quarks and the gluons g(E).Another frequently quoted transport coefficient of heavy quarks is their spatial diffusion coefficient D s , which can be related to q via D s (2πT ) = 8πT 3 /q.In Fig. 4, we present  this diffusion coefficient as a function of the temperature for heavy quark energy E = 10 GeV.Our constraints are consistent with other model results in the literature [19,20,24,25,99] as well as direct lattice QCD calculations [125][126][127].We have also verified that the q parameter we obtain is consistent with the constraint provided by the JET Collaboration study [96] using the inclusive hadron R AA .
With the effective coupling between quasi-particles extracted from the heavy flavor data, one may also calculate the temperature dependent shear viscosity of the quasi-particle system using the relaxation time approximation [108,128], as shown in Fig. 5. Our result is also compared with the Bayesian analysis from the Duke The specific shear viscosity (η/s) as a function of the medium temperature, compared between extractions using two different values of Tc, the Bayesian analysis from the Duke group [7] and the lattice result for pure glue QCD [129].
group [7] and the lattice result for pure glue QCD [129].We have verified that little difference is observed in the soft hadron yield and elliptic flow from CLVisc hydrodynamic simulations between using the shear viscosities obtained in the quasi-particle model and a constant value of η/s = 0.08 that was used for hydrodynamic evolution of QGP in this study.In principle, one can also use the quasi-particle model for shear and bulk viscosity in the hydrodynamic simulations and include experimental data on soft hadron spectra in future combined Bayesian analyses of both hard and soft probes.
Summary.Within the QLBT model for jet and heavy quark transport in QGP which evolves according to the CLVisc hydrodynamics and a hybrid fragmentationcoalescence model for hadronization, we have carried out a Bayesian analysis of the experimental data on D meson spectra and anisotropy v 2 at both RHIC and LHC.We realized a simultaneous constraint on the properties of the QGP and heavy quark probes.The QGP EOS we extract is consistent with the lattice QCD results, and the heavy quark diffusion coefficient we obtain agrees with results from other models and lattice calculations.
In future studies, one can incorporate a more extensive set of parameters, e.g.phase transition temperature T c , that characterize the medium properties, and involve a broader range of jet observables and soft hadron spectra in order to accomplish the goal of constraining the properties of nuclear matter using hard probes.

2 FIG. 1 :
FIG.1:(Color online) Calibration of the QLBT calculation (using Tc = 150 MeV) against the D meson RAA and v2 data at RHIC[114,115] and LHC[116,117], upper panel for calculations with the prior distributions of parameters and lower panel with the posterior distributions after the calibration.The inserted figures are RAA and v2 at the LHC in the same pT range as RHIC results (i.e., pT = 0 ∼ 12 GeV).The ALICE data[118] are also presented for comparison, though they have not been included in our model calibration.

FIG. 3 :
FIG.3:(Color online) Entropy density as a function of temperature, left for its prior range mapped from Tab.I and right for its posterior range constrained by the heavy flavor observables, compared between results obtained with two different values of Tc and their corresponding lattice data[120,121].

FIG. 4 :
FIG. 4: (Color online) The spatial diffusion coefficient of charm quarks Ds as a function of the medium temperature, compared between extractions using two different values of Tc and results from other model calculations.

T
FIG. 5: (Color online)The specific shear viscosity (η/s) as a function of the medium temperature, compared between extractions using two different values of Tc, the Bayesian analysis from the Duke group[7] and the lattice result for pure glue QCD[129].

TABLE I :
The ranges of model parameters used in the prior distributions, for two different values of Tc.

TABLE II :
The integrated autocorrelation time τ f of model parameters in the posterior distributions, for two different values of Tc.