Fluid-dynamics of charm quarks in the quark–gluon plasma

A ﬂuid-dynamic approach to charm-quark diﬀusion in the quark-gluon plasma (QGP) is developed for the ﬁrst time. Results for integrated yields and momentum distributions of charmed hadrons obtained with a ﬂuid-dynamic description for the dynamics of the QGP coupled to an additional heavy-quark-antiquark current are shown. In addition to the thermodynamic Equation of State (EoS), this description uses a heavy-quark diﬀusion constant which we take from Lattice QCD calculations. The results describe quantitatively experimental data measured at the LHC at the center-of-mass energy of √ s NN = 5.02 TeV up to p T ∼ 4-5 GeV/ c , showing that charm quarks undergo a very fast hydrodynamization in the medium created by ultrarelativistic heavy-ion collisions.

Introduction.Heavy quarks, i.e., charm and bottom, are produced in heavy-ion collisions (HICs) via hard partonic scattering processes.Due to their large mass and early production, they are suitable probes for studying the quark-gluon plasma (QGP).Their in-medium dynamics is tackled employing a Boltzmann/Langevin description in several transport models [1][2][3][4][5][6][7][8][9][10].However, recent experimental measurements [11,12] showed that J/ψ and D mesons display a positive elliptic flow, suggesting an early local thermalization or "hydrodynamization" of charm quarks within the QGP.The idea of charm thermalization was suggested by the Statistical Hadronization Model for charm [13] and supported by Lattice QCD (LQCD) calculations [14].In our recent work [15], the question of charm thermalization was addressed by studying the hydrodynamization time of charm quarks in the context of an expanding medium.It was shown that the time required for charm quark to hydrodynamize, and therefore to be included in the fluid-dynamic description of the QGP, is shorter than the typical expansion time scale of the medium.This result served as motivation to develop a fluid-dynamic description of charm quarks, which is the subject of the current work.We expect such a description to be relevant for the low transverse momentum (p T ) region, as it is for the light-flavor particles.At high momentum, the path-length-dependent energy loss mechanisms, are more important in defining the shape of the p T spectra.Fluid-dynamic equations.The fluid-dynamic equations to solve are mainly given by the system of equations which expresses the conservation of the energymomentum tensor T µν and of an additional conserved current N µ .The latter is associated with conserving * f.capellino@gsi.dethe number of charm-anticharm pairs [15].The Landau frame is chosen such that T µν and N µ can be decomposed as where ϵ, P , u µ , Π and π µν are the energy density, thermodynamic pressure, fluid four-velocity, bulk viscous pressure, and shear-stress tensor of the fluid, respectively.The charm-quark fields are the heavy-quark density n and the diffusion current ν µ .The local temperature T and the chemical potential to temperature ratio α are determined by the Landau matching conditions, ϵ(T ) ≡ ϵ equilibrium (T ) (5) n(T, α) ≡ n equilibrium (T, α) .(6) We assume that the energy density is approximately independent of the heavy-quark contribution, such that any energy density dependence on α is negligible.The thermal equilibrium heavy-quark density is taken to be one of the hadron-resonance gas, including all measured charm states (HRGc), where M i is the mass of each charm hadron, and q i is its charm charge.The HRGc is expected to give the correct limit for the thermodynamics of the charm density at temperatures close to the phase transition.This relation is assumed to also hold at high temperatures.In the temperature regime reached by the fireball in most central collisions, the HRGc yields larger values (of about a factor 5) than the density of the free charm quarks.Nevertheless, due to the absence of first principle calculations for the Equation of State of charm quarks at physical QGP temperatures, we assume this relation to hold also at high temperatures.In the future, a more realistic Equation of State will be developed.
The equations of motion for each of the dissipative currents in a second-order hydrodynamic formalism are solved, where one defines the projector and the vorticity tensor ω µν = (∇ µ u ν − ∇ ν u µ )/2.Here we introduced the transport coefficients for the bulk viscosity ζ, shear viscosity η, and the heavyquark diffusion coefficient κ n , with the corresponding relaxation times τ Π , τ π and τ n .The values of the viscosities are taken from Ref. [16], while the expression for the diffusion coefficient was derived in Ref. [15].We remark that κ n and τ n are proportional to the heavy-quark spatial diffusion coefficient D s .The equations are solved in Bjorken coordinates assuming boost and azimuthal rotation invariance, restricting effectively to 1+1 dimensions.We organize the fluid fields for the QGP into a Nambu spinor Φ = (T, u µ , π µν , Π), which satisfies the hyperbolic equation of motion.We assume that none of these fields or transport coefficients depend on the heavy-quark variables.Eqs. ( 1), ( 8) and ( 9), can be used to determine the time derivatives of the fluid fields explicitly.Let us now consider another Nambu spinor including also the heavy-quark fields Ξ = (T, u µ , π µν , Π, α, ν µ ).The new system of hyperbolic equations satisfied by Ξ can be numerically solved by setting the fluid fields contained in Φ on shell.This is equivalent to neglecting the back reaction of the heavyquark field on the fluid background evolution.More details regarding our approach to solve the fluid-dynamic equations, the numerical implementation and the validation of our fluid-dynamic framework can be found in the Appendix A, B and C. Initial conditions for the fluid fields.The initial condition for the entropy density is computed with T R ENTo [17] simulating Pb-Pb collisions at √ s NN = 5.02 TeV.The T R ENTo parameters are set based in Refs.[18,19]; the T R ENTo output is used as entropy density.The nucleonnucleon cross section for Pb-Pb collisions at √ s NN = 5.02 TeV is taken from [20], i.e. x = 67.6 mb.The nucleons in the Pb ion are sampled from a spherically symmetric Woods-Saxon distribution with radius R = 6.65 fm and surface thickness a = 0.54 fm.Using this set of parameters, the transverse density T R (x, y) is generated for 1.5 • 10 6 minimum-bias collisions, among which the ones belonging to the 10% most central are selected.The normalization of the T R ENTo profile is computed by fixing the multiplicity of protons to the measured one [21].The proton multiplicity in our calculation is obtained by employing a Cooper-Frye+FastReso approach at the end of the fluid-dynamic evolution as in Refs.[16,22,23].In the future, when performing a Bayesian analysis to fit the experimental measurements, the normalization will be left as a free parameter of our model as in Ref. [23].
The initial conditions for the temperature field are then obtained through the thermodynamics EoS described in Ref. [24].Radial fluid velocity, shear-stress tensor components, and bulk viscous pressure are initialized at zero.Initial conditions for charm fields.The midrapidity density of charm quarks at the initialization time of the hydrodynamic evolution τ 0 comes from the initial hard production, In the above expression, the QQ rapidity distribution in nucleus-nucleus collisions is set by the pQCD QQ crosssection where σ in is the inelastic proton-proton cross-section and σ QQ is the hard production cross-section.The average number of collisions N coll is computed with a Glauber model and depends on the impact parameter of the collision, providing: where n coll is assumed to be distributed according to the fluid energy density n coll ∝ T 4 .As a future development, one could evaluate the radial distribution of binary collisions directly from T R ENTo , not to neglect spacemomentum correlations that are important for flow observables.The integral of the density in the transverse plane provides the total number of heavy quarks to be conserved throughout the QGP evolution.As discussed in [15], we remark that the current associated with the number of heavy quark-antiquark pairs is accidentally conserved.The heavy-quark mass is too large for them to be produced thermally throughout the QGP evolution; moreover, the annihilation rate of a QQ pair is negligible within the lifetime of the plasma.
To fix at each point the initial value for α for the QQ pair, Taking the central prediction by FONLL [25] for collisions at √ s NN = 5.02 TeV, one gets, at y = 0, dσ QQ /dy = 0.463 mb, with σ in = 67.6 mb [20].At the beginning of the system evolution, the thermal distribution at zero chemical potential overshoots the density of charm quarks in the middle of the fireball.Therefore, α assumes negative values initially to match the hard production.This is not expected to happen at the fireball evolution's end, where the charm species' thermal abundance will be strongly suppressed.The total multiplicity of QQ pairs per unit of rapidity is given by the integrated density profile, e.g. at τ = τ 0 , In terms of fluid variables, due to the conservation of the charm current, the conserved charge is rewritten as, where |g| is the determinant of the metric.Besides the density, we can initialize the heavy-quark diffusion current.The assumed symmetries would allow a nonvanishing radial component, but we set it to zero in the absence of a more detailed initial state model.
Evolution of the fields.The initial conditions for the fields are set on a hypersurface at constant proper time τ 0 = 0.4 fm.In Fig. 1 (upper panel), the time evolution of the charm density times the longitudinal proper time as a function of the radial coordinate is reported for different values of τ .This is shown for a non-diffusive (D s = 0) and temperature-dependent D s case obtained by linearly fitting results from LQCD calculations [14].As expected, the density becomes more diluted when the temperature decreases.In the diffusive case, the density evolution is concurrent with developing the radial component of the diffusion current (Fig. 1, lower panel).Its values are always negative, thus negatively contributing to the conserved current N µ .This results in a higher density n in the diffusive case, as shown in Fig. 1.Comparing it to the equilibrium composition of the heavy-quark density n, one finds that the condition of |ν r | ≪ n is not satisfied in the entire radial region.This indicates that the out-of-equilibrium components of the heavy-quark distribution remain large throughout the evolution of the plasma.However, the magnitude of the diffusion current strongly depends on the spatial diffusion coefficient and its correspondent relaxation time.LQCD computations [14] favor a fast hydrodynamization of charm quarks and, thus, a reduction of the out-of-equilibrium correction.Around freeze-out we decompose the single-particle distribution functions, f i = f i,eq + δf i , where the equilibrium part f i,eq is given by the ideal gas distribution and δf i represents the out-of-equilibrium correction.In general, the δf i correction receives a contribution from all the dissipative stresses Π, π µν and ν µ , such that δf i = δf i,bulk + δf i,shear + δf i,diffusion .In our case, the open-charm hadrons distribution function includes both light and heavy components.To properly describe it, one should derive its expression in a multi-species fluid setup.
As for now, we neglect out-of-equilibrium corrections to the fluid variables at the freeze-out surface.In the future, we will address the inclusion of non-linear terms in the evolution equation for the dissipation current and the derivation of a more consistent expression of the total distribution function.Integrated yields.The charmed-hadron production is assumed to occur on a freeze-out hypersurface at a constant temperature.This chosen temperature is T fo = 156.5 MeV [13,26].The freeze-out hypersurface in the plane of Bjorken time τ and radius r is parametrized by a parameter γ ∈ (0, 1).According to the Cooper-Frye prescription, a sudden decoupling is assumed at the freeze-out temperature, and the thermal momentum distribution of the particles is computed according to where g hc accounts for the degeneracy of the produced charmed hadron and q accounts for the charm content of the hadron.The total integrated yield dN hc /dy per unit rapidity for charmed and anti-charmed hadrons is measured by integrating Eq. ( 17).The feed-down from resonance decays is calculated using the FastReso package [22].The list of resonances is taken from the PDG [27].In Fig. 2, the comparison between the obtained integrated yields and experimental measurements [28][29][30][31] is shown for the 0-10% centrality interval.The yields and the p T spectra correspond to the sum of particle and anti-particle divided by two, as reported by experiments.The p T integration range is from 0 to 10 GeV/c.These results are computed for D s = 0 since the integrated yield should not depend on the spatial diffusion coefficient.However, since out-of-equilibrium corrections to the single-particle distribution function at freeze-out are neglected, there can be a non-physical dependence of the yields on D s .While the relative abundance of each charmed-hadron species depends mainly on the mass of the hadron, the absolute value of the integrated yields strongly depends on the EoS for the charm density as a function of T and α.The HRGc as EoS is the most suitable choice to estimate the thermal production of the hadrons and resonances included in the HRGc.
The role played by the resonance decays is then to reshuffle the relative abundance of the hadrons while keeping the total number of charm quarks fixed.The agreement between the model and the measurements is quantified in the lower panel of Fig. 2. We observed that the mesons are compatible with the experimental uncertainties, computed as the sum in quadrature of the statistical and systematic uncertainties.A deviation of 2.4σ is observed for the Λ + c baryons.This deviation might be caused by missing higher resonance states in the PDG [13,32,33].Due to the resonances decay, the yield of the D 0 increases by a factor 4.3, while the one of the Λ + c of a factor 5. Another factor 2.3 would be needed to reproduce the experimentally measured yield.Estimates for the yields of the Ξ + c and Ω 0 c , whose experimental measurements are not yet available, are provided.Most likely, these values will underestimate the actual yields due to the lack of knowledge of higher resonance states.In other phenomenological models, the charmed-baryon enhancement is attributed to a recombination process between the heavy quark and light thermal partons [34][35][36][37].
Momentum distributions.In Fig. 3, the p T -differential spectra for the same hadron species are reported and compared with the experimental measurements [28][29][30][31].A ratio plot with the data to model comparison can be found in Appendix D. The bands correspond to a spread of the input value of the spatial diffusion coefficient D s going from a non-diffusive case (D s = 0) to a temperature-dependent 2πD s T [14].The fluid-dynamic description seems to capture the physics of D mesons up to p T ∼ 4-5 GeV/c.This implies that, even if the charm does not move collectively with the rest of the fluid in the early stage of the evolution, it relaxes to the same radial flow of the QGP before the freeze-out occurs.As observed for the integrated yield, the Λ + c calculation underestimates the experimental measurement.The J/ψ p T distribution describes the experimental measurements for p T < 3 GeV/c, while it overpredicts the yield for higher p T .This discrepancy for p T > 3 GeV/c might be attributed to the dominant contribution from primordial J/ψ, which is not accounted for in our model since it is not expected to reach thermal equilibrium [38][39][40], but is mainly sensitive to path-length-dependent effects, like survival probability and energy loss.It is also important to note that the experimental measurements consist of J/ψ directly produced in the collisions plus the contribution from beauty hadron decays.Including the out-of-equilibrium corrections in the model at the freeze-out surface will influence the shape of the momentum distributions.They would modify the spectra at intermediate/high p T .When adequately included, we do not expect such a strong dependence on D s in the spectra but rather only a tilt in the momentum distribution.A further remark regards the dependence of the final momentum distribution on the initial conditions for the charm fields.In particular, a broader initial distribution for the charm density results in a larger average p T at freeze-out.A more thorough study of the charm initial conditions will improve the description of the transverse momentum distribution of the charm hadrons, without of course impacting the results for the integrated yields.
Conclusions.A fluid-dynamic description of the charm quark is developed for the first time, unveiling that lowp T charm quarks undergo a very fast hydrodynamization in the QGP created during ultrarelativistic heavy-ion collisions.The developed model describes the charmedhadron yield and the p T -differential distribution up to p T ∼ 4-5 GeV/c.The calculations are carried out for a non-diffusive case and a temperature-dependent D s .The derivation of out-of-equilibrium corrections in a multispecies setup will be addressed.Additional constraints on the spatial diffusion coefficient will be set in future works via a Bayesian analysis using both p T -differential spectra and anisotropic flow coefficients.In addition, this study paves the way for a fluid-dynamic description of the beauty quark, which, despite its larger mass, might for Ξ 0 c and Ω 0 c baryon states, which have not been measured yet, are also shown.
still reach a partial local equilibrium allowing further constraining QGP parameters using heavy quarks as probes.
the evolution equation of the dissipative currents π µν , Π and ν µ .Schematically, the equations can be written as quasi-linear partial differential equations (PDEs).We will restrict ourselves to discussing the equations in one spatial dimension for simplicity.However, the extension to a higher number of dimensions is trivial.We consider a collection of independent variables called ϕ, whose equations of motion are where A(ϕ) is a matrix in field space that depends nonlinearly on the fields themselves, and S(ϕ) is a vector containing the source term in the equation.Usually, the numerical solutions of the fluid dynamic equations are discussed in a conservative form since the ideal limit of the equations is the divergence of a current -typically the energy-momentum and particle density current.Let be the conservation equation, where J µ represents generically the conserved current.However, including the dynamics of the dissipative current like diffusion and shear/bulk viscous stress spoils this property for Israel-Stewart-Müller theory [41].For this type of theory, the equations are non-conservative by construction, and it is impossible to cast them in a conservative form.In the relativistic viscous fluid dynamic literature, the equations are solved with a splitting algorithm: First, solve using a finite volume conservative scheme, then correct the intermediate solution using a central approximation of the dissipative equations, as in the so-called SHASTA algorithm [42], or some variations of it like KT [43].This type of algorithm performs well if the dissipative currents are minor corrections to the ideal step and do not modify the ideal evolution substantially.However, this is not always the case, especially when the system is far from the ideal approximation, meaning the non-equilibrium effects are important.In this work, we implement a different strategy.Instead of using the ideal-viscous splitting, we solve the equations together as a quasi-linear system of PDEs.The naive discretization of equations like Eq. (B1) can be obtained by replacing the first derivative with its central approximation.Denoting x i the central position of a cell of size ∆x, the the central derivative approximation is where ϕ i = ϕ(x i ).The semi-discretized version of the equations is This naive discretization, however, is unstable since there is no dissipation mechanism in the discretization to reduce the high-frequency mode of the discretized solution.
The physical motivation for this instability can be understood considering the nature of the PDE.The system of hyperbolic equations is a collection of propagating waves that interact non-linearly and with a non-constant velocity.The waves are usually (except in simple cases) a complicated combination of the primary variables ϕ, defined as the left eigenvector of the matrix A(ϕ).The eigenvalue is characteristic of the hyperbolic PDE and represents how fast the wave propagates.Each of the waves propagates at a different speed and direction.In a one-dimensional case, there will be right-and left-moving waves.To have a stable discretization, the numerical derivative should respect -up to some degree of accuracy -the direction of propagation of the different waves.If a wave is right-moving, the correct derivative discretization should involve only points in the past of the wavei.e. on its left -and vice versa.This mechanism is called upwinding [44].Therefore, the central approximation of the first derivative goes against this principle since it does not distinguish the direction of propagation of the waves.A natural solution is to separate right-moving and leftmoving waves and discretize them accordingly.By calling λ i the eigenvalues, one can separate them into positive and negative ones (λ + i and λ − i , respectively)1 , (B4) Each matrix has only information about the left and right propagating waves, respectively.With this construction, it is then easy (in principle) to construct an upwinding discretization as, ) where the derivatives are taken from the left or the right, respectively, The proposed discretization is sometimes called the fluxsplitting technique and was already introduced in [44][45][46].The drawback of this scheme is that it relies on the complete knowledge of the spectrum of the characteristic matrix.Only in a few cases is this achievable due to the complexity of the non-linearities of the characteristic matrix A.
The discretization reported in Eq. (B5) can be expressed in terms of the absolute value of the matrix A, such that Therefore, Eq. (B5) becomes The derivative operators now become leading to a discretized equation of the form The extra contribution introduced to upwind the derivative acts like a viscous terms into the equation, with an amplitude proportional to the lattice spacing ∆x.
A standard approximation for the absolute value of the matrix is |A| = λI where λ = max(|λ i |), which is the fastest characteristic speed in the system.Under this assumption, the scheme can be considered a nonconservative version of the Lax-Friedrichs scheme.However, this requires knowledge of the characteristic structure, which is possible only for exceptional cases.
An appealing alternative is to approximate |A| with a suitable expansion, as discussed in [46,47].Among the possible expansion choices that one can make, the simplest is a polynomial approximation around max( assuming that all the |λ i | < 1, such that the fastest wave speeds are modified correctly.Different and more performing possibilities are Chebyshev polynomials and rational functions.However, in this work, we restrict ourselves to the simplest choice and will study these possibilities in the future. For the evolution, we use the explicit Runge-Kutta with adaptive time-step as described in [48,49] and with the Proportional-Integral-Derivative (PID) controller as described in [50][51][52].For the integration on the freeze-out surface we used [53][54][55][56].

Appendix C: Validation against Gubser flow
Comparing it against a known analytic (or semianalytic) solution is useful to verify and validate the numerical scheme.For Israel-Stewart-type theories, such a solution with azimuthal rotation symmetry, longitudinal boost invariance, and an additional conformal symmetry has been found by Gubser [57].For symmetry reasons, the evolution of the diffusion current in this setup is trivial.So we will leave it out of the discussion in the rest of this section.The set of equations for the evolution of temperature, fluid velocity, shear stress, number density in the presence of a conformal symmetry reads, where θ = 2 tanh ρ is the scalar expansion rate for Gubser flow.In de Sitter space, by applying the Gubser flow profile ûµ = (1, 0, 0, 0), the equations read, where ρ is the Gubser conformal time variable and πµν = π µν /(ŝ T ).The transformation rules to obtain the fluid variables in Milne coordinates are given by T (τ, r) = T (ρ(τ, r))/τ , (C8) The conformal Equation of State at finite α can be written as where one defines the dimensionless coefficients Here we use p 0 = (16 + 10.5N f )π 2 /90 and the number of flavors N f = 2.5.In this setup, the equations for the charge current are decoupled from the rest of the system.In Fig. 4 the comparison between the semi-analytical solution by Gubser and the one obtained numerically is presented.The initialization time is τ 0 = 1 fm, the shear viscosity to entropy ratio is 0.2; the shear relaxation time is τ S = 5η/(sT ).The overall agreement is good for all fields in the full radial range.In Fig. 5, the percent deviation of the numerical solution for the temperature field with respect to Gubser's solution is shown for different numbers of discretization points at τ = 2 fm.As one can see, the finer the spatial grid is, the smaller the deviation.In particular, the deviation around 2 fm, corresponding to the maximum of the temperature profile, is progressively suppressed.
Appendix D: Details on the results of charm-hadron momentum distributions In Fig. 6 the results for the ratio between the experimental measurements of charm-hadron momentum dis- tributions and the results from our fluid-dynamic model are shown.The bands correspond to a spread of the input value of the spatial diffusion coefficient D s going from a non-diffusive case (D s = 0) to a temperature-dependent 2πD s T obtained by linearly fitting results from LQCD calculations [14].The fluid-dynamic descriptions captures the behavior to D 0 and J/Ψ up to p T ∼ 2 GeV.At intermediate transverse momentum, our calculation for the D mesons deviates of 25% from the experimental measurements for the D s = 0 case.A larger deviation is hereby observed for J/Ψ attributed to the dominant contribution from primordial J/ψ, which is not accounted for in our model since it is not expected to reach thermal equilibrium [38][39][40], but is mainly sensitive to path-length-dependent effects, like survival probability and energy loss.As observed for the integrated yield, the Λ + c calculation underestimates the experimental measurement.This deviation might be caused by missing higher resonance states in the PDG [13,32,33].At p T larger than 5 GeV, the fluid-dynamic model seems no longer applicable since it's not able to capture the behavior of the particle spectra.
A further remark regards the dependence of the final momentum distribution on the initial conditions for the charm fields.In particular, a broader initial distribution for the charm density results in a larger average p T at freeze-out.A more thorough study of the charm initial conditions will improve the description of the transverse momentum distribution of the charm hadrons, without of course impacting the results for the integrated yields., Λ + c , and J/ψ momentum distributions.Experimental measurements are taken from [28][29][30][31].

FIG. 1 .
FIG.1.Charm density times the longitudinal proper time (upper panel) and diffusion current (lower panel) as a function of radius for different times.Solid lines correspond to an ideal hydrodynamic evolution, with Ds = 0. Dashed lines correspond to a diffusive hydrodynamic evolution, with 2πDsT taken from LQCD[14].

D 0 DFIG. 2 .
FIG.2.Charmed-hadron integrated yields with and without feed-down contributions from resonance decays and comparison with experimental data from the ALICE collaboration.

FIG. 4 .
FIG.4.Temperature (upper panel) and chemical potential to temperature ratio α (lower panel) as a function of radius r at Bjorken times τ = 1.5 fm/c and τ = 2 fm/c.The solid lines correspond to the semianalytic Gubser solution, while the dashed lines are the numerical result with N = 200 discretization points.We have here chosen the maximal radius to be 10 fm.

FIG. 5 .
FIG. 5. Percent deviation of the absolute value of ∆T =T numerical /T semianalytical − 1 of the temperature at τ = fm as a function of radius and the number of discretization point N .The numerical solution converges to the semi-analytical one as the number of points increases.