# Shear and bulk viscosity of quark-gluon plasma with Gribov gluons and
quasiparticle quarks

Sadaf Madnia<EMAIL_ADDRESS>Arghya Mukherjeeb<EMAIL_ADDRESS>Amaresh Jaiswala<EMAIL_ADDRESS>Najmul Haquea<EMAIL_ADDRESS>aSchool of Physical Sciences, National Institute of Science Education and
Research, An OCC of Homi Bhabha National Institute, Jatni-752050, India
bRamakrishna Mission Residential College (Autonomous), Narendrapur,
Kolkata-700103, India

###### Abstract

In this paper, we investigate the transport coefficients, specifically the
bulk ($\zeta$) and shear ($\eta$) viscosity coefficients of hot QCD matter at
finite temperature and vanishing chemical potential. For gluonic plasma, we
consider the Gribov’s prescription for the gluon propagator, which has been
proven to be effective in exploring the non-perturbative region. For a system
of quark-gluon plasma (QGP), we adopt the quasiparticle prescription for
quarks along with the Gribov’s prescription for gluons. We determine the
Gribov parameter and quasi-quark mass by utilizing lattice QCD data for the
equation-of-state of pure glue and QGP, respectively. Our primary objective is
to examine the influence of quasi-quarks on the transport coefficients of QGP.
In the end, we observed a significant decrease in the scaled transport
coefficients with temperature across the entire relevant temperature range
($1\leq T/T_{c}\leq 3$).

##  I Introduction

The Large Hadron Collider (LHC) at CERN, Geneva, and the Relativistic Heavy
Ion Collider (RHIC) at the Brookhaven National Laboratory (BNL) in New York
are two remarkable scientific facilities that serve as vital platforms for
studying the primordial form of matter that existed in the universe few
microseconds after the Big Bang [1]. High-energy nucleus-nucleus collision
experiments at these facilities have the potential to convert a portion of the
kinetic energies of the colliding nuclei into heating the QCD vacuum within an
extremely small volume. This could result in achieving temperatures millions
of times hotter than the core of the sun. At temperatures around $\sim$150
MeV, the transition of nuclear matter from the hadronic phase, consisting of
protons and neutrons, to the deconfined phase, consisting of quarks and
gluons, occurs. This transition from the hadronic phase to the nearly
equilibrated deconfined phase leads to the emergence of a new state of matter
termed as the Quark-Gluon Plasma (QGP) [2].

Analysis of spectra and azimuthal anisotropy of particles produced in
relativistic heavy-ion collisions suggested that the QGP produced is strongly
coupled [3, 4, 5, 6]. Application of viscous hydrodynamics to understand the
evolution of the fireball produced in heavy-ion collisions gained widespread
interest ever since the ratio of shear viscosity to entropy density ($\eta/s$)
was estimated to be close to the absolute lower bound for all fluids, $1/4\pi$
[7]. This allows one to investigate the behavior of QGP in the light of
viscous hydrodynamics [8, 9, 10, 11]. Apart from the equation of state,
transport coefficients such as shear viscosity and bulk viscosity are major
supplements in the hydrodynamic equation, which governs the evolution of the
fireball. While the equation of state for QGP is well understood from lattice
QCD calculations, the transport coefficients are still not very well
constrained. Therefore, it is important to accurately and consistently
quantify the transport coefficients of QGP in order to understand the
evolution of strongly interacting QCD matter formed in high-energy nucleus-
nucleus collisions.

Considering these, there has been a notable endeavor in the past to determine
the ratio $\eta/s$ as well as the ratio of bulk viscosity to entropy density
($\zeta/s$) of the QGP as a function of the temperature ($T$). The transport
coefficients of the hot QCD matter have been an active area of research for
more than two decades and have been dealt with various approaches which
include the kinetic theory approach, quasiparticle model, Chapman-Enskog
method, Gribov plasma, Polyakov-quark-meson model [12, 13, 14, 15, 16, 17, 18,
19]. The transport properties of the QGP have also been extensively studied in
the high-temperature regime using resummed perturbation theory [20] as well as
lattice calculations [21, 22, 23]. In the phenomenologically relevant
temperature range, typically $T\sim(1-4)T_{c}$, where $T_{c}$ represents the
critical temperature, it is evident that the perturbative calculation alone is
not an efficient tool due to the limitations that arise when the coupling
strength of the system becomes of the order of unity [24]. Due to this, a
theoretical model that encapsulates the non-perturbative features is more
favored for exploring the medium QGP.

Lattice QCD is a numerical approach that is best suited for reliable
computations in the nonperturbative regime. However, lattice QCD is
computationally expensive, and it is challenging to calculate dynamical
quantities. Moreover, past studies based on lattice calculation have revealed
that the results for transport coefficients are plagued by large
inconsistencies, particularly in the phenomenologically relevant temperature
regime [22, 23] and are still inconclusive. Therefore, it is useful to have an
alternative approach to incorporate nonperturbative features in the theory. In
this work, we adopt the prescription suggested by Gribov in his seminal work
Ref. [25]. Following this idea, Zwnaziger, in the later year, formulated the
local and renormalizable action discussed in detail in Ref. [26].

The Gribov’s prescription has been explored quite extensively encompassing a
broad spectrum of topics [27, 28, 29, 30, 31, 32, 33, 34, 35]. In [36], the
connection between the Gribov parameter and the topological structure of the
QCD vacuum has been explored. The results of this study shed light on a
profound relationship between the Gribov-Zwanziger (GZ) quantization and the
confinement-deconfinement phase transition. Earlier works, such as [37, 38,
39, 40], have yielded promising agreement with lattice results, particularly
in the infrared regime.

These compelling consistencies with lattice findings have motivated us to
consider a Gribov-improved gluon propagator in investigating transport
coefficients of the hot QCD matter. To consider the effect of quarks on the
transport coefficients, we adopt a quasiparticle model [41]. A similar
approach has been explored previously in [42], where the influence of
dynamical quarks on QGP transport coefficients was explored by incorporating
the interactions within the effective mass through a temperature-dependent
running coupling derived from lattice QCD thermodynamics.

The structure of this article is as follows: In Sec. II, we outline the Gribov
prescription and the motivation behind it. In Sec. III, we delve into the
thermodynamics of gluonic plasma using the Gribov prescription within the
kinetic theory. Additionally, we employ a quasiparticle approach for quarks
within the framework of kinetic theory to analyze the thermodynamic properties
of the Quark-Gluon Plasma (QGP). We ensure the thermodynamic consistency is
maintained within the desired temperature range by fixing the Gribov parameter
and quasiparticle mass appropriately. Following this, we discussed the
velocity of sound($c_{s}^{2}$) in Sec. IV, which plays a major role in the
determination of bulk viscosity. In the next section V, we utilize the
theoretical framework established in [43] to derive the coefficient of shear
and bulk viscosity for the Gribov gluons and subsequently for QGP. We found
that the inclusion of quasi-quarks diminishes the transport coefficients for
the entire range of temperature($1\leq T/T_{c}\leq 3$). In Sec. VI, we explore
the behavior of the ratio $\zeta/\eta$ for Gribov-modified gluonic plasma as
well as for the QGP. In Sec. VII, we summarize our findings, provide
concluding remarks, and establish connections between our findings and related
studies conducted using various approaches.

Throughout this paper, we adopt natural units with $c=\hbar=k_{B}=1$. We
denote the magnitude of three vectors using bold and italic font (p), while
four vectors are represented in a standard font. The center-dot notation
signifies the scalar product of four vectors employing the metric
$g^{\mu\nu}={\rm diag}(+1,-1,-1,-1)$. The four-vector
$p=(p_{0},\vec{\textbf{{p}}})$, and the magnitude of $\vec{\textbf{{p}}}$ is
denoted by p.

##  II Gribov Gluons

In order to counter the divergence in the functional integral in the
quantization of the non-abelian gauge theory, gauge fixing seemed inevitable.
The gauge fixing method in the quantization of the non-abelian gauge field was
first proposed by Faddeev and Popov back in 1967, where the gauge fixing was
ensured by incorporating the delta function into the path integral which,
without the source term reads [44]

$\displaystyle\mathcal{Z}(0)=\int\mathcal{D}\alpha\int\mathcal{D}A\hskip
4.26773pt\delta(G(^{\alpha}\\!A))e^{-S_{YM}}\Delta_{FP}\leavevmode\nobreak\ ,$
(1)

where $\Delta_{FP}$ is the Faddeev-Popov determinant given as

$\displaystyle\Delta_{FP}=\bigg{|}\frac{1}{g}{\rm
det}\mathcal{M}_{ab}\bigg{|}=\bigg{|}{\rm
det}\frac{\delta\mathcal{G}^{a}(^{\alpha}\\!A))}{\delta\Lambda^{b}}\bigg{|}_{G(^{\Lambda}\\!A)=0}\leavevmode\nobreak\
,$ (2)

where $\mathcal{G}(^{\Lambda}\\!A)$ is the gauge constraint. $\Lambda^{b}$ is
the infinitesimal parameter of a particular gauge transformation
$U=1-\Lambda^{b}T^{b}+\mathcal{O}(\Lambda^{2})$. The Faddeev-Popov operator is
given as

$\displaystyle\mathcal{M}^{ab}=-\partial_{\mu}D_{\mu}^{ab}\leavevmode\nobreak\
.$ (3)

In Eq. (3), $D_{\mu}^{ab}$ is the covariant derivative. While the FP
quantization has shown considerable success across various aspects, it remains
incomplete. The root of the issue lies in our insistence on an ideal gauge
fixing condition by stipulating that the gauge orbit must intersect the gauge
fixing constraint surface precisely once and only once. However, in practical
scenarios, this is not true and is known as Gribov ambiguity. Gribov’s
approach, extensively reviewed in Refs [45] and [46], focuses on improving the
infrared dynamics of non-abelian gauge theory by fixing the residual gauge
transformations that persist even after applying Faddeev-Popov quantization.
Mathematically, it was achieved by restricting the path integral to a Gribov
region $\Omega$ in a way so chosen that for any gauge field configuration
$A^{\mu}_{a}$,

$\displaystyle\Omega=\\{A_{a}^{\mu},\partial_{\mu}A^{\mu}_{a}=0;\mathcal{M}_{ab}>0\\}\leavevmode\nobreak\
,$ (4)

where the gauge field $A^{\mu}_{a}$ obeys the Landau gauge condition for which
the Faddeev-Popov
operator$(\mathcal{M}_{ab}=-\partial_{\mu}\mathcal{D}^{\mu}_{ab})$ is strictly
greater than 0 to ensure the no-pole condition required the overcome the
Gribov ambiguity. Following this restriction on the gauge field configuration,
the partition function in the Euclidean space-time reads as [45]

$\displaystyle\mathcal{Z}$
$\displaystyle=\int[\mathcal{D}c]\int[\mathcal{D}\bar{c}]\int[\mathcal{D}A]\hskip
2.84544ptV(\Omega)\delta(\partial\cdot A)$ $\displaystyle\times{\rm exp}\hskip
2.84544pt[-S_{YM}-\int
d^{4}x\bar{c}^{a}(x)\partial_{\mu}D_{\mu}^{ab}c^{b}(x)]\leavevmode\nobreak\ ,$
(5)

where $c(x)$ and $\bar{c}(x)$ are the ghost and anti-ghost fields
respectively. The term $V(\Omega)=\theta[1-\sigma_{0}]$ is the step function
needed to take care of the no-pole condition. Also,
$[1-\sigma_{0}]=[1-\sigma(P=0)]$ is the inverse of ghost dressing function,
defined as $Z_{G}=[1-\sigma(P)]^{-1}$. In order to obtain the gluon two-point
function, the constraint on the gauge fields is inserted in the Eq. (5) by re-
expressing the parameter $V(\Omega)$ as

$\displaystyle V(\Omega)=\frac{1}{2\pi
i}\int_{-i\infty+0^{+}}^{i\infty+0^{+}}\frac{ds}{s}\hskip 2.84544pt{\rm
exp}[s(1-\sigma_{0})]\leavevmode\nobreak\ .$ (6)

The Gribov modified gluon propagator in Landau gauge, with the dynamical mass
term, $m_{g}$, is given as [47]

$\displaystyle\mathcal{D}^{ab}_{\mu\nu}(p)=\delta^{ab}\frac{p^{2}}{p^{4}+m_{g}^{2}p^{2}+\gamma_{G}^{4}}\left[\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{p^{2}}\right]\leavevmode\nobreak\
,$ (7)

where $p=(p_{0},\vec{\textbf{{p}}})$, and the thermal mass $m_{g}$ has the
form $m_{g}=gT$ and $\gamma_{G}$ is the Gribov parameter which can be fitted
from known lattice data of thermodynamics quantities or is solved self-
consistently from a gap equation that is defined to infinite loop orders. In
comparison to the normal gluon propagator, the Gribov-modified gluon
propagator is infrared suppressed and has complex poles as $p^{2}=\pm
i\gamma_{G}^{2}$. This configuration (with imaginary poles) prevents us from
ascribing the conventional particle interpretation to the gluons, undermining
the notion of gluons as excitations of the physical spectrum. Consequently,
gluons experience confinement within the Gribov approach.

At one-loop, the analytical solution at asymptotically high temperature gives
us [48]

$\displaystyle\gamma_{G}=\left(\frac{d-1}{d}\right)\frac{\rm
N_{c}}{4\sqrt{2}\pi}g^{2}T,$ (8)

where $d$ is the space-time dimension and $\rm N_{c}$ is the number of colors,
and $g$, the running coupling with $T$ being temperature. It is evident from
Eq. (8) that the Gribov parameter is proportional to the magnetic scale, i.e.
$\gamma_{G}\propto g^{2}T$. It is worth noting that the inclusion of the
magnetic scale inherently suggests a promising resolution to IR catastrophe
caused by the magnetic scale, called the Linde problem (see [49] for detail).
Therefore, the dynamical gluonic mass of the electric scale ($gT$) along with
the Gribov parameter of magnetic scale ($g^{2}T$) could be utilized to explore
the dynamics of the deconfined medium. A similar approach has been previously
studied in the Refs [27] where dilepton rate and quark number susceptibility
have been worked out and also in Ref [29] to examine the quark self-energy at
one-loop level at high temperature.

##  III Formalism

###  III.1 Thermodynamics

In our investigation of the system’s thermodynamics, we have employed the
well-established kinetic theory [50] approach. We have dedicated separate
sections to develop the thermodynamics of quasi-quarks and Gribov-modified
gluons, adopting distinct methodologies for each.

To characterize the system’s dynamics, the energy-momentum tensor takes on the
following form [51]

$\displaystyle T_{(0)}^{\mu\nu}=\int dp\cdot
p^{\mu}p^{\nu}f_{0}+B_{0}(T)g^{\mu\nu}\leavevmode\nobreak\ ,$ (9)

where $B_{0}(T)$ is bag pressure, which is added to take care of thermodynamic
consistency in equilibrium, $g^{\mu\nu}$ is the metric tensor where
$g^{\mu\nu}={\rm diag}(+1,-1,-1,-1)$ and $f_{0}$ is the equilibrium
distribution function.

The pressure and energy density are expressed as[12]

$\mathcal{P}=-\frac{1}{3}\Delta_{\mu\nu}T_{(0)}^{\mu\nu}=\mathcal{P}^{Kinetic}-B_{0}(T)\leavevmode\nobreak\
.$ (10)
$\mathcal{E}=u_{\mu}u_{\nu}T_{(0)}^{\mu\nu}=\mathcal{E}^{kinetic}+B_{0}(T)\leavevmode\nobreak\
.$ (11)

where $\Delta_{\mu\nu}=g_{\mu\nu}-u_{\mu}u_{\nu}$, $u_{\mu}$ is the fluid four
velocities satisfying $u^{\mu}u_{\mu}=1$. In fluid rest frame,
$u^{\mu}=(1,0,0,0)=(1,\vec{0})$. $\mathcal{P}^{Kinetic}$ and
$\mathcal{E}^{kinetic}$ are the particle contribution to pressure and energy
density, respectively.

####  III.1.1 Gribov gluons

In order to formulate the thermodynamics of the gluonic plasma, we have
considered Gribov’s idea discussed briefly in section II where the propagator
is given by Eq. (7). Also, gluons are spin-one particles that obey the Bose-
Einstein statistics, where the distribution function is

$\displaystyle f_{g}(\epsilon)=\left[{\rm
exp}\left(\frac{\epsilon}{T}\right)-1\right]^{-1}\leavevmode\nobreak\ ,$ (12)

where $\epsilon$ is the energy. The Lorentz invariant momentum integral is
given by

$\displaystyle\int dp=\frac{g_{gl}}{(2\pi)^{3}}\int d^{3}\textbf{{p}}\int
dp_{0}\hskip 2.84544pt2\times\Theta(p_{0})\hskip
2.84544pt\delta(\mathcal{K})\leavevmode\nobreak\ ,$ (13)

where $\mathcal{K}=(p^{2}-m_{g}^{2}+\gamma_{G}^{4}/p^{2})$ with
$m_{g}=m_{g}(T)$ is the thermal mass and $\gamma_{G}=\gamma_{G}(T)$ is the
Gribov parameter defined in Eq. (8). Also, $g_{gl}$ is the degeneracy factor
for gluonic plasma, which is $g_{gl}=2\times(\rm N_{c}^{2}-1)=16$ and
$\Theta(p_{0})$ is the Heaviside step function to ensure that only the energy
states with positive real part are considered.

Using Eq. (10), Eq. (11) and Eq. (13), for gluonic plasma, the pressure and
energy density are obtained as

$\begin{split}\mathcal{P}_{Gl}&=-\frac{1}{3}\Delta_{\mu\nu}T_{(0)}^{\mu\nu}=\mathcal{P}_{Gl}^{Kinetic}-B^{Gl}(T)\\\
&=\frac{2g_{gl}}{3\times(2\pi)^{3}}\int d^{3}\textbf{{p}}\hskip
2.84544pt\textbf{{p}}^{2}\int dp_{0}\hskip 2.84544pt\Theta(p_{0})\hskip
2.84544pt\delta(\mathcal{K})f_{g}-B^{Gl}(T)\leavevmode\nobreak\ .\end{split}$
(14)

and

$\begin{split}\mathcal{E}_{Gl}&=u_{\mu}u_{\nu}T_{(0)}^{\mu\nu}=\mathcal{E}_{Gl}^{kinetic}+B^{Gl}(T)\\\
&=\frac{2g_{gl}}{(2\pi)^{3}}\int d^{3}\textbf{{p}}\int dp_{0}\hskip
2.84544pt\Theta(p_{0})\hskip
2.84544pt\delta(\mathcal{K})p_{0}^{2}f_{g}+B^{Gl}(T)\leavevmode\nobreak\
.\end{split}$ (15)

Note that the bag pressure $B^{Gl}(T)$ is added to the energy density and
subtracted from the pressure to maintain the thermodynamics consistency.

The following identity of the delta function has been implemented to deal with
the delta function inside the integral

$\displaystyle\delta(f(x))=\sum_{i=1}^{n}\frac{\delta(x-x_{i})}{\mid{f^{\prime}(x_{i})}\mid},$
(16)

where $x_{i}$ are the roots of the function $f(x)$.

Using Eq. (9), Eq. (10), Eq. (11) and Eq. (16), the expressions for the
gluonic pressure and energy density is given by

$\mathcal{P}_{Gl}=\frac{g_{gl}}{3\pi^{2}}\int d\textbf{{p}}\hskip
1.42271pt\textbf{{p}}^{4}\hskip 1.42271pt\left[\frac{f_{g}(\epsilon_{1})}{\mid
f^{\prime}(\epsilon_{1})\mid}+\frac{f_{g}(\epsilon_{2})}{\mid
f^{\prime}(\epsilon_{2})\mid}\right]-B^{Gl}(T)\leavevmode\nobreak\ ,$ (17)

$\mathcal{E}_{Gl}=\frac{g_{gl}}{\pi^{2}}\int d\textbf{{p}}\hskip
1.42271pt\textbf{{p}}^{2}\left[\frac{\epsilon_{1}^{2}f_{g}(\epsilon_{1})}{\mid
f^{\prime}(\epsilon_{1})\mid}+\frac{\epsilon_{2}^{2}f_{g}(\epsilon_{2})}{\mid
f^{\prime}(\epsilon_{2})\mid}\right]+B^{Gl}(T)\leavevmode\nobreak\ ,$ (18)

respectively, where fluid velocity $u^{\mu}=(1,0,0,0)=(1,\vec{0})$. Notice
that in Eq. (17) and Eq. (18)

$\displaystyle\epsilon_{1}=\frac{1}{\sqrt{2}}\sqrt{m_{g}^{2}+2\textbf{{p}}^{2}-\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\leavevmode\nobreak\
,$ (19)

and

$\displaystyle\epsilon_{2}=\frac{1}{\sqrt{2}}\sqrt{m_{g}^{2}+2\textbf{{p}}^{2}+\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\leavevmode\nobreak\
.$ (20)

which is obtained by solving the dispersion relation, Eq. (7). Also, following
the identity given in Eq. (16), we get

$\displaystyle f^{\prime}(\epsilon_{1})=\frac{4\hskip
1.42271pt\epsilon_{1}\left(-m_{g}^{4}+4\gamma_{G}^{4}+m_{g}^{2}\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}\
\right)}{\left(m_{g}^{2}-\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}\
\right)^{2}}\leavevmode\nobreak\ ,$ (21)

and

$\displaystyle f^{\prime}(\epsilon_{2})=-\frac{4\hskip
1.42271pt\epsilon_{2}(m_{g}^{4}-4\gamma_{G}^{4}+m_{g}^{2}\sqrt{m_{g}^{4}-4\gamma_{G}^{4}})}{(m_{g}^{2}+\sqrt{m_{g}^{4}-4\gamma_{G}^{4}})^{2}}\leavevmode\nobreak\
.$ (22)

Being equipped with pressure, the energy density, and hence the entropy
density relations, we now fix the running coupling $g$, using lattice data.

In order to fit the scaled trace anomaly($\mathcal{I}/T^{4}$) with the lattice
findings, we have made use of the following fit function, also explored in
Ref. [52]

$\displaystyle\frac{\mathcal{I}}{T^{4}}=$ $\displaystyle{\rm
Exp}\left[-\left(\frac{h_{1}}{\hat{T}}+\frac{h_{2}}{\hat{T}^{2}}\right)\right]\times$
$\displaystyle\left(\frac{h_{0}}{1+h_{3}\hat{T}^{2}}+\frac{f_{0}({\rm
Tanh}[f_{1}\hat{T}+f_{2}]+1)}{1+g_{1}\hat{T}+g_{2}\hat{T}^{2}}\right),$ (23)

where $\hat{T}$ is the scaled temperature ($=T/T_{c}$).

$h_{0}$ | $h_{1}$ | $h_{2}$ | $h_{3}$ | $f_{0}$ | $f_{1}$ | $f_{2}$ | $g_{1}$ | $g_{2}$
---|---|---|---|---|---|---|---|---
$0.23$ | $-1.83$ | $2.92$ | $0.07$ | $0.32$ | $62.39$ | $-62.55$ | $-1.98$ | $1.08$
Table 1: The values of the fit parameters in the Eq. (23) as a result of
fitting with the lattice data of the scaled trace anomaly($\mathcal{I}/T^{4}$)
as a function of scaled temperature($T/T_{c}$) available in Ref [53].

For this set of parameters, the scaled pressure can be obtained by

$\displaystyle\frac{\mathcal{P}(T)}{T^{4}}=\frac{\mathcal{P}(T_{0})}{T^{4}}+\int_{T_{0}}^{T}\frac{d\hat{T}}{\hat{T}}\times\frac{\mathcal{I}}{\hat{T}^{4}}\leavevmode\nobreak\
,$ (24)

where ($\mathcal{P}(T_{0}=0.7T_{c})/T^{4}=0.0015$). The energy density is
another important parameter that needs to be known in order to have a complete
idea of the thermodynamic behavior of the medium under consideration.
Following the Eq. (23) and Eq. (24), $\mathcal{I}=\epsilon-3\mathcal{P}$
provides the variation of energy density $\epsilon$ as a function of scaled
temperature. In order to have the entropy density, we made use of
thermodynamic relation
$s_{0}=d\mathcal{P}_{0}/dT=(\mathcal{P}_{0}+\epsilon_{0})/T$. Finally, the
entropy density of the Gribov-modified gluon is been equated with that of the
entropy density from the lattice in order to fit the running coupling, $g$,
which is a common parameter in the dynamical mass $m_{g}(=gT)$ as well as in
the Gribov parameter $\gamma_{G}$, defined in Eq. (8).

Figure 1: Plot of scaled fitted entropy density$(s/T^{3})$ and running
coupling $g$ for Gribov-modified gluon as a function of scaled
temperature$(T/T_{c})$ in the temperature range ($0.7\leq T/T_{c}\leq 3.5$).
The Lattice data(Green dots) are taken from [53]. The red-dotted plot
corresponds to the entropy density, which is obtained using the thermodynamic
relation $s_{0}=d\mathcal{P}_{0}/dT=(\mathcal{P}_{0}+\epsilon_{0})/T$. Figure
2: The dynamical gluon mass, $m_{g}$, along with the Gribov mass parameter,
$\gamma_{G}$, both scaled with temperature $T$, plotted against scaled
temperature in the range($1\leq T/T_{c}\leq 3.5)$.

Taking the derivative of Eq. (17) and using the thermodynamic relation
$s=\frac{d\mathcal{P}_{Gl}}{dT}=\frac{\mathcal{P}_{Gl}+\mathcal{E}_{Gl}}{T}=\frac{\mathcal{P}_{Gl}^{Kinetic}+\mathcal{E}_{Gl}^{Kinetic}}{T}$(see
Eq. (10) and Eq. 11) in the above equation, bag pressure for gluonic plasma
takes the form

$\displaystyle\frac{dB^{Gl}(T)}{dT}=\frac{d\mathcal{P}_{Gl}^{kinetic}}{dT}-\frac{\mathcal{P}_{Gl}^{Kinetic}+\mathcal{E}_{Gl}^{Kinetic}}{T}$
(25)

Note that $``Gl"$ is specifically used here for gluons.

####  III.1.2 Quasi quarks

In the study of quark-gluon plasma, the quasi-particle concept describes the
collective behavior of quarks and gluons within the plasma. The quasi-
particles are effective descriptions and provide a useful way to understand
and model the complex dynamics of strongly interacting systems under extreme
conditions. In the quasi-particle approach, the energy of the particle is
taken as $\epsilon_{q}=(\textbf{{p}}^{2}+m_{q}^{2})^{1/2}$, where $m_{q}\equiv
m_{q}(T)$ denotes the temperature-dependent mass of the quasi-particles, here,
quarks. The distribution function of the quark, following the Fermi Dirac
statistics, is given as

$\displaystyle f_{q}(\epsilon_{q})=\left[{\rm
exp}\left(\frac{\epsilon_{q}}{T}\right)+1\right]^{-1}\leavevmode\nobreak\ .$
(26)

In order to formulate the thermodynamics of quasi-quarks, the pressure and
energy density are evaluated using Eq. (10) and Eq. (11). In this case, the
Lorentz invariant momentum integral takes the form

$\displaystyle\int dp=\frac{g_{q}}{(2\pi)^{3}}\int d^{3}\textbf{{p}}\int
dp_{0}\hskip 2.84544pt2\ \Theta(p_{0})\hskip
2.84544pt\delta(p^{2}-m_{q}^{2})\leavevmode\nobreak\ .$ (27)

Using the similar methods discussed in the previous subsection, we obtain
respectively the pressure ($\mathcal{P}_{Q}$) and energy density
($\mathcal{E}_{Q}$) for the quasi-quarks as

$\displaystyle\mathcal{P}_{Q}=\frac{g_{q}}{6\pi^{2}}\int d\textbf{{p}}\hskip
1.42271pt\textbf{{p}}^{4}\frac{f_{q}(\epsilon_{q})}{\epsilon_{q}}-B^{Q}(T)\leavevmode\nobreak\
,$ (28) $\displaystyle\mathcal{E}_{Q}=\frac{g_{q}}{2\pi^{2}}\int
d\textbf{{p}}\hskip
1.42271pt\textbf{{p}}^{2}f_{q}(\epsilon_{q})\epsilon_{q}+B^{Q}(T)\leavevmode\nobreak\
.$ (29)

The degeneracy factor, in this case, is $g_{q}=2\times 2\times{\rm N_{c}\times
N_{f}}$, where the factors correspond to the contributions from anti-quarks,
spin, number of colors, and flavors, respectively. Also, $B^{Q}(T)$ is the
bag-like term in quasi-quarks pressure and energy density. In order to ensure
thermodynamic consistency in the calculation, it is required that
$s=d\mathcal{P}_{Q}/dT=(\mathcal{P}_{Q}+\mathcal{E}_{Q})/T$. Therefore, taking
the derivative of the pressure from Eq. (28), it is obtained that the
thermodynamic consistency is maintained when

$\displaystyle\frac{dB^{Q}(T)}{dT}+\frac{g_{q}}{2\pi^{2}}m_{q}m_{q}^{\prime}\displaystyle\int
d\textbf{{p}}\hskip
3.41418pt\textbf{{p}}^{2}\frac{f_{q}(\epsilon_{q})}{\epsilon_{q}}=0\leavevmode\nobreak\
,$ (30)

where $m_{q}^{\prime}=dm_{q}/dT$.

It should be noted here that in the present framework, the quark-gluon
interaction is captured in the quasi-mass $m_{q}$. In order to obtain the
temperature dependence of the quasi mass $m_{q}$, we equate the entropy
density of the lattice finding to that of the calculated entropy density by
summing up the pressure and energy density using the Eq. (17), Eq. (18), Eq.
(28) and Eq. (29), i.e.

$\displaystyle\sum_{g,q,\bar{q}}s=s_{QGP}^{Lattice}\leavevmode\nobreak\ ,$
(31)

where $g$, $q$ and $\bar{q}$ refers to the gluons, quarks and anti-quarks
respectively and
$s=d\mathcal{P}/dT=(\mathcal{P}^{kinetic}+\mathcal{E}^{kinetic})/T$ is the
entropy density in the quasi-particle approach. The lattice data considering
the quark effect is taken from Ref. [52].

In QGP case, scaled trace anomaly ($\mathcal{I}/T^{4}$) were fitted using the
following fit function

$\displaystyle\frac{\mathcal{I}}{T^{4}}={\rm
Exp}\left[-\left(\frac{h_{1}}{\hat{T}}+\frac{h_{2}}{\hat{T}^{2}}\right)\right]\left(h_{0}+\frac{f_{3}({\rm
Tanh}[f_{1}\hat{T}+f_{2}]+1)}{1+g_{1}\hat{T}+g_{2}\hat{T}^{2}}\right)\leavevmode\nobreak\
.$ (32)

where $\hat{T}=T/200({\rm MeV})$. Those parameters, for (2+1) flavour QCD take
the values [52] $h_{0}=0.139$, $h_{1}=-0.18$, $h_{2}=0.035$, $f_{3}=2.76$,
$f_{1}=6.79$, $f_{2}=-5.29$, $g_{1}=-0.47$, $g_{2}=1.04$. For this set of
parameters, again, the scaled pressure can be obtained by the equation

$\displaystyle\frac{\mathcal{P}_{QGP}(T)}{T^{4}}=\frac{\mathcal{P}_{QGP}(T_{0})}{T^{4}}+\int_{T_{0}}^{T}\frac{d\hat{T}}{\hat{T}}\frac{\mathcal{I}}{\hat{T}^{4}}\leavevmode\nobreak\
.$ (33)

Equipped with the pressure and trace anomaly fit functions corresponding to
lattice data available in Ref. [52], the sum of entropy densities of quark
(with quasi mass $m_{q}$) and gluon has been equated with that of the entropy
density of the QGP, like Eq. (31). As a result, we got the temperature
dependence of scaled quasi mass $m_{q}(T)/T$, shown in Fig (3).

Figure 3: Quasi mass, $m_{q}$, scaled with temperature $T$, plotted as a
function of scaled temperature in the range $0.8\leq T/T_{c}\leq 3.1$.

Since $m_{q}$ has been obtained in a way that also encodes the information
about quark-gluon interaction, therefore it could be said that

$\displaystyle\mathcal{P}_{Q}+\mathcal{P}_{Gl}=\mathcal{P}_{QGP}\leavevmode\nobreak\
,$ (34)

where $\mathcal{P}_{QGP}$ is the pressure of medium QGP. Now, since
$\mathcal{P}_{Q}=\mathcal{P}_{Q}^{kinetic}-B^{Q}(T)$ and
$\mathcal{P}_{Gl}=\mathcal{P}_{Gl}^{kinetic}-B^{Gl}(T)$, therefore we get

$\displaystyle
B^{Q}(T)=\mathcal{P}_{Q}^{Kinetic}+\mathcal{P}_{Gl}^{Kinetic}-\mathcal{P}_{QGP}-B^{Gl}(T)\leavevmode\nobreak\
.$ (35)

In our investigation of transport coefficients incorporating quasi-quark, we
have made the significant assumption of treating light and strange quarks
equally, despite their distinct physical characteristics, and neglected any
possible dependence of transport coefficients on quark flavors.

##  IV Speed of sound

The speed of sound, $c_{s}^{2}$, plays a pivotal role in characterizing the
bulk viscosity and quantifying deviations from the conformal limit.
Specifically, the square of the speed of sound can be determined by evaluating
the entropy density as a function of temperature. This fundamental
relationship is essential in understanding the thermodynamic properties of the
system under consideration, especially in scenarios where the conformal
symmetry is perturbed. Consequently, this insight contributes significantly to
comprehending fluid dynamics and thermodynamic behavior in various physical
systems, making it paramount in heavy ion collision. Mathematically, the speed
of sound, $c_{s}^{2}$ is expressed as

$\displaystyle
c_{s}^{2}=\left(\frac{\partial\mathcal{P}}{\partial\epsilon}\right)=\frac{s}{T}\left(\frac{\partial
s}{\partial T}\right)^{-1}\leavevmode\nobreak\ ,$ (36)

where $\mathcal{P}$ is the pressure, $\epsilon$ is energy density and $s$ is
entropy density while $T$ being temperature.

In this study, we investigated two scenarios, one involving Gribov-modified
gluons and the other involving quasi-quarks along with Gribov-modified gluons,
with the primary objective to examine the impact of quasi-quarks on the
transport coefficients of the Quark-Gluon Plasma (QGP). As a result of our
analysis, we identified two distinct values for the velocity of sound in the
system: one corresponding to a purely gluonic medium($N_{f}=0$) and the other
representing the velocity of sound in the QGP medium($N_{f}=3$).

For medium composed of Gribov modified gluons, using Eq. (36), Eq. (17) and
Eq.(18), we get the velocity of sound, $c_{s}^{2}$, whose temperature
dependence is depicted in Fig(4).

Figure 4: Temperature dependence of the speed of sound, $c_{s}^{2}$ of Gribov-
modified gluonic plasma as well as QCD($N_{f}=3$) as a function of scaled
temperature ($T/T_{c}$). Lattice result for gluonic plasma has been obtained
from [53], while for $N_{f}=3$ has been taken from [52]. The horizontal solid
grey line corresponds to the speed of sound $(c_{s}^{2}=1/3)$, which is
attained in the conformal limit of the medium QGP.

Through our analysis, it is evident that the squared speed of sound in the
Gribov-modified gluons exhibits remarkable consistency, not only within the
confined phase but also in the deconfined phase when compared with the one
obtained from lattice data concerning the pressure and energy density in pure
Yang-Mills theory. Also, we witness a pronounced escalation in the squared
speed of sound ($c_{s}^{2}$) as the temperature reaches the transition point,
i.e., $T/T_{c}=1$. This notable increase can be attributed to a sudden
discontinuity in the entropy density, indicative of a first-order phase
transition in gluonic plasma.

Secondly, while adding the quasi-quark contribution, we have the speed of
sound ($c_{s}^{2}$) of the medium QGP, composed of Gribov-modified gluons and
quasi-quarks. In this case, in order to evaluate the speed of sound, we took
the pressure ($P$) and energy density ($\epsilon$) as

$\displaystyle\mathcal{P}$ $\displaystyle=$
$\displaystyle\mathcal{P}_{Gl}+\mathcal{P}_{Q}\leavevmode\nobreak\ ,$ (37)
$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\epsilon_{Gl}+\epsilon_{Q}\leavevmode\nobreak\ .$ (38)

Here we have used Eq. (17), Eq. (28), Eq. (18) and Eq. (29). The speed of
sound ($c_{s}^{2}$) of medium QGP is depicted in fig (4) again, along with the
lattice finding of the same.

It is clear from Fig (4) that our calculation matches the lattice in and out
of the deconfinement phase well. Compared with the speed of sound
($c_{s}^{2}$) of purely gluonic medium(in the same figure), we notice in the
present case that the variation with scaled temperature is smoother near the
transition region. Also, we see that as the temperature increases, the speed
of sound approaches the conformal limit, i.e., $c_{s}^{2}\rightarrow 1/3$ for
$N_{f}=0$ as well as for $N_{f}=3$.

Also, apart from being an important input parameter in the transport
coefficient(bulk viscosity), a decent match with the available lattice data of
the speed of sound to that of the calculated results solidifies the idea of
thermodynamics being well settled for the Gribov-modified gluons ($N_{f}=0$)
as well as for medium QGP ($N_{f}=3$). The next is to formulate the transport
coefficients for the two cases.

##  V Transport coefficients

In this section, we have presented the theoretical framework utilized for
calculating the transport coefficients, namely the Shear ($\eta$) and Bulk
($\zeta$) viscosities for the hot QGP medium. The centerpiece of our approach
lies in solving the Boltzmann equation, which serves as the fundamental
equation in kinetic theory and plays a crucial role in describing the
transport properties. As a result, the Boltzmann equation holds significant
importance in the field of heavy ion collisions. The Boltzmann equation is
mathematically expressed as

$\displaystyle\epsilon^{-1}p^{\mu}\partial_{\mu}f^{eq}=\frac{\partial
f^{eq}}{\partial t}+v\cdot\nabla f^{eq}=\textit{C}_{f}\leavevmode\nobreak\ ,$
(39)

where $f=f(\textbf{x},t,\textbf{p})$ is distribution function and
$\textit{C}_{f}$ is the collision kernel. In the relaxation time
approximation(RTA), the collision kernel $\textit{C}_{f}$ is

$\displaystyle\textit{C}_{f}=-\frac{u\cdot p}{\tau_{a}}\delta
f\leavevmode\nobreak\ ,$ (40)

where $\tau_{a}$ is the relaxation time and the deviation away from
equilibrium manifests itself via $\delta f$. Therefore, in RTA, the Boltzmann
equation takes the form

$\displaystyle\frac{\partial f^{eq}}{dt}+v\cdot\nabla
f=-\omega(\epsilon)\delta f\leavevmode\nobreak\ ,$ (41)

where $\omega(\epsilon)=1/\tau(\epsilon)$. Following the arguments of [43],
the expression for bulk ($\zeta$) and shear ($\eta$) viscosity takes the form

$\displaystyle\zeta=\frac{1}{3}\displaystyle\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{2}}{\epsilon}f^{eq}(\epsilon/T)\mathcal{A}(\epsilon)\leavevmode\nobreak\
.$ (42)
$\displaystyle\eta=\frac{2}{15}\displaystyle\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{4}}{\epsilon}f^{eq}(\epsilon/T)\mathcal{C}(\epsilon)\leavevmode\nobreak\
.$ (43)

As it is evident from the above expression of $\zeta$ (Eq. (42)) and $\eta$
(Eq .(43)), we need to find the unknown scalar parameters $\mathcal{A}$ and
$\mathcal{C}$ to find the temperature dependence of these two transport
coefficients. These parameters will be determined by solving the Boltzmann
equation using the local equilibrium distribution function in Gribov-modified
gluonic plasma and then separately for the quasi-quark approach. This is not
the end of it since we also have to see the consequence of the Landau-Lifshitz
condition on these parameters, which will define the final expression of bulk
($\zeta$) and shear ($\eta$) viscosity.

We confine our analysis to the context of zero chemical potentials. The
underlying assumption considers the system to be in a state of approximate
local equilibrium, characterized by local temperature distribution, $T(x)$,
along with a flow velocity, $u(x)$. Within the Landau-Lifshitz framework,
$u(x)$ signifies the energy transport velocity.

###  V.1 Shear and Bulk viscosity of Gribov-modified gluons

In this section, the transport coefficients of the Gribov-modified gluonic
plasma have been formulated in light of the theoretical setup discussed in
detail in Ref [43]. It is important to note that in Gribov-modified gluonic
plasma, there are two quasi-poles, given by Eq. (19) and Eq. (20). For this
set of energies, we will have two different forms of $\mathcal{A}$ and
$\mathcal{C}$ corresponding to each of them, determining the form of $\zeta$
and $\eta$, respectively.

While formulating the Bulk viscosity, the quantity of interest is
$\mathcal{A}$, which can be put back in Eq. (42) to get a complete form of
bulk viscosity, $\zeta$ for the Gribov-modified gluons. Here, we will have
$\mathcal{A}_{1}$ corresponding to the energy $\epsilon_{1}$(Eq. (19)), and in
the same way, we will have $\mathcal{A}_{2}$ corresponding to the energy
$\epsilon_{2}$(Eq. (20)). We first work for $\mathcal{A}_{1}$. Solving the
Boltzmann equation (LHS of Eq. (39)) we get

$\displaystyle\epsilon_{1}^{-1}p^{\mu}\partial_{\mu}f_{g}^{eq}$
$\displaystyle=-\frac{p^{\mu}f_{g}^{eq}\tilde{f}_{g}^{eq}}{\epsilon_{1}}$
$\displaystyle\\!\\!\\!\times\left(\frac{1}{T}p^{\nu}\partial_{\mu}u_{\nu}-\frac{u_{\nu}p^{\nu}}{T^{2}}\partial_{\mu}T+\frac{1}{T}u_{\nu}\partial_{\mu}p^{\nu}\right)\leavevmode\nobreak\
,$ (44)

where $\tilde{f}_{g}=1+f_{g}$ and $f_{g}^{eq}$ is the equilibrium distribution
function for gluons given by Eq. (12). Solving for the third term in Eq. (V.1)
in the local rest frame of the fluid $u=(1,0,0,0)$, for
$\epsilon=\epsilon_{1}$ (taken from Eq. (19)), we get

$\displaystyle\frac{p^{\mu}}{T}\partial_{\mu}\epsilon_{1}=-\frac{c_{s}^{2}}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{4\gamma_{G}^{3}\gamma_{G}^{\prime}-m_{g}^{3}m_{g}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\partial_{\alpha}u^{\alpha}\leavevmode\nobreak\
.$ (45)

where we used

$\displaystyle\partial_{\mu}T=T(u^{\alpha}\partial_{\alpha}u_{\mu}-c_{s}^{2}u_{\mu}\partial_{\alpha}u^{\alpha})\leavevmode\nobreak\
.$ (46)

Using Eq. (45) in Eq. (V.1) and doing further simplification, we get

$\displaystyle\epsilon^{-1}_{1}p^{\mu}\partial_{\mu}f_{g}^{eq}=-\frac{f_{g}^{eq}\tilde{f}_{g}^{eq}}{2T\epsilon_{1}}\left(p^{\mu}p^{\nu}\mathcal{Y}_{1}-\frac{2}{3}\partial_{\alpha}u^{\alpha}\mathcal{Y}_{2}\right)\leavevmode\nobreak\
,$ (47)

where

$\displaystyle\mathcal{Y}_{1}=D_{\mu}u_{\nu}+D_{\nu}u_{\mu}+\frac{2}{3}\Delta_{\mu\nu}\partial_{\rho}u^{\rho}\leavevmode\nobreak\
.$ (48)

$\displaystyle\mathcal{Y}_{2}=(1-3c_{s}^{2})\epsilon_{1}^{2}-p^{2}+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\leavevmode\nobreak\
.$ (49)

Now, following the argument of Ref [43], we see that in relaxation time
approximation,

$\begin{split}\mathcal{A}^{par}_{1}&=\frac{\tau(\epsilon_{1})}{3T\epsilon_{1}}\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{1}^{2}\times\\\
&+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}\leavevmode\nobreak\
.\end{split}$ (50)

In the same way, doing the calculation corresponding to $\epsilon_{2}$(Eq.
(20)), we get

$\begin{split}\mathcal{A}^{par}_{2}&=\frac{\tau(\epsilon_{2})}{3T\epsilon_{2}}\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{2}^{2}\\\
&+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{4\gamma_{G}^{3}\gamma_{G}^{\prime}-m_{g}^{3}m_{g}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}\leavevmode\nobreak\
.\end{split}$ (51)

Note that $\mathcal{A}^{par}_{1}$ and $\mathcal{A}^{par}_{2}$ are the
particular solution of $\mathcal{A}^{1}$ and $\mathcal{A}^{2}$ satisfying the
Landau-Lifshitz condition as discussed in detail in Appendix[A].

Now, considering the mean field effect along with Landau-Lifshitz condition(as
discussed in Appendix[A]), the final expression for Bulk viscosity($\zeta$)
for Gribov-modified gluons takes the form(corresponding to
$\mathcal{A}^{par}_{1}$), using Eq. (A) and Eq. (50)

$\displaystyle\zeta_{g}^{1}$
$\displaystyle=\frac{1}{9T}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\tau(\epsilon_{1})}{\epsilon_{1}^{2}}f_{g}^{eq}(\epsilon_{1}/T)\tilde{f}_{g}^{eq}(\epsilon_{1}/T)$
$\displaystyle\\!\\!\times\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{1}^{2}+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}^{2}\\!.$
(52)

Likewise, Bulk viscosity($\zeta$) expression, (corresponding to
$\mathcal{A}_{par}^{2}$), using Eq. (A) and Eq. (51) is

$\displaystyle\zeta_{g}^{2}$
$\displaystyle=\frac{1}{9T}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\tau(\epsilon_{2})}{\epsilon_{2}^{2}}f_{g}^{eq}(\epsilon_{2}/T)\tilde{f}_{g}^{eq}(\epsilon_{2}/T)$
$\displaystyle\\!\\!\times\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{2}^{2}+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{4\gamma_{G}^{3}\gamma_{G}^{\prime}-m_{g}^{3}m_{g}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}^{2}\\!.$
(53)

Hence, the final expression of the bulk viscosity of the Gribov-modified gluon
is

$\displaystyle\zeta_{g}=\zeta_{g}^{1}+\zeta_{g}^{2}$ (54)

In the above equation, $c_{s}^{2}$ is the speed of sound(for gluonic medium),
which is discussed in section (IV) in detail.

Now, we work to formulate the expression of shear viscosity ($\eta$) for the
Gribov-modified gluonic plasma. Following Eq. (43), we see that in order to
know about the shear viscosity, the parameter of interest is $\mathcal{C}$.
Here also we will have two separate forms of $\mathcal{C}$, i.e.
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ corresponding to two different sets of
energies Eq. (19) and Eq. (20). Using Eq. (47), Eq. (39) and Eq. (41), we see
that

$\displaystyle\mathcal{C}_{1}(\epsilon_{1})=\frac{\tau(\epsilon_{1})}{2T\epsilon_{1}}\hskip
7.11317pt{\rm and}\hskip
7.11317pt\mathcal{C}_{2}(\epsilon_{2})=\frac{\tau(\epsilon_{2})}{2T\epsilon_{2}}\leavevmode\nobreak\
.$ (55)

Therefore, using Eq. (43), the coefficient of shear viscosity ($\eta$) is
given by

$\displaystyle\eta_{g}^{1}=\frac{1}{15T}\displaystyle\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{4}}{\epsilon_{1}^{2}}\tau(\epsilon_{1})f_{g}^{eq}(\epsilon_{1}/T)\tilde{f}_{g}^{eq}(\epsilon_{1}/T)\leavevmode\nobreak\
.$ (56)

similarly,

$\displaystyle\eta_{g}^{2}=\frac{1}{15T}\displaystyle\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{4}}{\epsilon_{2}^{2}}\tau(\epsilon_{2})f_{g}^{eq}(\epsilon_{2}/T)\tilde{f}_{g}^{eq}(\epsilon_{2}/T)\leavevmode\nobreak\
.$ (57)

The final expression of the bulk viscosity of the Gribov gluon is

$\displaystyle\eta_{g}=\eta_{g}^{1}+\eta_{g}^{2}\leavevmode\nobreak\ .$ (58)

Figure 5: Temperature dependence of Shear ($\eta$, Eq. (58)) and Bulk($\zeta$,
Eq. (54)) viscosities, scaled with (“$s\tau_{R}T_{c}$”), where $s$ is the
entropy density of Gribov-modified gluonic plasma. $T_{c}$ for Gribov-Modified
gluonic plasma has been taken as 260 MeV.

It is interesting to see that mean field contribution has no effect on shear
viscosity $\eta$ since it does not explicitly depend on the quantities which
is modified by mean field effect, here, effective mass. The influence of
medium effects, specifically the equation of state (EoS), has a more
pronounced impact on the bulk viscosity than the shear viscosity. This
distinction arises from the fundamental nature of the bulk viscosity
($\zeta$), which is primarily relevant in processes that modify the volume of
a system, as opposed to its shape, in contrast to the shear viscosity
($\eta$), which in turn is prominent in the situation where there is a change
in shape. Therefore, the interplay between medium effects and transport
properties manifests more prominently through the bulk viscosity($\zeta$), as
it captures the response of the system to volume-changing processes.

The ratio $\zeta/\eta$ is also parameterized in terms of the speed of sound
($c_{s}^{2}$), allowing for a comprehensive investigation of the behavior of
the QCD system. Hence, a physical significance is associated with this ratio
($\zeta/\eta$). With shear and bulk viscosities in hand, we see the variation
of the ratio($\zeta/\eta$) as a function of scaled temperature($T/T_{c}$) in
the fig(6).

Figure 6: Temperature dependence of ($\zeta/\eta$) of Gribov-modified gluonic
plasma as a function of scaled temperature ($T/T_{c}$) in the range
($0.975\leq T/T_{c}\leq 3$). Lattice data has been taken from [22, 23], [54,
55] and [56, 57].

Near the transition temperature $T_{c}$, the ratio $\zeta/\eta$ obtained using
the Gribov-modified approach shows good agreement with the corresponding
quantity derived from existing lattice data from Refs [22, 23] and also from
Refs [54, 55]. Above $T/T_{c}=1.4$, the ratio aligns well with the lattice
finding based on perturbative Quantum Chromodynamics (pQCD) from Refs [56,
57]. A comparison between our findings and the results deduced in Ref [42],
which is based on the quasi-particle approach, reveals that our approach
provides a better match with the lattice results, particularly within the
range ($1\leq T/T_{c}\leq 2$). We also notice a bit of gradual restoration of
the conformal limit as the temperature goes up, in the Gribov-modified
approach.

###  V.2 Shear and Bulk viscosity of QGP

The inclusion of quasi-quarks plays a crucial role in shaping the
thermodynamic properties of the system under consideration. Likewise, the
transport coefficients are anticipated to exhibit substantial disparities
between Quantum Chromodynamics (QCD) and the pure gluon plasma. Investigating
the matter in which these distinctions manifest in the specific viscosity
coefficients would yield valuable insights. By exploring these variations, we
can gain a deeper understanding of the presence of quasi-quarks, which brings
to the system’s thermodynamics and transport properties.

In this section, we formulate the expressions of shear and bulk viscosities
for a system composed of quasi-quarks in relaxation time approximation. Since
the quasi mass $m_{q}$ is a function of temperature($T$), which in itself is a
function of space-time ($T=T(x,p)$), therefore considering the mean field
effects for the transport coefficients is required. For a system of quark with
energy $\epsilon_{q}=(\textbf{{p}}^{2}+m_{q}^{2})^{1/2}$, the distribution
function is as Eq. (26). Solving the Boltzmann equation for this system(LHS of
Eq. (39))

$\displaystyle p^{\mu}\partial_{\mu}f_{q}^{eq}$
$\displaystyle=\frac{-p^{\mu}f_{q}^{eq}\tilde{f}_{q}^{eq}}{\epsilon_{q}}$
$\displaystyle\times\Bigg{(}\frac{1}{T}p^{\nu}\partial_{\mu}u_{\nu}-\frac{u_{\nu}p^{\nu}}{T^{2}}\partial_{\mu}T+\frac{1}{T}u_{\nu}\partial_{\mu}p^{\nu}\Bigg{)}\leavevmode\nobreak\
.$ (59)

Simplifying the terms in the bracket and using [43]

$\displaystyle\frac{1}{T}\partial_{\mu}T=u^{\rho}\partial_{\rho}u_{\mu}-v_{s}^{2}u_{\mu}\partial_{\rho}u^{\rho}\leavevmode\nobreak\
,$ (60)

where $c_{s}^{2}=(d\mathcal{P}/d\epsilon)$, is the velocity of sound, we
finally get the expression as

$\displaystyle p^{\mu}\partial_{\mu}f_{q}^{eq}$
$\displaystyle=-\frac{p^{\mu}f_{q}^{eq}\tilde{f}_{q}^{eq}}{T\epsilon_{q}}\Big{[}p^{\mu}p^{\nu}(D_{\mu}u_{\nu}+v_{s}^{2}u_{\nu}u_{\mu}\partial_{\alpha}u^{\alpha})$
$\displaystyle+\frac{p^{\mu}}{\epsilon_{q}}Tm_{q}m_{q}^{\prime}(u^{\alpha}\partial_{\alpha}u_{\mu}-c_{s}^{2}u_{\mu}\partial_{\alpha}u^{\alpha})\Big{]}\leavevmode\nobreak\
,$ (61)

where $\tilde{f}^{eq}=(1-f^{eq})$. Now, since
$u^{\mu}=(p^{\mu}/\epsilon_{q})$, $(u\cdot p)=\epsilon$ and
$u^{\mu}u_{\mu}=1$, therefore we get

$\displaystyle p^{\mu}\partial_{\mu}f_{q}^{eq}$
$\displaystyle=\frac{-p^{\mu}f_{q}^{eq}\tilde{f}_{q}^{eq}}{2T\epsilon_{q}}\left[p^{\mu}p^{\nu}\left(D_{\mu}u_{\nu}+D_{\nu}u_{\mu}\right)\right.$
$\displaystyle\left.+2c_{s}^{2}\epsilon_{q}^{2}\partial_{\alpha}u^{\alpha}-2Tm_{q}m_{q}^{\prime}c_{s}^{2}\partial_{\alpha}u^{\alpha}\right]\leavevmode\nobreak\
.$ (62)

Now, further simplification to this equation like did in the Gribov-modified
gluonic section, we see that $\mathcal{A}$ finally takes the form

$\displaystyle\mathcal{A}(\epsilon_{q})=\frac{\tau_{a}(\epsilon_{q})}{3T\epsilon_{q}}[(1-3c_{s}^{2})\epsilon_{q}^{2}-m_{q}^{2}+3Tm_{q}m_{q}^{\prime}c_{s}^{2}]\leavevmode\nobreak\
.$ (63)

In the same way, $\mathcal{C}$ finally takes the form

$\displaystyle\mathcal{C}(\epsilon_{q})=\frac{\tau_{a}(\epsilon_{q})}{2T\epsilon_{q}}\leavevmode\nobreak\
.$ (64)

Using Eq.(42) and considering the mean field effect along with the Landau-
Lifshitz condition(as developed in [43] for quasi-particles), the final
expression for the bulk viscosity($\zeta$) takes the form,

$\begin{split}\zeta_{q/\bar{q}}=\frac{1}{9T}\int&\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\tau(\epsilon_{q})}{\epsilon_{q}^{2}}f_{q}^{eq}(\epsilon_{q}/T)\tilde{f}_{q}^{eq}(\epsilon_{q}/T)\\\
&\times\left[\textbf{{p}}^{2}-3c_{s}^{2}\left(\epsilon_{q}^{2}-T^{2}\frac{dm_{q}^{2}}{dT^{2}}\right)\right]^{2}\end{split}$
(65)

Similarly, substituting Eq.(64) in Eq.(43), shear viscosity($\eta$) takes the
form

$\displaystyle\eta_{q/\bar{q}}=\frac{1}{15T}\displaystyle\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{4}}{\epsilon_{q}^{2}}\tau_{a}(\epsilon_{q})f_{q}^{eq}(\epsilon_{q}/T)\tilde{f}_{q}^{eq}(\epsilon_{q}/T).\
\ $ (66)

Here, the subscript $``q"$ stands for a system composed of dynamical quarks.
Now, in order to see the temperature dependence of scaled transport
coefficients of the medium QGP, we sum up the expression of the individual
transport coefficients of Gribov-modified gluons and quasi-quark from Eq.
(54), Eq. (58), Eq. (65), Eq. (66), i.e.

$\displaystyle\zeta_{QGP}=\zeta_{g}+\zeta_{q/\bar{q}}\leavevmode\nobreak\ ,$
(67)

similarly,

$\displaystyle\eta_{QGP}=\eta_{g}+\eta_{q/\bar{q}}\leavevmode\nobreak\ .$ (68)

The variation of Eq. (67) as a function of scaled temperature($T/T_{c}$) has
been displayed in fig(7). While dealing with the transport coefficients of
medium QGP, particularly bulk viscosity($\zeta$), we make use of the velocity
of sound($c_{s}^{2}$) of the medium QGP, which is composed of quasi-quarks
along with Gribov-modified gluons, as discussed in detail in section IV).

Figure 7: Bulk viscosity($\zeta$) scaled with “$s\tau_{R}T_{c}$” plotted
against scaled temperature($T/T_{c}$). The blue curve corresponds to Eq. (54)
while the red one corresponds to Eq. (67). Figure 8: Shear viscosity($\eta$)
scaled with “$s\tau_{R}T_{c}$” plotted against scaled temperature($T/T_{c}$).
The blue curve corresponds to Eq. (58)) while the red one corresponds to Eq.
(68).

##  VI Bulk to Shear viscosity analysis

In this section, we present an analysis of the ratio bulk to shear
viscosity($\zeta/\eta$) for Gribov-modified gluonic plasma as well as the
plasma consisting of Gribov gluons along with quasi-quarks. The importance of
the ratio($\zeta/\eta$) lies in the fact that it carries information about the
breaking of conformal symmetry in the medium under consideration. The value of
the ($\zeta/\eta$) ratio can have significant consequences for the behavior of
the QGP. A low value of ($\zeta/\eta$) implies that shear viscosity dominates
over bulk viscosity, indicating that the QGP medium responds more readily to
changes in shape rather than volume. This can lead to a more fluid-like
behavior and a smaller tendency for the system to develop mechanical
instabilities. Conversely, a high value of ($\zeta/\eta$) suggests that bulk
viscosity is relatively more significant, implying that volume-changing
processes play a more prominent role, which could be seen near the transition
region. Physically, it leads to slower equilibration and an increased
likelihood of mechanical instabilities within the QGP.

Figure 9: The scaled temperature($T/T_{c}$) against the conformality
measure($\Delta c_{s}^{2}=(1/3-c_{s}^{2})$), normalized by its value at
$T_{c}$ which we considered 0.271 GeV for Gluonic plasma and 0.160 GeV for
medium QGP.

At high temperatures, the perturbative calculation predicts
$\zeta/\eta\propto\Delta c_{s}^{2}(=c_{s}^{2}-1/3)^{2}$ [51, 10]. The same
trend is also observed for an interacting photon gas [58] and in scalar field
theory[51]. Interestingly, near the transition region from hadronic to QGP
phase, calculation from gauge-gravity duality predicts the linear dependence
on $\Delta c_{s}^{2}$[59, 60].

We have compared our finding of ratio($\zeta/\eta$) to that of linear and
quadratic dependence on $\Delta c_{s}^{2}$, whose forms are given as

$\displaystyle{\rm Linear}:\frac{\zeta}{\eta}=\alpha(\Delta
c_{s}^{2})+\beta,{\rm and}$ (69) $\displaystyle{\rm
Quadratic}:\frac{\zeta}{\eta}=\gamma(\Delta
c_{s}^{2})^{2}+\delta\leavevmode\nobreak\ ,$ (70)

where $\alpha$, $\beta$, $\gamma$ and $\delta$ are the fit parameters. Similar
analysis has also been explored in[42] in pure Yang-Mills as well as medium
QGP in the flavor dependence of transport coefficients in a quasi-particle
approach. In our analysis, for $N_{f}=0$ and $N_{f}=3$, the speed of
sound($c_{s}^{2}$) has been discussed in detail in section (IV).

Figure 10: Bulk to Shear($\zeta/\eta$) viscosity ratio against the
conformality measure($\Delta c_{s}^{2}=(1/3-c_{s}^{2})$), for Gribov-modified
gluonic plasma. The green dashed line corresponds to linear dependence with
the fit parameters $\alpha=0.9757$ and $\beta=-0.1593$ while the red dashed
line is that for quadratic dependence with the fit parameters $\gamma=1.7616$
and $\delta=0.00608$.

In Fig(10), the variation of ratio($\zeta/\eta$) against $\Delta c_{s}^{2}$
normalized by its value at a transition temperature($T_{c}$) of Gribov-
modified gluonic plasma has been depicted. We notice a good agreement with the
quadratic scaling in the high-temperature domain, while for the intermediate
temperature range, linear scaling is a good match, which is in line with the
previous findings based on the behaviour of strongly interacting mediums. From
Fig(9), it can be seen that the shift from the linear scaling to quadratic
scaling(${\rm for}\hskip 2.84544ptN_{f}=0$) occurs around $T/T_{c}\sim 1.5$.

Another approach which describes the pure Yang-Mills plasma based on the
Gribov–Zwanziger quantization(based on a constant Gribov parameter,
$\gamma_{G}$), leads to the ratio $\zeta/\eta$ linearly proportional to the
quantity $\Delta c_{s}^{2}$[61]. In the high-temperature limit

$\displaystyle\zeta/\eta=\mathcal{X}_{GZ}\left(\frac{1}{3}-c_{s}^{2}\right)+....(T\gg\gamma_{G})\leavevmode\nobreak\
,$ (71)

where $\mathcal{X}_{GZ}=5/2$. However, in the present analysis, we have
considered the temperature dependence Gribov parameter $\gamma_{G}(T)$, whose
form is given by Eq. (8). For BE plasma(with effective mass, $m_{eff}$) at
high temperatures, again in Ref [61], it was found that, at high temperatures

$\displaystyle\frac{\zeta}{\eta}=\mathcal{X}_{BE}\left(\frac{1}{3}-c_{s}^{2}\right)^{3/2}+....(T\gg
m_{eff})\leavevmode\nobreak\ ,$ (72)

with $\mathcal{X}_{BE}\sim 5.81$. Also, at the low temperature expansion,
where $c_{s}^{2}\rightarrow 0$, corresponds to $\mathcal{X}_{GZ}=5$ while
$\mathcal{X}_{BE}=2$.

Figure 11: Bulk to Shear($\zeta/\eta$) viscosity of QGP(using Eq. (67) and Eq.
(68) ratio against conformality measure($\Delta c_{s}^{2}=(1/3-c_{s}^{2})$).
The fit parameters corresponding to linear dependence are $\alpha=0.1049$ and
$\beta=-0.01141$ while for quadratic dependence, $\gamma=0.47$ and
$\delta=0.0013$.

A similar analysis is done considering the quasi-quarks’ effects and Gribov-
modified gluons. Shown in Fig(11), the ratio ($\zeta/\eta$)(${\rm for}\hskip
2.84544ptN_{f}=3$). The transition from linear dependence to quadratic
dependence is quite evident and occurs near $T/T_{c}\sim 2.5$. Also, the
linear scaling shows an exact overlap with our findings till $T/T_{c}=1.5$ and
below, which implies the strong interaction among the constituents of the
medium(QGP).

##  VII Summary and Outlook

In this study, we have developed a framework to explore the transport
coefficients: shear viscosity ($\eta$) and bulk viscosity ($\zeta$) of the
medium QGP. While extensive literature exists on the investigation of these
coefficients in the perturbative domain, both for gluonic plasma [62, 28]
along with the inclusion of quark [52, 42], it is, however, crucial to
encompass the non-perturbative feature which has a sound effect near the phase
transition region(from hadronic to QGP medium) when examining the transport
coefficients of the QGP. Thus, our work aims to fill this gap by incorporating
non-perturbative effects and considering the influence of dynamical quarks and
Gribov-modified gluons to understand the transport properties near the
crossover region comprehensively.

The Gribov idea, known for its efficacy in capturing non-perturbative effects,
has been employed to analyze the thermodynamics of gluonic plasma. In the
context of quark systems, we have explored the quasi-particle approach, which
characterizes the energy as $\epsilon_{q}=(\textbf{{p}}^{2}+m_{q}^{2})^{1/2}$,
where the quasi-mass of quark, denoted as $m_{q}$, has been determined by
using lattice data of entropy density($s/T^{3}$) from [52]. It is worth noting
that we have captured the interaction of quarks and gluons in the quasi-mass
$m_{q}$ itself. The variation of the quasi-mass $m_{q}$ as a function of
scaled temperature has been depicted in Fig (3).

In Fig(6), we see the variation of the ratio of bulk to
shear($\zeta/\eta$)viscosity for the plasma consisting of Gribov-modified
gluons along with the lattice findings from the mentioned sources. We see a
decent match with the lattice findings, specifically near the transition
region ($T/T_{c}\sim 1$). However, it should also be noted that the transport
coefficients by lattice finding are plagued by large uncertainties near the
phase transition region. We find that the ratio $(\zeta/\eta)$ tends to zero
at high temperatures. This is also termed the Stefan-Boltzmann limit. In other
words, the system tends to regain its conformality.

The observations from Fig (7) indicate that in the context of strongly
interacting matter, the ratio of specific bulk viscosity to entropy density
$(\zeta/s\tau_{R}T_{c})$ exhibits a peak in the vicinity of the critical
temperature ($T_{c}$). This finding suggests the presence of a maximum value
for the ($\zeta/s\tau_{R}T_{c}$) near transition temperature, underscoring the
nontrivial behavior of the bulk viscosity in relation to temperature.
Interestingly, we observed a significant change(suppression) in the behavior
of the scaled bulk viscosity($\zeta/s\tau_{R}T_{c}$) throughout the entire
range of scaled temperature ($T/T_{c}$). The small value of the scaled bulk
viscosity across the whole temperature range serves as a mitigating factor
against the tendency of the fluid to undergo mechanical instability [63]. Near
$T\sim T_{c}$, the bulk viscosity of the QGP is relatively high, which implies
that the medium is less efficient at restoring equilibrium. However, as the
temperature increases beyond $T_{c}$, the specific bulk viscosity diminishes
rapidly, converging towards its conformal limit.

In the Fig(8), we also see in the proximity of the critical temperature
$T_{c}$, the scaled shear viscosity ($\eta/s\tau_{R}T_{c}$) demonstrates a
pronounced minimum. The minimum value of shear viscosity near the critical
temperature ($T_{c}$) in QGP is a consequence of the collective behavior of
the system. The subsequent increase in the temperature can be understood in
terms of the changing nature of particle interactions and the transition from
a strongly interacting, collective state to a weakly interacting state. The
rise in shear viscosity with increasing temperature is attributed to the
increasing dominance of particle collisions and scattering processes. These
interactions contribute to the transfer of momentum across different fluid
regions, resulting in a higher resistance to shear deformation and, therefore,
an elevated shear viscosity. It is worth noting that while incorporating the
effects of quasi-quarks along with Gribov-modified gluons, the scaled
transport coefficients($\zeta/s\tau_{R}T_{c}$ and $\eta/s\tau_{R}T_{c}$) got
suppressed for the entire range of temperature($1\leq T/T_{c}\leq 3$).

In Fig(10) and Fig(11), we have mapped the variation of the
ratio($\zeta/\eta$) to that on the linear and quadratic dependence on the
conformal parameter, $\Delta c_{s}^{2}=(1/3-c_{s}^{2})$ for $N_{f}=0$ and
$N_{f}=3$ case. The fit parameters are ($\alpha$ and $\beta$ for linear
dependence) and ($\gamma$ and $\delta$ for quadratic dependence). In our
findings, we notice that both fit well with the quadratic dependence at higher
temperatures, while near the transition region, we have a good fit with linear
scaling. The similar calculation in [42] shows a good fit with linear scaling
for $T/T_{c}\leq 1.3$ for $N_{f}=0$ and beyond $T/T_{c}\sim 1.4$, the
quadratic dependence is observed. A similar nature is observed for the
$N_{f}=3$ case; however, the shift from linear to quadratic dependence is now
observed near $T/T_{c}\sim 2$. In our case, for $N_{f}=0$, the shift from
linear to quadratic scaling occurs near $T/T_{c}\sim 1.5$, and for $N_{f}=3$,
the same shift occurs near $T/T_{c}\sim 2.5$.

The current approach in exploring the transport coefficients of QGP does not
aim to serve as a precise phenomenological model for a thorough investigation
of the transport coefficients ($\eta$ and $\zeta$). Instead, it could be
considered as a starting point of yet another means to incorporate non-
perturbative effects by Gribov prescription while studying the transport
coefficients of the hot QGP medium. Incorporating non-perturbative effects is
crucial, particularly in the vicinity of the transition region where the
deconfinement occurs. Such an investigation holds immense significance as it
deepens our understanding of the complex connection between non-perturbative
physics, transport coefficients, and the phenomenon of deconfinement. By
considering these non-perturbative effects, our approach contributes to a more
comprehensive analysis of the transport properties of the QGP medium,
especially in the critical region where significant changes in the system
occur.

In future research, the relaxation time ($\tau_{R}$) could be further
investigated using either momentum-dependent or momentum-independent
approaches. Additionally, to enhance the reliability of this framework, it
would be valuable to consider alternative means of determining the Gribov
parameter $\gamma_{G}$, which is not just the analytic form of it in the
asymptotic limit as presented in Eq. (8). One possibility is by fitting it
within a desired temperature range using available lattice data of
thermodynamic parameters, running coupling, or solving the one-loop gap
equation as demonstrated in [38]. Apart from these, the flavor dependence of
the system to recover the conformal limit could also be explored exclusively.
By incorporating these refinements, we can establish a more robust foundation
for assessing the validity of our findings, particularly by comparing them to
lattice calculations of scaled transport coefficients available in Refs. [64,
65, 54, 22]. This comparative analysis encompasses both the pure Yang-Mills
case and also the case with quark effects incorporated. Such endeavors will
contribute to a comprehensive understanding of the transport properties of the
system and provide valuable insights into the interplay between non-
perturbative effects, thermodynamics, and the role of quasi-quarks in the
transport properties of medium QGP.

##  VIII Acknowledgements

The author acknowledges the fruitful discussion with Ritesh Ghosh and Sumit
Rana. N.H. is supported in part by the SERB-MATRICS under Grant No.
MTR/2021/000939.

##  Appendix A Landau-Lifshitz condition: Gribov gluons

Here, we see the effect of the Landau-Lifshitz condition along with the mean-
field effect on the transport coefficients of Gribov-modified gluonic plasma.
This is achieved by following method [43].

For a small deviation away from the equilibrium, the change in distribution
function is

$\displaystyle f(x,p)=f^{eq}(\epsilon_{0})+\delta f(x,p)\leavevmode\nobreak\
,$ (73)

where $\epsilon_{0}$ correspond to the equilibrium energy. In order to capture
the non-equilibrium deviation in the energy, we define $\epsilon$, such that

$\displaystyle\epsilon=\epsilon_{0}+\delta\epsilon$ (74)

The distribution function considering the off-equilibrium shift is

$\displaystyle
f(x,p)=f^{eq}(\epsilon)+\delta\tilde{f}(x,p)\leavevmode\nobreak\ ,$ (75)

where

$\displaystyle\delta\tilde{f}(x,p)=\delta f(x,p)-\frac{\partial
f^{eq}(\epsilon)}{\partial\epsilon}\delta\epsilon\leavevmode\nobreak\ .$ (76)

Now, for out of equilibrium, we defined

$\displaystyle\Delta
T^{ij}=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\textbf{{p}}^{i}\textbf{{p}}^{j}}{\epsilon}\delta\tilde{f}\leavevmode\nobreak\
.$ (77)

Using Eq. (76) in Eq. (77) and a bit of simplification, the shift in energy
density is

$\displaystyle\Delta T^{00}$
$\displaystyle=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\left(\epsilon\delta\tilde{f}+\frac{\partial
f^{eq}(\epsilon)}{\partial\epsilon}\epsilon\delta\epsilon\right)$
$\displaystyle=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\left(\epsilon\delta\tilde{f}+\frac{f^{eq}(\epsilon)\tilde{f}^{eq}(\epsilon)}{T}\epsilon\delta\epsilon\right)$
(78)

Eq. (78) for $\epsilon=\epsilon_{1}$, given by Eq. (19), we get

$\displaystyle\Delta T^{00}$
$\displaystyle=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\Bigg{[}\epsilon_{1}\delta\tilde{f}$
$\displaystyle+\frac{f^{eq}\tilde{f}^{eq}}{2T}\left(m_{g}m_{g}^{\prime}-\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\right)\delta
T\Bigg{]}\leavevmode\nobreak\ ,$ (79)

where, for obvious reason,
$f^{eq}\tilde{f}^{eq}=f^{eq}(\epsilon_{1})\tilde{f}^{eq}(\epsilon_{1})$. Now,
for the deviation away from equilibrium, the distribution function takes the
form

$\displaystyle\delta
f=-f^{eq}(\epsilon_{1})\tilde{f}^{eq}(\epsilon_{1})\left(\frac{\delta\epsilon_{1}}{T}-\frac{\epsilon_{1}}{T^{2}}\delta
T\right)\leavevmode\nobreak\ .$ (80)

Using Eq. (80) in Eq. (76) and with a bit of simplification, we land to, for
$\epsilon=\epsilon_{1}$,

$\displaystyle\delta
T=\frac{T^{2}}{f^{eq}(\epsilon_{1})\tilde{f}^{eq}(\epsilon_{1})\epsilon_{1}}\delta\tilde{f}$
(81)

Now, we substitute Eq. (81) back in Eq. (79)

$\displaystyle\Delta T^{00}_{\epsilon=\epsilon_{1}}$
$\displaystyle=\Delta\epsilon_{1}=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{1}{\epsilon_{1}}$
$\displaystyle\times\left[\epsilon_{1}^{2}-\frac{T}{2}\left(m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\right)\right]\delta\tilde{f}$
(82)

Following the argument of [43], the more generalized form of out-of-
equilibrium energy-momentum tensor, for $\epsilon=\epsilon_{1}$, is

$\displaystyle\Delta T^{\mu\nu}$
$\displaystyle=\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{1}{\epsilon_{1}}\Bigg{[}p^{\mu}p^{\nu}$
$\displaystyle-u^{\mu}u^{\nu}\frac{T}{2}\left(m_{g}m_{g}\prime+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\right)\Bigg{]}\delta\tilde{f}$
(83)

In the context of hydrodynamics, the Landau-Lifshitz condition is often used
to analyze the behaviour of fluctuations in a fluid. When applied to the study
of viscosity, it specifically affects the behaviour of bulk viscosity($\zeta$)
and not shear viscosity($\eta$). As it is a well-established fact that the
shear viscosity($\eta$) is related to the dissipation of momentum due to the
relative motion of adjacent fluid layers, it does not directly involve the
dispersion relation or the speed of sound, $v_{s}^{2}$. Therefore, the Landau-
Lifshitz condition primarily influences the behavior of bulk viscosity, which
is associated with energy dissipation due to fluid volume changes.

One of the subtleties in $\mathcal{A}$ is an infinite number of solutions
could be generated by shifting

$\mathcal{A}(\epsilon_{1})\rightarrow\mathcal{A}^{\prime}(\epsilon_{1})=\mathcal{A}(\epsilon_{1})-a-b\epsilon_{1}\leavevmode\nobreak\
,$

where for zero chemical potential, $a=0$ since it is connected with the
particle number conservation.

For a particular solution, say $\mathcal{A}^{par}_{1}$, we take

$\displaystyle\mathcal{A}=\mathcal{A}^{par}_{1}-b\epsilon_{1}\leavevmode\nobreak\
.$ (84)

Using the Landau-Lifshitz condition, $u_{\mu}\Delta T^{\mu\nu}=0$, we get

$\displaystyle\int$
$\displaystyle\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{f_{g}^{eq}(\epsilon_{1})\tilde{f}_{g}^{eq}(\epsilon_{1})}{\epsilon_{1}}(\mathcal{A}^{par}_{1}-b\epsilon_{1})$
$\displaystyle\times\left[\epsilon_{1}^{2}-\frac{T}{2}\left(m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\right)\right]=0$
(85)

Making use of the thermodynamic relation
$(d\mathcal{P}/dT)=(d\mathcal{P}/d\epsilon)(d\epsilon/dT)=c_{s}^{2}(d\epsilon/dT)$,
we get

$\displaystyle b$
$\displaystyle=\frac{c_{s}^{2}}{T^{2}s}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{f_{g}^{eq}(\epsilon_{1})\tilde{f}_{g}^{eq}(\epsilon_{1})}{\epsilon_{1}}\times$
$\displaystyle\left[\epsilon_{1}^{2}-\frac{T}{2}\left(m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\right)\right]\mathcal{A}^{par}_{1}$
(86)

where,

$\displaystyle
s=\frac{1}{3T^{2}}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\textbf{{p}}^{2}f_{g}^{eq}(\epsilon_{1})\tilde{f}_{g}^{eq}(\epsilon_{1})\leavevmode\nobreak\
.$ (87)

From Eqs. (42), (84) and (86), we get

$\displaystyle\zeta_{gluon}^{1}$
$\displaystyle=\frac{1}{9T}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\tau(\epsilon_{1})}{\epsilon_{1}^{2}}f_{g}^{eq}\tilde{f}_{g}^{eq}\mathcal{A}_{par}^{1}$
$\displaystyle\times\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{1}^{2}+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{m_{g}^{3}m_{g}^{\prime}-4\gamma_{G}^{3}\gamma_{G}^{\prime}}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}\leavevmode\nobreak\
,$ (88)

where
$f_{g}^{eq}\tilde{f}_{g}^{eq}=f_{g}^{eq}(\epsilon_{1})\tilde{f}_{g}^{eq}(\epsilon_{1})$.
Doing the same calculation for $\epsilon=\epsilon_{2}$, we get

$\displaystyle\zeta_{gluon}^{2}$
$\displaystyle=\frac{1}{9T}\int\frac{d^{3}\textbf{{p}}}{(2\pi)^{3}}\frac{\tau(\epsilon_{2})}{\epsilon_{2}^{2}}f_{g}^{eq}\tilde{f}_{g}^{eq}\mathcal{A}_{par}^{2}\bigg{[}\textbf{{p}}^{2}-3c_{s}^{2}\epsilon_{2}^{2}$
$\displaystyle+\frac{3c_{s}^{2}T}{2}\bigg{(}m_{g}m_{g}^{\prime}+\frac{4\gamma_{G}^{3}\gamma_{G}^{\prime}-m_{g}^{3}m_{g}\prime}{\sqrt{m_{g}^{4}-4\gamma_{G}^{4}}}\bigg{)}\bigg{]}\leavevmode\nobreak\
,$ (89)

where
$f_{g}^{eq}\tilde{f}_{g}^{eq}=f_{g}^{eq}(\epsilon_{2})\tilde{f}_{g}^{eq}(\epsilon_{2})$.

## References

  * [1] C. H. Chen, Nucl. Phys. B Proc. Suppl. 246-247, 38-41 (2014)
  * [2] R. Pasechnik and M. Šumbera, Universe 3, 7 (2017) [arXiv:1611.01533 [hep-ph]].
  * [3] P. Braun-Munzinger and J. Stachel, Nature 448, 302-309 (2007)
  * [4] M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30-63 (2005) [arXiv:nucl-th/0405013 [nucl-th]].
  * [5] E. Shuryak, Prog. Part. Nucl. Phys. 62, 48-101 (2009) [arXiv:0807.3033 [hep-ph]].
  * [6] U. W. Heinz, [arXiv:nucl-th/0512051 [nucl-th]].
  * [7] P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005) [arXiv:hep-th/0405231 [hep-th]].
  * [8] A. Jaiswal and V. Roy, Adv. High Energy Phys. 2016, 9623034 (2016) [arXiv:1605.08694 [nucl-th]].
  * [9] P. Arnold, W. Florkowski, Z. Fodor, P. Foka, J. Harris, M. Lisa, H. Meyer, A. Milov, A. Mischke and B. Müller, et al. PoS ConfinementX, 030 (2012)
  * [10] W. Florkowski, A. Jaiswal, E. Maksymiuk, R. Ryblewski and M. Strickland, Phys. Rev. C 91, 054907 (2015) [arXiv:1503.03226 [nucl-th]].
  * [11] W. Florkowski, ‘Phenomenology of Ultra-Relativistic Heavy-Ion Collisions,” World Scientific, 2010.
  * [12] A. Jaiswal and N. Haque, Phys. Lett. B 811, 135936 (2020) [arXiv:2005.01303 [hep-ph]].
  * [13] C. Sasaki and K. Redlich, Phys. Rev. C 79, 055207 (2009) [arXiv:0806.4745 [hep-ph]].
  * [14] G. Torrieri, B. Tomasik and I. Mishustin, Phys. Rev. C 77, 034903 (2008) [arXiv:0707.4405 [nucl-th]].
  * [15] S. Plumari, W. M. Alberico, V. Greco and C. Ratti, Phys. Rev. D 84, 094004 (2011) [arXiv:1103.5611 [hep-ph]].
  * [16] O. Soloveva, P. Moreau and E. Bratkovskaya, Phys. Rev. C 101, 045203 (2020) [arXiv:1911.08547 [nucl-th]].
  * [17] H. Berrehrah, E. Bratkovskaya, T. Steinert and W. Cassing, Int. J. Mod. Phys. E 25, 1642003 (2016) [arXiv:1605.02371 [hep-ph]].
  * [18] S. Mitra and V. Chandra, Phys. Rev. D 96, 094003 (2017) [arXiv:1702.05728 [nucl-th]].
  * [19] P. Singha, A. Abhishek, G. Kadam, S. Ghosh and H. Mishra, J. Phys. G 46, no.1, 015201 (2019) [arXiv:1705.03084 [nucl-th]].
  * [20] J. O. Andersen and M. Strickland, Annals Phys. 317, 281-353 (2005) [arXiv:hep-ph/0404164 [hep-ph]].
  * [21] A. Bazavov et al. [HotQCD], Phys. Rev. D 90, 094503 (2014) [arXiv:1407.6387 [hep-lat]].
  * [22] H. B. Meyer, Phys. Rev. Lett. 100, 162001 (2008) [arXiv:0710.3717 [hep-lat]].
  * [23] H. B. Meyer, Phys. Rev. D 76, 101701 (2007) [arXiv:0704.1801 [hep-lat]].
  * [24] N. Su, Int. J. Mod. Phys. A 30, 1530025 (2015) [arXiv:1502.04589 [hep-ph]].
  * [25] V. N. Gribov, Nucl. Phys. B 139, 1 (1978)
  * [26] D. Zwanziger, Nucl. Phys. B 323, 513-544 (1989)
  * [27] A. Bandyopadhyay, N. Haque, M. G. Mustafa and M. Strickland, Phys. Rev. D 93, 065004 (2016) [arXiv:1508.06249 [hep-ph]].
  * [28] V. Begun, W. Florkowski and R. Ryblewski, Acta Phys. Polon. B 48, 125 (2017) [arXiv:1602.08308 [nucl-th]].
  * [29] N. Su and K. Tywoniuk, Phys. Rev. Lett. 114, no.16, 161601 (2015) [arXiv:1409.3203 [hep-ph]].
  * [30] D. Dudal, C. P. Felix, M. S. Guimaraes and S. P. Sorella, Phys. Rev. D 96, no.7, 074036 (2017) [arXiv:1704.06529 [hep-ph]].
  * [31] A. Bandyopadhyay, [arXiv:2307.09656 [hep-ph]].
  * [32] W. Wu, G. Huang, J. Zhao and P. Zhuang, Phys. Rev. D 107, 114033 (2023) [arXiv:2207.05402 [hep-ph]].
  * [33] F. E. Canfora, D. Dudal, I. F. Justo, P. Pais, L. Rosa and D. Vercauteren, Eur. Phys. J. C 75, no.7, 326 (2015) [arXiv:1505.02287 [hep-th]].
  * [34] E. Gotsman and E. Levin, Phys. Rev. D 103, no.1, 014020 (2021) [arXiv:2009.12218 [hep-ph]].
  * [35] M. Debnath, R. Ghosh and N. Haque, [arXiv:2305.16250 [hep-ph]].
  * [36] D. E. Kharzeev and E. M. Levin, Phys. Rev. Lett. 114, no.24, 242001 (2015) [arXiv:1501.04622 [hep-ph]].
  * [37] S. Madni, A. Mukherjee, A. Bandyopadhyay and N. Haque, Phys. Lett. B 838, 137714 (2023) [arXiv:2210.03076 [hep-ph]].
  * [38] K. Fukushima and N. Su, Phys. Rev. D 88, 076008 (2013) [arXiv:1304.8004 [hep-ph]].
  * [39] D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel and H. Verschelde, Phys. Rev. D 78, 065047 (2008) [arXiv:0806.4348 [hep-th]].
  * [40] G. Burgio, M. Quandt and H. Reinhardt, Phys. Rev. Lett. 102, 032002 (2009) [arXiv:0807.3291 [hep-lat]].
  * [41] M. I. Gorenstein and S. N. Yang, Phys. Rev. D 52, 5206-5212 (1995)
  * [42] V. Mykhaylova and C. Sasaki, Phys. Rev. D 103, no.1, 014007 (2021) [arXiv:2007.06846 [hep-ph]].
  * [43] P. Chakraborty and J. I. Kapusta, Phys. Rev. C 83, 014906 (2011) [arXiv:1006.0257 [nucl-th]].
  * [44] T. Sainapha, [arXiv:1910.11659 [hep-th]].
  * [45] N. Vandersickel and D. Zwanziger, Phys. Rept. 520, 175-251 (2012) [arXiv:1202.1491 [hep-th]].
  * [46] R. F. Sobreiro and S. P. Sorella, [arXiv:hep-th/0504095 [hep-th]].
  * [47] R. F. Sobreiro, S. P. Sorella, D. Dudal and H. Verschelde, Phys. Lett. B 590, 265-272 (2004) [arXiv:hep-th/0403135 [hep-th]].
  * [48] D. Zwanziger, Phys. Rev. D 76, 125014 (2007) [arXiv:hep-ph/0610021 [hep-ph]].
  * [49] A. D. Linde, Phys. Lett. B 96, 289-292 (1980)
  * [50] S. Jeon, Phys. Rev. D 52, 3591-3642 (1995) [arXiv:hep-ph/9409250 [hep-ph]].
  * [51] S. Jeon and L. G. Yaffe, Phys. Rev. D 53, 5799-5809 (1996) [arXiv:hep-ph/9512263 [hep-ph]].
  * [52] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP 11, 077 (2010) [arXiv:1007.2580 [hep-lat]].
  * [53] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP 07, 056 (2012) [arXiv:1204.6184 [hep-lat]].
  * [54] N. Y. Astrakhantsev, V. V. Braguta and A. Y. Kotov, Phys. Rev. D 98, 054515 (2018) [arXiv:1804.02382 [hep-lat]].
  * [55] N. Astrakhantsev, V. Braguta and A. Kotov, JHEP 04, 101 (2017) [arXiv:1701.02266 [hep-lat]].
  * [56] P. B. Arnold, C. Dogan and G. D. Moore, Phys. Rev. D 74, 085021 (2006) [arXiv:hep-ph/0608012 [hep-ph]].
  * [57] P. B. Arnold, G. D. Moore and L. G. Yaffe, JHEP 05, 051 (2003) [arXiv:hep-ph/0302165 [hep-ph]].
  * [58] S. Weinberg, Astrophys. J. 168, 175 (1971)
  * [59] A. Buchel, Phys. Rev. D 72, 106002 (2005) [arXiv:hep-th/0509083 [hep-th]].
  * [60] P. Benincasa, A. Buchel and A. O. Starinets, Nucl. Phys. B 733, 160-187 (2006) [arXiv:hep-th/0507026 [hep-th]].
  * [61] W. Florkowski, R. Ryblewski, N. Su and K. Tywoniuk, Phys. Rev. C 94, no.4, 044904 (2016) [arXiv:1509.01242 [hep-ph]].
  * [62] A. Czajka, S. Hauksson, C. Shen, S. Jeon and C. Gale, Phys. Rev. C 97, no.4, 044914 (2018) [arXiv:1712.05905 [nucl-th]].
  * [63] G. Torrieri and I. Mishustin, Phys. Rev. C 78, 021901 (2008) [arXiv:0805.0442 [hep-ph]].
  * [64] A. Nakamura and S. Sakai, Phys. Rev. Lett. 94, 072305 (2005) [arXiv:hep-lat/0406009 [hep-lat]].
  * [65] S. Sakai and A. Nakamura, PoS LATTICE2007, 221 (2007) [arXiv:0710.3625 [hep-lat]].