---
abstract: 'A large-$N_{\rm c}$ expansion is combined with the Kubo formalism to study the shear viscosity $\eta$ of strongly interacting matter in the two-flavor NJL model. We discuss analytical and numerical approaches to $\eta$ and investigate systematically its strong dependence on the spectral width and the momentum-space cutoff. Thermal effects on the constituent quark mass from spontaneous chiral symmetry breaking are included. The ratio $\eta/s$ and its thermal dependence are derived for different parameterizations of the spectral width and for an explicit one-loop calculation including mesonic modes within the NJL model.'
author:
- Robert Lang
- Wolfram Weise
date: 'March 26, 2014'
title: 'Shear viscosity from Kubo formalism: NJL-model study'
---

Introduction
============

Heavy-ion collisions at RHIC [@BRAHMSatRHIC05; @PHENIXatRHIC05; @STARatRHIC05; @PHOBOSatRHIC05] and at the LHC [@AamodtALICE10; @AamodtALICE11; @CaffarriALICE12] explore strongly interacting matter under extreme conditions. The quark-gluon matter created in such collisions at temperatures exceeding $T_{\rm c}\approx 0.2$ GeV leaves its indirect signatures in the produced particles at lower temperatures long after thermalization. Transport properties such as the shear viscosity $\eta$ of the highly excited matter are of prime interest in this context. Inertial anisotropies in the collision plane translate into non-trivial particle flow patterns [@DerendarzATLAS2013; @RoyMohantyChaudhuri2013; @DuslingTeaney2008; @LuzumGombeaudOllitrault2010]. In particular the elliptic-flow parameter, $v_2$, features a strong dependence on the ratio of shear viscosity to entropy density, $\eta/s$, of the dissipative quark-gluon matter formed in the collision. It is known that the temperature dependence of $\eta/s$ is crucial in order to describe the elliptic flow of hadrons in ultrarelativistic heavy-ion collisions at RHIC and LHC [@Niemi2011]. At the very high LHC energies, the temperature dependence of $\eta/s$ for $T>T_{\rm c}$ becomes dominant for the elliptic flow.

Two basic approaches are commonly used to deal with non-equilibrium systems and transport properties: the Boltzmann equation [@HuiDefuLic2006; @SasakiRedlich2010; @KhvorostukhinToneevVoskresensky2011] (most frequently applied in relaxation-time approximation) and the Kubo formalism using retarded correlators of the energy-momentum tensor [@Alberico2008; @HidakaKunihiro2011; @HidakaKunihiro2011NJL; @NamKao2013]. When examined in comparison, the kinetic approach seems generally to underestimate the shear viscosity [@PlumariEtAl2012; @PlumariEtAl2013]. The calculation in relaxation-time approximation makes use of the thermal cross sections of the colliding particles, whereas the Kubo formalism asks for the spectral functions of the basic degrees of freedom. These two approaches are connected via the optical theorem. In the present paper we choose the Kubo formalism, realizing at the same time that, in this approach, a perturbative treatment of transport coefficients is insufficient [@JeonSkeleton1995; @JeonYaffeSkeleton1996] and requires resummation techniques even in a weak-coupling situation.

We use the Nambu–Jona-Lasinio (NJL) model [@Nambu1961tp; @Nambu1961; @Klevansky1992; @KLVW1990; @KLVW1990:2; @VoglWeise1991; @BuballaHabil2005] as a schematic, non-perturbative approach to the thermodynamics of quark matter. Gluonic degrees of freedom are integrated out and hidden in point vertices of the effective interaction between quarks, while all relevant chiral and flavor symmetries and symmetry-breaking patterns of QCD are taken properly into account. The applicability of such a model is supposed to cover a temperature range $T_{\rm c} \lesssim T < \Lambda$, where $T_{\rm c}\approx 0.2$ GeV is the transition temperature from the hadronic to the quark phase and $\Lambda\approx 0.6\;{\rm GeV}$ is the characteristic NJL cutoff scale. We combine the NJL model with a large-${N_{\rm c}}$ expansion [@tHooft1974; @QuackKlevansky94; @BuballaMuellerWambach2010] to study the shear viscosity of quark matter.

In Section \[Sec:EtaLO\] the shear viscosity $\eta$ in leading order is deduced. We follow mainly the developments in [@Fukutome2006; @Fukutome2008Nucl; @Fukutome2008Prog] but do not restrict ourselves to the chiral limit and arrive at results under more general assumptions. Taking only the dominant scalar and pseudoscalar channels into account, an infinite number of ring diagrams reduces to just one single generic diagram. Corrections to correlation functions by ladder diagrams are suppressed in a large-${N_{\rm c}}$ expansion. However, for the shear viscosity itself a resummation of these subleading diagrams is potentially important, depending on the ${N_{\rm c}}$ scaling of the spectral width [@HidakaKunihiro2011; @HidakaKunihiro2011NJL]. This effect has been studied first in [@JeonSkeleton1995] for a bosonic field theory. In the present work resummations in the Kubo sector will not be included. A discussion concerning the conditions under which such resummations are necessary will however be given.

The general derivation is followed by a detailed parameter study in Section \[Sec:ParStud\]: for $\eta[\Gamma]$ as a functional of the quasiparticle spectral width $\Gamma(p)$ of the quarks. First an analytical result is derived assuming a constant spectral width to start with. Furthermore, the implementation of different parameterizations for $\Gamma(p)$ teaches us about the general dependence of the shear viscosity on the spectral width for a variety of examples. A strong dependence on the pertinent momentum-space cutoff, $\Lambda$, is found, reflecting the sensitivity to physical scales: the characteristic cutoff fixed by the NJL gap equation excludes up to $90\%$ of the mathematically accessible high-momentum contributions to $\eta$, a feature that actually turns out to be a prerequisite for achieving physically meaningful results within this framework. We also investigate to what extent the functional $\eta[\Gamma(p)]$ can be treated perturbatively by expanding in a Laurent series for a “small” spectral width and comparing with the full result. The impact of thermal constituent quark masses on the shear viscosity is investigated in Section \[Sec:ThermoMass\]. The thermal quark masses are generated dynamically by the spontaneous chiral symmetry breaking mechanism through the NJL gap equation. Not surprisingly, we find that thermal effects are crucial to obtain physically relevant results for the shear viscosity.

In Section \[Sec:GammaNJL\] an explicit calculation of the spectral width is performed within the NJL model, using the one-loop mesonic contributions to the quark self-energy at next-to-leading order in the large-${N_{\rm c}}$ expansion. Two different physical effects contribute to this width: Landau damping and mesonic recombination. For temperatures well above the critical/crossover temperature the resulting spectral width decreases, implying an increasing shear viscosity $\eta(T)$ in this temperature range.

Shear viscosity at leading order {#Sec:EtaLO}
================================

In this work we model quark matter starting from the two-flavor NJL Lagrangian $$\label{NJL2}
  \mathcal{L}={\bar{\psi}}\left({{\rm i}}{\slashed{\partial}}-\hat{m}\right)\psi+\frac G2\left[({\bar{\psi}}\psi)^2+({\bar{\psi}}{{\rm i}}\gamma_5{\mbox{\boldmath $\tau$}}\psi)^2\right],$$ where $\psi = (u,d)^{\rm T}$ is the isospin doublet quark field, $\hat{m} = {\rm diag}(m_{\rm u},m_{\rm d})$ is the current quark mass matrix (we work in the isospin limit, $m_{\rm u}= m_{\rm d} \equiv m$), $G$ denotes the scalar/pseudoscalar coupling, and ${\mbox{\boldmath $\tau$}}$ collects the three Pauli isospin matrices. Vector or axialvector terms are not considered in this work. The large masses of the corresponding quark-antiquark modes make their contributions to the relevant correlation functions far less important than those of pseudoscalar and scalar modes.

In the Kubo formalism [@Kubo57] transport coefficients are related to retarded correlators of energy-momentum tensors, i.e. to four-point functions in Matsubara space. The energy momentum tensor of the NJL model is simply $$\label{Tmunu}
T_{\mu\nu}={{\rm i}}{\bar{\psi}}\gamma_\mu\partial_\nu\psi-g_{\mu\nu}\mathcal{L}$$ in terms of the quark fields $\psi$. The Kubo formula for the shear viscosity reads $$\begin{aligned}
  \label{KuboEta}
\eta(\omega) =
  \frac{\beta}{15}\int_0^\infty{\text{d}}t\;{{\rm e}}^{{{\rm i}}\omega t}\int{\text{d}}^3r\; (T_{\mu\nu}({\mbox{\boldmath $r$}},t),T^{\mu\nu}(0))\;,
\end{aligned}$$ with $\beta = 1/T$ the inverse temperature. The correlator $(X,Y)$ is defined by[^1] $$(X,Y)=T\int_0^\beta{\text{d}}\xi\;\langle{{\rm e}}^{\xi H}X{{\rm e}}^{-\xi H}Y\rangle\;,$$ where $H$ is the NJL Hamiltonian and $\langle\mathcal{\cdot}\rangle={{\rm Tr}}\left(\cdot\,{{\rm e}}^{-\beta H}\right)$ denotes the thermal expectation value. An equivalent reduced expression for the shear viscosity is also frequently used in the literature [@Alberico2008]: $$\label{KuboEtaAlternative}
  \eta(\omega)=\beta\int_0^\infty{\text{d}}t\;{{\rm e}}^{{{\rm i}}\omega t}\int{\text{d}}^3r\; (T_{21}({\mbox{\boldmath $r$}},t),T_{21}(0))\;,$$ written in terms of only one component of the energy-momentum tensor. The relative factor $15$ in comparison with Eq. results from the following identity ($i,j\in\{1,2,3\}$ with $i\neq j$): $$\int{\text{d}}^3x\; x_i^2\,x_j^2\,f(x^2)=\frac{1}{15}\int{\text{d}}^3x\;x^4f(x^2)\;.$$ The (classical) components of the energy-momentum tensor are real quantities, $T_{\mu\nu}\in\mathds{R}$. It follows that the *static* shear viscosity $\eta(\omega = 0)$ is also real: $$\eta(\omega)^*=\eta(-\omega)\;\;\;\Rightarrow\;\;\;\eta:=\eta(0)\in\mathds{R}\;.$$ Neglecting surface terms at infinite time, one derives $$\label{DefRetGreenFctFromDifference}
  \eta(\omega)=\frac{{{\rm i}}}{\omega}\left[\Pi^{\rm R}(\omega)-\Pi^{\rm R}(0)\right],$$ with the retarded correlation function $$\label{DefPiROmega}
\Pi^{\rm R}(\omega)=-{{\rm i}}\int_0^\infty{\text{d}}t\;{{\rm e}}^{{{\rm i}}\omega t}\int{\text{d}}^3r\,\langle[T_{21}({\mbox{\boldmath $r$}},t),T_{21}(0)]\rangle\;.$$ The static shear viscosity ($\omega\to 0$) follows as $$\eta=-\left.\frac{{\text{d}}}{{\text{d}}\omega}\,{{\rm Im}\,}\Pi^{\rm R}(\omega)\right|_{\omega=0}\;.$$ The calculation of the retarded correlator can be performed switching to the Matsubara formalism and calculating $$\Pi(\omega_n)=\int_0^\beta{\text{d}}\tau\;{{\rm e}}^{{{\rm i}}\omega_n\tau}\int{\text{d}}^3 r\; \langle \mathcal{T}_\tau
\big(T_{21}({\mbox{\boldmath $r$}},\tau)T_{21}(0)\big)\rangle\;,$$ where have applied a Wick rotation $\tau={{\rm i}}t$ and introduced the time-ordering symbol in imaginary time, $\mathcal{T}_\tau$. Note that whereas the underlying (quark) degrees of freedom are fermionic, the Matsubara frequencies relevant for the correlator $\Pi$ are bosonic, $\omega_n=2\pi n T$, since the fermion fields under the integral group together to form quantities of bosonic character: ${\bar{\psi}}(\cdot)\psi$. The global sign of $\Pi(\omega_n)$ is fixed by the sign convention for analytical continuations: $$\left.\Pi(\omega_n)\right|_{i\omega_n=\omega\pm{{\rm i}}{\varepsilon}}= -\,\Pi^{\rm R/A}(\omega)\;,$$ where the upper and lower sign in $\pm{{\rm i}}{\varepsilon}$ corresponds to the retarded and advanced correlation function, respectively. The correlator $\Pi(\omega_n)$ is governed by non-perturbative physics resulting from the underlying interactions of the NJL model. We now apply a large-${N_{\rm c}}$ expansion and organize this correlator in ring diagrams, ladder diagrams and higher-order terms: $$\label{KuboNcExpansionFourPointFunction}
  \begin{minipage}{0.08\textwidth} \mbox{$\Pi(\omega_n)=$} \end{minipage}  
  \begin{minipage}{0.15\textwidth}
    \includegraphics[width=\textwidth]{NJLsk4pointfctFull.pdf}
  \end{minipage}
  \hspace{-0.2cm} \begin{minipage}{0.2\textwidth}
  \mbox{$= \mathcal{O}(N_{\rm c}^1)+\mathcal{O}(N_{\rm c}^0)+ \ldots $}
\end{minipage}$$ The four-point coupling of the NJL Lagrangian effectively incorporates gluonic degrees of freedom resulting in the scaling $G\sim 1/{N_{\rm c}}$. At leading order $\mathcal{O}(N_{\rm c}^1)$ there is just a one-loop diagram contributing to the four-point correlator, given that the NJL Lagrangian in its simplest form takes into account only scalar and pseudoscalar interactions: $\Gamma\in\{{{\mathds 1}},{{\rm i}}\gamma_5\}$. Iterating these interaction kernels in ring diagrams at leading order in $1/{N_{\rm c}}$ to $\Pi(\omega_n)$ does not affect the correlator: $$\label{Chain}
\hspace{0.5cm} \begin{minipage}{0.04\textwidth}
\end{minipage}
  \begin{minipage}{0.3\textwidth}
    \includegraphics[width=\textwidth]{NJLsk4pointfctABsnd.pdf}
  \end{minipage}
\begin{minipage}{0.05\textwidth} \!=\;0\;, \end{minipage}$$ because the trace (in momentum and Dirac space) in the first ring vanishes due to the orthogonal operator structure involving the combination of $\gamma_2$ and $\Gamma$: $$\begin{aligned}
  \label{KuboRingAZero1}
  &{T\sum_{n\in\mathds{Z}}\int\frac{{\text{d}}^3 p}{(2\pi)^3}}\; {{\rm Tr}}\left[\gamma_2\, G_\beta({\mbox{\boldmath $p$}},\nu_n)\, \Gamma\, G_\beta({\mbox{\boldmath $p$}},\nu_n)\right] =\\
  &  ={T\sum_{n\in\mathds{Z}}\int\frac{{\text{d}}^3 p}{(2\pi)^3}}\; \frac{1}{(\nu_n^2+{\mbox{\boldmath $p$}}^2+M^2)^2}\; \\
  &\hspace{0.5cm} \times\, {{\rm Tr}}\left[\gamma_2\Gamma M^2 +\gamma_2\slashed{p}\Gamma\slashed{p}+\gamma_2\slashed{p}\Gamma M+\gamma_2\Gamma\slashed{p} M\right]=0\;,
\end{aligned}$$ where we have used the notation $\slashed{p}=\nu_n\gamma_4-{\mbox{\boldmath $p$}}\cdot{\mbox{\boldmath $\gamma$}}$ and the full Matsubara propagator $$G_\beta({\mbox{\boldmath $p$}},\nu_n)=\frac{\slashed{p}+M}{\nu_n^2+{\mbox{\boldmath $p$}}^2+M^2}\;,$$ with frequencies $\nu_n=(2n+1)\pi T-{{\rm i}}\mu$. Exchange (ladder diagram) corrections to the chain in Eq.  are non-vanishing but of subleading order in $1/{N_{\rm c}}$, because each rank in the ladder gives rise to a suppression factor $G^2{N_{\rm c}}\sim 1/{N_{\rm c}}$. Note that adding one rank introduces two additional momentum integrations but only one additional color trace.

The shear viscosity in the NJL model has been deduced previously in Refs. [@Fukutome2006; @Fukutome2008Nucl; @Fukutome2008Prog] using the Kubo formula, but assuming the quarks to be in the chiral limit, $m=0$. We point out that this result can in fact be derived without assuming to work in the chiral limit. Setting the current quark masses to zero is *not necessary* to ensure the absence of iterated ring-diagram contributions when taking only scalar and pseudoscalar interactions of the NJL model into account. Iterated ring diagrams involving these interactions vanish naturally. (Note that even in the chiral limit and the Nambu-Goldstone phase, the second term of the trace in Eq. , ${{\rm Tr}}\left[\gamma_2\slashed{p}\Gamma\slashed{p}\right]$ would survive in the presence of vector interactions, but their contribution to the correlators would be small as mentioned before).

Collecting all arguments, we can summarize in general: for purely fermionic theories $\mathcal{L}=\mathcal{L}_{\rm kin}+\mathcal{L}_{\rm int}$ with momentum-independent pseudoscalar/scalar interactions and $2n$-vertices that scale as $G_{2n}\sim 1/{N_{\rm c}}^{n-1}$, the dominant contribution to the correlation function $\Pi^{\rm R}(\omega)$ in Matsubara space is: $$\label{PiRingLO}
  \begin{minipage}{0.08\textwidth} \mbox{$\Pi(\omega_n)=$} \end{minipage}  
  \begin{minipage}{0.125\textwidth}
    \includegraphics[width=\textwidth]{NJLsk4pointfctA0snd.pdf}
  \end{minipage}
  \begin{minipage}{0.1\textwidth} \mbox{$\;+\;\mathcal{O}({N_{\rm c}}^0)\;.$} \end{minipage}$$ With the definition of the spectral function, $$\rho(\omega,{\mbox{\boldmath $p$}})=-\frac{1}{\pi}{{\rm Im}\,}G^{\rm R}(\omega,{\mbox{\boldmath $p$}})=\frac{1}{2\pi{{\rm i}}}\left(G^{\rm A}(\omega,{\mbox{\boldmath $p$}})-G^{\rm R}(\omega,{\mbox{\boldmath $p$}})\right),$$ and using residue calculus one derives: $$\begin{aligned}
\label{EtaLeadingNcAsTwoSpectralDensityPap}
  \eta &={\pi\over T} \int_{-\infty}^\infty {\text{d}}{\varepsilon}\int\frac{{\text{d}}^3 p}{(2\pi)^3}\;p_x^2\, n_{\rm F}^+({\varepsilon})\big(1-n_{\rm F}^+({\varepsilon})\big) \\
  &\hspace{2.5cm} \times\,{{\rm Tr}}\left[\gamma_2\,\rho({\varepsilon},{\mbox{\boldmath $p$}})\,\gamma_2\,\rho({\varepsilon},{\mbox{\boldmath $p$}})\right],
\end{aligned}$$ with the Fermi-Dirac distribution $$n_{\rm F}^+(E)=\frac{1}{1+{{\rm e}}^{\beta (E-\mu)}}\;.$$ As in [@Fukutome2006] the dressed quark propagator is written as $$\label{KuboQuasiPartApprQuarkProp}
G^{\rm R/A}(p_0,{\mbox{\boldmath $p$}})=\frac{1}{\slashed{p}-M\pm{{\rm i}}\,{\rm sgn}(p_0)\Gamma(p)}\;,$$ with the quasiparticle mass $M$ and width $\Gamma(p)$. The next step is to relate this spectral width to the shear viscosity $\eta$. Even in the chiral limit the dynamical NJL mechanism of spontaneous chiral symmetry breaking generates a large constituent quark mass in the vacuum: $M\approx 0.3\;{\rm GeV}$, see the brief discussion in Section \[Sec:ThermoMass\]. Apart from this mechanism, the thermal environment at temperature $T$ and baryo-chemical potential $\mu$ of the quarks affects parameterically both the dynamical quark mass $M(T,\mu)$ and the spectral width $\Gamma(p;T,\mu)$.

The spectral function $\rho$ is represented in the standard form of a generalized Breit-Wigner shape, to be inserted in Eq.: $$\begin{aligned}
\label{EtaLeadingNcAsSpectralWidth}
  \eta[\Gamma(p)] &=\frac{16{N_{\rm c}}{N_{\rm f}}}{15\pi^3 T}\int_{-\infty}^\infty{\text{d}}{\varepsilon}\int_0^\infty{\text{d}}p\,p^6 \\
  &\times\,\frac{M^2\,\Gamma^2(p)\,n_{\rm F}^+({\varepsilon})(1-n_{\rm F}^+({\varepsilon}))}{\left[({\varepsilon}^2-p^2-M^2+\Gamma^2(p))^2+4M^2\Gamma^2(p)\right]^2}\;.
\end{aligned}$$ Even though we started from the NJL model in our derivation, the expression is generic for a system of strongly interacting Fermions, with real and imaginary parts of their self-energies encoded in $M$ and $\Gamma$, respectively. In general, the non-perturbative origin of $\Gamma(p)$ does not permit expanding the functional $\eta[\Gamma]$ in a Laurent series. In perturbative approaches (e.g. in chiral perturbation theory at low temperatures) such a treatment is possible: $\eta[\Gamma]\sim 1/\Gamma$ [@LangKaiserWeise2012].

In this context we comment briefly on the issue of ladder resummation and its correction to Eq. . If one assumes the spectral width $\Gamma\sim 1/{N_{\rm c}}$ to be suppressed for large ${N_{\rm c}}$ as suggested by hot-QCD calculations, then the superficial ${N_{\rm c}}$ counting of Eqs.  and , $\eta\sim{N_{\rm c}}$, is spoiled: in this case, the integrand becomes highly singular in the large-${N_{\rm c}}$ limit (“pinch poles” as described in [@JeonSkeleton1995; @HidakaKunihiro2011]), resulting in an additional factor ${N_{\rm c}}$, therefore $\eta\sim{N_{\rm c}}^2$. We find in the limit ${N_{\rm c}}\to\infty$, i.e. $\Gamma\to 0$: $$\label{EtaSmallGamma}
  \eta[\Gamma]\to\frac{2{N_{\rm c}}{N_{\rm f}}}{15\pi^2 T}\int_{|{\varepsilon}|>M}{\text{d}}{\varepsilon}\,\frac{({\varepsilon}^2-M^2)^{5/2}\,n_{\rm F}^+({\varepsilon})(1-n_{\rm F}^+({\varepsilon}))}{M\,\Gamma(\sqrt{{\varepsilon}^2-M^2})}\;.$$ In contrast to the non-perturbative result in Eq.  the ${\varepsilon}$-integration excludes the region $|{\varepsilon}|<M$. This is due to the delta functions appearing in the limit of small $\Gamma$. The momentum integration of the integrand involving $\delta({\varepsilon}^2-p^2-M^2)$ is readily carried out.

Exploratory studies of the shear viscosity {#Sec:ParStud}
==========================================

Analytical results
------------------

Consider now first the case of a constant spectral width, $\Gamma={\rm const}$. This rough schematic approximation allows for an analytical treatment of the momentum integral in Eq.  and one is left with the numerical ${\varepsilon}$-integration only. For the momentum integral we use the following identity: $$\begin{aligned}
\label{KuboAnalyitcalPIntegral}
  \int_0^\infty &{\text{d}}p\;\frac{p^6}{[(A-p^2)^2+B^2]^2}= \\
  &\hspace{-0.5cm}=\frac{\pi}{8\sqrt{2}}\frac{\sqrt{\sqrt{A^2+B^2}-A}}{B^4} \left[(2A^2+3B^2)\sqrt{A^2+B^2} \right.\\
  &\hspace{3.7cm} \left. +2A(A^2+2B^2)\right],
\end{aligned}$$ where we have introduced $A={\varepsilon}^2-M^2+\Gamma^2$ and $B=2M\Gamma$. This result is found by extending the integration region to negative $p$ (the integrand is an even function of $p$) and using residue calculus. Having performed the $p$-integration analytically reduces the computation time by roughly one order of magnitude. Furthermore, it helps finding an appropriate approximation scheme for the whole $({\varepsilon},p)$-integration when the spectral width is momentum dependent.

Fig. \[Fig1\] shows the results for $\eta$ assuming $\Gamma={\rm const.}$ For $\Gamma\to 0$ the shear viscosity diverges, as it follows from Eq. . This limit describes a system of free quarks for which the mean free path is infinite. With increasing temperature and chemical potential, the shear viscosity increases, but the dependence on temperature is more pronounced. Compare these figures to those in Ref. [@Fukutome2006], where $\eta(\Gamma)$ has been evaluated numerically without a momentum-space cutoff, equivalent to our analytical approach based on Eq. .

Inspecting the detailed behavior of the integrand in Eq. , a convergence criterion for the shear viscosity in the absence of a momentum-space cutoff can be derived:

> In order for the shear viscosity $\eta[\Gamma]$ as functional of $\Gamma(p)$ to be convergent, the asymptotic $\Gamma(p)$ should not converge too rapidly to zero: $$\label{CriterionEtaGamma}
>   \eta[\Gamma(p)]<\infty \;\;\; \Leftrightarrow \;\;\; p^3{{\rm e}}^{-\beta p/2}\in {\rm o}(\Gamma(p))\;,$$

where ${\rm o}(\cdot)$ denotes the little Landau symbol[^2]. Possible parameterizations of $\Gamma(p)$ satisfying this constraint are: $$\begin{aligned}
  \label{ModelsMomentumDepGamma}
  {\rm constant:}~~~~~\Gamma_{\rm const} &= 100\;{\rm MeV}\;,\\
  {\rm exponential:}~~~\Gamma_{\rm exp}(p) &= \Gamma_{\rm const}\,{{\rm e}}^{-\beta p/8}\;,\\
  {\rm Lorentzian:}~~~\Gamma_{\rm Lor}(p) &= \Gamma_{\rm const}\,\frac{\beta p}{1+(\beta p)^2}\;,\\
  {\rm divergent:}~~~\Gamma_{\rm div}(p) &= \Gamma_{\rm const}\,\sqrt{\beta p}\;.
\end{aligned}$$ Note that all these parameterizations lead to a finite shear viscosity and no mathematical regularization must be applied, compare the cutoff discussion in Section \[SectionCutOff\]. The particular shapes of these prototype widths have been chosen because of their different behavior at small and large momenta: vanishing or non-vanishing $\Gamma(p=0)$, convergent or divergent $\Gamma(p)$ for $p\to\infty$. These prototypes represent physical spectral widths in several theories [@LangKaiserWeise2012]: $\Gamma(p)$ in $\phi^4$ theory, for instance, is a monotonous function and converges to zero for large momenta. This can be described by the Lorentz parameterization for large momenta: $\lim_{p\to\infty}\Gamma_{\rm Lor}(p)\sim T/p$. In contrast, the spectral width of an interacting pion gas diverges for $p\to\infty$.

Numerical approximation scheme
------------------------------

Our numerical approximation of

![Accuracy of the numerical approximation scheme for the ${\varepsilon}$-integral, Eq. : to reach an accuracy of $10^{-4}$ it is sufficient to restrict its range to $|{\varepsilon}(p)|<1.3\,{\varepsilon}^*(p)$, see Eq. .[]{data-label="Fig2"}](Fig2.pdf){width="45.00000%"}

$\eta[\Gamma(p)]$ is based on the observation that its integrand typically ranges over many orders of magnitude. For every momentum $p$ there is a maximum of the integrand in Eq. , located at the denominator’s minimizer $$\label{EpsStarMaximizer}
  {\varepsilon}^*(p)=\sqrt{p^2+M^2(T,\mu)-\Gamma^2(p;T,\mu)}\;.$$ Adaptive methods do not work when facing a sharp peak structure: either the step size becomes too small for fast convergence (or convergence at all), or the most important contribution in the vicinity of the peak is not sampled by a step size that is too coarse. We overcome this numerical issue by cutting the ${\varepsilon}$-integration and allowing only $|{\varepsilon}(p)|<x{\varepsilon}^*(p)$ for some $x\gtrsim 1$. In comparison with the analytical result for $\Gamma(p)={\rm const.}$ we find that $x=1.3$ is sufficient to produce accurate results within a relative error of $10^{-4}$, see Fig. \[Fig2\].

For the momentum-dependent parameterizations of $\Gamma(p)$ the integrands for $\eta$ in Eq.  look qualitatively the same as for a constant spectral width. We therefore expect the described numerical scheme to work well also in these and more physical cases, where full momentum dependence and effects from the thermal environment are (parameterically) taken into account.

Cutoff dependence {#SectionCutOff}
-----------------

Generally, the shear viscosity increases when the spectral width decreases, compare Eq. . This behavior is also visible in Fig. \[Fig3\](a) when comparing our different parameterizations of $\Gamma(p)$: the “more divergent” the spectral width as $p\rightarrow \infty$, the smaller the corresponding shear viscosity: $$\eta_{\rm Lor}>\eta_{\rm exp}>\eta_{\rm const}>\eta_{\rm div}\,,$$ using notations as in Eq. . This sequence is implied by the corresponding (inverse) order for the spectral widths. These arguments hold also for non-vanishing chemical potentials. Assuming the spectral width to be independent of the chemical potential as in our parameterizations of $\Gamma(p)$ in Eq. , the shear viscosity increases for increasing $\mu$, but the qualitative shape of $\eta(T)$ does not change. We note that the results in Fig. \[Fig3\] have been derived using a constant constituent quark mass $M=325\;{\rm MeV}$, see the brief discussion in Section \[Sec:ThermoMass\].

The integrand of $\eta[\Gamma(p)]$, Eq. , is sizable for unphysically large momenta, so we expect a strong cutoff dependence. In the NJL model the quasiparticle interactions are restricted to quark momenta $p\le\Lambda=650\;{\rm MeV}$. Quarks with momenta $p>\Lambda$ do not interact and have infinite mean free paths. Retricting the momentum integration to the interval $p \le \Lambda$, we find a shear viscosity as shown in Fig. \[Fig3\](b). Excluding $p>\Lambda$ reduces the shear viscosity by one order of magnitude at low temperatures and even by two orders of magnitude at high $T$. As expected, this expresses a very strong cutoff dependence. In addition to these quantitative differences, the qualitative behavior of the shear viscosity also changes strongly and flattens for high temperatures.

This strong cutoff dependence is investigated in more detail in Fig. \[Fig4\]: the contributions taken into account (compared to the analytical result for $\eta$) depend strongly on temperature and just weakly on the chemical potential. At $T=200\,{\rm MeV}$ the momentum cutoff excludes about $90\%$ of the full integral extended to infinity, see Fig. \[Fig4\](a). As shown in Fig. \[Fig4\](b), varying the cutoff by up to $\pm 20\%$ implies for $\eta$ a change of up to $100\%$.

To assess the order of magnitude of the NJL shear viscosity, a comparison with $\eta(T)$ for other systems is instructive. For example, an interacting pion gas treated within the framework of chiral perturbation theory [@LangKaiserWeise2012] has a typical shear viscosity of order $\eta(T)\approx 40\,{\rm MeV}/{\rm fm}^2 \approx 1.6\cdot 10^{-3}\,{\rm GeV}^3$ at $T\approx 100\;{\rm MeV}$. This is a similar order of magnitude as the results shown in Fig. \[Fig3\](b) when applying the NJL cutoff $\Lambda=650\,{\rm MeV}$. We recall that this cutoff is fixed by reproducing physical observables such as the pion decay constant in vacuum. A physically meaningful order of magnitude for $\eta$ then follows naturally.

Perturbative aspects of [$\eta[\Gamma]$]{} {#Sec:PertEta}
------------------------------------------

![Scaling of $\eta\cdot\Gamma$ for different $T$ and $\mu$ as function of the inverse width expressed in units of the pion mass $m_\pi$. Solid horizontal lines correspond to the residues $A_{-1}$ of $\eta[\Gamma]$ in Eq. . A constituent quark mass $M=100\;{\rm MeV}$ has been used for convenience.[]{data-label="Fig5"}](Fig5.pdf){width="46.00000%"}

We have already mentioned that the shear viscosity $\eta$ diverges for non-interacting systems, i.e. for a vanishing spectral width, corresponding to infinite mean free path. Close to this limit $\eta$ can be expanded in a Laurent series (as realized for example analytically in ChPT and $\lambda\phi^4$ theory [@LangKaiserWeise2012]): $$\label{LaurentExpEtaGamma}
  \eta[\Gamma]=\frac{A_{-1}}{\Gamma}+A_0+A_1\Gamma+A_2\Gamma^2+\ldots$$ For small $\Gamma$, the combination $\eta\cdot\Gamma$ is just the residue $A_{-1}$. What does “small” mean in this context? In contrast to $\lambda\phi^4$ theory where $\Gamma\sim\lambda^2$, the NJL model is generically non-perturbative in its coupling, even though the scaling $G\sim 1/{N_{\rm c}}$ applies. The spectral width is therefore not expected to be sufficiently small in order to permit an expansion as in Eq. . Fig. \[Fig5\] shows results of the fully non-perturbative calculation of $\eta\cdot\Gamma$ as a function of the inverse width, conveniently written as $x=m_\pi/\Gamma$, at different $T$ and $\mu$ in comparison with the residue $A_{-1}$. As it can be seen from the figure, corrections to the leading term of the Laurent series are small for $x>1.5$ (demanding $10\%$ accuracy or better). From these considerations we conclude that a perturbative approach is justified only for spectral widths $\Gamma\ll\,m_\pi=140\;{\rm MeV}$.

The discussion of a perturbative treatment of $\eta[\Gamma(p)]$ is closely related to the resummation of ladder diagrams: if in the large-${N_{\rm c}}$ limit the spectral width decreases, i.e. $\Gamma\sim 1/{N_{\rm c}}$ as suggested by hot-QCD calcuations where the coupling $\alpha_{\rm s}\sim 1/{N_{\rm c}}$ becomes small, then the perturbative regime is reached in this limit and the Laurent series expansion in can be restricted to its leading-order term. As seen from Eq. , for a constant but small spectral width $\Gamma$ the residue $A_{-1}$ can be identified with the remaining ${\varepsilon}$-integral. While only this residue term of $\eta[\Gamma(p)]$ is relevant in this case, ladder diagrams now become sizable corrections and contribute also at leading order. Furthermore, the shear viscosity now scales as $\eta\sim{N_{\rm c}}^2$ and no longer linearly with ${N_{\rm c}}$ as Eqs.  and do for ${N_{\rm c}}$-independent spectral function $\
rho$ and width $\Gamma$, respectively.

We conclude that ladder diagram resummations are necessary in the perturbative regime of $\eta[\Gamma(p)]$ in Eq. , i.e. when the spectral width is small, $\Gamma\ll m_\pi$. In the NJL model with its genuine non-perturbative structure, the physical spectral width is large and outside the perturbative regime. Ladder diagram resummations are subleading corrections, while the shear viscosity functional is valid also for large spectral width when including all orders of the Laurent series expansion .

Effects of thermal quark masses and the ratio [$\eta/\lowercase{s}$]{} {#Sec:ThermoMass}
======================================================================

The constituent quark mass has so far been treated as a constant. We now proceed to incorporate its explicit $T$ and $\mu$ dependence. One of the key features of the NJL model is the spontaneous breaking of chiral symmetry: ${\rm SU}(2)_{\rm L}\times{\rm SU}(2)_{\rm R}\to{\rm SU}(2)_{\rm V}$. In addition, chiral symmetry is explicitly broken by the non-zero current quark mass $m$. Solving the NJL gap equation results in a dynamically generated constituent quark mass [@Nambu1961tp; @Nambu1961; @VoglWeise1991]: $$\begin{aligned}
\label{ThermalGapEquation}
  &M(T,\mu)=m-G\langle\bar{\psi}\psi\rangle =\,\\  
  &\!\!=m+\frac{2G{N_{\rm c}}M(T,\mu)}  {\pi^2}\int_0^\Lambda{\text{d}}p\,\frac{p^2}{E_p}\,\big[1-n_{\rm F}^+(E_p)-n_{\rm F}^-(E_p)\big],
\end{aligned}$$ with $E_p(T,\mu)=\sqrt{p^2 + M^2(T,\mu)}$ and the Fermi-Dirac distribution for quarks and antiquarks $$\label{DefFermiDistrMuT}
  n_{\rm F}^\pm(E)=n_{\rm F}(E\mp\mu)=\frac{1}{{{\rm e}}^{\beta (E\mp\mu)}+1}\;.$$ In the vacuum, $(T,\mu)=(0,0)$, the mass is determined to be about one third of the nucleon mass, $M=325\;{\rm MeV}$, where the input parameters of the NJL Lagrangian are chosen as $G=10.08\;{\rm GeV^{-2}}$, $m=5.5\;{\rm MeV}$ and $\Lambda=650\;{\rm MeV}$. This parameter set produces physical (vacuum) values of the pion mass, $m_\pi=140\;{\rm MeV}$, the pion decay constant, $f_\pi=94\;{\rm MeV}$, and the chiral condensate $\langle{\bar{\psi}}\psi\rangle=-(316.4\;{\rm MeV})^3$.

Fig. \[Fig6\](a) shows the shear viscosity $\eta$ for varying constituent quark mass $M$ treated as a parameter, assuming a constant spectral width $\Gamma_{\rm const}=100\;{\rm MeV}$. For $M\to 0$ the shear viscosity becomes divergent, again due to “pinch poles” appearing in Eq.  in this limit. In fact, the origin of this divergence is the same as for $\Gamma\to 0$, since $M$ and $\Gamma$ are formally (almost) interchangeable in the integrand of Eq. . For large constituent quark masses, two effects occur: first, the maximizer ${\varepsilon}^*(p)\sim M$ in Eq.  moves to larger values, and second, the integrand scales as $M^{-6}$. Both features result in a decreasing function $\eta(M)$.

Taking the full thermal dependence of the constituent quark mass into account has an essential influence on the shear viscosity, see Fig. \[Fig6\](b): for small $T$, a constant mass $M=325\;{\rm MeV}$ approximates the thermal constituent quark mass. In contrast, at large $T$, with a melting chiral condensate, the dropping dynamical quark mass implies a strongly increasing shear viscosity, qualitatively different from the case with constant quark mass. (We have chosen to compare thermal and non-thermal results for constant and exponential parameterizations of the spectral width. For $\Gamma_{\rm Lor}$ and $\Gamma_{\rm div}$ the results are qualitatively similar and therefore not shown.)

The shear viscosity itself is a dimensionful quantity. One usually compares $\eta$ to the entropy density $s$, in terms of the dimensionless ratio $\eta/s$. The corresponding ratio extracted from heavy-ion collisions is the smallest value of $\eta/s$ found so far in nature [@HeinzShenSong2012]. The entropy density of quark matter is given by [@Kapusta]: $$\begin{aligned}
  s(T,\mu)=\frac{{N_{\rm c}}{N_{\rm f}}}{\pi^2}\int_0^\infty{\text{d}}p\,p^2\, &\left[\ln\left(1+{{\rm e}}^{-\beta E_p^+}\right) \right.\\
  & \left. \hspace{-4cm}+\ln\left(1+{{\rm e}}^{-\beta E_p^-}\right)+\frac{\beta E_p^+}{1+{{\rm e}}^{\beta E_p^+}}+\frac{\beta E_p^-}{1+{{\rm e}}^{\beta E_p^-}}\right],
\end{aligned}$$ with $E_p^\pm=E_p\mp\mu$. The momentum integration ranges up to infinity. It can be performed without regularization. In the NJL model the thermal constituent quark mass, for momenta above the cutoff scale $\Lambda=650\;{\rm MeV}$, reduces to the bare current quark mass, $M\rightarrow m=5.5\;{\rm MeV}$. At low temperatures, $T\lesssim 150$ MeV, confinement implies that quarks are not the relevant physical degrees of freedom any more and the NJL model cannot be expected to give a realistic description of transport properties.

The evaluation of the ratio $\eta/s$ is shown in Fig. \[Fig7\] for ${N_{\rm f}}=2$ and ${N_{\rm c}}=3$, for the cases of a constant and an exponentially damped spectral width. The comparison between the panels (a) and (b) clearly demonstrates that implementing thermal masses is crucial in order to avoid an unphysical, continuously decreasing ratio $\eta/s(T)$. From experimental data [@Lacey07] and lattice calculations [@Nakamura05; @Meyer07] it is in fact suggested that this ratio increases for $T> 200\;{\rm MeV}$. We also compare to the benchmark $\eta/s=1/4\pi$ from AdS/CFT correspondence [@Maldacena99; @Kovtun05]. This value is known not to be a universal bound; it can be undershot in some field theories [@Cohen07; @Rebhan12; @Mamo12]. This is also found in Fig. \[Fig7\](a), additionally to the unphysical evolution of $\eta/s$ with increasing temperature. A constant quark mass $M=100\;{\rm MeV}$ has been chosen there for convenience. Its vacuum value, $M=325\;{\rm MeV}$, would reduce the scale of the $\eta/s$ ratio even more, compare Fig. \[Fig6\](a). However, taking the thermal constituent quark masses into account leads to an increasing function $\eta/s(T)$, see Fig. \[Fig7\](b). Despite the fact that $\eta$ itself rapidly increases at high $T$, the ratio $\eta/s$ flattens in that region. This flattening is expected to be shifted to higher temperatures for more realistic forms of the width $\Gamma(p)$ such the Lorentzian. In the considered temperature range $\eta/s$ stays above the AdS/CFT benchmark for all parameterizations of the width $\Gamma$.

Spectral width at one-loop level {#Sec:GammaNJL}
================================

So far we have discussed the impact of the shape of the spectral width on the shear viscosity of quark matter, its strong sensitivity to the NJL cutoff and to the thermal constituent quark masses. In this section an explicit calculation is performed including one-loop mesonic contributions at next-to-leading order in the large-${N_{\rm c}}$ expansion [@BuballaMuellerWambach2010]. The leading-order gap equation is modified by the mesonic insertions $$\label{NJLMesonInsertion}
\Sigma^{\rm S/P}_\beta({\mbox{\boldmath $p$}},\nu_n) = \includegraphics[width=0.15\textwidth]{NJLMesonInsertion.pdf}\;,$$ where $\nu_n=(2n+1)\pi T-{{\rm i}}\mu$ are the Matsubara frequencies for quarks. The corresponding expression for antiquarks is easily obtained by inserting $\bar{\nu}_n=\nu_n^*$. The spectral width is extracted from the imaginary part of this self-energy[^3]: $$\Gamma^{\rm S/P}_{\rm q}(p_0,{\mbox{\boldmath $p$}})=-\lim_{{\varepsilon}\to 0}\,{\rm Im}\,\Sigma^{\rm S/P}_\beta({\mbox{\boldmath $p$}},-{{\rm i}}p_0+{\varepsilon})\;.$$ It is of next-to-leading order using the mesonic modes generated by the Bethe-Salpeter equation (BSE) in the pertinent quark-antiquark channels (e.g. [@HatsudaKunihiro87; @KLVW1990; @BuballaMuellerWambach2010]). Whereas the gap equation is of leading order, $\mathcal{O}(1)$ in $N_c$, the BSE is of order $\mathcal{O}(1/{N_{\rm c}})$. In the two-flavor case the mesonic loop involves the three pions and the sigma mode which contribute with positive and negative signs, respectively, to the spectral width $\Gamma(p)$. This “antiscreening” by the sigma mode [@QuackKlevansky94] cancels roughly one of the three pionic contributions to the spectral width, as seen in Fig. \[Fig9\]a. The effective spectral width (for $N_f = 2$) reads: $$\label{QuarkSpectralWidthOneLoopAllMesons}
 \Gamma(p)=3\Gamma_{\rm q}^{\rm P}(p)+\Gamma_{\rm q}^{\rm S}(p)\;,$$ with the scalar and pseudoscalar contributions, $\widetilde{N}^{\rm S}=-2$ and $\widetilde{N}^{\rm P}=2$, respectively:

$$\label{QuarkSpectralWidthOneLoop}
  \Gamma^{\rm S/P}_{\rm q}(p)=\frac{M{g_{\rm Mqq}}^2 \widetilde{N}^{\rm S/P}}{4\pi |{\mbox{\boldmath $p$}}|}\int_{E_{\rm min}}^{E_{\rm max}} {\text{d}}E_f \left[n_{\rm B}(E_b)+n_{\rm F}^-(E_f)\right],$$

involving the Fermi distribution for antiquarks, $n_{\rm F}^{-}$ in Eq. , and the Bose distribution for mesons at zero (quark-)chemical potential: $$\label{BoseDistribution}
  n_{\rm B}(E)=\frac{1}{{{\rm e}}^{\beta E}-1}\;.$$ Energy-momentum conservation implies $E_b=E_f+p_0$ in Eq.  and leads to the restricted range of integration: $$\begin{aligned}
\label{QuarkSpectralWidthEminEmax}
  E_{{\rm min},{\rm max}}(p)&=\frac{1}{2M^2}\left[\left(m_{\rm M}^2-2M^2\right)\sqrt{M^2+p^2}\right. \\
  &\hspace{2cm}\left.\mp\, p\,m_{\rm M}\sqrt{m_{\rm M}^2-4M^2}\right].
\end{aligned}$$ Here $m_{\rm M}$ denotes the $T$- and $\mu$-dependent (pseuodoscalar or scalar) meson mass and $M(T,\mu)$ is the dynamical quark mass in the thermal medium, c.f. Eq. . The quark-meson vertex is of Yukawa type with a (thermal) coupling constant $g_{\rm Mqq}(T,\mu)$ arising from the renormalized meson propagator $D_{\rm M}(p_0,{\mbox{\boldmath $p$}})$ at $(p_0,{\mbox{\boldmath $p$}})=(m_{\rm M},{\mbox{\boldmath $0$}})$: $$g_{\rm Mqq}=\left({\rm Res}\, D_{\rm M}\right)^{-1/2}\;.$$ The details of the derivation of Eqs.  and are presented in the Appendix \[Sec:App\]. The spectral width is a function of the energy $p_0$ and the momentum $p=|{\mbox{\boldmath $p$}}|$ of the quark propagating in the (isotropic) thermal medium. In the on-shell case, $p_0=\sqrt{p^2+M^2}$, the width $\Gamma$ becomes a function of $p$ only.

Taking the imaginary part of the mesonic insertion means cutting the loop diagram and forcing it on-shell. Two different dissipative mechanisms contribute to the width, as sketched in Fig. \[Fig8\]. The term involving the Bose distribution describes Landau damping, a process also known from high-$T$ QCD calculations [@Boyanovsky1998], and recently applied to the calculation of the shear viscosity from kinetic theory [@Ghosh2013]. In the NJL model Landau damping arises from the scattering of quarks on the mesonic collective modes in the thermal medium. This scattering process dissipates energy from the quark sector and contributes to the shear viscosity. The second mechanism is a recombination process: a collective mesonic mode is created by quark-antiquark rescattering. This is described by the term involving the Fermi distribution of thermal antiquarks in Eq. .

Results for the calculated on-shell spectral width including all mesons, $\Gamma(p)$ in Eq. , are shown in Fig. \[Fig9\]b. Due to the “antiscreening” caused by the sigma mode, the total $\Gamma(p)$ is rougly twice the single pion contribution, c.f. Fig. \[Fig9\]a. For small quark momenta the spectral width can become quite large. At $p\approx 200\;{\rm MeV}$ it is of the order of $100\;{\rm MeV}$, roughly as large as the dynamical (constituent) quark mass. We recall from Fig. \[Fig5\] that the shear viscosity $\eta[\Gamma]$ can be treated perturbatively only if $\Gamma\ll m_\pi$. Hence for the one-loop spectral width calculated in the NJL model, the perturbative regime is not reached, given the large values of $\Gamma(p)$. With rising temperature the spectral width decreases leading to an increasing shear viscosity $\eta(T)$, as explored parametrically in Fig. \[Fig1\]. One important reason for this behavior is the range of integration, $E_f(p)\in[E_{\rm 
min}(p),E_{\rm max}(p)]$, which has its support for $m_\pi>2M$, the kinematic threshold condition for a pion to decay on-shell into two constituent quarks. As the temperature increases this range of integration is shifted to higher values of $E_f$ and $p$, but these are exponentially suppressed by the Bose and Fermi distributions in Eq. . For more details, see also Fig. \[Fig11\] in the Appendix.

![Shear viscosity per entropy density at zero chemical potential from the two-flavor NJL model at one-loop level in a large-${N_{\rm c}}$ expansion (solid line). Shown for comparison are results using the kinetic approach (open circles and triangles), those of a related calculation using the Kubo formalism (solid triangles) and pure-gauge lattice data (solid squares). The dashed horizontal line is the AdS/CFT benchmark $1/4\pi$.[]{data-label="Fig10"}](Fig10.pdf){width="49.00000%"}

Finally the resulting shear viscosity $\eta(T,\mu = 0)$ (its ratio $\eta/s$) is shown in Fig. \[Fig10\]. Only temperatures large enough to satisfy $m_\pi>2M$ give rise to a finite shear viscosity. Therefore the one-loop results are restricted to $T\gtrsim 210\;{\rm MeV}$. We compare our results to those from a kinetic approach using the two-favor NJL model [@SasakiRedlich2010] (open circles) and [@Ghosh2013] (open triangles). Our values based on the Kubo formalism are larger but feature the same order of magnitude and the same qualitative behavior. These findings are consistent with the general observation that the kinetic approach seems to underestimate the shear viscosity as pointed out in Refs. [@PlumariEtAl2012; @PlumariEtAl2013]. Their calculations have also been performed using the Kubo formalism (c.f. Fig. 8 in [@PlumariEtAl2012]) and lie slightly above our result for $\eta/s$ as shown in Fig. \[Fig10\] (filled triangles). A flattening of $\eta/s$ is observed at higher temperatures around $T\approx 300\;{\rm MeV}$, indicating that the shear viscosity scales roughly as $T^3$ in this range. The same behavior has been found in our preceding parameter study, c.f. Fig. \[Fig6\]. All results for $\eta/s$, both from kinetic and Kubo approaches, exceed those from pure-gauge lattice QCD [@Nakamura05; @Meyer07] that are close to the AdS/CFT bound as shown by the solid squares in Fig. \[Fig10\].

Summary and Conclusions
=======================

In the present work we have used the Kubo formalism to derive a general functional $\eta[\Gamma(p)]$ for a class of fermionic theories, in particular the two-flavor NJL model with scalar and pseudoscalar interactions. At leading order in $1/{N_{\rm c}}$ the retarded correlation function $\Pi^{\rm R}(\omega)$ reduces to a single generic diagram at $\mathcal{O}({N_{\rm c}}^1)$. We have found that it is not necessary to work in the chiral limit to obtain this result derived previously in [@Fukutome2006; @Fukutome2008Nucl; @Fukutome2008Prog] using stronger assumptions.

The detailed study of the functional $\eta[\Gamma(p)]$ leads to a convergence criterion, Eq. , to be fulfilled by the spectral width $\Gamma(p)$ in order for the shear viscosity $\eta$ to be finite without regularization. Four different schematic parameterizations have been chosen for $\Gamma(p)$ and a numerical approximation scheme suitable for arbitrary, momentum-dependent spectral widths has been introduced. The results for $\eta$ show a strong cutoff dependence: restricting the momentum integration to the typical NJL cutoff, $p<\Lambda=650\;{\rm MeV}$, changes the shear viscosity drastically as shown in Fig. \[Fig3\]. Such a cutoff places $\eta$ at physically meaningful values by excluding up to $90\%$ of the numerically accessible high-momentum region, as demonstrated in Fig. \[Fig4\].

An exploration has been performed determining the range of $\Gamma(p)$ for which a perturbative treatment of the (generally non-perturbative) shear viscosity is adequate. The functional $\eta[\Gamma(p)]$ as given in Eq.  is valid in particular for a large spectral width as it is expected in the NJL model. For a small spectral width, $\Gamma\ll m_\pi$, a Laurent series expansion of $\eta$ with restriction to the leading term $\eta \sim 1/\Gamma$ is justified. In this limit, subleading effects from resummations of ladder diagrams are expected to become non-negligible as known from studies in bosonic field theories [@JeonSkeleton1995; @HidakaKunihiro2011]. On the other hand, for the NJL model with its comparatively large spectral width, it is consistent to omit ladder diagram resummations while taking all orders of the Laurent expansion in the width $\Gamma(p)$ into account.

Including thermal quasiparticle masses (the temperature and density dependence of the dynamically generated constituent quark masses in NJL-like models) is of crucial importance. As experimental data indicate, the ratio $\eta/s$ is expected to increase for high temperatures $T > T_{\rm c}$ where $T_{\rm c}\approx 180-200\;{\rm MeV}$ is the typical temperature interval of the chiral crossover. A constant quark mass would instead lead to a continuously decreasing ratio $\eta/s$ with increasing temperature, as seen in Fig. \[Fig7\]. In comparison with other approaches to the shear viscosity using NJL-type models, such as the Boltzmann equation in relaxation time approximation, one can find different results: an increasing ratio $\eta/s(T)$ for two flavors and restricting to scalar and pseudoscalar interactions [@SasakiRedlich2010], but also a decreasing behavior [@MartyFrankfurt2013; @HidakaKunihiro2011NJL].

From AdS/CFT correspondence, the scaling of the ratio $\eta/s$ is expected to be of $\mathcal{O}({N_{\rm c}}^0)$ in the large-${N_{\rm c}}$ limit. With the entropy density of quarks and the functional $\eta[\Gamma(p)]\sim{N_{\rm c}}$ in Eq. , we find indeed $\eta/s\sim \mathcal{O}(1)$ in the NJL model. Its applicability is, however, limited to the temperature range $T_{\rm c} \lesssim T < \Lambda$. At much higher temperatures gluonic degrees of freedom are dominant and contribute to the entropy density as $s_{\rm gluon}\sim{N_{\rm c}}^2$ due to their adjoint representation instead of the fundamental representation for fermions. At the same time the spectral width is expected to became small ($\Gamma\to 0$) and scales as $1/N_c$, leading to “pinch poles” and a scaling of the shear viscosity as $\eta \sim N_c^2$, as seen from Eq. . Consequently the ratio still scales as $\eta/s\sim\mathcal{O}(1)$ at large $T$ in the limit ${N_{\rm c}}\to\infty$.

From the systematic parameter study the conclusion can be drawn that the momentum dependence of the spectral width $\Gamma(p)$ determines primarily the overall scale of the ratio $\eta/s$. Its actual behavior as a function of temperature is largely governed by the thermal properties of the quark propagator. A consistent one-loop NJL model calculation of the quark self-energy in the thermal medium has been performed at next-to-leading order in the large-${N_{\rm c}}$ expansion. At this order the gap equation for the thermal constituent quark mass is corrected by the mesonic insertions . Two generic dissipative contributions to the shear viscosity emerge: Landau damping and a quark-antiquark recombination process, resulting in the spectral width $\Gamma(p)$ given in Eq.. The restricted range of the energy integration leads to a decreasing function $\Gamma(T)$ as the temperature increases, which implies an increasing shear viscosity for high temperatures as seen in Fig. \[Fig10\].

Additional contributions to $\Gamma$ from quark-quark scattering via exchange of mesonic quark-antiquark modes still have to be included. Such processes appear at order ${N_{\rm c}}^{-1}$ in the large-${N_{\rm c}}$ expansion and are calculated in ongoing work. Results will be reported elsewhere.

Acknowledgments {#acknowledgments .unnumbered}
===============

This work is partially supported by BMBF and by the DFG Cluster of Excellence “Origin and Structure of the Universe”. R. L. thanks the ECT\* for kind hospitality. Valuable discussions with Y. Hidaka and T. Hatsuda are acknowledged. Thanks go to N. Kaiser for useful advices and reading of the manuscript. R.L. has been supported by the TUM Graduate School (TUM-GS), by the RIKEN IPA and the RIKEN iTHES projects.

Derivation of the NJL spectral width {#Sec:App}
====================================

In this appendix we focus on the analytical treatment of the mesonic insertion $\Sigma_\beta^{\rm S/P}$ in Eq.  and the extraction of the spectral width $\Gamma(p)$ as given in Eqs.  and . The Matsubara frequencies for a thermal quark are $\nu_n=(2n+1)\pi T-{{\rm i}}\mu$, whereas they read for an antiquark: $\bar{\nu}_n=(2n+1)\pi T+{{\rm i}}\mu=\nu_n^*$. The diagram for the mesonic insertion represents both pions (pseudoscalar channel) and the sigma meson (scalar channel) contributing the quark spectral width: $$\begin{aligned}
  \Sigma^{\rm S/P}_\beta({\mbox{\boldmath $p$}},\nu_n)\; &=-4M g_{\rm Mqq}^2 \widetilde{N}^{\rm S/P} \\
  & \hspace{-1.5cm}\times{T\sum_{m\in\mathds{Z}}\int\frac{{\text{d}}^3 q}{(2\pi)^3}}\,\frac{1}{\nu_m^2+E_f^2}\,\,\frac{1}{(\nu_m-\nu_n)^2+E_b^2}\;,
\end{aligned}$$

with $E_f^2={\mbox{\boldmath $q$}}^2+M^2$ and $E_b^2=({\mbox{\boldmath $q$}}-{\mbox{\boldmath $p$}})^2+m_{\rm M}^2$. We have introduced $\widetilde{N}^{\rm P}=2$ and $\widetilde{N}^{\rm S}=-{N_{\rm f}}=-2$ reflecting the opposite parities of pion and sigma boson, leading to screening and antiscreening of the quark mass, respectively. Carrying out the Matsubara sum introduces combinations of Fermi and Bose distribution functions, Eq.  and , respectively: $$\begin{aligned}
  Z_1 &=1+n_{\rm B}(E_b)-\frac 12\left(n_{\rm F}^+(E_f)+n_{\rm F}^-(E_f)\right), \\
  Z_2 &= n_{\rm B}(E_b)+\frac 12\left(n_{\rm F}^+(E_f)+n_{\rm F}^-(E_f)\right), \\ 
  Z_3 &=n_{\rm F}^+(E_f)-n_{\rm F}^-(E_f)\;.
\end{aligned}$$ The thermal self-energy of the quark is given by:

$$\label{WidthFirstItThermalSelfEnergy}
  \Sigma^{\rm S/P}_\beta\!=\!-4Mg_{\rm Mqq}^2\widetilde{N}^{\rm S/P}\!\!\int\frac{{\text{d}}^3q}{(2\pi)^3}\left[\frac{1}{2E_bE_f}\left[\frac{(E_f+E_b)Z_1}{(E_f+E_b)^2+\nu_n^2}+\frac{(E_f-E_b)Z_2}{(E_f-E_b)^2+\nu_n^2}\right]+\frac{{{\rm i}}\nu_n Z_3}{[(E_f+E_b)^2+\nu_n^2][(E_f-E_b)^2+\nu_n^2]} \right].$$

The analytical continuation of the thermal self-energy to the retarded self-energy in Minkowskian spacetime, $$\Sigma^{\rm S/P}_{\rm R}(p_0,{\mbox{\boldmath $p$}}):=\Sigma^{\rm S/P}_\beta({\mbox{\boldmath $p$}},-{{\rm i}}p_0+{\varepsilon})\;,$$ leads to the spectral width $$\label{IterationOneSpectralWidthPart1}
  \Gamma^{\rm S/P}_{\rm q}(p_0,{\mbox{\boldmath $p$}}) :=-\lim_{{\varepsilon}\to 0}{{\rm Im}\,}\Sigma^{\rm S/P}_{\rm R}(p_0,{\mbox{\boldmath $p$}})\;.$$ The non-vanishing spectral width is induced by the pole structure of the propagators: $$\begin{aligned}
\label{DeltaFctArrImPart}
  \lim_{{\varepsilon}\to 0}{{\rm Im}\,}\left.\frac{Z}{x^2+\nu_n^2}\right|_{\nu_n\mapsto -{{\rm i}}p_0+{\varepsilon}} &= Z\pi\delta(x^2-(p_0)^2) \\ & \hspace{-1.5cm}=\frac{\pi Z}{2p_0}\left[\delta(x-p_0)+\delta(x+p_0)\right].
\end{aligned}$$ This means for the $Z_1$ term: $E_f+E_b\pm p_0=0$, where only the minus sign can be realized. For the $Z_2$ term, $E_f-E_b\pm p_0=0$, both signs can be realized for the time being. We will see that only the plus sign contributes to the (on-shell) spectral width, so there is just one contribution from $Z_2$. The $Z_3$ term is considered later.

From the identity in Eq.  the following structure is found for the $Z_1$ and $Z_2$ terms: $$\label{CalcZ1Z2PhaseSpace}
\begin{aligned}
  &\int\frac{{\text{d}}^3q}{(2\pi)^3}\frac{\pi Z}{2p_0}\frac{1}{2E_bE_f}\,\delta(E_b-(*))= \\
  &=\int\frac{{\text{d}}^3q}{(2\pi)^3}\frac{\pi Z}{2p_0E_f}\,\delta(E_b^2-(*)^2)= \\
\end{aligned}$$ $$\begin{aligned}
  &=2\pi\int_{-1}^1{\text{d}}\xi\int_0^\infty \frac{{\text{d}}q\, q^2}{(2\pi)^3}\frac{\pi Z}{2p_0E_f }\delta(E_b^2(\xi)-(*)^2)=\\
  &=2\pi\int_{M}^\infty\frac{{\text{d}}E_f}{(2\pi)^3}\frac{\pi Z}{4p_0p}\,\Theta(1-\xi^2)\;,
\end{aligned}$$ with $\xi =\cos\theta$. In order to carry out this integral over the delta function we have used $$E_b^2 =m_{\rm M}^2+({\mbox{\boldmath $p$}}-{\mbox{\boldmath $q$}})^2=m_{\rm M}^2+p^2+q^2-2pq\xi\;,$$ from which follows $$\left|\frac{\partial E_b^2}{\partial\xi}\right|=2pq\;.$$ In addition the integral over momentum has been translated into an energy integral using $q\,{\text{d}}q=E_f\,{\text{d}}E_f$. The ill-conditioned $\Theta$ term can be removed by the following consideration: The condition $|\xi|\leq 1$ is equivalent to $$\begin{aligned}
\label{DefCapitalFfunctionPhaseSpace}
  -1\leq & \frac{E_b^2-m_{\rm M}^2-p^2-q^2}{2pq}\leq 1 \\
  & \Leftrightarrow F(E_f,p)\geq 0\;,
\end{aligned}$$ where we have defined $$F(E_f,p):=4p^2(E_f^2-M^2)-\left[E_b^2-m_{\rm M}^2-p^2+M^2-E_f^2\right]^2.$$

![Range of integration for $E_f\in[E_{\rm min},E_{\rm max}]$ as function of absolute momentum $p=|{\mbox{\boldmath $p$}}|$ for different temperatures. The solid line displays the orbit of minimal values.[]{data-label="Fig11"}](Fig11.pdf){width="49.00000%"}

At given momentum $p$ the roots of $F(\,\cdot\, ,p)$ read, for the plus-sign case of both the $Z_1$ and $Z_2$ terms, where $E_b^2=(E_f+p_0)^2$: $$\begin{aligned}
  E_{{\rm max},{\rm min}} &=\frac{1}{2M^2}\left[\left(m_{\rm M}^2-2M^2\right)\sqrt{M^2+p^2}\right. \\ 
   & \hspace{1.5cm} \left.\pm pm_{\rm M}\sqrt{m_{\rm M}^2-4M^2}\right].
\end{aligned}$$ From this we find $$\label{EmaxEminDifferenceZeroValue}
\begin{aligned}
  & E_{{\rm max}}-E_{{\rm min}}=\frac{p\,m_{\rm M}}{M^2}\sqrt{m_{\rm M}^2-4M^2}\sim p\;, \\
  & E_{{\rm max},{\rm min}}(p=0)=\frac{m_{\rm M}^2}{2M}-M>M\;.
\end{aligned}$$ The momentum-dependent phase space for $E_f$ is shown in Fig. \[Fig11\] for different temperatures. Due to the larger meson mass at large $T$, the curves $E_{\rm min,max}$ are shifted to higher energies and momenta when increasing the temperature. There is an exact linear dependence of the integration range, $E_{\rm max}-E_{\rm min}$, on the incoming momentum $p$. At zero momentum the phase space collapses to one single point. Under the condition $m_{\rm M}>2M$ the phase space is always non-empty and compact with $M$ as minimal value of $E_{\rm min}$: $\emptyset\neq [E_{\rm min},E_{\rm max}]\subseteq [M,\infty)$. We also have derived the following substitution rule $$\int_{M}^\infty\frac{{\text{d}}E_f}{(2\pi)^3} \big(\cdot\big) \Theta(1-\xi^2) = \int_{E_{\rm min}}^{E_{\rm max}} \frac{{\text{d}}E_f}{(2\pi)^3} \big(\cdot\big)\;,$$ which leads finally to a well-conditioned one-dimensional numerical integral.

Inserting the minus-sign case into Eq., $E_b^2=(E_f-p_0)^2$, the phase space vanishes for any incoming momentum, since both the minimal and maximal energy of the loop fermion need to be negative: $$E_{\rm min}'=-E_{\rm max}\, ,\;\;\; E_{\rm max}'=-E_{\rm min}\;.$$ We conclude that only the case $E_b=E_f+p_0$ leads to an on-shell contribution to the mesonic insertion. With this in mind the third term, involving $Z_3$, in Eq. becomes: $$\begin{aligned}
  \lim_{{\varepsilon}\to 0} &\;{{\rm Im}\,}\left.\frac{{{\rm i}}\nu_n Z_3}{[(E_f+E_b)^2+\nu_n^2][(E_f-E_b)^2+\nu_n^2]}\right|_{{{\rm i}}\nu_n\mapsto p_0+{{\rm i}}{\varepsilon}}= \\
   &=\frac{p_0\pi Z_3}{2p_0}\,\delta\big(\underbrace{[(E_f+E_b)^2-p_0^2]}_{=4E_fE_b}[E_b-E_f-p_0]\big)=\\
  &=\frac{\pi Z_3}{4E_f}\,\delta\left(E_b^2-(E_f+p_0)^2\right).
\end{aligned}$$ Note that due to the ${{\rm i}}\nu_n$ factor in the first line, the $p_0$ terms cancel in the final result. As done in the calculation the momentum integral can be performed: $$\begin{aligned}
  \int\frac{{\text{d}}^3q}{(2\pi)^3} \frac{\pi Z_3}{4E_f}\, &\delta(E_b^2-(E_f+p_0)^2)= \\
  &=2\pi\int_{M}^\infty\frac{{\text{d}}E_f}{(2\pi)^3}\frac{\pi Z_3}{8p}\,\Theta(1-\xi^2)\;.
\end{aligned}$$ Combining all contributions, we find for the spectral width of an incoming quark: $$\begin{aligned}
  &\Gamma^{\rm S/P}_{\rm q}(p) = \\
  &=-8\pi Mg_{\rm Mqq}^2\widetilde{N}^{\rm S/P}\!\int_{E_{\rm min}}^{E_{\rm max}}\frac{{\text{d}}E_f}{(2\pi)^3}\frac{\pi}{4|{\mbox{\boldmath $p$}}|}\left(\frac{E_f-E_b}{p_0}Z_2+\frac{Z_3}{2}\right)=\\
  &=\frac{Mg_{\rm Mqq}^2\widetilde{N}^{\rm S/P}}{4\pi|{\mbox{\boldmath $p$}}|}\int_{E_{\rm min}}^{E_{\rm max}}{\text{d}}E_f\left[n_{\rm B}(E_b)+n_{\rm F}^-(E_f)\right],
\end{aligned}$$ using $E_b=E_f+p_0$. The effective spectral width in Eq.  is determined by summing the pseudoscalar (pion) and scalar (sigma boson) channels weighted by their respective multiplicities, $3=N_{\rm f}^2-1$ and $1$.

[9]{} I. Arsene et al. (BRAHMS Collaboration), Nucl. Phys. A [**757**]{}, 1 (2005). K. Adcox et al. (PHENIX Collaboration), Nucl. Phys. A [**757**]{}, 184 (2005). B. B. Back et al. (PHOBOS Collaboration), Nucl. Phys. A [**757**]{}, 28 (2005). J. Adams et al. (STAR Collaboration), Nucl. Phys. A [**757**]{}, 102 (2005). K. Aamodt et al. (ALICE Collaboration), Phys. Rev. Lett. [**108**]{}, 252302 (2010). K. Aamodt et al. (ALICE Collaboration), Phys. Lett. B [**696**]{}, 30 (2011). D. Caffarri et al. (ALICE Collaboration), Nucl. Phys. A [**904–905**]{}, 643c (2013). D. Derendarz et al. (ATLAS Collaboration), Nucl. Phys. A [**904–905**]{}, 373c (2013). V. Roy, B. Mohanty, and A.K. Chaudhuri, J. Phys. G [**40**]{}, 065103 (2013). K. Dusling, and D. Teaney, Phys. Rev. C [**77**]{}, 034905 (2008). M. Luzum, C. Gombeaud, and J.Y.  Ollitrault, Phys. Rev. C [**81**]{}, 054910 (2010). H. Niemi, G.S. Denicol, P. Huovinen, E. Molnár, and D.H. Rischke, Phys. Rev. Lett. [**106**]{}, 212302 (2011). H. Liu, D. Hou, and J. Li, Eur. Phys. J. C [**45**]{}, 459 (2006). Ch. Sasaki and K. Redlich, Nucl. Phys. A [**832**]{}, 62 (2010). A.S. Khvorostukhin, V.D. Toneev, and D.N. Voskresensky, Phys. Rev. C [**84**]{}, 9 (2011). 

W.M. Alberico, S. Chiacchiera, H. Hansen, A. Molinari, and M. Nardi, Eur. Phys. J. A [**38**]{}, 97 (2008). Y. Hidaka and T. Kunihiro, Phys. Rev. D [**83**]{}, 076004 (2011). Y. Hidaka and T. Kunihiro, J. Phys.: Conf. Ser. [**270**]{}, 012050 (2011). S.-I. Nam and C.-W. Kao, Phys. Rev. D [**87**]{}, 114003 (2013). S. Plumari et al., Phys. Rev. C [**86**]{} 054902, (2012). 

S. Plumari et al., J. Phys.: Conf. Ser. [**420**]{}, 012029 (2013). S. Jeon, Phys. Rev. D [**52**]{}, 3591 (1995). S. Jeon and L.G. Yaffe, Phys. Rev. D [**53**]{}, 5799 (1996). Y. Nambu, G. Jona-Lasinio, Phys. Rev. [**122**]{}, 345 (1961). Y. Nambu, G. Jona-Lasinio, Phys. Rev. [**124**]{}, 246 (1961). S.P. Klevansky, Rev. Mod. Phys. [**64**]{}, 649 (1992). S. Klimt, M. Lutz, U. Vogl, and W. Weise, Nucl. Phys. A [**516**]{}, 429 (1990). S. Klimt, M. Lutz, U. Vogl, and W. Weise, Nucl. Phys. A [**516**]{}, 469 (1990). U. Vogl and W. Weise, Prog. Part. Nucl. Phys. [**27**]{}, 195 (1991). M. Buballa, Phys. Rept. [**407**]{}, 205 (2005). 

G. ’t Hooft, Nucl. Phys. B [**72**]{}, 461 (1974). E. Quack and S.P. Klevansky, Phys. Rev. C [**49**]{}(6) 3283 (1994). D. Müller, M. Buballa, and J. Wambach, Phys. Rev. D [**81**]{}, 1 (2010). M. Iwasaki, and H. Ohnishi, and T. Fukutome, arXiv:hep-ph/0606192 (2006). T. Fukutome and M. Iwasaki, Prog. Theor. Phys. [**119**]{}, 991 (2008). M. Iwasaki, H. Ohnishi, and T. Fukutome, arXiv:hep-ph/0703271 (2007). R. Kubo, J. Phys. Soc. Jpn. [**12**]{}, 570 (1957). R. Lang, N. Kaiser, and W. Weise, Eur. Phys. J. A [**48**]{}, 109 (2012). U. Heinz, Ch. Shen, and H. Song, AIP. Conf. Proc. [**1441**]{}, 766 (2012). J.I. Kapusta and Ch. Gale, *Finite-Temperature Field Theory* (Cambridge, 2006). R.A. Lacey et al., Phys. Rev. Lett. [**98**]{}, 092301 (2007). A. Nakamura and S. Sakai, Phys. Rev. Lett. [**94**]{} 072305 (2005). 

H.B. Meyer., Phys. Rev. D [**76**]{}, 101701(R) (2007). J. Maldacena, Int. J. Theor. Phys. [**38**]{}, 1113 (1999). P.K. Kovtun, D.T. Son, and A.O. Starinets, Phys. Rev. Lett. [**94**]{}, 111601 (2005). T.D. Cohen, Phys. Rev. Lett. [**99**]{}, 021602 (2007). A. Rebhan and D. Steineder, Phys. Rev. Lett. [**108**]{}, 021601 (2012). K.A. Mamo, JHEP [**1210**]{}, 070 (2012). 

T. Hatsuda and T.Kunihiro, Prog. Theor. Phys. [**91**]{}, 284 (1987). D. Boyanovsky et al., Phys. Rev. D [**58**]{} 125009 (1998). S. Ghosh et al., Phys. Rev. C [**88**]{}(6) 068201 (2013). R. Marty, E. Bratkovskaya, W. Cassing, J Aichelin, and H. Berrehrah, Phys. Rev. C [**88**]{}(4) 045204 (2013).

[^1]: Due to ${{\rm e}}^{-\beta H}X(t){{\rm e}}^{\beta H}=X(t+{{\rm i}}\beta)$, which implies $\langle X(t)Y(t'+{{\rm i}}\beta)\rangle=\langle Y(t')X(t)\rangle$, it follows that this correlator is symmetric in its arguments: $(X,Y)=(Y,X)$.

[^2]: The notation $f\in{\rm o}(g)$ is used to express accurately that “$f$ is growing less fast than $g$”, meaning that $f(x)/g(x)\to 0$ for $x\to\infty$. More intuitively, this also means that “$g$ grows much faster than $f$”.

[^3]: Note that in this definition as well as in previous expressions the spectral width $\Gamma$ corresponds to half the FWHM (full width at half maximum) of the resulting quark spectral function, $\pi\rho(p)=-{{\rm Im}\,}\,G^{\rm R}(p)$, where the retarded quark propagator has been defined in Eq. .
