Heavy quark diffusion in a Polyakov loop plasma

We calculate the transport coefficients, drag and momentum diffusion, of a heavy quark in a thermalized plasma of light quarks in the background of Polyakov loop. Quark thermal mass and the gluon Debye mass are calculated in a non-trivial Polyakov loop background. The constituent quark masses and the Polyakov loop is estimated within a Polyakov loop quark meson (PQM) model. The relevant scattering amplitudes for heavy quark and light partons in the background of Polyakov loop has been estimated within the matrix model. We have also compared the results with the Polyakov loop parameter estimated from lattice QCD simulations. We have studied the temperature and momentum dependence of heavy quark drag and diffusion coefficients. It is observed that the temperature dependence of the drag coefficient is quite weak which may play a key role to understand heavy quark observables at RHIC and LHC energies.

This paper is organized as follows, in section II we give the formalism for calculating drag and diffusion of heavy quarks by employing Boltzmann equation in soft momentum exchange between heavy quark and bulk medium [45]. In section III we recapitulate and summarize the calculation of the Debye mass and the quark thermal mass in a Polyakov loop background as has been outlined in Refs. [43,57]. In these calculations, we have also kept the effects of a possible finite quark mass. Such an effect can be important near the transition temperature where the light quark condensates could still be relevant. The drag and the diffusion coefficients are evaluated in section (IV) where we discuss their behavior as a function of temperature as well as momentum. Finally, in section (V) we summarise the results and present a possible outlook. We summarise the salient features of PQM model in Appendix(A). Further, in Appendix(B), we give some details of the calculation for the square of matrix elements for the relavant 2 → 2 processes.

II. FORMALISM
In the QGP phase, the Boltzmann equation for charm quark distribution function, neglecting any mean-field term, can be written as [45,46]: where f HQ represents the spatially integrated non-equilibrium distribution function for heavy quark. The right hand side of Eq.(1) is the collision integral where the phase-space distribution function of the bulk medium appears as an integrated quantity. If we define ω(p, k) as the transition rate of collisions of the heavy quark with the heat bath particles (light quarks/antiquarks and gluons) that change the heavy quark momentum from p to p − k, then we can write [45] ∂f HQ ∂t col The first term in the integrand represents a gain of probability through collisions which knock the charm quark into the volume element of momentum space at p and the second term represents the loss out of that volume element. ω(p, k) is the total contributions coming from heavy quark scattering from gluon and light quark/anti-quark. Furthermore, assuming the scattering processes to be dominated by small momentum transfer, we can expand ω(p+k, k)f HQ (p+k) around k, ω(p + k, k)f HQ (p + k) ≈ ω(p, k)f HQ (p) + k · ∂ ∂p (ωf HQ (p)) + 1 2 k i k j ∂ 2 ∂p i ∂p j (ωf HQ (p)).
The higher power of the momentum transfer, k i 's, are assumed to be small [47]. Keeping up to the second term and substituting in Eq.(2), we get: Now Eq.(1) is reduced to Fokker-Planck equation, where the kernels stand for the drag and the diffusion coefficients respectively. The function ω(p, k) is given by ω(p, k) = g q,g dq (2π) 3 f l (q)vσ p,q→p−k,q+k , where f l (q) is the thermal phase space distribution of the particles which constitute the heat bath which in the present case stands for light quarks/anti-quarks and gluons, v = |v p − v q | is the relative velocity between the two collision partners, σ denotes the interaction cross section and g q/g is the statistical degeneracy factor for light quarks/antiquarks and gluons.
In particular A i and B ij , for the (generic) process, HQ(p) + l(q) → HQ(p ′ ) + l(q ′ ) (l stands for light quarks and gluon), are given by [45,[48][49][50]: g HQ is the statistical degeneracy of the charm quark. The factor f l (q) denotes the thermal phase space factor for the gluons and light quarks/anti-quarks in the incident channel and 1 ± f l (q ′ ) is the final state Bose/Fermi enhanced/suppression phase space factor. The above expression indicates that the drag coefficient is a measure of the thermal average of the momentum transfer, p − p ′ , weighted by the elastic heavy quark-bulk interaction through the square of the invariant amplitude, | M | 2 . Similar, heavy quark diffusion coefficients can be defined as: From the above expression it is clear that the diffusion coefficient is a measure of the thermal average of the square of momentum transfer weighted by the elastic heavy quark-bulk interaction through the square of the invariant amplitude, |M | 2 . Since A i and B ij depend only on the vector p, we may write [45]: where, The integrals appearing in the above equations can be further simplified by solving the kinematics in the center of mass frame of the colliding particles and both the drag and diffusion coefficients can be defined from a single expression: with an appropriate choice of Γ(p ′ ). As in Ref. [45], we shall consider 2 → 2 processes which involve Coulomb scattering i.e., qQ → qQ through gluon exchange and Compton scattering of gluon and heavy quark i.e., gQ → gQ. In the present work, we shall estimate the scattering amplitudes in the background of Polyakov loop. This makes the square of the corresponding matrix element as well as distribution function dependent on the color indices (see e.g., Eqs. (19) and (20)). Therefore, the expression for Γ(p ′ ) becomes where |M C | 2 abef is matrix element squared for q a Q b → q e Q f with ab(ef ) as initial(final) quark color indices and |M Cm | 2 abef gh is matrix element squared for g ef Q b → g gh Q a scatterings with ef, a(gh, b) as initial(final) gluon and quark color indices. Here the color indices a, b, e, f, g, h = 1, 2, 3 are in fundamental representation. Furthermore, in Eq. (15), λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2zx is the triangular function. E p and m C are the heavy quark energy and mass, respectively. E q is the energy of the light quark/gluon. s is the Mandelstam variable. To compute the heavy quark transport coefficient, one needs, therefore, the heavy quark-light quark/gluon scattering matrix along with the thermal distribution functions, the mass of light quarks and gluons and Debye screening mass. The divergence in t-channel diagram here is regulated by a Debye mass [45].
In literature several attempts have been made, over the years, to compute the heavy quark drag and diffusion coefficients in QGP within different models. A recent study indicates that non-perturbative contributions are essential for the simultaneous description of heavy quarks R AA and v 2 [20]. Quasi-particles model is a way to take into account the non-perturbative effect. This can be done in a number of possible ways which differ in how the effects of QCD interactions are modeled. To study the heavy quark transport properties in QGP, the quasi-particle approaches [20,24] that have been recently used in literature include the interaction in the effective masses of the light quark and gluons. In these quasi-particle models strong coupling constant [51], g(T ), is the only free parameter which can be obtained by making a fit of the energy density obtained by lattice QCD calculations. The main feature of these quasi-particle approach is that the resulting coupling is significantly stronger than the one obtained from pQCD running coupling particularly near the quark-hadron transition temperature (T c ). In this present study we adopted a different model to include the non-perturbative effects. The statistical distribution function, thermal mass and Debye mass have been obtained in presence of a non-trivial Polyakov loop background. In the following section we attempt to estimate quark thermal mass and the Debye mass in Polyakov loop background.

III. THERMAL AND DEBYE MASSES IN POLYAKOV LOOP BACKGROUND
In this section, we shall estimate the non-perturbative Debye screening mass and quark thermal mass in a nontrivial Polyakov loop background to be used in the estimation of the drag and diffusion coefficients using Eqs. (11) and (12). Such a calculation has been performed in detail in Refs. [43,44,57,59] using a matrix model for semi-qgp and used for estimating shear viscosity to entropy ratio as well as to dilepton and photon production and energy loss of heavy quark in the medium. We recapituate the salient features of such a calculation including also the possible effects from a finite mass of the light quarks which can arise from a nonvanishing scalar quark-antiquark condensate.
Polyakov loop is a particular case of the Wilson loop where the gluon field is time-like. The background gauge field can be taken as a constant diagonal matrix A ab µ = δ µ0 δ ab Q a /g, where, the color index a is not summed and g is the gauge field coupling constant. The Wilson line in the temporal direction is given by where, P denotes path ordering in the imaginary time, with τ being the imaginary time τ : 0 → β. In the mean field level, neglecting the fluctuations and with the choice of the time-independent constant background field, the path ordering becomes irrelevant and one can perform the integration over the imaginary time leading to P = exp(igβA 0 ). The trace of the Wilson line is the Polyakov loop φ given as In a SU(N) gauge group the vector potential A 0 is traceless so the sum over all the Q's vanishes i.e., a Q a = 0, for SU (3), one can parameterize Q a = 2πT (−q, 0, q) , where we have introduced a dimensionless Polyakov loop dependent parameter "q" [44], so that.
Physically, such a nontrivial background field A 0 can be thought of as an imaginary chemical potential [54]. The thermal distribution function for the quarks/anti-quarks and the gluons are given respectively by [43] f a (E) = 1 Let us note that the quark distribution function involves only one color index because these are represented in fundamental representation. For gluons, the adjoint representation leads to two fundamental indices. For three colors, the color averaged statistical distribution function of the gluons becomes Similarly, for three colors, the color averaged distribution functions of the quark/anti-quark is It may be noted that for pure gluon case, φ = 1 in the confined phase and φ = 0 in the deconfined phase. This leads to the gluon distribution function in the confined phase and in the deconfined phase. In the presence of quarks, one does not have a rigorous order parameter for deconfinement, however in φ = 0 case the color averaged quark/anti-quark distribution reduces to f q/q (E) = 1 e 3βE + 1 (25) so that quark are suppressed statistically. In the perturbative limit i.e., φ = 1 it becomes The color averaged distribution function of quark/anti-quark as given in Eq. (22) is exactly the same as that in PQM model within mean field approximation [42]. For the computation of Debye and thermal mass, we use double line notation [55,56] which is convenient for large N c calculations. For SU (N ) gauge group, the generators λ A satisfy the following relation [57] T where A and B are adjoint indices and takes the values A, B = 1, 2, 3, .., N 2 − 1. Each adjoint indices can be denoted by a pair of fundamental indices. For double line notation, the quantity that we need here is the projection operator, with adjoint indices it is written as In the calculation of quark and gluon self energies, one needs the vertices for quark-antiquark-gluon(qqg) interaction, which is proportional to the generators. In the double line notation, the generators in the fundamental representation are written as Here upper pair ab denotes the adjoint index while the lower pair cd denotes the components of this matrix in the fundamental representation. Similarly, the triple-gluon vertex is proportional to structure constants which in the double line notation can be written as A. Quark loop contribution to Debye mass Generally, Debye mass (m D ) is defined through the pole of effective propagator in the static limit i.e., ω = 0, p → 0 and is related to the time-like component of gluon self-energy Π 44 (ω = 0, p → 0) [61]. It turns out that, in the presence of a static background field, apart from the usual T 2 dependent term similar to as in perturbative HTL calculations, there is an additional T 3 dependent contribution to the gluon self energy. The later component arises because the background field induces a color current which couples to the gluon. While the T 2 dependent term in Π µν is transverse (i.e., P µ Π µν (P ) = 0), the T 3 dependent term is not and spoils the transversality relation which is required for the gauge invariance. Therefore, one needs an additional contribution which may be of non-perturbative origin to the gluon self energy to cancel such a term. Similar to Ref. [59], we assume that such a term exists and cancels this undesirable T 3 term. Under these assumptions, the Polyakov loop dependent resummed propagator can be written as [59] Dµν; abcd = P L µν where µν respectively are the longitudinal and the transverse projection operators and are defined as where with x = k0 k and m 2 = (m 2 ) abcd is the thermal mass of the gluon. Under the assumptions taken here, it is clear that the pole (F) of the longitudinal propagator can be related to Π 44 component of gluon self energy. Furthermore, in the static limit, this term can be defined as Debye mass [57].
In this work, we shall focus only on the time like component of the gluon self energy with the assumption that T 3 dependent term is cancelled. For massless quarks, Debye mass has already been computed in Ref. [57]. We include here the effect of finite constituent quark mass in the quark loop contribution to the Debye mass. We work in the imaginary time formalism of thermal field theory for evaluating the corresponding diagrams. In this formalism, because of the boundary conditions of imaginary time, the energy of a fermion p 4 is an odd multiple of πT while that for a boson is an even multiple of πT . For calculating the Debye mass, we first evaluate the quark loop in the gluon self-energy for which the corresponding diagram is shown in Fig.(1), where the loop momentum four vector is written asK e µ = (K +Q e ) µ = (ω n +Q e , k) withQ e = Q e + πT . In t'hooft double line notation, the polarization tensor can be written as where aa ′ , bb ′ (e, e ′ ) are color indices of gluons (quark/antiquark), N f is quark flavor number and ∆(K , ω n = (2n + 1)πT and P 4 = ω. T r D is trace in Dirac space and Q i is the diagonal matrix in color space which is given as Q a = (−2πT q, 0, 2πT q) and q is related to the Polyakov loop expectation value as given in Eq. (18). Here, we take hard thermal loop (HTL) approximation and also assume that m ≪ T . Thus, taking HTL limit and the trace over Dirac space, Eq.(36) reduces to As we are interested in calculating Debye mass for which we need time-like component (Π 44 ) of the gluon self-energy. So from here onwards, we shall proceed with this term. For this purpose, we write the integration in Eq.(37) as The frequency sums in Eq.(38) over discrete Matsubara frequencies are somewhat involved but can be performed routinely leading to where Q2 = Q bb ′ −Q e , Q1 =Q e and f (E ± iQ) is Bose-Einstein distribution function. In Eqs. (39) and (40), the term which is independent of distribution function is the vacuum contribution which can be dropped when one considers the medium dependent terms only. First and third term in Eq. (39) contribute to the T 3 dependent term. Such a term exists only in the presence of a background gauge field in the HTL approximations [57]. As mentioned earlier, this term spoils the transversality condition and we shall not consider this undesirable contribution. Furthermore, the T 2 dependent contributions are given by second and fourth term of Eq.(39) as well as by the medium dependent term in Eq. (40). In the static limit, the time like component of the gluon self-energy can be written as where we have defined the dimensionless variables x = βk, y = βm and q1 = βQ1. Further, f (x, y, iq)'s are the Bose distribution functions in terms of these dimensionless variables as e.g., Also note that although distribution function is a complex quantity, the functions With further simplification, Π 44 can be written as where the dimensionless real functions D, F and B are In the limiting case of vanishing quark masses i.e. y = 0, the function B(q, y) do not contribute to Π q 44 as it is multiplied by a y 2 term while the functions D(q, y = 0) and F(q, y = 0) become equal and can be written in terms of Polylog functions Li 2 (z) as The Polylog function Li 2 (z) can also be written in terms of Clausen functions Cl 2 (z) e.g.
that has been used in Ref. [57]. In the present investigation, however, we will keep the effect of masses in Eqs (47), (48), (49) and integrate it numericaly to estimate the Debye mass. Generators appearing in the right side of Eq. (46) can be simplified by using projection operators, so that the product of two generators becomes Note that Π 44 depends on the color of quark and gluon and has a, b, a ′ , b ′ as free color indices. So we need to sum over other repeated color indices (i.e., e, e ′ ) which can be done by contracting color indices of Eq.(52) with that of Eq.(46). Using Eq.(52) along with Eq.(46) and summing over contracted color indices, gluon self energy can be written as

B. Gluon contribution to Debye mass
Gluon loop contribution to the gluon self energy has already been evaluated in Ref. [57]. For the sake of completeness, we recapitulate the results here. Gluon loop diagram with tri-gluon vertex is shown in Fig.(2).
Gluon loop in gluon self energy in double line notation approximation, the sum of gluon loop, four gluon vertex and ghost loop contribution to the gluon self energy can be written as As explained earlier, the time like component of the self energy is needed for the Debye mass which can be written as where Similar to quark loop, we shall not consider T 3 dependent term here and the summation over discrete Matsubara frequencies are same as in Eqs. (39) and (40). Using these summations and taking static limit, the T 2 dependent contribution to gluon self energy can be written as where Same as in the case of quark loop, for gluon loops, gluon self energy depends on the color of the gluon, and these color indices are free. Other repeated color indices can be summed by using Eq.(30) for structure constant. Thus Eq.(56) becomes To get the total Debye mass we need to add both the contribution which are given in Eqs. (53) and (58). Taking both the contributions into account, Debye mass can be given as leading to As Debye mass is color dependent and therefore, one need to sum the contributions from all the colors and then average over the number of colors to get the total Debye mass i.e., In the large N limit (i.e., neglecting 1/N terms in Eq. (60)), the Debye mass is diagonal and its components can be written in the limit quark mass m = 0 as which is same as was derived in Ref. [57] It is easy to check that, in the limit Q = 0 and m = 0, the Debye mass as written in Eq.(60) reduces to its familiar HTL limit given as In our calculation for the heavy quark transport coefficients, however, we will use the color averaged Debye mass as given in Eq.(61).

C. Light Quark thermal mass
In the double line notation, the standard diagram of one loop quark self energy is shown in Fig.(3) where a and a ′ respectively are the color indices for incoming and outgoing quark. It is expected that similar to the gluon self energy, the quark self energy also depends on the colors of incoming and outgoing quark and in the presence of a background gauge field the same can be written as where g is coupling constant, To solve the integration in Eq. Eq.(39) with product of two propagators ∆(K)∆(P − K) (arising from the term proportional to m) and another is ω n ∆(K)∆(P − K) (arising from the / K b term). The later one can be written as We take HTL approximation and evaluate only T 2 dependent term in quark-self energy. We note here that, unlike gluon self energy, one does not get any extra term different in structure as compared to the usual perturbative HTL approximation for the quark self energy. The leading contribution arises from the terms having E q − E k in the denominators of Matsubara frequency sums and in Eq.(66) comes from the first and the third terms. Simplifying Eq. (65) with Eqs. (39) and (66), quark self energy becomes In the above equation, we have used HTL approximation so that . After simplifying Eq.(67) further, it can be written as where as before, y = βm, q1 = βQ1; F(q) is same as given in Eq. (49) and J is given as It is easy to see that to estimate the quark thermal mass from its self energy, one need to sum over colors in Eq.(68) keeping a and a ′ open indices. After performing this color sum, quark self energy reduces to In the HTL approximation, the effective fermion mass (thermal mass) can be written as [60] 4m 2 th = T r( / P Σ(P )).
From Eqs. (71) and (70), the color dependent quark thermal mass a function of Polyakov loop parameter q can be written as In the limit of vanishing quark mass, using Eq.(50), it is easy to show that In the subsequent calculationS that follow, we however, keep the quark mass dependence as in Eq.(72). Similar to Eq.(61), one can define a color averaged quark thermal mass as so the total quark mass becomes Thus the color averaged Debye mass for the gluons and color averaged thermal mass for quarks as given by Eqs. (61) and (74) depend upon the Polyakov loop parameter. For the Polyakov loop parameter, we adopt here two approaches. Firstly, we estimate the same from a phenomenological 2 flavor PQM model [37,42]. The salient features of the model and the parameters taken in the model is discussed in Appendix A. With this parameterization the critical temperature for the crossover transition T c ≈ 176 MeV. We also take the Polyakov loop parameter from lattice simulations as in Ref. [58]. The variation of the Polyakov loop with temperature (T) is shown on the left of Fig(4). Clearly, compared to the lattice simulations, the Polyakov  [42]. The blue curve is from the lattice results of Ref. [58]. Right panel: Debye mass (mD) as a function of temperature. The black curve corresponds to pQCD hard thermal loop calculations [61]. The blue curve corresponds to large N limit for mD as given in Eq. (60). The green curve correspond to taking all the terms in Eq.(60) for N = 3. Here Polyakov loop is taken from lattice data [58].
loop parameter φ in PQM model shows a sharper rise and reaches its asymptotic value φ = 1 at the temperature around 320 MeV. On the other hand, in the lattice simulations, this happens at a much higher temperature. This means that the non-perturbative effects are significant up to temperature as high as 400 MeV in lattice. However in PQM these effects are significant only temperatures upto around 320 MeV. On the right side of Fig(4), Debye mass as a function of temperature is shown. Here the black curve corresponds to the Debye mass in pQCD, while the blue and the green curves are in the presence of Polyakov loop. The blue curve corresponds to the large N limit (i.e., dropping 1/N terms in Eq. (53)). On the other hand, the Green curve corresponds to including the 1/N terms in Eq. (53). Clearly, the large N limit approaches the perturbative limit faster compared to the one including 1/N terms for the Debye mass.
For the light quarks, different contributions to the masses as a function of temperature (T ) are shown on the left side of Fig. (5). Red and blue curves correspond to quark thermal masses (m th ) as given in Eq.(74) evaluated in the HTL approximation in the presence of a background gauge field. The red curve corresponds to Polyakov loop value taken from PQM model while the blue curve corresponds to the same taken from lattice simulations. The HTL perturbative QCD thermal masse as in Ref. [61] is shown by the black curve. Clearly, with the lattice value of the Polyakov loop, thermal masses approach the perturbative results at a much higher temperature while with values taken from PQM, the perturbative limit reaches at a relatively lower temperature around 320 MeV. It ought to be mentioned that beyond 330 MeV φ value is larger than one in which case q becomes imaginary. We have taken here the real part of q for estimating the thermal masses. Beyond temperature 330 MeV the real part of q vanishes which leads to the perturbative limit. As compared to PQM model, the color averaged thermal mass is smaller for Polyakov loop expectation value taken from lattice simulation. This is because, with the smaller value of φ, statistical distribution functions are suppressed more. The magenta curve is the constituent quark mass estimated in PQM model. The right side of Fig.(5), shows the behavior of color averaged Debye mass in the large N limit of Eq (60). The red and the blue curves correspond to the masses with Polyakov loop value taken from PQM and lattice simulations respectively. Debye mass is smaller as compared to the perturbative QCD Debye mass and this suppression is more when φ is taken from the lattice simulations. The reason for this is the same as that for the case of quark thermal mass. In the estimation of the transport coefficients, we shall use the Debye mass and thermal masses of quarks as in Eq.(75). It is clear that the non-perturbative effects which are in the distribution function and the masses of quarks and gluons can significantly affect these transport coefficient as compared to the perturbative QCD.

IV. RESULTS AND DISCUSSIONS
With the thermal mass of the quarks and the Debye mass as computed in the background of a nontrivial Polyakov loop, we next numerically compute the drag and diffusion coefficients using Eq. (15). For the heavy quark elastic interaction with the light quarks and gluons, qQ → qQ and gQ → gQ scattering processes are considered where Q stands for heavy quark, q stands for light quarks and g stands for the gluon. In the case of massless light quark and gluon, the leading-order (LO) matrix elements for qQ → qQ and gQ → gQ scattering have been calculated in Ref. [45,52]. These pQCD cross sections have to be supplemented by the value of the coupling constant and the Debye screening mass which is needed to shield the divergence associated with the t-channel diagrams to compute the heavy quark transport coefficients. For massive light quark and gluon, the calculation of the scattering matrix, M (q,g)+Q→(q,g)+Q , is performed considering the leading-order (LO) diagram with massive quark and gluon propagators for gQ → gQ and a massive gluon propagator for qQ → qQ scatterings [21,23]. Within the matrix model, the scattering amplitudes are summarised in Appendix(B). Similar to previous work [21,23], massive gluon propagator for qQ → qQ and t-channel of gQ → gQ is used. We estimate the transport coefficients for the charm quark whose mass is taken as m C = 1.27 GeV. Here we use the two loop running coupling constant given as [62] α s = 1 4π with Λ = 260 MeV and N f = 2.
We evaluate the drag and diffusion coefficients of heavy quark in QGP with Polyakov loop value from two different models. In one case the Polyakov loop value, hence the Debye mass and thermal masses, has been taken from PQM calculation as inputs to compute the heavy quark transport and we label it as PQM. In the other case, Polyakov loop value has been taken from the lattice simulatons and hence, we label it as lattice in the following discussions. The temperature variation of the drag coefficient has been shown in Fig.(6) for charm quark interaction with light quarks and gluon for a given momentum (p=0.1 GeV) obtained for both PQM and lattice Polyakov loop values. We obtain quite a mild temperature dependence of heavy quark drag coefficient for the case of PQM. However, with lattice, we obtained a quite stronger temperature dependence of heavy quark drag coefficient than the one with PQM. We notice that the drag coefficient obtained with PQM input is larger at low temperature than the one obtained with lattice inputs whereas the trend is opposite at high temperature. This is mainly because of the interplay between the Debye mass and Polyakov loop value obtained within both the models.
A smaller value of the Polyakov loop, as shown on the left side of Fig.(4), in case of lattice reduces the magnitude of the drag coefficients at low temperature. However, at high temperature, with smaller Debye mass, as shown in Fig.(5), obtained with lattice input enhances the magnitude of heavy quark drag coefficients. Hence, at low temperature Polyakov loop value plays the dominant role ( e.g., at T=180 MeV, the Polyakov loop value obtained within both the models differ by a factor about 2) whereas at high temperature the Debye mass plays the dominant role ( e.g., T=300 MeV, the differences between the Polyakov loop value obtained within both cases reduced significantly) for the behavior of the drag coefficient. We observed temperature dependence of heavy quark drag coefficient obtained with PQM Polyakov loop value is quite consistent with the results obtained with other quasi-particle models [20,23] and T-matrix approach [18]. It is important to mention that the temperature dependence of the drag coefficient plays a significant role [20] to describe heavy quark R AA and v 2 simultaneously, which is a challenge to almost all the models on heavy quark dynamics. A constant or weak temperature dependence of the drag coefficient is an essential ingredient to reproduce the heavy quarks R AA and v 2 simultaneously, whereas in pQCD the drag coefficient increases with temperature.
The momentum variation of the drag coefficient has been shown in the right panel of Fig.(6) for charm quark interaction with light quarks and gluon obtained with PQM and lattice Polyakov loop value. We observe a strong momentum dependence of heavy quark drag coefficient as compared to the same estimated within pQCD [7]. This is mainly due to the inclusion of non-perturbative effects through the Polyakov loop background. At T=300 MeV the drag obtained with the PQM Polyakov loop (at p=0.1 GeV) is marginally larger than the drag obtained with lattice Polyakov value. Hence, the momentum variation of drag coefficients obtained with inputs from PQM is marginally larger than the one obtained with inputs from lattice simulation in the entire momentum range considered here.
In Fig.(7) heavy quark diffusion coefficient B 0 has been displayed as a function of temperature obtained with input parameter from PQM and lattice. The diffusion coefficients increases with temperature for both the cases as it involves the square of the momentum transfer. In terms of magnitude the diffusion coefficient obtained within both the cases follow similar trend of drag coefficient due to the same reason (i.e., interplay between Debye mass and Polyakov loop value). The momentum variation of the diffusion coefficient has been shown in Fig.(7) for charm quark interaction with light quarks and gluons for the same values of Polyakov loop. Similar to the drag coefficient, the diffusion coefficient also shows the same trend with PQM having larger value then that from the lattice as a function of momentum. Stronger suppression of distribution function at high momentum in lattice Polyakov loop than that of from PQM also play a marginal role in the momentum variation of heavy quark drag and diffusion coefficients obtained.
To understand the temperature dependence of the transport coefficients, we plot the temperature variation of the drag coefficient in Fig. 8 at different momentum obtained with Polyakov loop value from lattice simulations. We obtain almost similar temperature dependence of heavy quark drag coefficient at both the momentum having larger magnitude at p=2 GeV than at p=5 GeV. In Fig. 9 we have depicted the temperature variation of diffusion coefficient at different momentum for the same values of Polyakov loop. As expected, the magnitude of the diffusion coefficient is large at p=5 GeV than p=2 GeV having similar temperature variation for both the momenta.
In Fig. 8 we have shown the variation of drag coefficient with momentum at different temperature obtained with the lattice inputs. We observe a larger magnitude of the drag coefficient at T=320 MeV than T=200 MeV but the momentum variation is similar at both temperature. Momentum variation of the diffusion coefficient has been depicted in Fig. 9 at difference temperature. At both the momenta the diffusion increase with temperature having larger magnitude at T=320 MeV than T=200 MeV.
It is worth mentioning here that, non-perturbative effects from a different perspective has been investigated recently in Ref. [63][64][65] and employed to calculate the transport coefficients [65]. The method here consisted of using T-matrix with an in-medium potential for the heavy quarks. This potential is constrained by the heavy quark free energy from the lattice data. The lattice heavy quark free energy is directly related to the Polyakov loop and hence is correlated with the strength of the confining potential. Therefore it is nice to see that the behavior of drag coefficient being rather flat with regards to temperature dependence whereas the diffusion coefficient having a strong temperature dependence as observed here was also observed in Ref. [65]. This consistency suggest of having a possible existence of model independent correlation between Polyakov loop and the heavy quark transport coefficients.

V. SUMMARY
In this work, we have computed the heavy quark drag and diffusion coefficients in QGP including non-perturbative effects via a Polyakov loop background. In order to incorporate these effects we first calculate quark and gluon thermal masses also taking the quark constituent mass into account. We found that for temperatures below 300 MeV quark thermal mass and gluon Debye mass starts deviating from its perturbative value this effect significant for even higher temperatures when Polyakov values are taken from the lattice simulations. This decrease in the Debye mass of gluon and the thermal mass of light quarks is due to color suppression manifested in the quark and gluon distribution functions in the presence of a background Polyakov loop field. In the calculation of HQ diffusion coefficient the distribution function of the light quark and the Debye mass play complimentary roles. While the distribution function with Polyakov loop tend to decrease the HQ transport coefficient the Debye mass has the effect of increasing these transport coefficients. We have found a weak temperature dependence of the heavy quark drag coefficient with Polyakov loop value taken from PQM which is consistent with other models like T-matrix and quasi particle model which also take into account the non-perturbative effects in a different manner. This consistency suggests existence of possible model independent correlations between the results obtained with the Polyakov loop and other non-perturbative models and reaffirm the temperature and momentum dependence of heavy quark transport coefficients. In the present investigation, we have aconfined our attention to the elastic 2 → 2 processes within the matrix model. Inclusion of other effects arising from 2 → 3 processes, LPM effects are expected to be sub-dominant due to the large mass of the heavy quark [66] but, none the less, can be important at high parton density. We plan to explore different possible phenomenological implications of the present investigation in future.

Acknowledgment
In the above, the first term is the kinetic and interaction term for the quark doublet ψ = (u, d) interacting with the scalar (σ) and the isovector pseudoscalar pion (π) field. The scalar field σ and the pion field π together form a SU (2) isovector field. The quark field is also coupled to a spatially constant temporal gauge field A 0 through the covariant The mesonic potential U χ (σ, π) essentially describes the chiral symmetry breaking pattern in strong interaction and is given by The last term in the Lagrangian in Eq.(A1) is responsible for including the physics of color confinement in terms of a potential energy for the expectation value of the Polyakov loop φ andφ which are defined in terms of the Polyakov loop operator which is a Wilson loop in the temporal direction In the Polyakov gauge A 0 is time independent and is in the Cartan subalgebra i.e. A a 0 = A 3 0 λ 3 + A 8 0 λ 8 . One can perform the integration over the time variable trivially as path ordering becomes irrelevant so that P(x) = exp(βA 0 ). The Polyakov loop variable φ and its hermitian conjugateφ are defined as In the limit of heavy quark mass, the confining phase is center symmetric and therefore φ = 0 while for deconfined phase φ = 0. Finite quark masses break this symmetry explicitly. The explicit form of the potential U p (φ,φ) is not known from first principle calculations. The common strategy is to choose a functional form of the potential that reproduces the pure gauge lattice simulation thermodynamic results. Several forms of this potential has been suggested in literature. We shall use here the following polynomial parameterization [37] with the temperature dependent coefficient b 2 given as The numerical values of the parameters are a 0 = 6.75, a 1 = −1.95, a 2 = 2.625, a 3 = −7.44 b 3 = 0.75, b 4 = 7.5 (A7) (A8) The parameter T 0 corresponds to the transition temperature of Yang-Mills theory. However, for the full dynamical QCD, there is a flavor dependence on T 0 (N f ). For two flavors we take it to be T 0 (2) = 192 MeV as in Ref. [37]. The Lagrangian in Eq.(A1) is invariant under SU (2) L ×SU (2) R transformation when the explicit symmetry breaking term cσ vanishes in the potential U χ in Eq.(A2). The parameters of the potential U χ are chosen such that the chiral symmetry is spontaneously broken in the vacuum. The expectation values of the meson fields in vacuum are σ = f π and π = 0. Here f π = 93 MeV is the pion decay constant. The coefficient of the symmetry breaking linear term is decided from the partial conservation of axial vector current (PCAC) as c = f π m 2 π , m π = 138 MeV, being the pion mass. Then minimizing the potential one has v 2 = f 2 π − m 2 π /λ. The quartic coupling for the meson, λ is determined from the mass of the sigma meson given as m 2 σ = m 2 π + 2λf 2 π . In the present work we take m σ = 600MeV which gives λ=19.7. The coupling g σ is fixed here from the constituent quark mass in vacuum M q = g q f π which has to be about (1/3)rd of nucleon mass that leads to g σ = 3.3 [53].
To calculate the bulk thermodynamical properties of the system we use a mean field approximation for the meson and the Polyakov fields while retaining the quantum and thermal fluctuations of the quark fields. The thermodynamic potential can then be written as The fermionic part of the thermodynamic potential is given as modulo a divergent vacuum part. In the above, ω ∓ = E p ∓ µ, with the single particle quark/anti-quark energy E p = p 2 + M 2 . The constituent quark/anti-quark mass is defined to be The divergent vacuum part arises from the negative energy states of the Dirac sea. Using standard renormalisation, it can be partly absorbed in the coupling λ and v 2 . However, a logarithmic correction from the renormalisation scale remains which we neglect in the calculations that follow [53]. The mean fields are obtained by minimizing Ω with respect to σ, φ,φ, and π. Extremising the effective potential with respect to σ field leads to where, the scalar density ρ s = − ψ ψ is given by In the above, f ∓ (p) are the distribution functions for the quarks and anti quarks given as and, The condition ∂Ω ∂φ = 0 leads to where , Similarly, ∂Ω ∂φ = 0 leads to with, Iφ = ∂Ωq q ∂φ = −6N f T dp (2π) 3 e −2βω− 1 + 3φe −βω− + 3φe −2βω− + e −3βω− + e −βω+ 1 + 3φe −βω+ + 3φe −2βω+ + e −3βω+ , (A19) By solving Eqs.(A12),(A16) and (A18) self consistently one can get the values of constituent quark mass, Polyakov loop variable and the conjugate Polyakov loop variable as a function of temperature.

Appendix B: Scattering amplitudes
There are two types of scatterings that contribute to the drag and the diffusion coefficients namely Coulomb scattering i.e., scattering off of HQ from light quark and Compton scattering i.e., scattering off of gluon from HQ [45]. The dominant contribution for these scatterings arise from the gluon exchange in the t-channel which is infrared divergent [16,46]. This is regularised by introducing the Debye screening [16,45] which we have evaluated in the HTL limit in the background of Polyakov loop. In s and u channel, however, there is no such infrared divergences. Note that in the matrix model, m D , in Eq.(60) is color dependent so the propagator is also color dependent. For N f flavor of light quark, the spin averaged matrix element squared for Coulomb scattering as shown on the left side of Fig.(10), can be written as . (B1) where M is HQ mass. For the q a Q b → q e Q f scattering, the product of distribution function and matrix element squared that appears in Eq. (15) can be simplified by summing over color of initial and final light/heavy quarks. Note that for light quarks, the colors appearing in Eq.(B3) has to be summed with the distribution function and can be written as where f (q) q is the average distribution function of quark as defined in Eq. (22). Similarly for the t channel Compton scattering shown on the right side of Fig.(10), one can write .
( Eq.(11) can be written as Note here that the propagator has no color dependent term. Matrix element squared for s channel Compton scattering as shown on the left side of Fig.(11) is |M s | 2 = 8g 4 8(N 2 − 1) P ef bc P ef bc ′ P gh ca P gh c ′ a M 4 − us + M 2 (u + 3s) (s − M 2 ) 2 . (B7) There are interferences between different scatterings contributing to g ef Q b → g gh Q a that can be written as 2(N 2 − 1) P ef bc P gh ca P lm ab (if dc,f e,hg ) −8(M 4 − 2M 2 s + us) (s − M 2 )(t + (m 2 D ) mlcd ) .

(B10)
Total matrix element squared that contribute to Compton scattering i.e., gQ → gQ is |M Cm | 2 abef gh = |M s | 2 +|M u | 2 These matrix elements are used in Eq.(15) to estimate the drag and the diffusion coefficient.