Transport coefficients for multi-component gas of hadrons using Chapman Enskog method

TL;DR

This paper calculates the transport coefficients of a multi-component hadronic gas at different chemical potentials using the Chapman-Enskog method within an S-matrix based model, providing temperature-dependent results and comparisons with prior studies.

Contribution

It introduces a novel calculation of transport coefficients for hadronic gases using an S-matrix based approach combined with Chapman-Enskog method, including temperature dependence and chemical potential effects.

Findings

Transport coefficients depend on temperature and chemical potential.

Results are consistent with previous literature.

Provides detailed shear viscosity, bulk viscosity, heat conductivity, and diffusion coefficients.

Abstract

The transport coefficients of a multi-component hadronic gas at zero and non-zero baryon chemical potential are calculated using the Chapman-Enskog method. The calculations are done within the framework of an $S$-matrix based interacting hadron resonance gas model. In this model, the phase-shifts and cross-sections are calculated using $K$-matrix formalism and where required, by parameterizing the experimental phase-shifts. Using the energy dependence of cross-section, we find the temperature dependence of various transport coefficients such as shear viscosity, bulk viscosity, heat conductivity and diffusion coefficient. We finally compare our results regarding various transport coefficients with previous results in the literature.

Equations

$$S_{ab\to cd}=⟨cd|S|ab⟩,$$

$$f(\sqrt{s},\theta )=\frac{1}{q_{ab}}\sum _l(2l+1)T^l{P}_l(\cos \theta ),$$

$$\sigma (\sqrt{s},\theta )={|f(\sqrt{s},\theta )|}^2.$$

$$S=I+2iT,$$

$${(T^{-1}+iI)}^{†}=T^{-1}+iI.$$

$$K^{-1}=T^{-1}+iI.$$

$$\frac{\text{Im}T}{\text{Re}T}=K.$$

$$K_{ab\to cd}^{I,l}=\sum _R\frac{g_{R\to ab}(\sqrt{s})g_{R\to cd}(\sqrt{s})}{m_R^2-s},$$

$$g_{R\to ab}^2(\sqrt{s})=m_R{\Gamma }_{R\to ab}(\sqrt{s}),$$

$${\Gamma }_{R\to ab}(\sqrt{s})={\Gamma }_{R\to ab}^0\frac{m_R}{\sqrt{s}}\frac{q_{ab}}{q_{ab0}}{\left(B^l(q_{ab},q_{ab0})\right)}^2,$$

$$q_{ab}(\sqrt{s})=\frac{1}{2\sqrt{s}}\sqrt{\left(s-{(m_a+m_b)}^2\right)\left(s-{(m_a-m_b)}^2\right)},$$

$$K=\tan {\delta }_l^I.$$

$$\ln Z=\ln Z_0+\ln Z_{int}.$$

$$\ln Z_0=\sum _h\frac{V{g}_h{m}_h^{2}T}{2{\pi }^2}\sum _{j=1}^{\infty }(\pm 1{)}^{j-1}(z_h^j/j^2)K_2(j{m}_h\beta ),$$

$$\ln Z_{int}=\sum _{i_1,i_2}{z}_1^{i_1}{z}_2^{i_2}b(i_1,i_2),$$

$$b(i_1,i_2)=$$

$${\left[A\left\{S^{-1}\frac{\partial S}{\partial \epsilon }-\frac{\partial S^{-1}}{\partial \epsilon }S\right\}\right]}_c.$$

$$S(\epsilon )=\sum _{l.I}(2l+1)(2I+1)\exp (2i{\delta }_l^I),$$

$$b_2=\frac{1}{2{\pi }^3\beta }{\int }_M^{\infty }d\epsilon {\epsilon }^2{K}_2(\beta \epsilon )\sum _{l,I}{}^{{}^{\prime }}{g}_h\frac{\partial {\delta }_l^I(\epsilon )}{\partial \epsilon }.$$

$$p^{\mu }{\partial }_{\mu }{f}_1=C[f,f].$$

$$C[f,f]=\frac{1}{2}\int \frac{d^3{p}_2}{p_2^0}\frac{d^3{p}_3}{p_3^0}\frac{d^3{p}_4}{p_4^0}[f_3{f}_4(1+\theta f_1)(1+\theta f_2)$$

$$-f_1{f}_2(1+\theta f_3)(1+\theta f_4)]W(p_3,p_4|p_1,p_2),$$

$${\partial }^{\mu }\to U^{\mu }D+{\nabla }^{\mu },$$

$$f=f^{(0)}+ϵf^{(1)}.$$

$$p^{\mu }{U}_{\mu }D{f}_1^{(0)}+p^{\mu }{\nabla }_{\mu }{f}_1^{(0)}=-f_1^{(0)}\mathcal{L}[\varphi ],$$

$$\mathcal{L}[\varphi ]$$

$$W(p_3,p_4|p_1,p_2).$$

$$f_i^{(0)}=\exp \left(\frac{{\mu }_i(x)-p_i^{\nu }{U}_{\nu }(x)}{T(x)}\right)$$

$$n=\int \frac{d^{3}p}{{(2\pi )}^3{p}^0}(p^{\mu }{U}_{\mu })f^{(0)},$$

$$en=\int \frac{d^{3}p}{{(2\pi )}^3{p}^0}{(p^{\mu }{U}_{\mu })}^2{f}^{(0)}.$$

$$D{f}^{(0)}$$

$${\nabla }^{\alpha }{f}^{(0)}=\left[T{\nabla }^{\alpha }\left(\frac{\mu }{T}\right)+p^{\mu }{U}_{\mu }{\nabla }^{\alpha }\log T-p^{\mu }{\nabla }^{\alpha }{U}_{\mu }\right]\frac{f^{(0)}}{T},$$

Full text

I Introduction

One of the important discoveries from experiments at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) in search of the quark-gluon plasma (QGP) is the fact that the deconfined quark- gluon matter behaves as an almost-perfect fluid Adams et al. (2005); Gyulassy and McLerran (2005); Hirano and Gyulassy (2006); Luzum and Romatschke (2008); Song et al. (2011); Schenke et al. (2011); Jacak and Steinberg (2010); Roy and Chaudhuri (2011); Heinz and Snellings (2013). The property that quantifies the ’liquidness’ of a liquid are its transport coefficients, for e.g. the ratio of shear viscosity to its entropy density $\eta_{s}/s$ or the ratio of bulk viscosity to its entropy density $\eta_{v}/s$. Experimentally measured elliptic flow, through the azimuthal correlation of produced particles with respect to the reaction plane is a sensitive probe for obtaining the transport coefficients. For example, the magnitude of elliptic flow depends quite sensitively on the shear viscosity of the QGP fluid, which is estimated to be around $\eta_{s}/s\approx$ $0.08-0.20$ during the hydrodynamic evolution Drescher et al. (2007); Chaudhuri (2009); Roy et al. (2013). These small values of $\eta_{s}/s$ make the system (QGP), a near perfect fluid.

Theoretically, the value of the shear viscosity depends on the model under consideration. For example, at weak coupling the dimensionless ratio $\eta_{s}/s$ is proportional to the ratio of mean free path to the mean spacing between the particles and weaker coupling means a larger value of this ratio Arnold et al. (2000, 2003); York and Moore (2009). For example, for a weakly coupled gas of gluons the $\eta_{s}/s\sim{(\alpha_{s}^{4}\log(1/\alpha_{s}}))^{-1}$, where $\alpha_{s}$ is the QCD coupling constant. For $\alpha_{s}=0.1$, the value of $\eta_{s}/s\sim 4\times 10^{3}$. On the other hand, in strongly coupled system their is efficient momentum transfer and the ratio $\eta_{s}/s$ is significantly smaller Busza et al. (2018). For example, a lower bound of the ratio $\eta_{s}/s=1/4\pi$ for strongly coupled field theories using the anti-de Sitter/conformal field theory (AdS/CFT) correspondence has been conjectured in Ref. Policastro et al. (2001). In this work, we shall find that the ratio $\eta_{s}/s$ values varies from weak coupling to strong coupling regime as one goes from lower to higher temperatures.

There are two other important reasons for studying the temperature dependence of transport coefficients. First, experimentally it has been observed Lacey et al. (2007); Csernai et al. (2006) that $\eta/s$ shows a minimum near the liquid-gas phase transition for different substances, this might help in studying QCD phase diagram. Such a minimum has also been observed in model calculation. For example in Ref. Prakash et al. (1993), it was shown for massless pions the ratio $\eta_{s}/s$ diverges as temperature $T\rightarrow 0$, whereas in the quark gluon phase, at one loop order in $\alpha_{s}$, the $\eta_{s}/s$ is an increasing function of $T$. Second, it has been predicted that the ratio $\eta_{v}/s$ should show a maximum near the phase transition Kharzeev and Tuchin (2008); Chakraborty and Kapusta (2011); Dobado and Torres-Rincon (2012). For example in Ref. Chen and Wang (2009) the ratio $\eta_{v}/s$ for massless pions goes to zero in the $T\rightarrow 0$ limit and also to zero in the quark gluon phase for asymptotically high $T$ Arnold et al. (2006).

In this work, we calculate various transport coefficients of a hadronic gas consisting of baryon and meson octets namely $\pi,\eta,K,N,\Lambda,\Sigma,\Xi$. The corresponding resonances which appear as interactions among these hadronic constituents are handled using the $K$- matrix formalism Wiranata et al. (2013); Dash et al. (2018). The formalism that we use for calculating the equilibrium thermodynamic quantities like entropy density, enthalpy density, number density etc. is through the $S$-matrix based Hadron Resonance Gas (HRG) model Dash et al. (2018, 2019). The cross-sections that are used in the calculation of transport coefficients are calculated in the $K$-matrix formalism for all hadrons except for the nucleons, where we directly use the experimental phase-shifts Workman et al. (2016). The calculations are done for a system with vanishing baryon chemical potential ($\mu_{B}$) as well as for finite $\mu_{B}=100$ MeV. The transport coefficients are obtained using Chapman-Enskog (CE) method developed in Refs. van Leeuwen et al. (1973, 1974); van Leeuwen (1975). In this method the solution of the transport equation, i.e. the distribution function to be determined, is first written as an infinite series of Laguerre polynomials. With the help of this expansion, the transport equation could be transformed into an infinite set of linear, algebraic equations. From this infinite set of equations, a finite number of equations are taken and solved to get an approximate solution for the distribution function. This solution is used to compute the transport coefficients.

Pioneering work on transport coefficients in the CE method has been done in Ref. von Oertzen (1992) for quark and gluon system and in Ref. Prakash et al. (1993), for various binary combinations in a system consisting of $\pi-K-N$ using experimental cross-sections. Similarly, in Ref. Wiranata et al. (2013) the calculations of $\eta/s$ for a multi-component system consisting of $\pi-K-N-\eta$ at vanishing $\mu_{B}$ in the $K$-matrix formalism, has been carried out but without including $NN$ interaction. In Ref. Moroz (2013), CE method was used for calculating $\eta_{s}/s$ and $\eta_{v}/s$ using UrQMD cross-sections. The current work improves upon all these previous work to incorporate a larger spectrum of interacting hadronic states (7 stable hadrons + 112 resonances) and also extends to finite chemical potential. By adding more stable hadrons into the mixture, one hopes that new channels of interaction (through resonance formation) could open up, which would relax the system to equilibrium quicker, than with fewer hadrons considered in earlier works. Also, the degeneracy of the system changes, which affects equilibrium quantities like entropy density and number density etc., and in turn also affects the dimensionless transport coefficient ratios. It is also interesting to compute transport coefficients at non-zero chemical potential, since finite baryon density, affects the concentration of various species interacting in the mixture and thus the overall weight coming from different channels, on the final value of transport coefficient. In regards to other formalism for e.g in Refs. Gavin (1985); Chakraborty and Kapusta (2011); Kadam et al. (2019), which uses relaxation time approximation (RTA), the present formalism is better in the sense that small angle scattering is taken care of naturally, where as RTA uses thermal averaged cross-sections. Similarly, compared to models like ideal hadron resonance gas, excluded volume approach Gorenstein et al. (2008); Denicol et al. (2013); Kadam and Mishra (2015); Ghosh et al. (2016); Dobado and Llanes-Estrada (2004) which uses constant values of cross-section, the present formalism utilizes the energy dependence of cross-sections to calculate the temperature dependence of transport coefficients. Calculations of shear viscosity has also been done using the Kubo formalism in transport models Demir and Bass (2009); Rose et al. (2018). Our results on transport coefficients are in reasonable agreement with that from the transport models in the temperature range of $T=80-110$ MeV.

The paper is organised as follows. In Sec. II, we describe the $K$-matrix formalism for calculating the scattering phase-shifts and cross-section. In Sec. III we have discussed the thermodynamics of interacting hadron gas using the formalism of Sec. II. Then in Sec. IV, we describe the CE method for calculating the transport coefficients for single, binary and multi-component system of hadrons at zero and finite $\mu_{B}$. Finally, in Sec. V, we summarize our findings.

II $K$-matrix Formalism

A theoretical way of calculating the attractive part of the phase shifts is to use the $K$-matrix formalism. In this section, we briefly discuss the $K$-matrix formalism.

For the process $ab\rightarrow cd$

[TABLE]

where $S$ is the scattering matrix operator. The scattering amplitude for the process can be expressed in terms of interaction matrix $T^{l}$ as,

[TABLE]

where $P_{l}(\cos\theta)$ are the Legendre polynomials for the angular momentum $l$ and $\theta$ is the scattering angle in the center of mass frame. The cross section for the process can be given in terms of terms of scattering amplitude,

[TABLE]

The $T$ matrix is related to the $S$-matrix by the following equation

[TABLE]

where $I$ is the unit matrix. Using the unitarity of $S$-matrix one can show

[TABLE]

Therefore, one can define a Hermitian K matrix through

[TABLE]

The $K$-matrix formalism preserves the unitarity of $S$-matrix and neatly handles multiple resonances Chung et al. (1995). In addition to that, widths of the resonances are handled naturally in the above formalism. For overlapping resonances the $K$-matrix gives good description of the phase shifts.

One can write real and imaginary part of the $T$-matrix in terms of $K$ matrix as

[TABLE]

In Ref. Wiranata et al. (2013) the $K$-matrix formalism was used to study the shear viscosity of an interacting gas of hadrons. Recently in Ref. Dash et al. (2018, 2019), $K$-matrix formalism is used to study the equation of state and susceptibilities pertaining to conserved charges.

For the process $ab\rightarrow R\rightarrow cd$, resonances appear as sum of poles in the $K$-matrix as $K_{ab\rightarrow cd}^{I,l}$

[TABLE]

where $a$, $b$ and $c$, $d$ are hadrons, $R$ is the resonance with mass $m_{R}$. The sum over $R$ is restricted to the addition of resonances which have the same spin $l$ and isospin $I$. The residue functions are given by

[TABLE]

where $\sqrt{s}$ is the energy in the center of mass frame and $\Gamma_{R\rightarrow ab}(\sqrt{s})$ is the energy dependent partial decay widths. For the process $R\rightarrow ab$, decay width can be written as Chung et al. (1995)

[TABLE]

where $\Gamma^{0}_{R\rightarrow ab}$ is the partial width of the pole of the resonance at half maximum for the channel ${R\rightarrow ab}$, $B^{l}(q_{ab},q_{ab0})$ are the Blatt-Weisskopf barrier factors Chung et al. (1995) which can be expressed in terms of daughter momentum $q_{ab}$ and resonance momentum $q_{ab0}$ for the orbital angular momentum $l$. The momentum $q_{ab}$ in the last expression is given as

[TABLE]

where $m_{a}$ and $m_{b}$ being the mass of hadrons $a$ and $b$ decaying from resonance $R$. In Eq. (10), $q_{ab0}=q_{ab}(m_{R})$ is the resonance momentum at $\sqrt{s}=m_{R}$.

Therefore, to compute the $K$-matrix one needs the relevant masses and widths of resonances. Further, in partial wave decomposition, the $K$-matrix can be written in terms of phase shift $\delta_{l}^{I}$ as  Chung et al. (1995),

[TABLE]

Since the $K$-matrix formalism is applicable only for attractive interaction, for the repulsive as well as $NN$ interaction we have used the experimental data of phase shift Workman et al. (2016).

III Thermodynamics of interacting hadrons

The most natural way to incorporate interaction among a gas of hadrons is to use relativistic virial expansion Dashen et al. (1969); Venugopalan and Prakash (1992) where the logarithm of the partition function can be separated into two parts, non-interacting ($Z_{0}$) and interacting ($Z_{int}$) parts i.e.,

[TABLE]

The non-interacting part of the partition function can be written as

[TABLE]

where $h$ denotes the index of stable hadron, $V$ is the volume of the system, $g_{h}$ is the degeneracy, $m_{h}$ is the mass of the hadron, $z_{h}=\exp(\mu_{h}\beta)$ is the fugacity, $\beta$ is the inverse of the temperature (T) and $K_{2}$ is the modified Bessel function of the second kind. For the conserved quantities like baryon number, strangeness and electric charge, $\mu_{h}$ can be written as $\mu_{h}=B_{h}\mu_{B}+S_{h}\mu_{S}+Q_{h}\mu_{Q}$. Here $B_{h},S_{h},Q_{h}$ are respectively the baryon number, strangeness and electric charge of the hadron and the $\mu^{,}s$ are the chemical potentials of the corresponding conserved charges. In Eq. (14), $+(-)$ sign refer to bosons (fermions) and $j=1$ term corresponds to the classical ideal gas.

The interacting part of Eq. (13) can be written as

[TABLE]

where $b(i_{1},i_{2})$ is the virial coefficient defined as,

[TABLE]

In the above equation, the labels $i_{1}$ and $i_{2}$ refer to channel of the $S$-matrix which has initial state containing $i_{1}+i_{2}$ particles. The symbol $A$ denotes the symmetrization (anti-symmetrization) operator for a system of bosons (fermions), the subscript $c$ refers to trace over all the linked diagrams.

The lowest virial coefficient i.e., the second virial coefficient, $b_{2}=b(i_{1},i_{2})/V$ as $V\rightarrow\infty$, corresponds to interaction between two hadrons ($i_{1}=i_{2}=1$). The higher order virial coefficients give interaction among many hadrons. In the present, work we will consider only up to the second virial coefficient.

The $S$-matrix can be expressed in terms of phase shifts ${\delta}_{l}^{I}$ as  Dashen et al. (1969)

[TABLE]

where $l$ and $I$ denote angular momentum and isospin, respectively. Integrating Eq. (III) over the total momentum we get

[TABLE]

The factor $g_{h}=(2I+1)(2l+1)$ is the degeneracy factor, $M$ is the invariant mass of the interacting pair at zero centre of mass momentum. The prime over the summation sign denotes that for given $l$ the sum over $I$ is restricted to values consistent with statistics.

The Eq. (18) shows that the contribution arising from interaction to thermodynamic variable depends on the derivative of the phase shift. The positive values of derivative of phase shifts (attractive interactions) give positive contributions to the thermodynamical variables and negative value of derivative of phase shifts (repulsive interactions) give negative contributions.

Once we know the partition function (Eq. 13) we can calculate various thermodynamic quantities like pressure, energy density, entropy density, number density etc. In Fig. 1, we show the scaled number density ($n/T^{3}$) as a function of temperature, using the interacting model described above, for mesons ($\pi,K,\eta$) and baryons ($N,\Lambda,\Sigma,\Xi$) at $\mu_{B}=0$ MeV and at $\mu_{B}=100$ MeV.

IV Transport coefficients

The relativistic Boltzmann equation, describing the space-time evolution of the phase space density $f=f(x,p)$, where $x$ is position and $p$ is momentum, is given by De Groot (1980),

[TABLE]

The collision term $C[f,f]$, in the Boltzmann approximation, is given by,

[TABLE]

where $p_{1},p_{2}$ are momenta of incomimg and $p_{3},p_{4}$ are momenta of outgoing particles respectively. $W(p_{3},p_{4}|p_{1},p_{2})$ is the transition rate in the collision process $p_{1}+p_{2}\leftrightarrow p_{3}+p_{4}$. The constant $\theta=\pm 1$ for bosons or fermions and [math] for classical Maxwellian particles. We shall employ the Chapman Enskog method as discussed in Refs. van Leeuwen et al. (1973, 1974); van Leeuwen (1975) to linearize and solve the kinetic equation Eq. (20). We split the derivative operator $\partial^{\mu}$ into a time-like and space-like part

[TABLE]

where $D=U^{\nu}\partial_{\nu}$ and $\nabla^{\mu}=\Delta^{\mu\nu}\partial_{\nu}$ and $\Delta^{\mu\nu}=g^{\mu\nu}-U^{\mu}U^{\nu}$ is the projection operator. Here, $U^{\mu}$ is the hydrodynamic four velocity, as discussed in Ref. van Leeuwen et al. (1973). Taking $\theta=0$, i.e., assuming the particles to be classical, we expand the distribution function $f$ into an equilibrium part $f^{(0)}$ and a deviation $\epsilon f^{(1)}$, i.e.,

[TABLE]

To order $\epsilon$, substituting Eq. (22) into the the transport equation Eq. (19) gives

[TABLE]

where $\mathcal{L}[\phi]$ is the linearized collision operator (found from Eq. (20), using Eq. (22) and invoking the principle of detailed balance given as, $f_{1}^{(0)}f_{2}^{(0)}=f_{3}^{(0)}f_{4}^{(0)}$).

Hence,

[TABLE]

The $\phi_{i}$ is the ratio $f^{(1)}_{i}/f^{(0)}_{1}$. The equilibrium distribution functions $f^{(0)}_{i}$ are assumed to be Maxwell Boltzmann type

[TABLE]

To identify the functions $\mu(x),U^{\mu}(x)$ and $T(x)$ with the usual definitions of chemical potential, hydrodynamic velocity and temperature of the system, we demand that the particle density $n$ and energy density $en$ be determined solely by the local equilibrium distribution function in Eq. (25) as,

[TABLE]

[TABLE]

The choice of distribution function given in Eq. (25) along with condition given in Eqs. (26) and (27) determines the set of independent variables $T,\mu,U^{\nu}$. The derivative of the distribution function $f^{(0)}$, then depends only on the above set of independent variables. Then one can express $Df^{(0)}$ as,

[TABLE]

and $\nabla^{\alpha}f^{(0)}$ as,

[TABLE]

expressed in terms of temperature $T$, density $n$, hydrodynamic four-velocity $U^{\mu}$ and the chemical potential $\mu$.

Multiplying Eq (23) with $\int d^{3}pp^{\mu}/p^{0}$ and contracting with $U^{\mu}$, gives De Groot (1980),

[TABLE]

where $\gamma=c_{P}/c_{V}$ is the ratio of heat capacities at constant pressure $c_{P}$ and constant volume $c_{V}$. Similarly, on multiplying Eq (23) with $\int d^{3}pp^{\mu}/p^{0}$ and contracting with projection operator $\Delta^{\mu\nu}$, gives the equation of motion

[TABLE]

where, $wn$ is the enthalpy density, $wn=en+P$ and $P$ is the pressure van Leeuwen et al. (1973). Also, the continuity equation, for e.g. given in Ref. De Groot (1980)

[TABLE]

can be used to express the time derivative of number density in terms of gradients of hydrodynamic velocity. Eqs. (30-32) are used in Eqs. (28) and (29), to express time derivative of $T,n$ and $U^{\mu}$ in terms of gradients of $U^{\mu}$ and $P$ respectively.

The expressions of $Df^{(0)}$ and $\nabla^{\mu}f^{(0)}$ given in Eqs. (28) and (29) can be substituted in the linearized transport equation Eq. (23). Thus, one can express the transport equation in terms of thermodynamics forces, whose components include scalar force, vectorial force and tensorial force respectively. The scalar force can be expressed as the divergence of hydrodynamic velocity

[TABLE]

the vectorial force, due to temperature gradient and pressure gradient is given as

[TABLE]

and tensorial forces (traceless indicated by “$\langle~{}\rangle$”, due to gradient of hydrodynamic velocity is given as

[TABLE]

In terms of these forces, the transport equation is then given as

[TABLE]

The quantity $Q$ is defined as

[TABLE]

where the relativistic version of Gibbs-Duhem relation De Groot (1980) $T\nabla^{\mu}(\mu/T)=-w(\nabla^{\mu}T/T-\nabla^{\mu}P/wn)$ was used for the derivation of Eq. (36).

An equation similar to Eq. (36) can also be derived for a two component mixture with components labeled by subscripts 1 and 2. Here, we indicate the few differences pertaining to extension of derivation of Eq. (36) for binary mixtures. Interested readers may refer van Leeuwen et al. (1974) for the complete derivation.

The analogous linearized transport equation for mixtures can be written as,

[TABLE]

An equation similar to the Eq. (38) also holds for component 2. The right hand side takes collisions of the form $1(2)+1(2)\rightarrow 1(2)+1(2)$ and $1+2\rightarrow 1+2$ into account. The linearized operator is given by

[TABLE]

The factor $(1-\delta_{1k}/2)$ takes into account the correct weighting for same or different species which interact in the scattering process.

However, an extra thermodynamics force called the diffusion force van Leeuwen et al. (1974), given by

[TABLE]

needs to be introduced when dealing with mixtures. Here, $n$ is the particle density and $w_{i}$ is the enthalpy per particle of component $i$. Further derivation of the transport equation in terms of thermodynamic forces proceed along lines similar to single component system and can be found in Ref. van Leeuwen et al. (1974). Here we state the final result analogous to Eq. (36) for component 1 as

[TABLE]

where, $x_{i}=n_{i}/(n_{1}+n_{2})$ being the particle number density fraction. An equation similar to above, also holds for component 2. The linear equations given in, Eqs. (36) and (41) are used in the later sections to derive explicit expressions for the transport coefficients.

IV.1 Single component system

In the present section we derive the transport coefficients for a single component system as described by the transport equation given in Eq. (36).

The observation that thermodynamic forces $X$, $Y^{\mu}$ and $\langle Z^{\mu\nu}\rangle$ appear as linearly independent quantities in Eq. (36), enables us to write the function $\phi$ of Eq. (IV) as

[TABLE]

where the unknown coefficients $A$, $B_{\mu}$ and $C_{\mu\nu}$ are still to be determined. The sign of $B_{\mu}$ is chosen in accordance with the sign of the vector force in Eq. (36). Inserting Eq. (42) into Eq. (36), the transport equation can be separated into three independent equations, given as

[TABLE]

where $\mathcal{L}[\phi]$ is the linearised collision operator, as defined in Eq. (IV).

We next define the macroscopic dissipative quantities, such as the viscous pressure and the heat flow which are functions of $\phi$. The viscous pressure is defined as van Leeuwen et al. (1973)

[TABLE]

the heat flow is defined as

[TABLE]

and the traceless viscous pressure is defined as

[TABLE]

The dissipative quantities can be written in a more transparent way using the following dimensionless inner product bracket notation

[TABLE]

where the quantities $z=m/T$ and $\tau=p^{\mu}U_{\mu}/T$ have been used. Inserting the expression for function $\phi$, given in Eq. (42) into the definitions of dissipative quantities defined in Eqs. (46-48), expresses these dissipative quantities in terms of bracket notation. Hence, the bulk viscous pressure is given as,

[TABLE]

such that $\pi^{\mu}=\Delta^{\mu\nu}p_{\nu}/T$. The heat flow is given as

[TABLE]

and the shear viscous flow as

[TABLE]

The quantities $\eta_{v}$, $\lambda=\Delta_{\mu\nu}\lambda^{\mu\nu}/3$ and $\eta_{s}$ stand for the bulk (volume) viscosity, heat conductivity and shear viscosity coefficients that appear as a constant of proportionality between thermodynamic forces and the dissipative quantities.

The technical details needed to compute the still unknown quantities $A$, $B^{\mu}$ and $C^{\mu\nu}$ into a tractable form, using collision integrals is given in Ref. van Leeuwen et al. (1973). Here, we simply write the expressions that can be used for computational purposes. The bulk viscosity is given by

[TABLE]

the heat conductivity is given by

[TABLE]

and the shear viscosity is given by

[TABLE]

The definitions of symbols $\alpha_{2}$, $\beta_{1}$ and $\gamma_{0}$ and the expression for the quantities $a_{22}$, $b_{11}$ and $c_{00}$ are given in Appendix A.

In Fig. 2 we use the relations as given in Eqs. (53)-(55) to calculate various transport coefficients for single component gas of baryons and mesons. The differential cross-sections that go into the expression of $a_{22}$, $b_{11}$ and $c_{00}$ are calculated using $K$-matrix formalism described in Sec. II for $\pi$, $K$ and $\eta$ while for nucleons ($N$) differential cross-section, we use the experimental phase-shift data from Ref. Workman et al. (2016). One must note that the temperature dependence of transport coefficients are highly dependent on the energy dependence of differential cross-sections. This is because the transport coefficients, through the quantities $a_{22}$, $b_{11}$ and $c_{00}$ (as given in Appendix A) are inversely proportional to the interaction cross-sections.

Fig. 3 shows, the role of cross-sections on the temperature dependence of transport coefficients. It can be seen from the comparision of the current algebra (CA) cross-sections of massive pions Weinberg (1966) which increase with the centre of mass energy to that of $K$-matrix cross-section which shows various peaks corresponding to various resonances that occur in $\pi\pi$ interaction throughout the energy spectrum. This makes the transport coefficients as shown in Fig. 2(2a-2c) to decrease with $T$ for current algebra and increase with $T$ for $K$-matrix. Similarly, for $\eta\eta$ interaction which has only a few resonances, the temperature dependence of transport coefficients show a dip at some given range of temperature, which can be alluded to the sharp rise in the cross-sections at corresponding energies (shown in Fig. 3). Comparing the transport coefficients among various mesons, we find that transport coefficients of a gas of $\eta>K>\pi$. This is because the total cross-section of $\pi>K>\eta$. For nucleons, the elastic cross-section decreases with the centre of mass energy, the same is reflected in the transport coefficients of nucleons at low $T$, where it drops even lower than for $\pi$s, but with increasing $T$, increases faster than for $\pi$s.

IV.2 Binary component system

The equation needed to obtain the transport coefficients for a mixture of two component gas is given in Eq. (41). The trial function is a linear combination of thermodynamic forces i.e.

[TABLE]

The only differences between the trial function for single component system Eq. (42) and $\phi_{k}$ of binary-component system is the diffusion force $Y_{1}^{\mu}$. Substituting function $\phi_{k}$ in Eq. (41) gives us

[TABLE]

where the factors $A_{1}$, $B_{1}^{\mu}$, $B^{\mu}_{1k}$ and $C^{\mu\nu}_{1}$ are unknown functions that are determined later. The law relating the traceless viscous pressure tensor to the hydrodynamic velocity and the law relating the viscous pressure to the divergence of hydrodynamic velocity as in Eqs. (50) and (55) do not change for mixtures. However, the law relating the heat flow to the temperature and pressure gradient, as in Eq. (51) needs to modified as,

[TABLE]

where $X^{\mu}_{q}$ is the generalized driving force of heat flow and $X_{1}^{\mu}$ is the diffusion driving force, which accounts for the flow due to gradients of different constituents of the system. The transport coefficients are defined as

[TABLE]

for the thermal conductivity and

[TABLE]

for the Dufour coefficients. This coefficients accounts for the heat flow in the presence of density gradients in a mixture. The other new coefficient for a mixture is the diffusion flow given by

[TABLE]

where the coefficient $l_{1q}$ is equal to the Dufour coefficient $l_{q1}$. The second coefficient $l_{11}$ is related to diffusion coefficient through the relation van Leeuwen et al. (1974), $D_{d}=\frac{l_{11}T}{nx_{1}x_{2}}$. This is given as

[TABLE]

As in a single component system, the transport coefficients can be written in collision bracket form, the details of which can be found in Ref. van Leeuwen et al. (1974) and in the Appendix A. Here we write the expression which can be used for computational purposes. The bulk viscosity is given as

[TABLE]

the shear viscosity is given as

[TABLE]

and the diffusion coefficient is given as

[TABLE]

The symbols and their relations to collision brackets are explained in the Appendix A.

One should note that the expressions given in Eqs. (53-55) for single component system and Eqs. (66-68) for binary component system corresponds to the first non-vanishing approximation of the transport coefficients (by approximation, we mean that the unknown coefficients $B^{\mu},C^{\mu\nu}$, etc. are expanded using a infinite series of Laguerre polynomials truncated at some order). Except for bulk viscosity, the first approximation corresponds to first non-vanishing value. For bulk viscosity, the non-vanishing value happens to be the third order approximation for single component system and second order approximation for binary component system. Thus, bulk viscosity for binary mixtures in the second order approximation calculated in this work depends only on the interaction among dissimilar species. The coefficient of shear viscosity, on the other hand, depends on $c_{12}$, $c_{11}$ and $c_{22}$ where, $c_{12}$ describes the interaction between dissimilar species and $c_{11},c_{22}$, describe the interaction among similar species (see Eqs. (86-88)).

The resulting transport coefficients for various binary mixtures are shown in Fig. 4. We have found both shear and bulk viscosities of the mixtures of two species lie in between the transport coefficients of the individual species. The dip seen in the shear viscosities of $\pi N$ and $KN$ can be attributed to resonances that appear in $\pi N$ and $KN$ interaction at the relevant energies which leads to an increase in the cross section and thus lowering the value of shear viscosity. Similarly, we show the diffusion coefficient of various binary components in Fig. 4(c) which depends on the density gradients in a mixture. We find that $KN$ system has largest diffusion coefficient at smaller temperatures and $\pi N$ system the lowest, but with increasing temperature, the coefficient for $KN$ system, shows a sharp decrease in its value and the $\pi K$ system shows a minimum. The open symbols in Fig. 4 correspond to transport coefficients at $\mu_{B}=100$ MeV. In the CE approximation $\mu_{B}$ enters implicitly in the expressions of transport coefficients via concentration or number densities of various reacting mixtures. The number densities were calculated using virial expansion that was described in Sec. II and are themselves function of temperature. We find that values of bulk viscosities are larger for large $\mu_{B}$ but gradually asymptotes towards $\mu_{B}=0$ MeV value, while shear viscosities values are smaller for large $\mu_{B}$ and gradually asymptotes towards the $\mu_{B}=0$ MeV values. The diffusion coefficient are mostly unaffected by the value of $\mu_{B}$ considered in the work.

IV.3 Multi component system

The derivation of transport coefficients for multi-component system follows the same line of reasoning as in case of the single and binary component system. The transport coefficient can be expressed transparently using the bracket notation which can be found in Refs. van Leeuwen et al. (1974); van Leeuwen (1975). Here, we only give the final expressions which can used for computational purposes. The bulk viscosity of a $N$ component gas can be written as

[TABLE]

while the coefficients $a_{k}$ satisfy the linear equations

[TABLE]

and the shear viscosity can be written as

[TABLE]

and the coefficients $c_{l}$ are solutions of

[TABLE]

In this work the Eq. (72) for the multi-component system can be written as

[TABLE]

and similarly for Eq. (70). The coefficients $c_{kl}$ and $a_{kl}$ depend on the scattering cross-section of the given channel $k$ and $l$ and the expressions in terms of collision integrals are given in the Appendix A (see Eqs. (91-93)). The zeros in $c_{kl}$ occur, when we do not have a resonance decaying in a channel $kl$.

The result of transport coefficients ($\eta_{v},\eta_{s}$), for various multi-channel processes is shown in Fig. (5a,5c) at $\mu_{B}=0$ MeV and Fig. (5b,5d) at $\mu_{B}=100$ MeV. We find that bulk viscosity turns out to be additive for a mixture of hadrons, in contrast to the shear viscosity, which decreases with the increase in number of components. This also explains why in RTA, for shear viscosities one should not add the relaxation time but the inverse of relaxation time for a multi-component system. The decrease in shear viscosities due to the increase in the components of the reacting mixture is evident, since it opens additional channels for reactions to occur and thus the overall cross-section of the system. Comparing the result of $\eta_{v}$ at $\mu_{B}=0$ MeV with that at $\mu_{B}=100$ MeV, we find that the values of $\eta_{v}$ are larger at large $\mu_{B}$. Similarly, we notice that at low $T$ the shear viscosities at finite $\mu_{B}$ is slightly lower than at zero $\mu_{B}$. However, with increasing temperature, the value of shear viscosity at finite $\mu_{B}$ overshoots that at zero $\mu_{B}$. This can be understood, since at large $T$ contributions from heavier baryonic states which have smaller cross-section increases and thus increases the viscosity. At lower $\mu_{B}$ their concentration is smaller, hence their effect is not noticeable but increasing $\mu_{B}$ increases their concentration (the cross-section remains the same) and hence their effect on viscosity also increases.

The variation of $\eta_{v}/s$ and $\eta_{s}/s$ as a function of temperature is shown in Fig. (6a-b). Our results of $\eta_{v}/s$ is an increasing function of $T$ for $T<150$ MeV and decreasing for $T>150$ MeV. At $\mu_{B}=100$ MeV, we find the magnitude of the peak, seen in $\eta_{v}/s$ is larger than at $\mu_{B}=0$ MeV. Similarly, we find that $\eta_{s}/s$ decreases with temperature consistent with previous results in this regard Gorenstein et al. (2008); Kadam and Mishra (2015); Moroz (2013); Wiranata et al. (2013); Weil et al. (2016). However, we find that the result of $\eta_{s}/s$ at $\mu_{B}=0$ violates the AdS/CFT bound around a temperature of $T=160$ MeV, while the result of $\mu_{B}=100$ MeV remains above the bound and asymptotically approaches it at higher temperatures. Of course, the temperature where the violation of the AdS/CFT bound occurs, is in between the deconfinement temperature which is around $T\approx 155-165$ MeV Aoki et al. (2009); Bazavov et al. (2012), where our model should break-down. It is also interesting to note that peak in the ratio $\eta_{v}/s$ happens to be around the same temperature where the ratio $\eta_{s}/s$ violates the AdS/CFT bound.

Let us now discuss the comparison of our result with the calculations that has been previously reported in the literature at $\mu_{B}=0$ MeV. In EVHRG (Excluded volume HRG model, $\eta_{v}/s$ monotonically decreases as a function of temperature $T$ in contrast to our results which shows a peak structure and further one can note that magnitude of $\eta_{s}/s$ in the EVHRG model is a factor of ten more than our results. The first reason for this is that, the calculation of $\eta_{v}/s$ is carried out using RTA Kadam and Mishra (2015), in the EVHRG model using momentum independent relaxation times which is quantitatively different from that of CE method used in the current work. The difference in temperature variation can be attributed to use of constant cross-section in the EVHRG model calculations compared to energy dependent cross-sections used in our work. Moroz Moroz (2013) uses cross sections from the UrQMD model, including elastic plus resonance processes calculated in the CE approximation. The $\eta_{s}/s$ result from Moroz calculation is qualitatively and quantitatively similar to our calculations. Some discrepancies are still present because of the use of some constant cross-sections to describe non-resonant interaction in Moroz’s calculation.

The $\eta_{s}/s$ calculation in EVHRG model Gorenstein et al. (2008) is done assuming all hadrons have the same hard-core radius $r=0.5$ fm. Apart from the fact that the value of $r$ used is model dependent, one must note that, they also assume that the shear viscosity is additive for a mixture of hadrons, contrary to our results. Although $\eta_{s}/s$ decreases with temperature, but the slope is less steeper than our calculation. This is because in Ref. Gorenstein et al. (2008) both $\eta_{s}$ and $s$ increase, as degeneracies increase. However, in our case $\eta$ decreases and $s$ increases as degeneracies increase. Both this feature make the slope of $\eta_{s}/s$ steeper than Ref. Gorenstein et al. (2008). Wiranata et al. Wiranata et al. (2013) used $K$-matrix formalism for calculating $\eta_{s}/s$ in a hadronic gas consisting of $\pi-K-N-\eta$. Their result is around six times larger than ours at low $T$ and about two times larger in high $T$. The discrepancies between the two results are first, due to the fact that we have used a larger spectrum of interacting hadrons and resonances. Secondly, and an important difference is that Ref. Wiranata et al. (2013) did not include the $NN$ mutual interaction, since their cross-section were solely using $K$- matrix formalism, where as we parameterize $NN$ experimental phase shifts to calculate the differential cross-section. Owing, to the fact that $NN$ cross-section are larger at small $\sqrt{s}$ as has been previously discussed, their contribution to transport coefficients is quite different and dramatic than other resonant interaction present in $K$-matrix formalism. SMASH (Simulating Many Accelerated Strongly-interacting Hadrons) Weil et al. (2016), which is a transport code, uses Green-Kubo formalism to calculate $\eta_{s}/s$ for hadronic gas mixture. One of the common feature between our model and SMASH is the treatment of interactions through resonances, which have a non-zero lifetime. Our result of $\eta/s$ is in good agreement with SMASH within temperature range of $80-110$ MeV. But after $T\sim 120$ MeV, we find that the SMASH result saturates and forms a plateau at higher temperature. The same trend is also seen in other transport codes for e.g UrQMD Bleicher et al. (1999). The crucial difference between our approach and SMASH is that, SMASH utilises a feedback loop between the relaxation time and resonance lifetimes whereas our approach does not Rose et al. (2018).

V Summary

In summary, we have calculated the transport coefficients for a multi-component hadronic gas. The thermodynamic quantities are calculated using the $S$-matrix based hadron resonance gas model. The phase-shifts required for the calculation of $S$-matrix was calculated using the $K$-matrix formalism for all hadrons except for nucleons, for which we directly parameterize the experimental phase-shifts. The transport coefficients were calculated using the Chapman-Enskog method. Such a method utilizes the energy dependence of cross-sections to calculate the temperature dependence of transport coefficients.

We start with various single component gas systems and gradually add different species of hadrons to finally form a multi-component gas mixture. We found that adding new species into the mixture, opens up new channels for interaction to occur, which leads to an increase in cross-section and thus reducing the shear viscosity. Similarly, we calculate the transport coefficients at zero and non-zero $\mu_{B}$. We found that increasing $\mu_{B}$ increase the contribution of higher mass baryons, which have smaller cross-section, in the transport coefficients. This leads to the increase in the value of $\eta/s$ at higher temperatures. Interestingly, we found that at the temperature around $T\approx 160$ MeV, the ratio $\eta_{v}/s$ shows a maximum and around the same temperature, the $\eta_{s}/s$ starts violating the AdS/CFT bound. A maximum in $\eta_{v}/s$ is a signature of crossover transition, that has been seen previously in molecular gases Meier et al. (2005). Similarly, the violation of AdS/CFT bound may signal the breakdown of a simple model like the HRG and that the non-perturbative nature of physics in this regime. However, increasing $\mu_{B}$, evades such a violation of the bound to larger temperature. Finally, we compute and compare the ratio of $\eta_{s}/s$ and $\eta_{v}/s$, with other models in the literature. Our calculation show qualitatively similar features, with models that use energy dependent cross-section in the relevant temperature ranges. It is also interesting to see that a model which assumes the hadrons to be gases with interaction governed by $S$-matrix elements, which are basically resonances, capture the essential physics of transport coefficients in the $T-\mu_{B}$ plane.

A few future directions for this work. The crucial assumption that has been done in this work, is the use of Maxwell-Boltzmann (MB) distribution function, which may not be valid for large chemical potentials. Then, one needs to solve the full quantum Boltzmann equation in the CE method. In that case, the polynomials (the Laguerre polynomials in the case of MB) satisfying them are not known, and we have to find them order by order, as has been done in Ref. Itakura et al. (2008). Another important direction would be to include a feedback mechanism between the relaxation time and resonance lifetimes as is done in transport codes Weil et al. (2016). For example, in this work we have considered resonances like $\rho$, $\Delta$, etc. as unstable particles, which although contribute to the cross-section, by themselves are not a part of the mixture. This is only valid, if the lifetime of resonance is shorter than mean free time of the system. However, if the resonance lifetime is comparable or larger to the mean free time of the system , interaction can only occur until the resonance decays. Thus, the relaxation time in such cases is limited by the lifetime of the resonance and not by the mean free time.

Appendix A A

In the following appendix, we define the various symbols and expressions that were used in the main text.

For single component system the symbols $\alpha_{2}$, $\beta_{1}$ and $\gamma_{0}$ are defined as

[TABLE]

where $\gamma=c_{p}/c_{v}$. The quantities $a_{22}$, $b_{11}$ and $c_{00}$ are defined in terms of relativistic omega integrals, $\omega_{i}^{(j)}$

[TABLE]

where the definitions of relativistic omega integrals are given in van Leeuwen et al. (1973) and can be written as

[TABLE]

[TABLE]

where $\sigma(\psi,\theta)$ is the differential cross-section for interaction between two identical particles, expressed through the quantities $\psi$ and angle $\theta$ between the initial and final hadrons defined as

[TABLE]

where $p_{1}$ and $p_{2}$ are the initial four-momenta of the two colliding hadrons.

For binary component system the symbols $\alpha_{i}$, $\delta_{i}$ and $\gamma_{i}$, where $i=1,2$ are defined as

[TABLE]

where $h_{i}=K_{3}(z_{i})/K_{2}(z_{i})$ is the specific enthalpy of species $i$, $c_{i}=\rho_{i}/\rho$ is the mass fraction of species $i$. $\rho_{i}$ is the mass density, which is mass times the number density of species $i$ and $\rho$ is the total mass density. Similarly $x_{i}=n_{i}/n$ is the number density fractions of species $i$, where $n_{i}$ is particle number density of species $i$ and $n$ is the total number density. The quantity $\gamma_{(i)}=c_{p,i}/c_{v,i}$ is the ratio of specific heats of species $i$. The quantities $a_{ii}$, $c_{ii}$, $c_{ij}$, $b_{ii}$ and $\Delta_{c}$ are defined in terms of relativistic omega integrals, $\omega_{ijkl}^{(m)}(\sigma_{uv})$

[TABLE]

where $\sigma_{uv}$ is the cross-section between particles $u$ and $v$. The coefficients $c_{00}(z_{k})$ accounts for contribution from interaction between identical species of type $k$ as given in Eq. (79). The reduced mass $\mu$ is given as $\mu=m_{1}m_{2}/(m_{1}+m_{2})$. The abbreviations $Z$ and $\zeta$ are given as $Z=M/T$ and $\zeta=2\mu/T$, where $M=m_{1}+m_{2}$ is the total mass. The definitions of relativistic omega integrals are given in Refs. van Leeuwen et al. (1974); van Leeuwen (1975) and we do not write them here.

For multi-component system, the coefficients $a_{kl}$, $c_{kl}$ are given as

[TABLE]

where $a_{22}(kl)$, $c_{22}(kl)$ and $c_{12}(kl)$ are the expressions given in Eqs. (85-88), with subscripts 1 and 2 replaced by $k$ and $l$.

Acknowledgement

BM acknowledges financial support from J C Bose National Fellowship of DST, Government of India. SS and AD acknowledge financial support from DAE, Government of India.