Theoretical Perspectives on Viscous Nature of Strongly Interacting Systems

Abstract

Matter prevailing during the early stages of the Universe or under extreme conditions in high-energy heavy-ion experiments supposedly possesses a rich phase structure. During the evolution of such a system, the complicated pictures of transitions among various phases are studied as part of hydrodynamics. This system, on most occasions, is considered to be non-viscous. However, various theoretical studies reveal the importance of incorporating viscous effects into the analysis. Here, the paper discusses the behavioral patterns of transport coefficients with varying temperatures and chemical potentials to obtain a qualitative, if not quantitative, picture of the same. Discussions are also shared regarding their impacts on such an exotic system for different energies, as explored in the experimental domain. This theoretical analysis, made using the structure of the Polyakov–Nambu–Jona-Lasinio (PNJL) model with a 2+1-flavor quark–antiquark system reveals important aspects of the inclusion of viscous effects in the hydrodynamic studies of QGP.

1. Introduction

Strongly interacting matter is believed to have pervaded the Universe after a few microseconds of the Big Bang as well as in the interior of neutron stars. Relativistic heavy-ion collision experiments give us the opportunity to study the corresponding physics. In these experiments, the collision of two heavy ions takes place at relativistic energies to form a deconfined state of matter known as quark–gluon plasma (QGP). Various heavy-ion collision experiments have been conducted over several years in the quest for QGP. In this regard, the role of certain observables, like radial, azimuthal and longitudinal flow, are significant [1,2]. In general, they are computed from azimuthal Fourier coefficients $v_n$ of the triple differential distributions of produced hadrons depending on the range of the impact parameter. For non-central heavy-ion collisions, elliptic flow is the most important. The occurrence of elliptic flow takes place with the plasma’s response to initial pressure gradients. Hydrodynamical evolution leads to the conversion of the initial pressure gradients to final state velocity gradients. The deformation of the initial state cannot be controlled in a heavy-ion collision. Plasma deformation is obtained by the shape of the overlap** region of the colliding nuclei. This is governed by the impact parameter b. Measurement of the impact parameter takes place on an event-by-event basis using azimuthal dependence of the produced particle spectra. Once the direction of the impact parameter is known, particle distribution is expanded in terms of the Fourier components of the azimuthal angle. The deformation of the final state is given by the Fourier coefficients ($v_2$, $v_4$, etc.) [3]. A positive value of $v_2$ indicates the preferential emission of particles in the short direction, thus pointing toward the occurrence of elliptic flow. The presence of shear viscous coefficients is believed to oppose the elliptic flow, which makes it necessary to incorporate it. Bulk viscosity also has a significant effect on elliptic flow and is found to oppose it. Thus, for a proper characterization of the behavior in high-energy collisions, it is extremely essential to quantify both shear and bulk viscous coefficients [4,5].

In this context, it is noteworthy that although most studies were carried out for infinite systems [5], it is very important to consider the finiteness of systems, which is believed to have an impact on the viscous effects [6,7,8,9,10,11,12,13,14,15,16,17,18]. The fireball formed in high-energy heavy-ion collisions is conjectured to have a finite spatial extent depending on certain factors, like colliding nuclei, the center of mass energy and collision centrality. Thus, these effects are believed to have an impact on the transport coefficients of a system.

From a theoretical front, various QCD-inspired models, like the Nambu–Jona-Lasinio (NJL), Polyakov-loop-enhanced Quark Meson (PQM), etc., provide a solid theoretical baseline for these studies, barring the complexities of the lattice at finite chemical potentials. Here, this paper outlines one such QCD-inspired model, namely, the Polyakov–Nambu–Jona-Lasinio (PNJL) model [19,20,21,22,23].

2. The Model Framework

The Polyakov–Nambu–Jona-Lasinio model, also called the PNJL model, was made by adding the Polyakov or P-loop to the NJL model, thus suitably incorporating both chiral and deconfinement transitions within a single framework. In the NJL model, the multi-quark interaction terms account for the interaction between quarks, obeying the global symmetries of QCD [24,25,26,27,28,29,30,31,32]. Dynamical fermion mass generation gives rise to the spontaneous breaking of chiral symmetry. With the gluon degrees of freedom integrated out, this model fails to encapsulate the very important crucial feature of deconfinement. The PNJL model, by appending the P-loop and thus bringing in a background temporal gluonic field, appends the important feature of the confinement–deconfinement transition. Considering the ${SU\left(3\right)}_f$ version of the model kee** up to eight-quark interactions, the Lagrangian takes the form:

$$\begin{gathered} \mathcal{L}=\sum _{f=u, d, s}{\overline{\Psi }}_f{\gamma }_{\mathrm{\mu }}i{D}^{\mathrm{\mu }}{\Psi }_f-\sum _{f=u, d, s}{m}_{\mathrm{f}}{\overline{\Psi }}_f{\Psi }_f+\sum _{f=u, d, s}\begin{array}{c}{\mu \gamma }_0{\overline{\Psi }}_f{\Psi }_f\\ \end{array}+\frac{g_S}{2} \sum _{a=0,...8}\left[{\left(\overline{\Psi }{\lambda }^a\Psi \right)}^2+{\left(\overline{\Psi }{i{\gamma }_5\lambda }^a\Psi \right)}^2\right] \\ -g_D\left[\det \left({\overline{\Psi }}_f{P}_{\mathrm{L}}{\Psi }_{f^{\prime }}\right)+\det \left({\overline{\Psi }}_f{P}_{\mathrm{R}}{\Psi }_{f^{\prime }}\right)\right]+8g_1{\left[\right({\overline{\Psi }}_i{P}_{\mathrm{R}}{\Psi }_m\left)\right({\overline{\Psi }}_m{P}_{\mathrm{L}}{\Psi }_i\left)\right]}^2 \\ +16g_2\left[\left({\overline{\Psi }}_i{P}_{\mathrm{R}}{\Psi }_m\right)\left({\overline{\Psi }}_m{P}_{\mathrm{L}}{\Psi }_j\right)\left({\overline{\Psi }}_j{P}_{\mathrm{R}}{\Psi }_k\right)\left({\overline{\Psi }}_k{P}_{\mathrm{L}}{\Psi }_i\right)\right]-{\mathcal{U}}^{\prime }\left(\mathsf{\Phi }\left[\mathrm{A}\right],\overline{\mathsf{\Phi }}\left[\mathrm{A}\right],\mathrm{T}\right) \end{gathered}$$

(1)

where $P_{\mathrm{L},\mathrm{R}}$ are the chiral projector matrices given by ${(1\pm \gamma }_5)/2$. The first term in Equation (1) denotes the Dirac term with the gauge field interactions $D^{\mathrm{\mu }}={\partial }^{\mathrm{\mu }}-i{A}_4{\delta }_{\mathrm{\mu }4}$, with $A_4$ being the time component of the Euclidean gauge field. This further helps to define the Polyakov line as $L\left(\overrightarrow{x}\right)=\mathcal{P}\mathrm{e}\mathrm{x}\mathrm{p}\left[i{\int }_0^{1/T}d\tau A_4 (\overrightarrow{x},t)\right]$, with $\mathcal{P}$ denoting the path ordering. The corresponding Polyakov loop field $\mathsf{\varphi }$ and its conjugate $\overline{\varphi }$ are thereby defined as:

$$\mathsf{\varphi }=\frac{{Tr}_{c}L}{N_c}, \overline{\varphi }=\frac{{Tr}_c{L}^{†}}{N_c}$$

The second and third term in Equation (1) incorporates the mass effect breaking the symmetry explicitly. Subscript ‘f’ denotes the flavors up (u), down (d) or strange (s). The fourth term represents the four-quark interaction with coupling $g_S$. The floowing one is for six-quark interactions, which remains invariant under ${SU\left(3\right)}_{L}X {SU\left(3\right)}_R$ but breaks ${U\left(1\right)}_A$, mimicking the chiral anomaly. The next two terms account for the spin zero eight-quark interactions, with $g_1$ and $g_2$ being the couplings. $g_1$ and $g_2$ are set to zero in the absence of such interactions. Local interactions are approximated here. As the model is non-renormalizable, the couplings are dimensionful. To take care, a three-momentum cut-off is used over the quark loops to make them finite.

${\mathcal{U}}^{\prime },$ the Polyakov loop potential, is expressed as [33]:

$$\frac{{\mathcal{U}}^{\prime }[\phi ,\overline{\phi }, T]}{T^4}=\frac{\mathcal{U}[\phi ,\overline{\phi }, T]}{T^4}-\kappa \ln \left(J\right[\phi ,\overline{\phi }\left]\right)$$

(2)

The second term on the right side is the Vandermonde term [23], which reflects the effect of the SU(3) Haar measure. Here, $\mathcal{U}$ is the Landau–Ginzburg type potential, which is chosen to be:

$$\frac{\mathcal{U}[\phi ,\overline{\phi }, T]}{T^4}=-\frac{b_2\left(T\right)}{2} \overline{\phi }\phi -\frac{b_3}{6} \left({\phi }^3+{\overline{\phi }}^3\right)+\frac{b_4}{4} {\left(\phi \overline{\phi }\right)}^2$$

(3)

Here, the coefficients $b_3$ and $b_4$ are constants, whereas b₂ shows a temperature dependence of the form:

$$b_2\left(T\right)=a_0+a_1\left(\frac{T_0}{T}\right)+a_2{\left(\frac{T_0}{T}\right)}^3$$

(4)

The associated parameters were set with physical constraints, and some were set by fitting the available lattice QCD results [33]. The set of parameters can be found in [34].

Using dynamical fermion mass generation, the current quark mass gains a constituent value given by:

$$M_f=m_f-g_S{\sigma }_f+\frac{g_D}{2}{\sigma }_{f+1}{\sigma }_{f+2}-2g_1{\sigma }_f\left({\sigma }_u^2+{\sigma }_d^2+{\sigma }_s^2\right)-4g_2{\sigma }_f^3$$

(5)

where σ_f is the chiral condensate given by ${\sigma }_f=<{\overline{\Psi }}_f{\Psi }_f>$ and ${\sigma }_f={\sigma }_u, {\sigma }_d, {\sigma }_s$ in a cyclic order.

In order to move to the regime of finite system sizes, it is essential to choose the correct boundary conditions, i.e., periodic for bosons and antiperiodic for fermions. This, in turn, leads to a sum of infinite extent over the discretized momentum values of $p_i=\frac{\pi n_i}{R}$. R is the dimension of cubical volume, and n_i are positive integers. Most of the studies in the realm of QCD-inspired models were undertaken with infinite system extents. In this study, as one of the first of its kind regarding the finite extents of systems, a few simplifications were considered. Surface and curvature effects were neglected. The infinite sum over the discrete momentum values was replaced by integration over the continuum momentum variation with an infrared cut-off. Also, there was no amendment in the mean-field values due to the finite sizes. The volume of the finite system, considered V, was on the same footing as T and the chemical potential µ.

The thermodynamic potential in the realm of the finite sizes is then obtained as:

$$\begin{array}{ll}\mathsf{\Omega }={\mathcal{U}}^{\prime }\left(\mathsf{\Phi }\left[\mathrm{A}\right],\overline{\mathsf{\Phi }}\left[\mathrm{A}\right],\mathrm{T}\right)& +2{\mathrm{g}}_{\mathrm{S}}{\displaystyle \sum _{\mathrm{f}=\mathrm{u},\mathrm{d},\mathrm{s}}}{{\mathsf{\sigma }}_{\mathrm{f}}}^2-\frac{{\mathrm{g}}_{\mathrm{D}}}{2}{\mathsf{\sigma }}_{\mathrm{u}}{\mathsf{\sigma }}_{\mathrm{d}}{\mathsf{\sigma }}_{\mathrm{s}}+3\frac{g_1}{2}{\left({{\mathsf{\sigma }}_f}^2\right)}^2+3g_2{{\mathsf{\sigma }}_f}^4\\ & -6{\displaystyle \sum _{\mathrm{f}}}{\int }_{\mathsf{\lambda }}^{\mathrm{ᴧ}}\frac{{\mathrm{d}}^3\mathrm{p}}{{\left(2\mathsf{\pi }\right)}^3}{\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}\mathsf{\Theta }\left(\mathsf{\Lambda }-\left|\overrightarrow{\mathrm{p}}\right|\right)\\ & -2{\displaystyle \sum _{\mathrm{f}}}\mathrm{T}{\int }_{\mathsf{\lambda }}^{\mathrm{\infty }}\frac{{\mathrm{d}}^3\mathrm{p}}{{\left(2\mathsf{\pi }\right)}^3}\mathrm{l}\mathrm{n}\left[1+3\left(\mathsf{\Phi }+\overline{\mathsf{\Phi }}{\mathrm{e}}^{-\frac{\left({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}-{\mathsf{\mu }}_{\mathrm{f}}\right)}{\mathrm{T}}}\right){\mathrm{e}}^{\frac{-({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}-{\mathsf{\mu }}_{\mathrm{f}})}{\mathrm{T}}}+{\mathrm{e}}^{\frac{-3({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}-{\mathsf{\mu }}_{\mathrm{f}})}{\mathrm{T}}}\right]\\ & -2{\displaystyle \sum _{\mathrm{f}}}\mathrm{T}{\int }_{\mathsf{\lambda }}^{\mathrm{\infty }}\frac{{\mathrm{d}}^3\mathrm{p}}{{\left(2\mathsf{\pi }\right)}^3}\mathrm{l}\mathrm{n}\left[1+3\left(\overline{\mathsf{\Phi }}+\mathsf{\Phi }{\mathrm{e}}^{-\frac{\left({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}+{\mathsf{\mu }}_{\mathrm{f}}\right)}{\mathrm{T}}}\right){\mathrm{e}}^{\frac{-({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}+{\mathsf{\mu }}_{\mathrm{f}})}{\mathrm{T}}}+{\mathrm{e}}^{\frac{-3({\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}+{\mathsf{\mu }}_{\mathrm{f}})}{\mathrm{T}}}\right]\end{array}$$

(6)

The single quasiparticle energy ${\mathrm{E}}_{{\mathrm{p}}_{\mathrm{f}}}$ is given by $\sqrt{p^2+{M_f}^2}$. The three-momentum cut-off Λ was used for the vacuum integral to remove discrepancies over the non-renormalizability of the model. The eight-quark interaction terms were incorporated to stabilize the vacuum [35,36].

3. Transport Coefficients

Transport coefficients can be obtained using different methods. Here, different methods are reviewed, and comparative studies are provided alongside as necessary. One popular method is the use of Kubo formalism [37]. This method relates the viscous coefficients to the correlation functions of the energy–momentum tensor. The shear viscous coefficient in this method is given by:

$$\eta \left(\omega \right)=\frac{1}{15 T} {\int }_0^{\infty }dt e^{i\omega t}\int d\overrightarrow{r}(T_{\mu \nu }\left(\overrightarrow{r},t\right),T^{\mu \nu }(\mathrm{0,0}))$$

(7)

where $\mathcal{L}$ is denoted as the PNJL Lagrangian and $T_{\mu \nu }=i\overline{\Psi } {\gamma }_{\mu }{\partial }^{\nu }\Psi -g_{\mu \nu }\mathcal{L}$, as the (µ,ν) component of the energy–momentum tensor, one component form of which can also be used to write the Kubo formula. Thus, neglecting the surface terms at infinity, the following form can be obtained using the retarded correlators as:

$$\eta \left(\omega \right)=\frac{i}{\omega } [{П}^R(\omega )-{П}^R(0)]$$

(8)

Considering the retarded correlator with one component energy–momentum tensor, the static shear viscosity can be expressed as [38]:

$$\eta =-\frac{d}{d\omega } Im {П}^R(\omega ){|}_{\omega =0}$$

(9)

In order to calculate the value of the retarded correlator, Matsubara formalism can be used. The interaction kernel can then be organized with n-loop ring diagrams considering large N_c-limit and scalar-pseudoscalar channels as in the PNJL model. However, the correlator in Matsubara space can still be restricted to a one-loop diagram effectively using one component form of the energy–momentum tensor. The frequency summation is evaluated using spectral representation of the full propagator, where the effect of the background Polyakov loop is experienced through the PNJL distribution function as:

$$f_{\Phi }^{+}\left(E_p\right)=\frac{\left(\overline{\Phi }+2\Phi e^{-\beta (E_p+\mu )}\right)e^{-\beta (E_p+\mu )}+e^{-3\beta (E_p+\mu )}}{1+3\left(\overline{\Phi }+\Phi e^{-\beta \left(E_p+\mu \right)}\right)e^{-\beta \left(E_p+\mu \right)}+e^{-3\beta (E_p+\mu )}}$$

(10)

$$f_{\Phi }^{-}\left(E_p\right)=\frac{\left(\Phi +2\overline{\Phi }{e}^{-\beta (E_p-\mu )}\right)e^{-\beta (E_p-\mu )}+e^{-3\beta (E_p-\mu )}}{1+3\left(\Phi +\overline{\Phi }{e}^{-\beta \left(E_p-\mu \right)}\right)e^{-\beta \left(E_p-\mu \right)}+e^{-3\beta (E_p-\mu )}}$$

(11)

where ‘+’ and ‘−’ refer to particle and anti-particle, respectively. Expressing the spectral function in terms of the advanced and retarded Green function [39,40], the trace can then be evaluated to obtain [41]:

$$\eta \left[\Gamma \left(p\right)\right]=\frac{16 N_c{N}_f}{15 {\pi }^{3}T} {\int }_{-\infty }^{\infty }dϵ {\int }_0^{\infty }dp p^6 \frac{M^2{\Gamma }^2\left(p\right)f_{\Phi }\left(\epsilon \right)(1-f_{\Phi }(\epsilon \left)\right)}{{\left(\right({\epsilon }^2-p^2-M^2+{\Gamma }^2\left(p\right))}^2+4M^2{\Gamma }^2\left(p\right){)}^2}$$

(12)

Here, N_c and N_f are the number of colors and flavors, respectively. M is the thermal constituent quark mass, as given by Equation (5).

Using the formalism of Kubo, one can also obtain the bulk viscous coefficient in terms of the correlation functions of the energy–momentum tensor [42,43]. Bulk viscosity reads as:

$$\zeta =\underset{\omega \to 0}{\lim }\frac{-Im G^R(\omega ,0)}{9\omega }=\underset{\omega \to 0}{\lim }\frac{\pi \rho (\omega ,0)}{9\omega }$$

(13)

The imaginary part of the retarded Green function is Im G^R, and the associated spectral density ρ is connected as [42,43,44,45]:

$$-Im G^R\left(\omega ,0\right)=\pi \rho \left(\omega ,0\right)={\int }_0^{\infty }dt\int d^{3}r e^{i\omega t}<\left[{\theta }_{\mu }^{\mu }\left(x\right),{\theta }_{\mu }^{\mu }\left(0\right)\right]>$$

(14)

Following the structures of the QCD low-energy theorem and thus the traces of energy–momentum tensor from [45] and using the Kramers–Kronig relation and basic thermodynamical relations, one can construct the retarded Green function neglecting the divergent contributions [46,47]. Associated ansatz for spectral density is then used as in [44] and similarly excludes the high-frequency perturbative continuum:

$$\frac{\rho \left(\omega ,0\right)}{\omega }=\frac{9\zeta }{\pi } \frac{{\omega }_0^2}{{\omega }_0^2+{\omega }^2}$$

(15)

with ${\omega }_0$ being the mass scale at par to the region of validity for the corresponding perturbation theory. Thus, the bulk viscosity coefficient is obtained as [48]:

$$\zeta =\frac{1}{9{\omega }_0} [Ts\left(\frac{1}{c_s^2}-3\right)+\left(\mu \frac{\partial }{\partial \mu }-4\right)T^5\frac{\partial \left(\frac{P}{T^4}\right)}{\partial T}+\left(T\frac{\partial }{\partial T}+\mu \frac{\partial }{\partial \mu }-2\right){<m\overline{q}q>}_T+6\left(f_{\pi }^2 M_{\pi }^2+f_K^2{M}_K^2\right)+16 |{ϵ}_v|]$$

(16)

Alongside, the expression of bulk viscosity using the kinetic theory approach based on relaxation time approximation (RTA) can also be obtained. The quasi-particle approach of Kubo formalism, like before, can be deployed for the same, expressing the spectral density with one-loop diagrams [49,50,51,52,53]. The finite thermal width $\Gamma =1/\tau$ is then employed in the internal line of the diagram. For zero chemical potentials, this expression reads:

$$\zeta =\frac{12}{T} \int \frac{d^{3}k}{{\left(2\pi \right)}^3} \frac{f_{\Phi }[1-f_{\Phi }]}{{\left(E_k\right)}^2\Gamma } {\left[\left(\frac{1}{3}-c_s^2\right)k^2-c_s^2\frac{d}{d{\beta }^{2 }}({\beta }^2{m}^2)\right]}^2$$

(17)

τ is the in-medium relaxation time, which is inversely determined by $\Gamma$. Here, the effect of the Polyakov loop is absorbed in the PNJL distribution function f_Φ, which is the modified Fermi–Dirac distribution function, where the effects of the background Polyakov-loop fields are taken [54], as in Equations (10) and (11). Further studies in this direction have been made by several groups, as can be found in [55,56,57,58,59,60,61,62].

It is significant enough to see these viscosities in the system at a finite volume [63]. Alongside shear and bulk viscosities, the electrical conductivity σ also plays a major role in the system’s behavior, thus making it important to study. One obtains the standard expression of σ, as in [63]:

$$\sigma =\frac{6{\breve{e}}_Q^2}{3T} \int \frac{d^3\overrightarrow{k}}{{\left(2\pi \right)}^3} {\tau }_Q{\left(\frac{\overrightarrow{k}}{{\omega }_Q}\right)}^2\left[f_Q^{+}\left(1-f_Q^{+}\right)+f_Q^{-}(1-f_Q^{-})\right]$$

(18)

4. Results

4.1. Shear Viscosity

Figure 1 shows the behavior of η with temperature at the zero chemical potential. With exact time evolution absent in this framework, T and µ serve as the two essential variables, and changes in these should mimic the evolution of the system. Variation can be seen with various choices of Г maintaining the convergence of η [64,65]. The forms thus chosen are of constant, exponential, Lorentzian and divergent natures under investigation:

$${\Gamma }_{const}=100 \mathrm{M}\mathrm{e}\mathrm{V}$$

$${\Gamma }_{exp}(p)={\Gamma }_{const}{e}^{-\frac{\beta p}{8}}$$

$${\Gamma }_{lor}\left(p\right)={\Gamma }_{const}\frac{\beta p}{1+{\left(\beta p\right)}^2}$$

$${\Gamma }_{div}\left(p\right)={\Gamma }_{const}\sqrt{\beta p}$$

It is seen that the shear viscous effect increases with temperature T. This behavior is found to be similar to a gaseous system discussed in [41]. The inset portrays the behavior of a two-flavor system kee** Г fixed at 100 MeV.

Variation with the quark chemical potential µ_q is shown in Figure 2. It is enticing to see the change at three different temperatures corresponding to the first-order phase transition, crossover and beyond crossover regions. In the case T = 100 MeV, if one moves along the µ-axis, the first-order phase transition line is expected to be encountered. This is reflected in the results, where η shows a jump at µ_q = 280 MeV. In juxtaposition, as expected, for T = 150 and 200 MeV, a smooth variation is seen along the µ_q direction.

The difference is more visible if the shear viscous coefficient is plotted as a function of the quark chemical potential with different temperatures as shown in Figure 3. One thus obtains a glimpse of the phase diagram from the behavioral change in η.

The value of specific shear viscosity (η/s), i.e., the shear viscosity to entropy ratio, has a very special significance, which calls for the need to study its behavior. For relativistic fluids, the Reynolds number is defined in terms of η/s, where s is entropy density. The variation in η/s with temperature and chemical potentials gives us an insight into the phase diagram and possible location of the critical endpoint. It is seen that η/s at first decreases with increasing temperature and reaches a minimum near the transition point for the corresponding value of µ_q. It then increases with temperature slowly after transition, achieving a stable value or slowly decreasing henceforth. The minimum value occurs close to the Kovtun–Son–Starinet (KSS) bound for a value of a quark chemical potential within 100 to 150 MeV, which was predicted to be the lower bound for $\frac{\mathsf{\eta }}{s}(\ge \frac{\mathrm{h}}{4\pi k_B})$ by Kovtun et al. [5]. For µ_q greater than or equal to 200 MeV, the system goes through a transition at lower temperatures, and η/s is shown to jump, thus pointing toward a first-order phase transition. Thus, it may be conjectured that this gives a probable location to the existence of a Critical Endpoint (CEP). Figure 4 shows the behavior of η/s as a function of temperature at the vanishing chemical potential and for different constant values of Г. The minima occur at the crossover temperature with the overall qualitative nature similar for all Г. The inset signifies the region of the crossover transition. Fixing Г at 100 MeV, the specific shear viscosity is then analyzed in Figure 5 for various chemical potentials covering the phase diagram. Figure 5a shows continuous behavior, thus indicating the crossover region, whereas, in Figure 5b, the two chemical potentials fall in the domain of the first-order transition.

The locations of minima and/or discontinuities thus help to obtain an idea of the phase diagram for the matter under concern, as shown in Figure 6. The location of the Critical Endpoint (CEP) is conjectured as the point from where discontinuities start appearing in the pattern of η/s.

To comprehend the behavior of the strongly interacting system near the CEP, studying observables, like specific heat $C_V$, can give extra impetus. $C_V$ is defined as:

$$C_V=\frac{\partial ϵ}{\partial T}=T\frac{{\partial }^{2}P}{\partial T^2}=T\frac{\partial s}{\partial T}$$

$C_V$, as the second-order derivative of pressure, is supposed to show a distinctive nature near the second-order phase transition. This could significantly contribute toward understanding the CEP region, where the transition is also supposedly of the second order. Its nature has been verified in [41]. The scaled specific heat, when studied in regions of crossover, the CEP, first order and beyond, is seen to show varied behavior, thus pointing toward the existence of the CEP.

To keep a juxtaposition with the experimental scenario, it becomes important to study observables as a function of center-of-mass collision energy ($\sqrt{s}$). Using freeze-out parametrization, the shear-viscosity-to-entropy ratio was studied with collision energy. As shown in Figure 7, a higher value of Г aids in reaching the KSS bound more easily.

4.2. Bulk Viscosity

The behavior of $\frac{\zeta }{s}$ can be obtained using two different approaches from the QCD sum rule and the quasi-particle approach, as discussed in [48]. From both approaches, it is seen that $\frac{\zeta }{s}$ decreases when T > T_c. The behavior in the temperature domain lower than the transition region is seen to quantitatively vary in the two approaches. In the high-temperature regime, $\frac{\zeta }{s}$ from both approaches is seen to approach zero. This gives the indication that at high-temperature QCD, the massless limit is attained from conformal symmetric behavior.

Its nature has been thoroughly discussed in [48]. In Figure 8, the primary reasons for peak-like behavior in both cases are attributed to the dependencies on quantities like (ε-3P) and $<m\overline{q}q>$. The first quantity depends on c_s², the square of the speed of sound, a deviation of which, from the classical limit (~1/3), gives the violation of conformality. The corresponding nature and influence of constituent quark masses are also discussed there. In both the approaches, the parameters were chosen kee** in parity the numerical strength of the bulk viscous coefficient. Changes in the behavioral patterns for ζ and ζ/s (specific bulk viscosity) for non-zero chemical potentials are also observed, as shown in Figure 9a,b. The peaks in ζ occur close to the transition points along the phase diagram shown above. ζ/s, on the other hand, decreases monotonically before approaching the classical limit of zero for all the cases under scrutiny.

The corresponding nature from various other calculations was compared vis-à-vis to obtain at least a qualitative idea in [48]. This also helps to characterize the importance of including the bulk viscous effects in the hydrodynamic evolution of strongly interacting matter. Apart from the mismatches in the two approaches, both of them indicate the significance of ζ in characterizing the phase diagram, where it increases quantitatively.

Moving on, the nature of shear and bulk viscosity changes significantly under the finite size approximation of the system under concern. The primary role is played by the constituent masses of up (or down) and strange quarks, as shown in Figure 10.

∆M was calculated using the difference of masses in the thermodynamic limit and that for finite system sizes. Changes in entropy under similar circumstances were also observed and are shown in Figure 11a. This led to the observation of η and η/s for the finite system volumes, as portrayed in Figure 11b. The negative values clearly indicate the increase in η for finite system sizes below a certain temperature.

Bulk viscosity, on the other hand, depends additionally on the slope of M_u,d and M_s vs. the T graph. The second quantity responsible for the observed behavior of ζ is the conformal breaking quantity involving c_s². They cumulatively influence bulk viscosity, which takes the shape shown in Figure 12a,b. The dominance of entropy density in the high-temperature domain is of crucial importance here because the second peak in ζ almost dilutes away in ζ/s.

As a significant part of the family of transport coefficients, thermal conductivity σ was also explored in [63]. As observed for low temperatures, with a decrease in the system size, shear viscosity and electrical conductivity increase. Bulk viscosity, on the other hand, measuring the scale violation of the medium, has a non-trivial connection with system size. The effect of the Polyakov loop in this scenario makes a significant contribution. This acts behind the confinement–deconfinement symmetry and thus alters the results toward more accuracy. This was also verified in Figure 13 of [63] when comparing similar quantities under similar conditions as obtained from the NJL and PNJL models at a finite system size and the thermodynamic result. It is seen that the transport coefficients obtained are significantly enhanced in the PNJL model compared to the NJL model and from infinite to finite system sizes. This increase occurs mainly in low-temperature regimes and almost vanishes at high temperatures.

5. Conclusions

This article collectively addresses the physics for the transport coefficients of the exotic matter of QGP, as obtained under the framework of a PNJL model. The shear and bulk viscous coefficients, in this regard, might play crucial roles in understanding the flow patterns and hydrodynamic behavior of such systems. Kubo formalism is one very effective way to address them theoretically alongside the RTA mechanism. Based on this analysis, it is interesting to see that $\frac{\eta }{s}$ calculated in the background of the PNJL model, when studied along the freeze-out curve, agrees with the RHIC and LHC flow analysis. It is also significant to see how the results reflect the conjectured position of the CEP. As argued, this can also be confirmed using the analysis of the corresponding observables showing the transition of the second order near the critical point. The discontinuities in the associated quantity have been used to capture the transition locations providing a unique method of study. The finite density results of bulk viscosity also reflect the fact of the QCD phase diagram that the transition occurs at low temperatures as the quark chemical potential increases when the peak position in ζ shifts accordingly with changes in the chemical potential. The entire picture changes meaningfully when investigating the regime of finite system sizes. The constituent strange quark mass, however, dilutes more slowly than expected and thereby leads to a possible model artifact, which will be addressed elsewhere. Also, it is noteworthy that the absence of the hadronic degrees of freedom explicitly in the low-temperature domain makes the results qualitative to some extent. Proper incorporation of them in finite temperature and chemical potential windows of interest will set quantitative accuracy more profoundly and will be performed in the future.