Estimation of the diffusion coefficient of Heavy Quarks in light of Gribov-Zwanziger action

The heavy quark momentum diffusion coefficient ($\kappa$) is one of the most essential ingredients for the Langevin description of heavy quark dynamics. In the temperature regime relevant to the heavy ion collision phenomenology, a substantial difference exists between the lattice estimations of $\kappa$ and the corresponding leading order (LO) result from the hard thermal loop (HTL) perturbation theory. Moreover, the indication of poor convergence in the next-to-leading order (NLO) perturbative analysis has motivated the development of several approaches to incorporate the non-perturbative effects in the heavy quark phenomenology. In this work, we estimate the heavy quark diffusion coefficient based on the Gribov-Zwanziger prescription. In this framework, the gluon propagator depends on the temperature-dependent Gribov mass parameter, which has been obtained self-consistently from the one-loop gap equation. Incorporating this modified gluon propagator in the analysis, we find a reasonable agreement with the existing lattice estimations of $\kappa$ within the model uncertainties.


I. INTRODUCTION
Since the seminal work by Svetitsky [1], there has been a concerted effort to understand the heavy quark dynamics in the presence of the evolving background medium produced in the relativistic heavy-ion collision experiments [2][3][4][5]. An intriguing feature associated with the heavy quarks is the inherent distinction of their dynamics compared to the light quark sector. With orders of magnitude differences in their masses, the heavy quarks can only achieve a slower equilibration rate relative to their lighter counterparts, rendering their thermalization time comparable to or even larger than the lifetime of the fireball. As a consequence, an imprint of the interaction history is expected to be retained in the final state heavy quark spectra [6], which makes them an important probe for the characterization of the strongly interacting background. On the other hand, the typical momenta carried by the heavy quarks being much larger than the ambient temperature scale, a single collision with the medium constituents can hardly change the momentum significantly. Hence, a Brownian motion with successive uncorrelated momentum kicks has been a well-accepted description of the dynamics of the heavy quarks. The corresponding Langevin equation governing the heavy quarks' momentum evolution possesses a drag force whose strength is determined by the associated drag coefficient ( ). Whereas the random momentum kicks are incorporated in the momentum evolution through a stochastic force term ( ) having a correlation ( ) ( ) proportional to ( − ). The proportionality constant is related to the mean squared momentum transfer per unit time and is known as the momentum diffusion coefficient [7]. In equilibrium, the drag coefficient is related to the momentum diffusion coefficient through the Einstein relation = /(2 ) where corresponds to the mass of the propagating heavy quark and is the temperature of the medium. At vanishing momentum transfer, the two coefficients can further be related to the spatial diffusion coefficient as = /( ) = 2 2 / . This simplification allows one to conveniently encode the medium characteristics in a single transport coefficient conventionally defined as (2 ) , where the spatial diffusion coefficient is scaled with the thermal wavelength. Within the perturbative approach, the heavy quark diffusion coefficient has been obtained in the LO [1,7] as well as in the NLO [8] with the strong coupling. From those studies, it has been observed that for realistic coupling strength, the NLO corrections become several times larger than the LO estimates, reflecting the poor convergence of the perturbation series. Hence a non-perturbative estimation is crucial to constrain the temperature dependence of the diffusion coefficient, which has important implications in the heavy ion phenomenology [9]. Since the first attempt [10] based on the definition obtained in Ref. [11], there have been several non-perturbative estimations for the diffusion coefficient using the lattice discretized QCD approach (see for example Refs. [12][13][14][15][16] ) which significantly enriched our knowledge of the temperature dependence of the diffusion coefficient. Though large uncertainties exist in the lattice estimates, their consistency within the error reconfirms the poor convergence of perturbative estimation. Although it should be noted here that these non-perturbative estimations of are restricted to the quenched approximation in the lattice framework, which assumes a gluonic plasma and hence does not represent the realistic scenario expected in the heavy ion collision experiments. However, the quenched approximated results are of similar magnitude as expected in the temperature regime relevant in the phenomenology (see, for example, Refs. [3,17] for a comprehensive comparison between the lattice results and the estimations from different phenomenological models as well as data-driven approaches [18]). At this point, it is interesting to note that the estimation of the heavy quark transport coefficient in the phenomenologically relevant regime of temperatures is not the only scenario where the perturbative estimation differs from the lattice results obtained with pure gluonic plasma. For example, in the temperature regime close to ∼ 270 MeV, a significant deviation can be observed when the free energy data obtained from the pure Yang-Mills lattice simulation [19] is compared with the resummed perturbative theory estimates [20]. In this context, a remarkable improvement over the resummed perturbative approach has been achieved in Ref. [21] implementing the Gribov-Zwanziger (GZ) prescription [22,23] in the evaluation of the free energy.
The Gribov-Zwanziger formalism [24,25] improves the infrared behaviour of QCD by fixing the residual gauge transformations that remain after the implementation of the Faddeev-Popov quantization procedure. A review on the Gribov prescription can be found in Refs. [26,27] (for recent progress on the refined Gribov approach, see also [28][29][30][31][32][33] and the references therein). Within the GZ framework, the gluon propagator in the Landau gauge is expressed as [21] where is the Gribov parameter. This parameter in the denominator shifts the poles of the gluon propagator off the energy axis to an unphysical location 2 = ± 2 , suggesting that the gluons are not physical excitations. In the context of finite temperature Yang-Mills (YM) theory, the temperature dependence of the Gribov parameter can be determined by solving a gap equation self-consistently. As a result, one obtains a nonperturbative IR cut-off in the gluon propagator, which is absent in the traditional YM approach where the finite temperature perturbative expansion breaks down at the (chromo)magnetic scale 2 , known as the Linde problem [35]. It was first pointed out in Ref. [36] that, in the limit of asymptotically high temperatures, becomes proportional precisely to the (chromo)magnetic scale 2 indicating a promising resolution to such IR-catastrophe. Recently, the Gribov modified gluon propagator with this asymptotic temperature dependence has been implemented in the hard-thermal-loop (HTL) effective theory to study the quark dispersion [34] relations, the dilepton rate, and the quark number susceptibility [37], which revealed some intriguing physical properties of the said observables. Thus, it is interesting to investigate whether the effects from the magnetic scale incorporated via the Gribov-Zwanziger prescription can improve the leading order perturbative estimation of the heavy quark diffusion coefficient in the phenomenologically relevant temperature regime. In the present work, we address this issue considering two scenarios: first, we consider pure gluonic plasma as the background medium for the heavy quarks. In this case, the temperature dependence of the Gribov parameter is extracted using the one-loop gap equation. Next, we consider the heavy quarks in the QGP background with only the asymptotic temperature dependence of . It turns out that the corresponding momentum diffusion coefficient in both scenarios shows significantly different temperature dependence compared to the leading order perturbative result.
The letter is organized as follows. At first in section II, the momentum diffusion coefficient is obtained from the elastic scattering amplitude in the non-relativistic limit implementing the Gribov-modified gluon propagator. The obtained diffusion coefficient depends explicitly on the temperature-dependent Gribov parameter ( ), which is subsequently determined in section III. Next, in section IV, the numerical estimation of the scaled diffusion coefficient is compared with the available lattice data and the perturbative results. Finally, we summarize the main results in section V and conclude with a brief discussion of the possible future directions.

II. FORMALISM
The Langevin equation corresponding to the momentum evolution of the non-relativistic heavy quarks in a background medium is given by [7] = − + ( ) , Here the momentum diffusion coefficient can be obtained from the mean squared momentum transfer per unit time, which is 3 and which can be computed from the expression [7]: where |M | 2 quark and |M | 2 gluon represent the squared 2 ↔ 2 scattering matrix elements corresponding to → and → processes respectively. Note that, here, the and respectively correspond to the incoming or outgoing massless light quarks (as well as the anti-quarks), and gluons of the background medium with four-momentum denoted as ≡ ( 0 , ) or . Unlike the heavy quarks (represented by ), these medium constituents possess the corresponding statistical weight factors, which are expressed in terms of the Fermi-Dirac ( ) and the Bose-Einstein ( ) distribution functions, respectively, which, in the massless limit, depending on the magnitude of the three momentum vectors denoted as or . Moreover, in the above expression, the heavy quarks are considered non-relativistic, which essentially makes the energies of the incoming ( ) and the outgoing ( ) heavy quarks identical to their rest mass energy . Now, in the leading order, the squared scattering matrix element for the → process via the -channel gluon exchange, together with the anti-quark contribution, can be expressed as where is the color Casimir constant, is the gluon propagator that depends on the momentum transfer = − = − . Note that the squared matrix amplitude has been summed and averaged over the final and the initial quantum numbers (spins and colors), respectively, and multiplied with the degeneracy factor of the initial light quark (2 ) and the number of flavors ( ) along with a factor of 2 for the antiquark contribution. The Dirac traces over the quark external lines along with the metric and in the gluon propagators provide Now, considering the heavy quark approximately at rest, i.e. = ( , 0), becomes purely spatial, and we have, = − , which gives the following simplified relations between the momenta: Note that, as the heavy quarks are considered to be approximately at rest, we essentially need to incorporate only the temporal part of the gauge boson propagator, i.e., 00 ( ) in the squared matrix element. Now, in the standard approach, to incorporate the screening in the t-channel amplitude, the free gluon propagator is replaced with the HTL resummed one loop effective propagator, which, in the limit of small energy transfer, is equivalent to introducing the Debye mass ( ) in the denominator [7]. However, in the present case, we consider the 00 ( ) component of the modified gluon propagator within GZ approach, which, along with the fact that the energy transfer is negligible (i.e., is purely spatial), gives Incorporating this modified propagator in Eq. (4) along with the simplified relationships given in Eq. (6) and (7), we obtain the final expression for the squared matrix element for the → channel as The scattering matrix for the → process involves the three gluon vertex Γ ( , , ) in the t-channel. The complete expression of the squared matrix element considering the s, t, and u channel processes can be found in Refs. [38,39]. As argued in [7], in the rest frame of the heavy quark, the matrix element is dominated by the t-channel gluon exchange, and the Compton amplitude can be safely ignored. Incorporating the modified gluon propagator in the t-channel, we obtain where an additional degeneracy factor 2( 2 −1) corresponding to the initial gluon has been multiplied with the spin and coloraveraged (summed) squared amplitude. Now, the obtained squared matrix elements from Eqs. (9) and (10) can be inserted within the integrals given in Eq. (3). It is convenient to shift the integral over to an integral over the momentum transfer = − and express the diffusion coefficient as Here we note that the squared matrix amplitudes depend only on the magnitude of the integration variables and as their angular dependence can be re-expressed using Thus, the integral over the three-momentum delta function can be easily performed, and we obtain Next, the energy delta function can be integrated over the variable cos by using the property which also restricts the upper limit of the integral to max = 2 . After this angular integral, we obtain where the dependence only appears in the factor where¯= ( )/ and the integral is represented in terms of the dimensionless variable = / . Note that, similar to the Debye mass, the Gribov parameter itself is a function of temperature, and it plays a crucial role in determining the overall temperature dependence of the diffusion coefficient, as we will see later. Let us first obtain the diffusion coefficient in the limit of small¯. In that case, it can be compared with the standard / expanded results, which, including the NLO corrections, are given by [8] where the NLO correction is incorporated in the last term. The LO result in Eq. (18) can be obtained from Eq. (15) by simply replacing the factor ( , ) with the Debye screened expression ( , ) given by and performing the subsequent and integration. We can further establish the correlation between the LO Debye screened case explored earlier [7] and our results by denoting the LO expansions as (3 ) LO and (3 ) respectively for the Debye screened case and the Gribov case. A straightforward calculation shows : where, the perturbative definition of 2 = (1/3) 2 2 ( + /2) is used in the last term, which essentially arises from the difference in the momentum dependence of the gluon propagator compared to the Debye screened propagator in the standard approach. At asymptotically high temperatures, this additional contribution becomes ∼ 2 (as corresponds to the (chromo)magnetic scale 2 ) whereas the leading log dependence i.e. ∼ 4 log( / ) now becomes ∼ 4 log( / ), thus, remains ∼ 4 log(1/ ). However, for phenomenologically relevant temperatures, one may require a non-perturbative estimation of the temperature dependence of the Gribov parameter, as will be discussed in the following.

III. FIXING THE GRIBOV PARAMETER
As mentioned earlier, the Gribov prescription improves the Faddeev-Popov quantization procedure by modifying the YM partition function in the Euclidean space-time as [27] = ∫

[D ] [D¯] [D ] (Ω) ( )
where, is the covariant derivative, and¯are the ghost, and anti-ghost fields, and the no pole condition is implemented by the step function (Ω) = (1 − 0 ) that constraints the path integral over the field configurations within the Gribov region Ω defined as The function 0 = ( = 0) appears in the inverse of the ghost dressing function as −1 = 1 − ( ) corresponding to the ghost field propagator given by The standard procedure to obtain the momentum space gluon propagator ( ) ( ) from the path integral as given in Eq. (21), is to re-express the constraint as and consider the steepest descent approximation for the integral over the variable . Introducing a mass parameter corresponding to the specific value = 0 that minimizes the entire exponent of the integral, one obtains the gap equation for the Gribov parameter (see Ref. [27] for details). In the case of finite temperature, the gap equation is given by where is the space-time dimension, and the temporal component of the Euclidean four-momentum = ( , 4 ) possesses the discreet Matsubara frequencies 4 = 2 . It is straightforward to perform the frequency sum and evaluate the momentum integral using dimensional regularization. In the MS renormalization scheme with = 4, one obtains [40] where 0 represents the regularization scale and the argument in the distribution function ( ) = ( / − 1) −1 corresponds to the dispersion relations, which in this case are, ± = √︃ 2 ± 2 . On the other hand, the sum integral corresponding to the ( ) in the ghost dressing function is given by [21], To obtain the temperature dependence of the Gribov parameter from the gap equation given in Eq. (26), it is essential first to fix the scale 0 that appears in the temperature-independent part. Following Ref. [21] for this purpose, we first consider the = 0 contribution of the gap equation given by with 0 as the strong coupling at vanishing temperature. The above expression of 0 is then substituted in the temperatureindependent contribution of the inverse ghost dressing function, which can be obtained from Eq. (27) which provides the 0 as a function of the coupling strength 0 . Considering 0 = 3.13, one obtains 0 = 1.69 GeV, which will be used in the finite temperature gap equation along with the temperature-dependent running coupling ( ). In Ref. [21], the authors have parameterized the non-perturbative running coupling with a single parameter (say ) by utilizing the functional form of the one-loop perturbative coupling with = 0. It is shown that this simple, functional form.
can fit the lattice-measured coupling data [41] reasonably well for a wide range of temperatures. The fitted parameter values corresponding to the coupling data extracted from the large distance (IR) and the short distance (UV) behaviour of the heavy quark free energy are given by IR = 1.43 and UV = 2.97 respectively. In the present work, we consider this parametrized running coupling 2 ( / ) = 4 par ( / ) in the gap equation to numerically extract the temperature dependence of the Gribov parameter for the gluonic plasma ( = 0). The obtained temperature dependence of the scaled Gribov parameter corresponding to the IR and the UV parametrization is shown in Fig. 1 [21]. The uncertainty band, in this case, corresponds to the variation of the parameter from IR to UV . Note that the temperature dependence of the scaled Gribov parameter in the IR case is similar to the decreasing trend obtained using the kinetic theory framework considering the quasiparticle picture of the Gribov plasma [42][43][44]. However, unlike the IR scenario, the obtained¯for = UV remains small throughout the considered temperature range. Moreover, as discussed in Ref. [21], the fitted coupling, along with the obtained 0 , provide an almost temperature-independent ghost dressing function (as obtained from Eq. (27) numerically) which is consistent with the lattice results [45,46]. However, it should be mentioned that in the original GZ approach, the dressing function shows an IR-enhancement at → 0 whereas the recent lattice measurements at vanishing momentum support a finite decoupling behaviour which is in better agreement with the refined GZ approach [28,29]. Nevertheless, we adopt the original GZ procedure here due to its considerable simplicity compared to the refined GZ approach, especially for finite temperature applications.
It is also interesting to compare the temperature dependence of the Gribov-modified diffusion coefficient with the perturbative result, including the quark contribution. For that purpose, we utilize the asymptotic temperature dependence of the Gribov parameter that can be obtained from the gap equation considering the → ∞ limit as Thus, the temperature dependence of¯, in this case, is completely determined by the running coupling. Keeping in mind the functional form of the parametrized coupling for = 0, in this case, we use the one-loop running coupling with = 3, which is given by where, Λ MS = 0.176 GeV is obtained from the lattice measurement (1.5GeV, = 3) = 0.326 [47] and we consider = 0.16 GeV to express the running coupling in terms of the scaled temperature / . The uncertainty band, in this  Figure 2. The spatial diffusion coefficient (2 ) as a function of the scaled temperature is compared with the lattice data from Refs. [10,[12][13][14]. At = 1.5 , a horizontal shift is introduced in the data for visual clarity. The bands refer to the variation of the scale parameter from = 1 to 4 and IR to UV for the LO/NLO and the Gribov case, respectively. case, is obtained by varying the energy scale symmetrically around 2 by a factor of two (i.e. 1 ≤ ≤ 4). As can be observed from Fig. 1, with = 3, where the asymptotic form of the Gribov parameter is considered, the temperature dependence becomes similar to the UV case obtained previously with = 0. In the following, we consider the two scenarios with = 0 and three and discuss the influence of the temperature dependence of the Gribov parameter on the estimation of the heavy quark diffusion coefficient.

IV. NUMERICAL ESTIMATION OF THE DIFFUSION COEFFICIENT
To estimate the temperature dependence of the diffusion coefficient in the GZ framework, one has to perform the momentum integral in Eq. (17) numerically. Let us first consider the pure gluonic background. In that case, the non-trivial temperature dependence of the Gribov mass parameter (see Fig.  1) is obtained by numerically solving Eq. (26) along with the parametrized running coupling given in Eq. (30). The scaled spatial diffusion coefficient, (2 ) , can be obtained from the momentum diffusion coefficient ( ) by using the Einstein relation (2 ) = (4 3 )/ . The variation of the obtained (2 ) as a function of the scaled temperature / is shown in Fig. 2 along with the lattice estimations from Refs. [10,[12][13][14]. The uncertainty band, in this case, refers to the variation of the scale ( ) obtained by varying the fit parameter in the range IR to UV . The corresponding values of the Gribov parameter (as shown in Fig.1) in each of the cases are obtained by solving the gap equation given in Eq. (26). It should be mentioned here that because of this variation of the Gribov parameter with , the error band in the diffusion coefficient deviates from the trivial uncertainty expected from the overall multiplicative factor (see, for example, Eq. (17) with the overall 4 factor). Similar reasoning also ap- = 0 and the lattice data from [10] and [13] are also shown for comparison.
plies to the known LO and NLO expressions where the Debye mass varies with the parameter . As can be noticed from Fig. 2, the temperature dependence of the Gribov-modified diffusion coefficient shows a reasonable agreement with the lattice estimations within the uncertainties. Moreover, the linearly increasing nature of the scaled spatial diffusion coefficient is consistent with the data-driven approach [18,48] where a linearly increasing parametrization is usually implemented for the non-perturbative soft-part to compare the model output to the experimental data of the heavy-meson nuclear modification factor and the elliptic flow 2 . As shown in Fig. 2, the Gribov-modified approach significantly corrected the LO result. However, the behaviour of the Gribov improved spatial diffusion coefficient is qualitatively similar to the result obtained incorporating the NLO correction, which is also shown in Fig. 2 for comparison. In this case, the approximate expression for the two-loop perturbative running [49] is used in Eq. (18) along with /Λ MS = 1.15 [50] and the uncertainty band is obtained by varying the energy scale symmetrically around 2 by a factor of two. An overlap between the uncertainty bands of the NLO result and the Gribov-modified diffusion coefficient estimation can be observed throughout the considered range of temperatures. Nevertheless, one can notice a relatively narrower uncertainty band in the GZ case, which essentially stems from the uncertainties in the lattice determination of the running coupling at finite temperatures [41].
It is also interesting to compare the NLO result and the Gribov modified estimation, including the quark sector contribution. For this purpose, we consider the perturbative one loop running coupling given in Eq. (32) (with the energy scale ( ) varied considering 1 ≤ ≤ 4) along with the simple asymptotic temperature dependence of the Gribov parameter (proportional to the chromo-magnetic scale) as given in Eq. (31). Note that the asymptotic temperature dependence of the Gribov parameter has been considered earlier in the studies of the fermion dispersion relation [34] and the dilep-ton production rate [37] in the GZ framework. In the present case, we incorporate this asymptotic form of the Gribov parameter in Eq. (17) to obtain the corresponding temperature dependence of the scaled momentum diffusion coefficient as shown in Fig 3. For comparison, the Gribov modified result with = 0 is also shown along with the pure glue lattice data from Refs. [10,13]. It can be observed that the asymptotic form of the Gribov parameter with = 3 gives rise to an almost similar temperature dependence as obtained by solving the gap equation with the parametrized coupling in the = 0 case. Here also, an overlap between the uncertainty bands corresponding to the NLO correction and the Gribovmodified results can be observed throughout the considered temperature range. However, concentrating on a fixed energy scale 2 for both cases, one can infer that, close to the critical temperature, the Gribov modified analysis indicates a smaller momentum diffusion coefficient value than the NLO estimation.

V. SUMMARY AND OUTLOOK
In this work, the heavy quark diffusion coefficient up to the leading order in 1/ expansion (where represents the mass of the heavy quark) has been studied in the phenomenologically relevant temperature regime using the Gribov-Zwanziger prescription. The relevant t-channel matrix amplitude at low momentum transfer has been obtained with the IR-suppressed gluon propagator that depends on the Gribov parameter. For pure gluonic medium, the temperature dependence of this parameter has been obtained by solving the gap equation with a parametrized running coupling adopted from Ref. [21]. On the other hand, for = 3, we consider the asymptotic temperature dependence of the Gribov parameter, which is proportional to the chromo-magnetic scale. The temperature dependencies so obtained have been implemented in the estimation of the momentum diffusion coefficient based on the kinetic theory framework, and the corresponding spatial diffusion coefficient is obtained using the Einstein relation. We find that the Gribov-modified prescription significantly improves the LO estimation of the spatial diffusion coefficient resulting in a linearly increasing temperature dependence consistent with the lattice estimations as well as the data-driven approaches. For both = 0 and = 3 cases, the estimated diffusion coefficients in the Gribov approach and the NLO perturbative approach overlap within the uncertainties (though with reduced uncertainties in the former approach). Considering the median of the uncertainty bands with = 3, we find that, near the transition temperature, the GZ approach results in a smaller momentum diffusion coefficient (closer to the value obtained in the pure glue scenario) compared to the NLO estimation.
Note that the present study implements the original GZ approach for the simplicity of incorporating finite temperature modifications. A similar analysis based on the refined GZ approach [51] serves as an interesting future direction. Nevertheless, the significant improvement over the LO estimation of the diffusion coefficient obtained here (even with the asymptotic temperature dependence of the Gribov parameter) definitely encourages one to investigate the importance of the (chromo)magnetic scale in the estimation of other relevant transport properties of the medium [42,44] as well as in the transport studies in the presence of non-trivial background [52,53].