Theory of heavy-quarks contribution to the quark-gluon plasma viscosity

The shear viscosity of quark gluon plasma is customarily estimated in the literature using kinetic theory, which, however, is well known to break down for dense interacting systems. Here we propose an alternative theoretical approach based on recent advances in the physics of dense interacting liquid-like systems, which is valid for strongly-interacting and arbitrarily dense relativistic systems. With this approach, the viscosity of strongly interacting dense heavy-quarks plasma is evaluated analytically, at the level of special relativity. For QGP well above the confinement temperature, the theory predicts that the viscosity increases with the cube of temperature, in agreement with evidence.

Understanding shear viscosity at a microscopic level is vital to solving outstanding mysteries in physics, such as the so-called viscosity minimum [1] and the existence of nearly non-dissipative fluids described by the Euler hydrodynamics such as quark-gluon plasmas [2][3][4] and the microscopic mechanism of superfluidity in helium-4 [5].
Relativistic Heavy Ion Collider (RHIC) experiments at Brookhaven National Laboratory studying Au + Au collisions provided evidence of a new state of matter, the quark-gluon plasma (QGP), at extremely high temperature and density, where hadrons melt into a soup of their de-confined constituents, i.e. quarks and gluons [6].Surprisingly, the experimental data on QGP were initially found to be well fitted by hydrodynamic models assuming Euler inviscid flow.Such nearly dissipationless hydrodynamics is only possible when the shear viscosity η is very low.Finding a mechanistic explanation to this surprising result has been a grand challenge to theoretical physics for the last 30 years [7].
In parallel, calculation of the dimensionless ratio of shear viscosity η to entropy density s by Kovtun, Son and Starinets (KSS) within anti-de Sitter/conformal field theory (AdS/CFT) gave a lower bound on the viscosity η/s ≥ ℏ/4πk B .This has been believed to be a universal lower bound to be obeyed by all physical systems, although violations of this bound have been reported [8][9][10], including, possibly, in the QGP viscosity [11].Recent Bayesian analysis of experimental data support values of the QGP viscosity that are comparable to the KSS bound [12].
Much effort has been devoted to computing the viscosity of QGP numerically and from first-principles theoretical approaches [13].Most recent theoretical models based on AdS/CFT [14] confirm the existence of a nearly dissipationless perfect fluid regime at high densities.However, effective field theory approaches, such as those based on AdS/CFT, cannot provide a microscopic mechanistic understanding of how viscosity vanishes based on microscopic interactions.Other approaches widely used in the literature are based on kinetic theory, i.e. the Boltzmann equation with the Chapman-Enskog expansion in the relaxation time approximation [15,16].As is well known in statistical physics [17], kinetic theory based on the Chapman-Enskog framework can, at best, provide a semi-quantitative description of the viscosity of gases, and even with the most modern corrections the theory is not fully quantitative [18], let alone for liquids.In practice, the kinetic theory (Chapman-Enskog) is applicable only when the dimensionless parameter r 3 c n is small, where r c is the interaction range between the particles and n is the number density [17].Considering the typical mass of a quark (3.56 -8.9 ×10 −30 kg) and the typical interaction range (1 fm), the dimensionless parameter is in the range r 3 c n ∼ 10 3 −10 4 , which is certainly not small and thus makes the kinetic theory inapplicable.
Alternatively, the Green-Kubo formalism based on time correlation functions of the stress-energy tensor [19,20], appears to be the best suited approach, although for calculations of real systems it involves the fitting of long-time tails of autocorrelation functions, which is an art in its own.Also, it does not give access to simple analytical formulae.Among recent approaches that focused on the diffusivity of heavy quarks, it is worth mentioning the nonperturbative treatment of Ref. [21].Recent state-of-art calculations of QGP viscosity using the Green-Kubo formalism within lattice QCD can be found in Ref. [22].Previous calculations [23] based on measuring correlations of the energy-momentum tensor within lattice gluodynamics found viscosity to entropy-density ratios close to the KSS bound.
Here, we present a new method of computing the viscosity of dense strongly-correlated matter based on combining the concept of nonaffine displacements in dense disordered systems [24,25] with the Langevin transport model of heavy quarks [26,27].
We start from basic viscoelastic theory (which applies to both solids and liquids [28]), by introducing the complex shear modulus G * [28][29][30].The real and imaginary part, G ′ and G ′′ , of G * correspond to the dissipationless and to the dissipative part of the response, respectively, and they are also known as the storage (elastic) modu- Since the deformation of a disordered dense system is always nonaffine due to the lack of centrosymmetry around a tagged particle, these particles are displaced with respect to the thick dashed lines in the sheared configuration.The distance between these actual positions and the affine positions on the thick dashed lines represent the nonaffine displacements.The lower panel illustrates the fact that in a centrosymmetric system such as a crystal lattice, all forces acting on a tagged particle in the affine position would cancel each other out by centrosymmetry (left panel) whereas in a dense disordered system such as a QGP the tagged particle in the affine position is not occupying a center of inversion symmetry, which means that the net force transmitted from its nearest neighbours is not zero.This net force has to be relaxed via an additional (nonaffine) displacement (see top panel).
lus and the loss (viscous) modulus.For a generic stress wave σ 0 e iωt , the strain wave lagging behind would be γ 0 e (iωt+ π 2 ) , hence the factor i ≡ e i π 2 in front of G ′′ , leading to the complex shear modulus defined as In the linear response regime, with intrinsic material properties that do not vary with time and with causality being obeyed, G ′ and G ′′ are related to each other by the Kramers-Kronig relations.The above relation applies equally to non-relativistic as well as to relativistic systems, as discussed below.Furthermore, one obtains the following splitting of the shear stress into elastic and viscous components where ϵ xy denotes the macroscopic shear strain.For a relativistic off-diagonal shear deformation µν = ik with i ̸ = k the second term in the above equation represents the dissipative stress [28] σ and the off-diagonal relativistic strain rate: In Eq. ( 2), the first term refers to the elastic part of the response, which takes place on an infinitely long timescale where relativistic effects are absent, by construction.The second term in Eq. ( 2) is the viscous stress described by Eqs. ( 3)-( 4) at the level of special relativity (see [28]).In particular, upon neglecting the elastic part (as appropriate for a nearly perfect fluid), the total stress reduces to the viscous stress, σ ≈ σ ′ .This, in turn, corresponds to an off-diagonal space-like component of the stress-energy tensor, for example σ = σ ′ = T 12 where T 12 is the µ = 1, ν = 2 component of the stress-energy tensor T µν with µ, ν = 0, 1, 2, 3.This component describes shear deformation in the spatial 12 plane of the Minkowski space.
The above equations yield the following identification [28][29][30]: between the viscosity η of the system and the loss modulus G ′′ .We now follow the Langevin transport model for heavy quarks [26,27], and we generalize it by making it relativistic.In the Langevin transport framework, the nonrelativistic Langevin equation describes the motion of heavy quarks, including charm.One uses the fact that the heavy quark mass is much larger than the temperature scale, and that most of the heavy quarks in the QGP will have small momentum in the frame of the bath.The formalism is best suited for the bottom quark, mass ∼ 5 GeV, and is also used for charm, mass ∼ 1.5 GeV.Typical estimate of the temperature of QGP gives ≤ 0.3 GeV.A systematic expansion in 1/M , the heavy quark mass, of the heavy quark motion has been worked out in Ref. [31].The size of the 1/M correction, following this formalism, was estimated in [32] and was found to be moderate also for charm.Note, however, that this analysis is valid for heavy quarks moving slowly with respect to the heatbath frame (heavy quark momentum ≪ mass).The QGP itself is not non-relativistic, its main constituents are gluons and light quarks with mass ≪ temperature.So, in principle, the viscosity of QGP has to be evaluated in a relativistic framework, and the hydrodynamic formalism used for the medium flow has to be relativistic.
In the following we will estimate the relativistic viscosity contribution of the heavy quarks using a method valid for dense strongly-correlated matter, which provides a direct link between the shear viscosity and the frictional damping.
Introducing the mass-scaled tagged-particle displacement s = Q √ m, with Q the canonical coordinates, the resulting non-relativistic generalized Langevin equation of motion for the displacement of the tagged particle becomes (cfr.[33] for the full derivation from a microscopically reversible Caldeira-Leggett particle-bath Hamiltonian): where F p (t) is the thermal stochastic noise with zero average, U is a local interaction potential (e.g. with the nearest neighbours), and ν is the friction resulting from many long-range interactions with all other particles in the system, imposed by the dynamical bi-linear coupling in the Caldeira-Leggett Hamiltonian [33][34][35][36].
As rigorously demonstrated in recent work [35,36], the above Langevin equation can be written in manifestly covariant form for a relativistic system in terms of the four-momentum: where is the rank-two friction tensor [36], and we have reintroduced the particle mass m.Clearly, the nonrelativistic and the relativistic Langevin equations, Eq. ( 6) and Eq. ( 7) respectively, have, mutatis mutandis, the same form.Hence, from now on we shall work using the coordinates measured in the lab frame and employ suitable transformation laws when necessary.The above relativistic Langevin equation can also be re-written in terms of the Cartesian 3-vector displacements of particle i (measured in the deformed lab frame) [36]: where γ is the Lorentz factor and we take an isotropic frictional memory-function, K µ ν (t) ≡ ν(t).The relativistic velocity transformation between the lab frame and the co-moving frame of the macroscopic deformation is given by: ṙi = 1 where the dot indicates a time derivative (with respect to time t measured in the lab frame) while the circle indicates quantities measured in the undeformed rest (lab) frame.Furthermore, where u = Ḟr i represents the local velocity of the deformed frame.
Since we are interested in the zero-shear rate limit of the shear viscosity (static shear viscosity), which is the parameter that enters the Navier-Stokes equation and the subsequent reduction to Euler equations for the QGP, then both Ḟ and u are small and certainly much smaller than the speed of light c.Hence, with u ≪ c, Eq. ( 9) reduces to the (Galilean) transformation: We can thus rewrite Eq. ( 8) for a tagged particle in d-dimensions, which moves with an affine velocity prescribed by the strain-rate tensor Ḟ as follows: We now can expand the time derivative on the l.h.s. to obtain: For heavy quarks, which are often treated nonrelativistically, as e.g. in Refs.[26,27,37], it is totally acceptable to assume that γ 3 ( ṙi) → 0, and therefore approximate Eq. ( 13) as: We then arrive at: where f µ i = −∂U/∂r µ i generalises the −U ′ (s) to the tagged particle.Furthermore, we used the following consideration.
This notation is consistent with the use of the circle on the particle position variables to signify that they are measured with respect to the reference rest frame (undeformed).In terms of this original rest frame {r i }, the equation of motion can be written, for the particle position averaged over several oscillations, as where ⟨F p ⟩ = 0 is used upon averaging over many particles.
We work in the linear regime of small strain ∥ F − 1 ∥≪ 1 by making a perturbative expansion in the small displacement {s i (t) = ri (t) −r i } around a known rest frame ri .That is, we take F = 1 + δF + ... where δF ≈ F − 1 is the small parameter.Replacing this back into Eq.( 16) gives For the term δf i , imposing mechanical equilibrium again, which is f i = 0, implies: where the double dot product is used because η is a ranktwo tensor.In the first term on the r.h.s. of the above expression Eq. ( 19) we recognise where H ij represents the Hessian matrix, defined as since r(γ)| γ→0 = r 0 , where U denotes the total energy of the system.For the second term in Eq. ( 19) we have the following identification: and the limit η κχ → 0 is implied.Here η κχ are the components of the Cauchy-Green strain tensor defined as η = 1 2 F T F − 1 .This is a second-rank tensor, and should not be confused with the fluid viscosity (a scalar), also denoted as η.From a physical point of view, Ξ i,κχ represents the force acting on a heavy quark in the affine position.
With these identifications, we can write Eq. ( 18), to first order: The above equation can be solved by Fourier transformation followed by a normal mode decomposition, as we shall see next.If we specialise on time-dependent uniaxial strain along the x direction, η xx (t), then the vector Ξ i,xx η xx (t) represents the force acting on particle i due to the motion of its nearest-neighbors which are moving towards their respective affine positions (see e.g.[38] for a more detailed discussion).Hence, all terms in the above Eq.( 23) are vectors in R 3 and the equation is in manifestly covariant form.
To make it convenient for further manipulation, we extend all matrices and vectors to N d × N d and N ddimensional, respectively, and we will select d = 3.After applying Fourier transformation to Eq. ( 23), we obtain where ν(ω) is the Fourier transform of ν(t) etc (we use the tilde consistently throughout to denote Fouriertransformed quantities).In the above equation, all the terms are are now vectors in R 3N space (we work in the heavy M → ∞ mass limit, where N , the number of heavy quarks, is conserved).Next, we apply normal mode decomposition in R 3N using the 3N -dimensional eigenvectors of the Hessian as the basis set for the decomposition.This is equivalent to diagonalize the Hessian matrix H.
Proceeding in the same way as in [39], we have that the m-th mode of displacement can be written as: ) It was shown [38], by means of MD simulations that Ξm,κχ = v m • Ξ κχ is self-averaging (even in the glassy state), and one might therefore introduce the smooth correlator function on eigenfrequency shells on frequency shells.Following the general procedure of [38] to find the oscillatory stress for a dynamic nonaffine deformation, the stress is obtained to first order in strain amplitude as a function of ω (note that the summation convention is active for repeated indices): In the thermodynamic limit and assuming a continuous vibrational spectrum, we can replace the discrete sum over 3N degrees of freedom with an integral over vibrational frequencies up to an ultraviolet cutoff frequency ω D (in the kinetic theory of matter, this is known as the Debye frequency and is an increasing function of density ω D ∼ (N/ V ) 1/3 ).In this case, we need to replace the discrete sum over the 3N degrees of freedom (eigenmodes) with an integral, 3N m=1 ... → ω D 0 g(ω p )...dω p , where g(ω p ) is the density of states (DOS) or energy spectrum of the system, which can be normalized to unity or to the total number of degrees of freedom (e.g.3N ) [40].
Choosing the latter normalization of the DOS, the complex elastic constants tensor can be obtained as: where N/ V denotes the density of the system in the initial undeformed state [39].In the above we used massrescaled variables throughout, so that the particle mass is not present explicitly in the final expression.If one, instead, uses non-rescaled variables and specializing to off-diagonal shear deformations, µνκχ = xyxy, we obtain the following relativistic expression for the complex shear modulus G * : The first term on the r.h.s. is the affine shear modulus G A , which is independent of ω.The low-frequency behaviour of Γ(ω p ) can be estimated analytically using the following result: derived originally in Ref. [24], which gives ⟨ Ξ2 p,xy ⟩ ∝ λ p , thus implying (from its definition above): Γ(ω p ) ∝ ω 2 p .This analytical estimate appears to work reasonably well in the low-eigenfrequency part of the Γ(ω p ) spectrum of amorphous solids [38,41,42].Furthermore, α B α,µνκχ is a geometric coefficient that depends only on the geometry of macroscopic deformation and its values can be found tabulated in Ref. [24].
By separating real and imaginary part of the above expression, we then get to the storage and loss moduli as [42]: It is easy to check that the storage modulus G ′ (ω) reduces to the affine Born modulus G A ≡ G ∞ in the limit of infinite oscillation frequency, ω → 0. From the point of view of practical computation, the DOS g(ω p ) can be obtained numerically via direct diagonalization of the Hessian matrix H ij , since its eigenvalues are related to the eigenfrequencies via λ p = mω 2 p .Similarly, the affine-force correlator Γ(ω p ) can also be computed from its definition by knowing the positions of all the particles, their interactions and forces (so that the affine force fields Ξ can be computed) as well as the eigenvectors of the Hessian v p .
We note that, in the zero-frequency limit for the quasistatic elastic shear modulus we have: which agrees with previous derivations, cfr.[38].For a fluid or liquid at thermal equilibrium, it has been rigorously demonstrated that the two terms in the above equation (i.e. the positive affine term and the negative non-affine term) are identical, with opposite sign, see [25,43], resulting in G ′ (0) = 0, in agreement with the well known fact that fluids and liquids possess zero shear rigidity and flow under an external mechanical perturbation.This also applies to the QGP which is well described by hydrodynamic equations of fluids, such as the Euler equation.
We recall that the viscosity can be obtained from the loss viscoelastic modulus G ′′ using nonaffine response theory by means of Eq. ( 5).The nonaffine response theory developed from first principles above thus provides the following form for the loss modulus G ′′ (cfr.Eq. ( 31)): ) For shear deformation (κχ = xy), Eq. ( 27) becomes: Since G * = G ′ + iG ′′ and there is an overall factor ω in the above expression for G ′′ in Eq. ( 33), the theory correctly recovers, for the dissipative part of the stress, σ ′ xy , the relativistic viscous flow given by Eq. ( 3).This, indeed, has a zero-frequency shear viscosity given by η = lim ω→0 G ′′ /ω, leading to: As a sanity check, this formula correctly recovers the formula for the viscosity of non-relativistic fluids [25] upon taking the non-relativistic limit, i.e. as γ → 1.
By comparing with the Kubo relationship used for QCD matter: η = lim ω→0 ρ(ω)/ω, with ρ(ω) the spectral function or spectral density [44], we obtain the following identification: Based on the most recent literature studies [37,45], the low-energy part of the spectral function for the deconfined QGP can be written as [46]: , where κ d is the momentum diffusion coefficient.
Since in the low-frequency sector the main excitations are acoustic waves with a linear dispersion, we take a Debye-type approximation for the DOS spectrum, , where c s is the speed of sound.Upon recalling Eq. ( 30) and λ p = mω 2 p , the correlator Γ(ω p ) is given by Γ(ω p ) = 1  5 mκR 2 0 ω 2 p [24].By further including the degeneracy factor of 12 for the 3 colors, 2 spin states and particle-antiparticle, we thus get: where the effective spring constant κ has dimensions of force per unit length, and in the high-T QGP, dominated by collisional physics, can be expressed in terms of the relevant energy scale divided by the separation distance squared, κ ∼ k B T R 2 0 . The rational behind this is that the interaction between heavy quarks will generate excitations of field quanta (gluons, or, for that matter, phonons), which eventually thermalize too.For the same reason, in the above relation, ω D ∼ T (vibrational temperature, ℏω = k B T ), and γ ∼ T for an ultrarelativistic fluid (since γmc 2 ≈ k B T ).Furthermore, in the regime T > 2T c we can assume the speed of sound c s of the QGP to be approximately T -independent [47].Therefore, putting all these things together, we obtain the following estimate for the temperature dependence of the viscosity: for the QGP viscosity in the high energy regime T > 2T c , which agrees very well with the available data and earlier calculations [44,48,49].
In conclusion, we have presented a microscopic derivation of the viscosity coefficient of strongly correlated dense QCD matter and we applied it to the case of heavy quarks in the QGP within a new relativistic Langevin transport framework.The derivation leads to a fundamental expression for the viscosity contribution of the heavy quarks which turns out to be proportional to the cube of temperature, in agreement with previous studies.This provides an analytical microscopic derivation of QGP viscosity treated as a dense relativistic liquid.Furthermore, it provides a new fundamental link between the viscosity, the Lorentz factor, and the speed of sound.In future work, this model can be extended to include further microscopic details e.g.along the lines of [50].

Figure 1 .
Figure1.The top panel shows particles in a strongly interacting dense fluid subjected to a shear deformation.Particles that, in the undeformed configuration (left) sit on the thick dashed vertical lines would land on the thick dashed lines of the deformed configuration (right) if the deformation were completely affine.Since the deformation of a disordered dense system is always nonaffine due to the lack of centrosymmetry around a tagged particle, these particles are displaced with respect to the thick dashed lines in the sheared configuration.The distance between these actual positions and the affine positions on the thick dashed lines represent the nonaffine displacements.The lower panel illustrates the fact that in a centrosymmetric system such as a crystal lattice, all forces acting on a tagged particle in the affine position would cancel each other out by centrosymmetry (left panel) whereas in a dense disordered system such as a QGP the tagged particle in the affine position is not occupying a center of inversion symmetry, which means that the net force transmitted from its nearest neighbours is not zero.This net force has to be relaxed via an additional (nonaffine) displacement (see top panel).