Kubo formulae for the shear and bulk viscosity relaxation times and the scalar field theory shear $\tau_\pi$ calculation

TL;DR

This paper derives Kubo formulas for shear and bulk relaxation times using quantum field theory, and computes the shear relaxation time in massless bb theory via diagrammatic methods, linking it to thermal width.

Contribution

It introduces Kubo formulas for relaxation times and calculates shear relaxation time in bb theory using a diagrammatic approach, connecting microscopic dynamics to hydrodynamic parameters.

Findings

Kubo formulas for shear and bulk relaxation times derived from stress-energy tensor response functions.

Shear relaxation time computed in massless bb theory using integral equations.

Shear relaxation time inversely related to thermal width, consistent with theoretical expectations.

Abstract

In this paper we provide a quantum field theoretical study on the shear and bulk relaxation times. First, we find Kubo formulas for the shear and the bulk relaxation times, respectively. They are found by examining response functions of the stress-energy tensor. We use general properties of correlation functions and the gravitational Ward identity to parametrize analytical structures of the Green functions describing both sound and diffusion mode. We find that the hydrodynamic limits of the real parts of the respective energy-momentum tensor correlation functions provide us with the method of computing both the shear and bulk viscosity relaxation times. Next, we calculate the shear viscosity relaxation time using the diagrammatic approach in the Keldysh basis for the massless $\lambda\phi^4$ theory. We derive a respective integral equation which enables us to compute $\eta\tau_\pi$ and then we extract the shear relaxation time. The relaxation time is shown to be inversely related to the thermal width as it should be.

Figures

Equations

$${\partial }_{\mu }{T}^{\mu \nu }$$

$$T^{\mu \nu }=ϵu^{\mu }{u}^{\nu }-{\Delta }^{\mu \nu }(P+\Pi )+{\pi }^{\mu \nu },$$

$${\Pi }_{\text{NS}}$$

$${\pi }_{\text{NS}}^{\mu \nu }$$

$${\Pi }_{\text{NS}}$$

$${\pi }_{\text{NS}}^{ij}$$

$$\displaystyle\delta T^{ij}=\pi^{ij}_{\textrm{NS}}+\Pi_{\textrm{NS}}=D_{T}\bigg{(}\partial^{i}T^{j0}+\partial^{j}T^{i0}-\frac{2}{3}g^{ij}\partial_{l}T^{l0}\bigg{)}+g^{ij}\gamma\partial_{l}T^{l0}.$$

$${\pi }^{\mu \nu }$$

$$\Pi$$

$${\partial }_t\Pi$$

$${\partial }_t{\pi }^{ij}$$

$${\partial }_t{T}^{0k}=-{\partial }_l{T}^{lk}.$$

$$0$$

$${\partial }_t^2ϵ$$

$$0$$

$$\delta ⟨\hat{A}(t,\mathbf{x})⟩=\int d^4{x}^{\prime }{G}_R(t-t^{\prime },\mathbf{x}-{\mathbf{x}}^{\prime })\theta (-t^{\prime })e^{\epsilon t^{\prime }}f(\mathbf{x}),$$

$$G_R(t-t^{\prime },\mathbf{x}-{\mathbf{x}}^{\prime })=-i\theta (t-t^{\prime })⟨[{\hat{A}}_H(t,\mathbf{x}),{\hat{A}}_H(t^{\prime },{\mathbf{x}}^{\prime })]⟩,$$

$$G_A(t-t^{\prime },\mathbf{x}-{\mathbf{x}}^{\prime })=i\theta (t^{\prime }-t)⟨[{\hat{A}}_H(t,\mathbf{x}),{\hat{A}}_H(t^{\prime },{\mathbf{x}}^{\prime })]⟩.$$

$$D_A{G}_R(t-t^{\prime },\mathbf{x}-{\mathbf{x}}^{\prime })=d_A\delta (t-t^{\prime })\delta (\mathbf{x}-{\mathbf{x}}^{\prime }),$$

$$D_A\delta ⟨\hat{A}(t,\mathbf{x})⟩=0.$$

$${\rho }^{AA}(k)=\int d^{4}x{e}^{ikx}⟨[{\hat{A}}_H(x),{\hat{A}}_H(0)]⟩,$$

$$\displaystyle\rho^{AA}(k)=\frac{1}{Z_{0}}\sum_{m,n}\big{(}e^{-\beta E_{n}}-e^{-\beta E_{m}}\big{)}(2\pi)^{4}\delta(k-p_{m}+p_{n})\big{|}\langle p_{n}|\hat{A}|p_{m}\rangle\big{|}^{2},$$

$${\rho }^{AA}(-\omega ,-\mathbf{k})=-{\rho }^{AA}(\omega ,\mathbf{k}).$$

$$G_{R/A}(\omega ,\mathbf{k})=\int \frac{d{\omega }^{\prime }}{2\pi }\frac{{\rho }^{AA}({\omega }^{\prime },\mathbf{k})}{{\omega }^{\prime }-\omega \mp iϵ},$$

$$\text{Re}{G}_R(\omega ,\mathbf{k})$$

$$\text{Im}{G}_R(\omega ,\mathbf{k})$$

$$\displaystyle\partial_{\alpha}\big{[}\bar{G}^{\alpha\beta,\mu\nu}(x,x^{\prime})-\delta^{(4)}(x-x^{\prime})\big{(}g^{\beta\mu}\langle\hat{T}^{\alpha\nu}(x^{\prime})\rangle+g^{\beta\nu}\langle\hat{T}^{\alpha\mu}(x^{\prime})\rangle-g^{\alpha\beta}\langle\hat{T}^{\mu\nu}(x^{\prime})\rangle\big{)}\big{]}=0,$$

$$\displaystyle k_{\alpha}\big{(}\bar{G}^{\alpha\beta,\mu\nu}(k)-g^{\beta\mu}\langle\hat{T}^{\alpha\nu}\rangle-g^{\beta\nu}\langle\hat{T}^{\alpha\mu}\rangle+g^{\alpha\beta}\langle\hat{T}^{\mu\nu}\rangle\big{)}=0,$$

$$G_R^{\mu \nu ,\alpha \beta }(x,x^{\prime })=-i\theta (x_0-x_0^{\prime })⟨[{\hat{T}}^{\mu \nu }(x),{\hat{T}}^{\alpha \beta }(x^{\prime })]⟩$$

$$\delta \hat{H}(t)=-\int d^{3}x\theta (-t)e^{\epsilon t}{\hat{T}}^{x0}(t,\mathbf{x}){\beta }_x(y)$$

$$\delta ⟨{\hat{T}}^{x0}(t,k_y)⟩={\beta }_x(k_y){\int }_{-\infty }^{\infty }d{t}^{\prime }\theta (-t^{\prime })e^{\epsilon t^{\prime }}{\bar{G}}_R^{x0,x0}(t-t^{\prime },k_y)$$

Full text

I Introduction

Relativistic viscous hydrodynamics seems to be perfectly suited to investigate and understand collective phenomena characteristic of strongly interacting matter produced in heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), see Refs. Gale:2013da ; Heinz:2013th and references therein. Transport coefficients are inherent ingredients of the hydrodynamic description. They control the dynamics of a fluid and to find any of the coefficients we need to employ a microscopic theory. Due to the multiple-scale nature of the problem, determination of the full set of the transport coefficients is not a trivial task. Since the shear and the bulk viscosity have already been examined in many papers, our particular interest is in the kinetic coefficients of the second-order hydrodynamics, namely the shear and the bulk relaxation times. These relaxation times fix the characteristic time scales at which the dissipative currents relax to their first-order solutions. In this study, we use field theoretical approach to investigate the relaxation times.

To understand the microscopic dynamics of plasma constituents and determine the values of any transport coefficient one needs to use either kinetic theory or quantum field theory in the weakly coupled limit. Nevertheless, the quark-gluon plasma, studied experimentally at RHIC and LHC, is believed to achieve this limit only at sufficiently high temperatures. It was just the investigation of the shear and bulk viscosities that established numerous methods of evaluation of transport coefficients in general. Within the imaginary-time formalism of the scalar field theory, shear and bulk viscosities were calculated Jeon:1992kk ; Jeon:1994if which showed how to properly handle the contributing diagrams to obtain the leading order contribution. In short, in any diagrammatic approach, the ladder diagrams have to be resummed to get the leading order transport coefficients due to the pinching pole effect arising when small frequency and long-wavelength limits of the correlation function are taken.

With this knowledge transport coefficients of QED Gagnon:2006hi ; Gagnon:2007qt and of QCD at leading log order ValleBasagoiti:2002ir were obtained. The kinetic theory approach has appeared, in turn, to be very effective, as in the case of scalar field theory Jeon:1995zm , and also successful in determining transport coefficients of QCD medium in the leading order Arnold:2000dr ; Arnold:2003zc . The latter study shows, in particular, the effectiveness of the kinetic theory to find leading order results and, at the same time, some difficulties to go beyond it Arnold:2002zm .

To date, many approaches to the relaxation times have been developed within kinetic theory. In York:2008rr it was shown that the ratio of the shear viscosity over its relaxation time is proportional to the enthalpy density and the proportionality factor is slightly different in QCD than in the $\phi^{4}$ theory. A similar relation for the ratio was established within the so-called 14-moment approximation to the Boltzmann equation. And generally the method of moments, first proposed by Grad and then further developed by Israel and Stewart, has been examined comprehensively in Denicol:2010xn ; Denicol:2012es ; Denicol:2012cn ; Denicol:2014vaa . The moment approaches show, in particular, that in order to compute the viscosity coefficients it is necessary to invert the collision operator and to determine the relaxation times one has to find the eigenvalues and eigenvectors of the operator.

There exist also field-theoretical studies on the shear relaxation time Denicol:2011fa ; Denicol:2011ef , where a general form of a retarded Green function is considered. This work actually shows the microscopic origin of the shear relaxation time, which is found to be inversely proportional to the imaginary part of the pole of the particle propagator. The Kubo formula engaging shear relaxation time was found for conformal systems Baier:2007ix ; Moore:2010bu via the response of a system to small and smooth perturbations of a background metric. The projection operator method was used in Koide:2009sy to obtain the shear relaxation time.

As for the bulk relaxation time, the projection operator method was used in Huang:2010sa to obtain it. Kubo formula for the product of the bulk viscosity and the bulk relaxation time, $\zeta\tau_{\Pi}$, can be also deduced from response functions studied in Ref. Hong:2010at using a mixture of the effective kinetic theory and metric perturbations.

Our goal here is to figure out, through a standard formulation of quantum field theory, how the shear and the bulk relaxation times are related to microscopic quantities. In order to obtain any transport coefficient on this ground, one needs to know the corresponding Kubo-type relation. In this study, by studying both the sound and the shear modes we are able to find a set of new Kubo-type relations. In particular, we find the Kubo formula which relates the product $\zeta\tau_{\Pi}$ to the second derivative of the real part of the pressure-pressure response function with respect to the frequency. It is worth emphasizing that the obtained Kubo relation provides us with an explicit prescription on how to examine the bulk relaxation time from the quantum field theory perspective. All formulas studied here were obtained by making use of hydrodynamic limits, general properties of a response function and the Ward identity. The method we employ was introduced in Jeon:2015dfa and here we provide its extension. In the low-frequency and long-wavelength limits, each of the correlation functions are related to some set of kinetic coefficients.

In this paper we focus only on the leading order of the response function, which is sufficient to study shear effects and, in particular, find the shear relaxation time. In order to analyze bulk effects one needs to go to the next-to-leading order calculation, which will be the purpose of a future work.

To find the value of the shear relaxation time we use diagrammatic methods of the closed-time path (Keldysh-Schwinger) formalism. As advocated in Wang:2002nba , the $(r,a)$ or Keldysh basis serves very convenient framework for such considerations, mostly due to the vanishing $aa$ propagator component. We work with the massless $\lambda\phi^{4}$ theory in the weakly coupled limit. We start with the one loop case and then perform the resummation over ladder diagrams, which contribute at the same order. The one-loop approximation allows us to determine what is the typical scale at which $\eta\tau_{\pi}$ appears. This is found to be of the order $1/\Gamma_{p}^{2}$, where $\Gamma_{p}$ is the thermal width and it is directly related to the mean free path of the system constituents. With the knowledge on the one-loop result for $\eta$ we are able to extract $\tau_{\pi}$, which scales as $1/\Gamma_{p}$. The summation of ladder diagrams leads us to manipulation on the four-point Green functions which couple to each other through the Bethe-Salpeter equation. In the $(r,a)$ basis, however, only the $G_{aarr}$ component matters and consequently the Bethe-Salpeter equation decouples. To compute the shear viscosity the integral equation for $\text{Im}G_{aarr}$ needs to be solved. We show that to compute $\eta\tau_{\pi}$ one needs to solve the integral equation for $\partial_{\omega}\text{Re}G_{aarr}$, which requires us to introduce a new type of the effective vertex. Both $\eta$ and $\eta\tau_{\pi}$ are evaluated numerically to extract the shear relaxation time.

The remaining part of the paper is organized as follows. In Sec. II we provide a brief introduction to viscous hydrodynamics, mainly to write the dispersion relations of the two hydrodynamic modes arising from the momentum and energy dissipation. We assume that there are no other currents coupled to the energy-momentum tensor. In Sec. III general properties and constraints of the response functions are discussed. Section IV presents the way on how to parametrize the response functions to the longitudinal and transverse hydrodynamic fluctuations so that to reproduce the corresponding dispersion relations. Subsequently, we find Kubo-type relations. In Sec. V the expression which allows us to extract the relaxation time for shear viscosity is derived within diagrammatic methods in the real-time formalism. We perform full leading order analysis providing summation over multi-loop diagrams. We also introduce a new effective vertex, which enters the formula for evaluation of $\eta\tau_{\pi}$. In Sec. VI we evaluate $\tau_{\pi}$ and $\langle\epsilon+P\rangle\tau_{\pi}/\eta$ numerically as a function of the constant coupling. We also discuss our results in the context of kinetic theory findings. We conclude in Sec. VII.

II Hydrodynamic modes

Here, we very briefly introduce basic equations of hydrodynamics, mostly to fix a starting point for the further discussions. For more comprehensive analysis we refer the reader to, for example, Denicol:2012es ; Jeon:2015dfa ; Betz:2009zz ; Kovtun:2012rj . The discussions in this section and Secs. III and IV closely follow Jeon:2015dfa .

Hydrodynamics as a long-wavelength and low-frequency many-body effective theory provides a macroscopic description of a system which is close to thermal equilibrium. It governs the evolution of any fluid in terms of flows of its conserved quantities, such as the energy, momentum or baryon current. Here, we study phenomena associated only with the energy-momentum conservation. The corresponding conservation law is the continuity equation of the energy-momentum tensor $T^{\mu\nu}$

[TABLE]

When the system is approaching local thermal equilibrium, its relevant behavior is fairly well described by the viscous hydrodynamics, for which the energy-momentum tensor takes the form

[TABLE]

where $\epsilon$ is the energy density, $P$ is the thermodynamic pressure, $u^{\mu}$ are the components of the flow velocity with the normalization condition $u^{\mu}u_{\mu}=1$, $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ is the projection operator with $u_{\mu}\Delta^{\mu\nu}=0$, and the Minkowski metric is $g^{\mu\nu}=(1,-1,-1,-1)$. The terms $\Pi$ and $\pi^{\mu\nu}$ are the bulk viscous pressure and the shear stress tensor, respectively. They are viscous corrections which contain the dynamics of the dissipative medium approaching the equilibrium state. The shear tensor is symmetric, traceless, $\pi^{\mu}_{\mu}=0$, and transverse, $u_{\mu}\pi^{\mu\nu}=0$, and the bulk pressure $\Pi$ is a correction to the thermodynamic pressure. These corrections are assumed to be small for the nearly equilibrium state of a system. If the static equilibrium limit is achieved, the dissipative corrections in Eq. (2) vanish and one reproduces the ideal hydrodynamic stress-energy tensor. Within the Navier-Stokes approach, the viscous corrections are obtained from the gradient expansion of the energy density and flow velocity. Then, up to the linear terms, only the corrections proportional to $\partial^{\mu}u^{\nu}$ matter. The dissipative currents $\Pi$ and $\pi^{\mu\nu}$ take the following forms

[TABLE]

where $\zeta$ and $\eta$ are the bulk and shear viscosities which determine transport phenomena of the energy and momentum, and $\Delta^{\mu\nu}_{\alpha\beta}\equiv(\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\mu}_{\beta}\Delta^{\nu}_{\alpha}-2/3\Delta^{\mu\nu}\Delta_{\alpha\beta})/2$ is the traceless and transverse projection operator. In the fluid cell rest frame the dissipative currents may be written as

[TABLE]

where $\gamma=\zeta/(\epsilon+P)$ and $D_{T}=\eta/(\epsilon+P)$. Accordingly, the full spatial viscous correction to the energy-momentum tensor of a viscous fluid is

[TABLE]

in the rest frame of the fluid cell at the position $x$.

In fact, the Navier-Stokes theory of a relativistic fluid is acausal and unstable. Generally speaking, the viscous currents must be allowed to take some time when responding to changes in the thermodynamic forces. This requires taking the next terms in gradient expansion, and the corresponding hydrodynamics is the second-order Israel-Stewart theory. The bulk viscous pressure and the shear stress tensor are then subject to the relaxation equations. In a general frame, these are given by

[TABLE]

where we ignored the non-linear terms as they are not relevant for our study. We have also used the notation $A^{\langle\mu\nu\rangle}\equiv\Delta^{\mu\nu}_{\alpha\beta}A^{\alpha\beta}$ for the spin-2 component of a rank-2 tensor. The new transport coefficients, $\tau_{\Pi}$ and $\tau_{\pi}$, are relaxation times for the bulk and shear viscosities, respectively. They determine how fast the bulk pressure and the shear tensor relax to the respective Navier-Stokes forms given by Eqs. (3) and (4), respectively. As a result, causality of the theory is maintained if the relaxation times satisfy certain restrictions Pu:2009fj . The relaxation equations take the following forms in the local rest frame

[TABLE]

If there are no other currents coupled to the energy-momentum tensor, there are two hydrodynamic modes that determine the behavior of the system. These are the diffusion and sound modes. The diffusion mode describes fluid flow in the direction transverse to the flow velocity. It appears as a consequence of the momentum conservation, which is

[TABLE]

When the relaxation equation (11) and the Navier-Stokes form of the shear tensor (6) are implemented to the momentum conservation law (12), we get the corresponding hydrodynamic equation, which is the equation of motion of the transverse part of stress tensor

[TABLE]

where $\pi^{i}_{T}=\epsilon_{ijk}\partial_{j}T^{k0}$. The corresponding dispersion relation is then found to be

[TABLE]

with $\omega$ and ${\bf k}$ being the frequency and wavevector of the momentum diffusion excitation.

The other mode is associated with small disturbances in dynamic variables propagating longitudinally in the medium. The conservation law in the local rest frame then is

[TABLE]

By multiplying Eq. (15) by $(\tau_{\pi}\partial_{t}+1)(\tau_{\Pi}\partial_{t}+1)$, making use of the relaxation equations (10) and (11) and using the Navier-Stokes forms of the stress tensor (6) and bulk pressure (5), we get the equation of motion for the energy density deviation $\delta\epsilon$,

[TABLE]

where $v_{s}^{2}=\partial P/\partial\epsilon$ is the speed of sound. The solution to Eq. (16) is provided by the following dispersion relation:

[TABLE]

The dispersion relations (14) and (II) play an essential role in further analysis as they encode full information on the relaxation times for viscosities. This information should be also contained in the pole structure of the respective retarded Green function and this is the subject of the next sections.

III Response functions

Linear response theory is a natural quantum-mechanical framework to examine systems exhibiting small deviations from equilibrium. Within the linear response theory one is able to express quantities characteristic of the nonequilibrium state of a fluid in terms of time dependent correlation functions of the equilibrium state. The linear response theory is explained in many textbooks, see for example Kapusta:2006pm , and here we restrict ourselves to discuss only those properties of response functions relevant to our study of transport coefficients.

When a system undergoes small perturbations, the deviation of an observable $A$ from equilibrium is encoded in equilibrium response function as

[TABLE]

where $t>0$, $\theta(t)$ is the standard step function, $f({\bf x})$ is an external perturbing force coupled to $\langle\hat{A}\rangle$ and acts on the system with infinitesimally slow rate $\varepsilon$, and $\langle\cdots\rangle$ means the thermal expectation value. $G_{R}$ is the retarded response function corresponding to the Hermitian operator $\hat{A}$,

[TABLE]

where $\hat{A}_{H}$ stands for the operator in the Heisenberg picture. For further analysis it is also convenient to introduce the advanced correlation function, which is

[TABLE]

Suppose the retarded Green function $G_{R}$ satisfies the following equation of motion:

[TABLE]

where $D_{A}$ is some operator such that $G_{R}$ is its generalized Green function and $d_{A}$ may contain a finite number of derivatives. For positive values of $t$, $t\neq t^{\prime}$ since $t^{\prime}$ is restricted by $\theta(-t^{\prime})$. Therefore, for positive $t$, $\delta\langle\hat{A}\rangle$ satisfies the following evolution equation:

[TABLE]

Accordingly, the evolution equation is known whenever one finds the pole structure of the response function Kadanoff .

The formula (18) shows explicitly that the linear response of the system is expressed in terms of a retarded Green function of Heisenberg operators. To study the retarded Green functions, it is convenient to introduce the spectral density defined by the thermal expectation value of the commutator

[TABLE]

where $k=(\omega,{\bf k})$, which may be expressed as

[TABLE]

when $\hat{A}$ is Hermitian. Here $|m\rangle$ is the simultaneous eigenstate of the system’s total Hamiltonian $\hat{H}$ and the total momentum $\hat{\bf P}$ with the eigenvalue $p_{m}=(E_{m},{\bf p}_{m})$. Relying on the fact that for any observable the corresponding operator must be Hermitian, $\hat{A}^{\dagger}=\hat{A}$, one can derive

[TABLE]

For an equilibrium system, which is isotropic, $\rho^{AA}(\omega,{\bf k})$ must preserve a rotational invariance so that it depends on momentum only through its absolute value $|{\bf k}|$. Therefore, the spectral density is an odd function of $\omega$, that is, $\rho^{AA}(-\omega,{\bf k})=-\rho^{AA}(\omega,{\bf k})$.

In the spectral representation the retarded and advanced Green functions are given by

[TABLE]

where the upper sign $(-)$ corresponds to the retarded function and the lower sign $(+)$ to the advanced one. By extracting the principal value of the integral in Eq. (26) from the imaginary part one obtains the following relations:

[TABLE]

where $\mathcal{P}$ stands for the principal value. By changing the sign of $\omega^{\prime}$ in the formula (27) and using the fact that the spectral function is an odd function of the frequency one observes that $\textrm{Re}\,G_{R}(\omega,{\bf k})=\textrm{Re}\,G_{R}(-\omega,{\bf k})$, that is, the real part of the retarded and the advanced Green function is an even function of frequency. Moreover, the imaginary part of the retarded response function, since related directly to the spectral function, is an odd function of $\omega$. These facts will be frequently used in the next parts of this paper.

Due to the fact that the stress-energy tensor represents both the conserved current as well as the generators of the space time evolution, the correlation functions of $T^{\mu\nu}$ are not so simple. To determine them correctly, one must first start with the following gravitational Ward identity Deser:1967zzf :

[TABLE]

which becomes in the momentum space

[TABLE]

where $k=(\omega,{\bf k})$. The identity (29) is most conveniently derived in the imaginary-time metric. The two-point functions of $T^{\mu\nu}$ is then obtained by taking the second functional derivative of the partition function with respect to the imaginary-time metric. Going to the flat space and then analytic continuing to the real space give (29). Let us stress that $\bar{G}_{R}^{\mu\nu,\alpha\beta}(x,x^{\prime})$ is the response function, which is not, in general, the same as

[TABLE]

due to the presence of the single stress-energy tensor average terms in Eq. (30). $\bar{G}_{R}$ differs from $G_{R}$ by terms containing $\delta(x-x^{\prime})$. Let us add that the formula (30), when combined with the continuity equations, fixes actually a set of constraints that the response functions, corresponding to different components of the energy-momentum tensor, must maintain.

IV Analytic structure of response functions and Kubo formulas

The general properties of a correlation function supported by the Ward identity constrain its analytical form enough so that one is able to parametrize it for both the propagating and diffusive mode. Here, we parametrize the correlation functions in terms of the linear response method.

IV.1 Response function to transverse fluctuations

The perturbing Hamiltonian for the shear flow is

[TABLE]

Note that the external force $\beta_{x}(y)$ is related to the flow velocity component in the $x$-direction which only varies in the perpendicular $y$-direction. Hence, at $t=0$, this sets up a system with a non-zero shear flow. The corresponding linear response is

[TABLE]

for $t>0$. From the Ward identity (30), one finds

[TABLE]

When these two equations are combined, one gets

[TABLE]

In the $\omega\rightarrow 0$ limit, $\bar{G}_{R}^{x0,x0}(\omega,k_{y})$ must have a well defined limit since it is a thermodynamic quantity. Moreover, both correlation functions must be well behaved in the $k_{y}\rightarrow 0$ limit. Using these arguments and the fact that the imaginary part of the retarded Green function must be an odd function of $\omega$, one can parametrize $\bar{G}_{R}^{xy,xy}(\omega,k_{y})$ as

[TABLE]

and

[TABLE]

where $g_{T}(k_{y})=\bar{G}^{x0,x0}(0,k_{y})=P+g_{\pi\pi}(k_{y})$ comes from Eq. (209) in Appendix A. The functions $A$, $B$, and $D$ have the form

[TABLE]

where $D_{R}(\omega,k_{y})$ and $D_{I}(\omega,k_{y})$ are real-valued even functions of $\omega$ and $k_{y}$. The real parts $D_{R}$ and $B_{R}$ must have a non-zero limit as $\omega\rightarrow 0$ and $k_{y}\rightarrow 0$. All other parts of $A$, $B$, and $D$ must have finite limits as $\omega\rightarrow 0$ and $k_{y}\rightarrow 0$. Dynamical information in the hydrodynamic limit is contained in the constants $D_{R}(0,0)$, $D_{I}(0,0)$, $B_{R}(0,0)$, etc.

The pole structure of a Green function determines the corresponding dispersion relation of a given hydrodynamic mode. Therefore, the dispersion relation of the diffusive excitation is dictated by the form of the denominator of the function (37), that is

[TABLE]

By comparing the pole structure (40) to the dispersion relation obtained from the conservation law (14), we find the following relations

[TABLE]

From this, one finds

[TABLE]

Hence, when the small $\omega$ and $k_{y}$ limits of the function (37) are taken, we find

[TABLE]

which is the Kubo relation for the shear viscosity. The hydrodynamic limits also enable us to find

[TABLE]

where we have used $P=g_{T}(0,0)$.

In the relation (45), there appear the constants $A_{R}(0,0)$ and $D_{I}(0,0)$, which we are not able to identify within this approach. However, other studies involving a slightly different perturbing Hamiltonian were able to identify the second term in (45) as a thermodynamic quantity. The Kubo relations for the second-order hydrodynamics coefficients were examined in Baier:2007ix ; Moore:2010bu ; Moore:2012tc , where they are provided by studying the response of a fluid to small and smooth metric perturbations. If one takes into account only linearized equations, the following relations are found:111In the formulas (46)–(49) we applied the same sign convention as used in the entire paper, that is, with the metric being mostly negative. In the original papers Baier:2007ix ; Moore:2010bu ; Moore:2012tc these formulas are given with the opposite sign convention since they are studied in the flat space which is convenient when one examines transport properties of a medium via background geometry perturbations.

[TABLE]

where $\kappa$ is an additional coefficient. Note that these formulas involve non-zero $k_{z}$ while our formulas involve $k_{y}$. When expanded around small $\omega$ and $k_{z}=0$, the correlation function $\bar{G}_{R}^{xy,xy}(\omega,k_{z})$ becomes

[TABLE]

Since the two-point functions (37) and (49) have different momentum arguments, their analytical structures are slightly different. See Appendix B for details. As pointed out in Baier:2007ix , the $\kappa$ term in the dissipative part of the stress-momentum tensor is proportional to $u^{\mu}$ so that these additional terms in the correlation function do not come from the contact term in the coordinate space. Nevertheless, providing $\bar{G}_{R}^{xy,xy}(\omega,k_{z})$ and $\bar{G}_{R}^{xy,xy}(\omega,k_{y})$ share the same diffusion pole structure, their small frequency and vanishing momentum limits should be consistent with each other.

In the small $\omega$ and vanishing $k_{y}$ limits, the function (37) is

[TABLE]

By comparing the function (IV.1) to (49) we can obtain the condition on the unknown functions $A_{R}(0,0)$ and $D_{I}(0,0)$ and clarify the relation (45). So we identify

[TABLE]

Then the condition on the contribution from $D_{I}(0,0)$ and $A_{R}(0,0)$ is

[TABLE]

Finally, we can write the Kubo relation for $\eta\tau_{\pi}$ as

[TABLE]

It is known Romatschke:2009ng ; Moore:2012tc that $\kappa=\mathcal{O}(\lambda^{0}T^{2})$ in the weak coupling limit. On the other hand, both $\eta$ and $\tau_{\pi}$ behave like the mean free path which depends inversely on the cross section. Hence, Eq. (53) can be still used to study the leading order shear relaxation time.

IV.2 Response function to longitudinal fluctuations

The perturbing Hamiltonian for the bulk flow is

[TABLE]

where $\hat{T}^{00}$ is the operator of energy density and $\beta_{0}({\bf x})$ is an space-dependent external force, which has driven the system off equilibrium. The response of the medium then is

[TABLE]

By applying the continuity equation to each index of $\bar{G}_{R}^{\alpha\beta,\mu\nu}$ in the Ward identity $(\ref{Ward-id-real})$, we get

[TABLE]

where $\bar{G}_{L}(\omega,{\bf k})$ is the response function to the longitudinal fluctuations and, through the Ward identity, it is related to the spatial stress-stress function $\bar{G}_{R}^{ij,mn}(\omega,{\bf k})$ as

[TABLE]

What is more, the Ward identity enables one to express a response function associated with an arbitrary energy-momentum tensor component via the spatial stress-stress response function $\bar{G}_{R}^{ij,mn}(\omega,{\bf k})$.

From the response function $\bar{G}_{L}$ (56), one finds

[TABLE]

We know that $\bar{G}_{R}^{00,00}(0,{\bf k})=Tc_{v}+O(k^{2})$ where $c_{v}$ is the specific heat per unit volume from Eq. (214) in Appendix A. It is also related to the speed of sound $Tc_{v}=(\epsilon+P)/v_{s}^{2}$. We also take into account that the imaginary part of $\bar{G}_{L}$ must be an odd function of $\omega$ and the full function should behave well in the ${\bf k}\to 0$ limit. All these arguments allow one to parametrize the most general form of the function $\bar{G}_{L}$ as Jeon:2015dfa 222 In Jeon:2015dfa , the numerator had $i\omega^{3}Q(\omega,{\bf k})$ instead of $\omega^{2}Q(\omega,{\bf k})$. Equation (59) is the correct form for the most general parametrization.

[TABLE]

The functions $Z(\omega,{\bf k})$, $R(\omega,{\bf k})$, and $Q(\omega,{\bf k})$ are all of the form

[TABLE]

where $Z_{R}(\omega,{\bf k})$ and $Z_{I}(\omega,{\bf k})$ are real-valued even functions of $\omega$ and ${\bf k}$. The real parts $Z_{R}$ and $R_{R}$ must have a non-zero limit as $\omega\rightarrow 0$ and ${\bf k}\rightarrow 0$. All other parts of $Z$, $Q$, and $R$ must have finite limits as $\omega\rightarrow 0$ and ${\bf k}\rightarrow 0$.

The pole structure of the correlation function (59) provides us with the dispersion relation

[TABLE]

Comparing the dispersion relation (61) to (II) one observes that in order to reproduce terms $\sim\omega^{2}{\bf k}^{2}$ and other terms of higher powers, it is enough to expand the real part of $Z(\omega,{\bf k})$ up to $\omega^{2}$ so that

[TABLE]

Then, the expression (61) takes the form

[TABLE]

Comparing this pole structure to the dispersion relation provided by a purely hydrodynamic framework (II) one finds the following relations

[TABLE]

The imaginary and real parts of the Green function $\bar{G}_{L}(\omega,{\bf k})$ in the thermodynamic limit becomes

[TABLE]

where all the constants $Z_{I}$, $Z_{R1}$, $Z_{R2}$, $R_{I}$, $R_{R}$, and $Q_{R}$ should be understood as $Z_{I}\equiv Z_{I}(0,0)$, etc. Using the relations (64)–(68), we find the following Kubo formulas

[TABLE]

Here as in the shear case, one constant, $Q_{R}$, is left undetermined through this analysis. In Hong:2010at , it was shown through the curved metric analysis that

[TABLE]

Hence, we may identify

[TABLE]

More convenient Kubo formulas can be obtained if one combines the Kubo formulas (71) and (72) with those for the shear viscosity and the relaxation time so that $\eta$, $\tau_{\pi}$, and $\kappa$ are eliminated. As shown in Appendix B, the ${\bf k}\to 0$ limit of the pressure-pressure correlation function accomplishes just that,

[TABLE]

Here the pressure operator is defined as $\hat{P}=\delta_{ij}\hat{T}^{ij}/3$. Hence our final Kubo formulas for the bulk viscosity and the relaxation time are

[TABLE]

The formula (78) is especially important, as it consists of the bulk relaxation time and encodes the prescription on how to compute it.

While the Kubo formulas for the linear second-order viscous hydrodynamics have been consistently derived here, it is worth mentioning that the stress-energy correlation functions were also examined in Young:2013fka . In the current paper we focus on derivation of the second-order fluid dynamics from the general analytic properties of the correlation functions while in Ref. Young:2013fka the correlation functions were obtained as the Green functions of the Israel-Stewart type second-order hydrodynamics.

V Shear relaxation time in the scalar field theory

We perform here the perturbative analysis of the stress-energy tensor response functions. The study is done in the leading order for the massless real scalar quantum field theory333 Our analysis can be easily generalized to the massive case by just substituting $m_{\rm th}^{2}\to m_{\rm phys}^{2}+m_{\rm th}^{2}$ where $m_{\rm phys}$ is the physical mass.

with the Lagrangian

[TABLE]

where $\lambda$ is the coupling constant, which is assumed to be small. Since the scalar field is real, there are no conserved number or charge operators coupled to the chemical potential.

In leading order the scalar field dynamics is governed by $2\leftrightarrow 2$ scatterings which give rise only to the shear viscosity effects being of the order $\mathcal{O}(\lambda^{-2}T^{3})$ where $T$ is the temperature Jeon:1995zm . The bulk viscosity strictly vanishes in the conformal limit. Since the conformal symmetry in the scalar theory is broken by the nonzero $\beta$ function, the bulk viscosity is much smaller than the shear viscosity $\zeta=\mathcal{O}(\lambda T^{3})$ Jeon:1995zm . It also requires inclusion of number changing inelastic processes at higher orders in the coupling constant.

Here we work in the leading order of the expansion of the response function and provide a systematic analysis to compute the shear relaxation time. To find it we make use of the formula (53), where we ignore the coefficient $\kappa$ since it only scales as $\mathcal{O}(\lambda^{0}T^{2})$. By evaluating the real and imaginary parts of $\bar{G}_{R}^{xy,xy}$ we are able to get $\eta\tau_{\pi}$ and $\eta$, respectively, and then extract the shear relaxation time. In the forthcoming derivation we employ the closed time path (Keldysh-Schwinger) formalism which is briefly summarized in Appendix C. The analysis of the four-point functions in Sec. V.5 closely follows Wang:2002nba . We also adopt the notations and sign conventions for the real-time $n$-point functions from the same reference throughout this section.

V.1 Definition of the retarded response function

Since any response function $\bar{G}_{R}$ differs from the standard retarded Green function due to the Ward identity let us start with the definition of $\bar{G}_{R}^{ij,mn}$. Making allowance for the Ward identity (29), $\bar{G}_{R}^{ij,mn}$ is defined by

[TABLE]

In the equilibrium state, the rotational invariance provides that $\langle\hat{T}^{ij}\rangle=\delta^{ij}P$, where $P$ is the thermodynamic pressure, so that the retarded Green function becomes

[TABLE]

Analogously, the advanced Green function is

[TABLE]

The commutator in Eqs. (V.1) and (V.1) consists of the Wightman functions and should be understood as

[TABLE]

where the reading from right to left inside the brackets should be understood as the evolution from the initial to final state; the indices 1 and 2 locate the operators on the upper (earlier) or lower (later) branch of the Keldysh time contour, respectively. The stress tensor operator is defined by

[TABLE]

with the Lagrangian given by (79). Rewriting the Lagrangian as

[TABLE]

where

[TABLE]

is the equation of motion for the field operator, we see that within the thermal average, $\langle{\cal L}\rangle=O(\lambda)$. Hence the leading order stress-energy tensor is dominated by the kinetic term

[TABLE]

Accordingly, the Wightman functions read

[TABLE]

which may be expressed in terms of the four-point Green functions as

[TABLE]

The four-point Green functions are defined as

[TABLE]

Then the retarded and advanced Green functions of the stress-energy tensors take the forms

[TABLE]

[TABLE]

where the limits $x_{1}\to x_{2}$ and $y_{1}\to y_{2}$ must be taken after the derivatives.

V.2 Real part of the retarded Green function in the free field theory

To find the real part of $\bar{G}^{ij,mn}$ we continue our considerations to study the one-loop diagram which appears in the free scalar quantum field theory; this is done mostly to work out the details.

This diagram is shown in Fig. 1 and it comes from the disconnected parts of four-point Green functions which are nothing but the products of the Wightman functions,

[TABLE]

When the points are joined so that $x_{1}=x_{2}=x$ and $y_{1}=y_{2}=y$, we get

[TABLE]

where we defined $\Delta_{12}(x,y)=\langle\phi_{1}(x)\phi_{2}(y)\rangle$ and $\Delta_{21}(x,y)=\langle\phi_{2}(x)\phi_{1}(y)\rangle$ also used the fact that $\Delta_{12}(x,y)=\Delta_{21}(y,x)$. Since the system considered is homogeneous, the two-point functions depend on $x$ and $y$ only through their difference $x-y$. Putting all these facts together, we get the retarded function $\bar{G}_{R}^{ij,mn}$ in the form

[TABLE]

At this point, it is more convenient to change the basis from $\phi_{1,2}$ to $\phi_{r,a}$ defined by

[TABLE]

Using the relations between the Green functions in the (1,2) and $(r,a)$ bases expressed by (229)–(232) and (234)–(237), we find

[TABLE]

Performing the Fourier transform we next get

[TABLE]

where $k\equiv(k_{0},{\bf k})\equiv(\omega,{\bf k})$. And analogously $\bar{G}^{ij,mn}_{A}(k)$ is

[TABLE]

The functions $\Delta_{ra}(p)$ and $\Delta_{ar}(p)$ are the usual retarded and advanced two-point Green functions which are of the following forms:

[TABLE]

and they satisfy

[TABLE]

$\Delta_{rr}(p)$ is the autocorrelation function and it is the only function where a distribution function explicitly enters. In thermal equilibrium all three of these functions are related via the fluctuation-dissipation theorem,

[TABLE]

where $N(p^{0})=1+2n(p_{0})$ and $n(p^{0})=1/(e^{\beta p^{0}}-1)$ is the Bose distribution function with $\beta$ being the inverse of temperature $T$.

The functions $\bar{G}^{xy,xy}_{R}(k)$ and $\bar{G}^{xy,xy}_{A}(k)$ are then obtained by setting $i=m=x$ and $j=n=y$ in Eqs. (V.2) and (V.2). We also choose ${\bf k}=(0,k_{y},0)$ to use our analysis in Sec. IV.1. The real and imaginary parts of $\bar{G}^{xy,xy}_{R}(k)$ are then obtained by the sum and the difference of $\bar{G}^{xy,xy}_{R}(k)$ and $\bar{G}^{xy,xy}_{A}(k)$, respectively. In the vanishing momentum limit they are as follows:

[TABLE]

where in the right hand side $k=(\omega,\mathbf{0})$.

V.3 Pinching poles

Provided that the fluctuation-dissipation theorem of the form (109) is applied to the formulas (110) and (111), there appear, in particular, terms of products of propagators with four poles lying symmetrically on both sides of the real axis in the complex $p_{0}$ plane if the small $\omega$ limit is used. These poles give rise to the pinching of the integration contour by the poles lying on opposite sides of the real energy axis, that is, the pinching pole effect. These terms indeed may be evaluated as

[TABLE]

so that they produce a singularity as $\epsilon\rightarrow 0$. In a non-interacting theory such an effect is natural and it means that since the emerged excitation is not subject to collisions it can propagate indefinitely long. It is reflected by $\delta$ functions carried by the spectral density. The width of such a peak, which is inversely proportional to the lifetime of the excitation, is vanishingly small. This implies that in the free theory there is no transport of conserved quantities and consequently transport coefficients cannot be defined. In an interacting system, transport coefficients are finite due to the finite mean free path (or lifetime) of a propagating excitation until it suffers from scatterings with constituents of the thermal bath. Thus, in thermal weakly interacting medium the spectral density can be approximated by Lorentzians Jeon:1994if

[TABLE]

where $E_{p}$ is the quasi-particle excitation energy and $\Gamma_{p}$ is the thermal width. The origin of such a form of the spectral density may be also understood if one uses the resummed propagators. These propagators carry information on the interaction of a given particle with the medium in terms of the self-energy $\Sigma=\textrm{Re}\,\Sigma+i\textrm{Im}\,\Sigma$. They are defined as $\Delta_{ra}(p)=[p^{2}-m^{2}-\Sigma(p)]^{-1}$ and $\Delta_{ar}(p)=\Delta^{*}_{ra}(p)$, where $m$ is the mass.

In this paper, we study massless theory, but the real part of the self-energy in the lowest order does not vanish. The leading order diagram is the tadpole diagram, which is momentum independent, and may be identified as the thermal mass ($m_{\rm th}$) squared. The spectral density is then given in terms of the resummed retarded and advanced propagators as

[TABLE]

When the spectral density has sharp peaks near $p_{0}=\pm E_{p}$, we can then say that the dispersion relation is $E_{p}^{2}={\bf p}^{2}+m_{\rm th}^{2}$, and the thermal width is related to the imaginary part of the self-energy as $\Gamma_{p}=\frac{\textrm{Im}\,\Sigma(E_{p},|{\bf p}|)}{2E_{p}}$. Consequently, the resummed retarded and advanced propagators can be approximated as

[TABLE]

Appearance of $\Gamma_{p}$ in the propagators shifts the pinching poles away from the real axis in the complex $p_{0}$ plane, which regulates the singularities in (110) and (111) making the integral finite. The pinching pole contribution is then of the order $\mathcal{O}(1/\Gamma_{p})$. What is more, the terms of the type $\Delta_{ra}(p)\Delta_{ra}(p)$ and $\Delta_{ar}(p)\Delta_{ar}(p)$ have poles on the same side of the real energy axis and thus they give much smaller contribution to the expressions (110) and (111) than the pinching poles, and may be safely ignored in further computations. The omission of these terms constitutes the pinching pole approximation.

Replacement of bare propagators by dressed ones means that we need to deal with the skeleton expansion where propagators are dressed and vertices remain bare. Here, we are to study the first loop of this expansion. However, since the thermal width is related to the imaginary part of a self-energy, some complications arise. In the weakly coupled $\lambda\phi^{4}$ theory the lowest contribution to $\textrm{Im}\,\Sigma$ comes from a two-loop diagram which is of the order $\mathcal{O}(\lambda^{2})$ and since the pinching pole contribution dominates, the one-loop diagram is of the order $\mathcal{O}(1/\lambda^{2})$ Jeon:1992kk . However, one realizes that there may be momentum exchange between the side rails of the loop. This is represented by the one-loop rungs connecting the two side rails as shown in Fig. 2.

Each rung introduces a factor of $\lambda^{2}$ coming from the vertices and a factor of the order $\mathcal{O}(1/\lambda^{2})$ coming from the pinching poles introduced by the additional pair of propagators. Therefore, all such multi-loop ladder diagrams contribute at the leading order. They must be resummed to give the full result in the leading order.

The situation described above holds when the single transport coefficient, such as the shear viscosity, is analyzed. In case of the combination $\eta\tau_{\pi}$, it gets more involved and it will be discussed in the next part of this work.

V.4 Evaluation of $\eta$ and $\eta\tau_{\pi}$ in the one-loop limit

Before we include all ladder diagrams let us consider first only the one-loop diagram with the resummed propagators. This is illuminating as we can find the typical scales of $\eta$ and $\eta\tau_{\pi}$.

The shear viscosity is related to the imaginary part of the relevant retarded Green function. It can be calculated in a few ways which are related to different choices of the correlation function. Usually, it is examined from the imaginary part of the Green function of the traceless, spatial part of the stress-energy tensor $\pi^{ij}$. It may be also computed in terms of the stress-stress function $\bar{G}_{R}^{xy,xy}$ as shown by the Kubo formula (44). We employ here the latter choice. Applying the fluctuation-dissipation theorem (109) and the pinching pole approximation to Eq. (111), the imaginary part is then given by

[TABLE]

with the propagators $\Delta_{ra}$ and $\Delta_{ar}$ being dressed and given by (115) and (116) and $k=(\omega,\mathbf{0})$ in the right hand side. The real part of the retarded Green function (110) is

[TABLE]

In pursuit of $\eta=\frac{1}{\omega}\text{Im}\bar{G}^{xy,xy}_{R}(\omega,k_{y})|_{\omega,k_{y}\to 0}=\partial_{\omega}\text{Im}\bar{G}^{xy,xy}_{R}(\omega,k_{y})|_{\omega,k_{y}\to 0}$, we find

[TABLE]

Realizing that

[TABLE]

the one-loop shear viscosity is

[TABLE]

The action of the second-order derivative and inclusion of the factor $-1/2$ to the real part given by the formula (118), as dictated by the Kubo formula (53), leads us to the following equation:

[TABLE]

Taking Eq. (120) into account and

[TABLE]

the expression (V.4) becomes

[TABLE]

The frequency integrals to perform in (121) and (124) are

[TABLE]

[TABLE]

The integrands in Eqs. (125) and (126) have four poles at $p_{1}=i\Gamma_{p}+E_{p}$, $p_{2}=i\Gamma_{p}-E_{p}$, $p_{3}=-i\Gamma_{p}+E_{p}$, and $p_{4}=-i\Gamma_{p}-E_{p}$. In Eq. (125) the poles are simple poles while in Eq. (126) they are double poles. By using the residue theorem and closing the contour in the upper-half plane, the sum of the residua in Eq. (125) is found in the leading order of $\Gamma_{p}/E_{p}$ as

[TABLE]

In Eq. (126) we handle the second-order poles. We recall that the residua of a function with second-order poles contain the derivative with respect to the complex argument. Thus, when the contour integration is carried out, the expression (126) becomes

[TABLE]

Finally, the formula (121) for the shear viscosity is

[TABLE]

and the product of the shear viscosity and its relaxation time is

[TABLE]

At this level one immediately notices that the shear relaxation time scales as $1/\Gamma_{p}$. The computation of its value is given in Sec. VI.

V.5 Summation over multiloop diagrams

The one-loop limit is, however, not sufficient and we need to resum ladder diagrams which requires us to manipulate the connected 4-point Green functions as well. To do so we employ the definitions (C) and (C) to get the retarded and advanced four-point Green functions as

[TABLE]

[TABLE]

The subscripts $R$ and $A$ in $\tilde{G}_{R}$ and $\tilde{G}_{A}$ do not mean that they themselves are four-point retarded and advanced functions. The subscripts just indicates that these functions will become the two-point retarded and advanced functions when $x_{1}$ is identified with $x_{2}$ and $y_{1}$ is identified with $y_{2}$. The real part is

[TABLE]

where we have used the relation (226), and the imaginary part is

[TABLE]

The four-point functions in $(1,2)$ basis may be transformed to the $(r,a)$ basis using the relation (238). Then, one finds

[TABLE]

Despite the fact that the shear viscosity, related to the imaginary part of a Green function, has been studied in literature many times, it may be illuminating to see some analogies and differences between the real and imaginary parts. Thus, when deriving the real part of the Green function we will be referring to $\textrm{Im}\tilde{G}_{R}$ quite frequently, as well. In particular, we will be quoting the results from Wang:2002nba , where $\eta$ was derived in the real-time formalism.

Each of the four-point Green functions in Eqs. (135) and (136) couples to the remaining ones. This is shown in Fig. 3 where the disconnected parts represent the one-loop diagram, the rung is a kernel and the shaded box is a sum of all possible combinations of four-point Green functions in the ($r,a$) basis. Both the kernel and the shaded box contribute to an effective vertex. Therefore, Fig. 3 presents an infinite series of diagrams which can be written as the Bethe-Salpeter equation (BSE). For an arbitrary function $G_{\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}}$, the BSE reads

[TABLE]

The analytic solution of the BSE is, in general, not readily available. However, as already mentioned, the use of the Keldysh basis accompanied by the pinching pole approximation makes it much simpler.

In Wang:2002nba it has been shown that by means of the fluctuation-dissipation theorems (FDT), developed in Wang:1998wg ; Carrington:2006xj , one may show that the only contribution to $\textrm{Im}\tilde{G}_{R}$ comes from $\textrm{Im}G_{aarr}$. Even without referring to FDT one may quickly observe that the functions $G_{aaar}$, $G_{aara}$, $G_{araa}$, and $G_{raaa}$ in Eqs. (135) and (136) do not contribute because these functions must contain at least one $\Delta_{aa}$.

In this way we are left only with $G_{rrra}$, $G_{rrar}$, $G_{rarr}$, and $G_{arrr}$. As shown in Appendix D, these functions are related to $G_{aarr}$ and $G_{rraa}$ by

[TABLE]

where $N_{p}=N(p_{0})$. The functions $G_{aarr}$ and $G_{rraa}$ are also related through FDT. The relevant relation is Carrington:2006xj

[TABLE]

So, ignoring again the contributions from the functions with three $a$ indices, we obtain $\textrm{Re}\tilde{G}_{R}$ and $\textrm{Im}\tilde{G}_{R}$ as follows:

[TABLE]

V.6 Bethe-Salpeter equation for $G_{aarr}$ and $G^{*}_{aarr}$

The Bethe-Salpeter equation for $G_{aarr}$ is analyzed in detail in Wang:2002nba . Nonetheless, we repeat here the main steps because unlike the viscosity calculations, the $\eta\tau_{\pi}$ calculation requires that the external frequency $k_{0}=\omega$ be kept until the derivatives are taken. Fully written out, the Bethe-Salpeter equation for $G_{aarr}$ is

[TABLE]

The external momentum $k$ can flow in an arbitrary way along a diagram, that is, it can enter the rungs but it does not have to. Since all opportunities are equally possible and lead to the same final result, we have a freedom to choose what is convenient for computations. Therefore, in this analysis the external momentum is flowing along the external lower side rail only, as indicated in Fig. 3. Then, the kernel couples to four-point Green functions as

[TABLE]

where we have taken into account that $K_{rrrr}=0$, which is the analog of $G_{aaaa}=0$ amputated of the external legs. By truncating external legs from $G_{rarr}$ and $G_{arrr}$, that is, by using the formulas (D) and (D), we find that the expression (144) becomes

[TABLE]

The Bethe-Salpeter equation for $G_{aarr}$ becomes now

[TABLE]

To find $G^{*}_{aarr}$ one just makes a complex conjugate of the formula (146). Taking into account the analysis of the complex conjugate procedure of the kernel rungs, shown in Appendix D, one gets

[TABLE]

Note that the combinations of the kernel functions in (146) and (147) become identical in the $k\to 0$ limit.

V.7 Evaluation of $\eta$ and $\eta\tau_{\pi}$

To calculate the shear viscosity and the shear relaxation time, one needs to evaluate Eq. (142) with the appropriate derivatives

[TABLE]

and Eq. (141),

[TABLE]

Accordingly, only $G_{aarr}$ and its complex conjugate matter when the transport coefficients are needed.

To obtain $\eta$ and $\eta\tau_{\pi}$ we move forward to find $\partial_{\omega}\text{Im}\bar{G}^{xy,xy}_{R}$ and $\partial_{\omega}^{2}\textrm{Re}\bar{G}^{xy,xy}_{R}$. Simultaneously, we will be applying the limit $\omega\rightarrow 0$. So, we get

[TABLE]

The first line in (151), however, gives vanishing contribution to the calculation of $\eta\tau_{\pi}$ since it is an odd function of $p_{0}$. It is seen when the following arguments are taken into account: $n(-p_{0})=-n(p_{0})-1$ and $\Delta_{ar}(-p_{0})=\Delta_{ra}(p_{0})$. Including the relation (120), Eqs. (150) and (151) are rewritten as

[TABLE]

where $\text{Im}G_{aarr}$ and $\partial_{\omega}\text{Re}G_{aarr}$ are to be found from the Bethe-Salpeter equations for $G_{aarr}$ and $G^{*}_{aarr}$. First, using the BSEs given by (146) and (147), $\text{Im}G_{aarr}$ and $\text{Re}G_{aarr}$ are given by

[TABLE]

[TABLE]

where, for clarity, we have introduced the following symbolic notations

[TABLE]

and

[TABLE]

The imaginary part $\text{Im}G_{aarr}$ has been expressed in the vanishing $k$ limit. It is also the case for the kernel $\mathcal{K}_{pl}$ of the real part of $G_{aarr}$. This is justified due to the following reason. In the next step we need to apply the derivative with respect to frequency. There would appear terms consisting of $\partial_{\omega}\mathcal{K}_{pl}(k)$ which are, however, of the order of $1/\Gamma_{p}$ or less and they give much smaller contribution to the final formula than the remaining ones, which are of the order $1/\Gamma_{p}^{2}$. Therefore, the hydrodynamic limits could have been applied at this stage, which has simplified a lot the notation of formula (155). The action of the derivative on Eq. (155) in the hydrodynamic limits produces

[TABLE]

To assess which terms contribute further let us check the small frequency behavior of $A_{p}$ and $B_{p}$ functions given by (159) and (160) and their derivatives,

[TABLE]

Due to Eq. (163), the second term in Eq. (V.7) vanishes. To see the behavior of the third term we need to include the integrals as given by the formula (V.7). Then the following expression needs to be considered:

[TABLE]

where $I_{p}=I(p)=p_{x}p_{y}$ and $F_{p}(0)$ is

[TABLE]

The frequency integral part of $I_{3}$ is

[TABLE]

and

[TABLE]

The function $B^{\prime}_{p}(0)$ has four poles at $p_{1}=i\Gamma_{p}+E_{p}$, $p_{2}=i\Gamma_{p}-E_{p}$, $p_{3}=-i\Gamma_{p}+E_{p}$, and $p_{4}=-i\Gamma_{p}-E_{p}$ and they are of the second order. We calculate them by using the residue theorem and closing the contour in the upper-half plane. We recall that the residua of a function with second-order poles contain the derivative with respect to the complex argument. Thus, upon carrying out the contour integration and summing up the residua, the expression (169) becomes

[TABLE]

If we now group the terms in Eq. (V.7) according to the derivative with respect to $p_{0}$, we can write

[TABLE]

The expressions in the first and the second lines of Eq. (172) are equal to 0. The other terms produce

[TABLE]

which is of the order of $1/\Gamma_{p}$. The derivative $\partial_{p_{0}}$ in (173) acts only on the kernel $K_{pl}$ included in $F_{p}(0)$ as shown by (168). The kernel, however, does not contribute any terms of the order of $1/\Gamma$ so neither can do the action of $\partial_{p_{0}}$. The contribution of $F_{p}(0)$ is $\mathcal{K}_{pl}(0)\,\textrm{Re}G_{lq}(0)\sim\mathcal{O}(1)$. Consequently, the third term in Eq. (V.7) is only $O(1/\Gamma_{p})$ and may be ignored as we know from the one-loop analysis that $\eta\tau_{\pi}=O(1/\Gamma_{p}^{2})$.

In this way we are left with the first and the fourth term of Eq. (V.7), that is

[TABLE]

From now on the vanishing $\omega$ and $k_{y}$ limits are implicit in all expressions.

To reduce the BSE to a more manageable form, it is convenient to define the effective vertex

[TABLE]

where again $I_{q}=I(q)=q_{x}q_{y}$. The effective vertex $D_{p}$ satisfies the following integral equation

[TABLE]

which is schematically shown in Fig. 4.

Then formula for $\eta$ becomes

[TABLE]

For the real part, including the integral over $q$ to the equation (V.7) and the expression for $D_{p}$, given by (176) we can define a new effective vertex $R_{p}$,

[TABLE]

which satisfies the following integral equation

[TABLE]

and the formula for $\eta\tau_{\pi}$ can be written in a compact form as

[TABLE]

Note that $R_{p}=O(1/\Gamma_{p})$.

If we insert $D_{p}$, given by (176), and $B_{p}$, given by Eq. (164), into Eq. (177) and perform the contour integration, then the shear viscosity is given by

[TABLE]

where the effective vertex satisfies

[TABLE]

When deriving this we have used $D(E_{l},{\bf l})=D(-E_{l},{\bf l})$, $\mathcal{K}(E_{p},E_{l})=\mathcal{K}(-E_{p},-E_{l})$, and $\mathcal{K}(-E_{p},E_{l})=\mathcal{K}(E_{p},-E_{l})$. The relations for the kernel are shown in Appendix D.

Then, using $R_{p}$, given by (179), and $B_{p}$ and $A^{\prime}_{p}$, given by (164) and (165), respectively, to Eq. (180), we get

[TABLE]

with the effective vertex $R(E_{p},{\bf p})$ satisfying

[TABLE]

VI Shear relaxation time

First of all, by comparing the one-loop results (130) and (129), one clearly sees that both the shear viscosity and the shear relaxation time are controlled by the thermal width, $\Gamma_{p}$. The thermal width is defined by the imaginary part of the self-energy and therefore is momentum dependent.

Let us first, however, consider a simple example by assuming that the thermal width is constant. Then by comparing the formulas (130) and (129) in the one-loop limit, it is found as

[TABLE]

Thus one can claim that the thermal width $\Gamma$, as directly related to the lifetime or the mean free path of a thermal excitation, introduces the only time scale into the system of interacting particles and the shear relaxation time is directly related to that scale. Hence the ratio of the one-loop results is

[TABLE]

With the Bose-Einstein momentum distribution $n(E_{p})=1/(e^{\beta E_{p}}-1)$, we get the value

[TABLE]

With the definition of the energy density for the one-component field,

[TABLE]

where $n(E_{p})$ is the Bose-Einstein statistics and the pressure is $\langle P\rangle=\frac{1}{3}\langle\epsilon\rangle$, the formula (187) may be rewritten in the form

[TABLE]

The similar relation was found within 14-moment approximation to the Boltzmann equation for the classical massless gas, studied in Ref. Denicol:2014vaa , but the energy density and thermodynamic pressure are defined there through the Boltzmann statistics. If we use the Boltzmann statistics in Eqs. (186) and (188) then the relation (189) holds approximately.

In general, one needs to maintain the momentum dependence of the thermal width and solve the integral equations numerically. By analyzing the integral equations for the effective vertices $D(E_{p},{\bf p})$ and $R(E_{p},{\bf p})$ given by Eqs. (182) and (184), one sees that they are of the same type but the inhomogeneous terms are different. This indicates that the order of these two equations is different since $D(E_{p},{\bf p})\sim\mathcal{O}(1)$ and $R(E_{p},{\bf p})\sim\mathcal{O}(1/\Gamma_{p})$. Also, to find the solution to $R(E_{p},{\bf p})$ one needs to first obtain the solution to $D(E_{p},{\bf p})$. The complication of solving double integral equations can be, however, avoided. Define the inner product of two functions as

[TABLE]

Then the viscosity is

[TABLE]

The integral equation for $D_{p}$ can be symbolically written as

[TABLE]

where

[TABLE]

and for $R_{p}$,

[TABLE]

Formally, the solutions are

[TABLE]

The product $\eta\tau_{\pi}$ is then

[TABLE]

provided that the kernel operator ${\cal C}$ is real and symmetric. The same formula was found in York:2008rr through the effective kinetic theory approach. The fact that ${\cal C}$ is real and symmetric is shown in Appendix D.

The form of the effective vertices actually reflects the fact that when the ladder diagrams are summed over to get $\eta\tau_{\pi}$, one out of all pairs of propagators in each diagram contributes one more factor of $1/2\Gamma_{p}$ when compared to equivalent resummation of diagrams corresponding to $\eta$ calculation. Since every loop may be cut so as to represent an elastic scattering process, one can say that each of these processes can contribute the lifetime associated with the momentum of incoming or outgoing particles, or also that of mediating particles, when loops with at least two rungs are considered. In the end the shear relaxation time is obtained when all distinguishable possibilities are included and they all give rise to a balanced relaxation process.

The shear relaxation time has been then obtained by evaluating the integrals in (181) and (183) numerically. The result is shown in Fig. 5, where the relaxation time is displayed as a function of $m_{\textrm{th}}/T$ with $m_{\textrm{th}}$ being the thermal mass, which introduces the natural cut off to infrared divergences. The ratio $m_{\textrm{th}}/T$ ought to be identified with $\sqrt{\lambda}$ since the thermal mass behaves as $m^{2}_{\textrm{th}}=\lambda T^{2}$ in the leading order and the range of the plot has been chosen in such a way to show a general behavior of the shear relaxation time. Yet, one needs to keep in mind that when the coupling constant is increasing, nonperturbative effects start to play more and more a role and the value of the relaxation time becomes less and less realistic. On the other hand, our analysis equally well applies even if the physical mass is non-zero. Hence, the large $m_{\rm th}/T$ part of the results shown in this section can be interpreted as those for the the massive scalar field but still weakly coupled.

The solid (blue) curve in Fig. 5 corresponds to the shear relaxation time obtained as the ratio of (183) over (181). The dashed (red) curve represents the shear relaxation time evaluated directly from one-loop expressions (130) and (129). As can be seen, the difference between these two treatments is noticeable. One also immediately observes that the bigger the coupling constant is, the shorter the shear relaxation time becomes, as expected. What is more, the shear relaxation time is around three times bigger than the corresponding one-loop finding so it is not justified to claim that the one-loop result dominates the behavior of the relaxation time.

We also present the ratio $\langle\epsilon+P\rangle\tau_{\pi}/\eta$ as a function of $m_{\textrm{th}}/T$, which is shown in Fig. 6. One can see that the ratio is decreasing when coupling constant decreases and it varies between 6.11 up to 6.55 in the range shown. Furthermore, the tendency of the one-loop result is just opposite.

The ratio $\langle\epsilon+P\rangle\tau_{\pi}/\eta$ as a function of $m_{\textrm{th}}/T$ was also obtained in York:2008rr within the effective kinetic theory for QCD, QED, and scalar $\lambda\phi^{4}$ theory. In the case of the scalar theory only the value at $m_{\textrm{th}}=0$ was presented in York:2008rr . As $m_{\textrm{th}}/T\to 0$, they observed this ratio to be 6.11 for the full theory and exactly 6 for the massless Boltzmann gas. As can be seen, the quantum field theoretical findings obtained here are in line with these of the effective kinetic theory in the regime of a very small coupling constant and one may expect that this equivalence holds in the whole range of $m_{\textrm{th}}/T$.

The ratio $\langle\epsilon+P\rangle\tau_{\pi}/\eta$ obtained here for the scalar theory behaves in opposite way than the one obtained in York:2008rr for QCD in the range of $m_{\textrm{th}}/T$ changing from 0 to 1. When the coupling constant increases the nonperturbative effects become, as mentioned, more and more essential but even for very small values of the coupling different nature of the interaction governing both theories should be taken into consideration. In QCD Coulomb-like interaction dominates the physics of the system which results in appearance of $\textrm{ln}(1/g)$, with $g$ being the coupling constant of strong interaction, in the parametric form of the transport coefficients. What is more, the cross section is strongly angle dependent which gives rise to the soft and near collinear singularities. As argued in York:2008rr the collinear splittings become more and more important when the coupling constant increases and they make the ratio $\langle\epsilon+P\rangle\tau_{\pi}/\eta$ go down. In the scalar $\lambda\phi^{4}$ theory the contact interaction causes the cross-section to be isotropic therefore only ladder diagrams constitute the leading order behavior of the shear viscosity and its relaxation time. In Ref. Jeon:1994if it was shown that indeed the soft and near collinear singularities are of smaller size than the ladder diagrams and they have not had to be discussed in detail here. A broader analysis of the next-to-leading order behavior of shear viscosity is given in Ref. Moore:2007ib , where the soft physics (quasiparticle momenta $k$ being of the order $m_{\textrm{th}}$) is shown to determine the NLO correction. Since the leading order equation does not differ in the form from the one that includes subleading corrections, we expect that higher order contributions will behave similarly, but we leave that for future analysis.

It is also illuminating to notice that the leading order behavior of different transport coefficients manifests different susceptibility to the thermal mass and to the way it emerges in the computational analysis. For instance, the thermal mass plays an essential role in the case of bulk viscosity where not only the leading order term but also higher order corrections in the coupling constant must be properly incorporated. Then, apart from the loops with an arbitrary number of rungs inside them, there appear so-called chain diagrams made of an arbitrary number of loops. Addition of a subsequent loop introduces an additional coupling constant. However, although chain diagrams emerge at every possible order in the coupling constant their resummation is of the size $\mathcal{O}(\lambda T^{2})$, that is, the same as the thermal mass squared. Therefore, the chain diagrams give a significant contribution to the bulk effects analysis. In the case of the shear viscosity or shear relaxation time calculation this effect is absent since rotational symmetry makes the chain diagrams vanish, as discussed in Jeon:1994if . Accordingly, these higher order corrections do not have to be considered here.

VII Conclusions

The first goal in this paper was to work out Kubo formulas for the shear and the bulk relaxation times. Since the Kubo formula for the shear relaxation time was studied before, our focus was on the Kubo formula which allows us to compute the bulk relaxation time. Our Kubo formula is different than the one obtained in Huang:2010sa using the projection operator method, but it is consistent with the one found in Hong:2010at using the metric perturbation method.

Although our ultimate goal is to compute the bulk relaxation time, in this paper we have concentrated on simpler task of computing the shear relaxation time. This is primarily to check the soundness of our overall formulations and refine the field theoretical techniques we will need for the much more involved bulk relaxation time computation. Calculation of the bulk relaxation time is in progress.

The behavior of the shear relaxation time obtained in this study is consistent with previous studies within the kinetic theory setting as well as within the field theory setting. Our conclusion that

[TABLE]

seems to be robust across different theories and also consistent with the values used in hydrodynamic calculations Vujanovic:2016anq .

Acknowledgments

A. Czajka acknowledges support from the program Mobility Plus of the Polish Ministry of Science and Higher Education. S.J. was supported in part by the Natural Sciences and Engineering Research Council of Canada. Discussions with C. Gale and G.D. Moore are very much appreciated.

Appendix A Ward identity analysis

The stress-energy tensor Ward identity (30) is

[TABLE]

In the local rest frame of the medium, the above expression can be decomposed as

[TABLE]

In the ${\bf k}\to 0$ limit, the right hand sides of Eqs. (200)–(205) must vanish. Hence,

[TABLE]

In the $\omega\to 0$ limit, the left hand side must vanish. Hence,

[TABLE]

where $\hat{\delta}_{ij}=\delta_{ij}-\hat{k}_{i}\hat{k}_{j}$ with $\hat{k}_{i}=k_{i}/|{\bf k}|$ is the transverse projection and $g_{ijlm}$ is transverse with respect to all its indices. Each of the $g_{\cdots}({\bf k})$ functions must be at least $O({\bf k}^{2})$ in the small ${\bf k}$ limit so that $\hat{k}_{i}\hat{k}_{j}$ in $\hat{\delta}_{ij}$ is well defined.

We also need $\bar{G}^{\epsilon\epsilon}(0,{\bf k})$. In the zero frequency limit, the retarded correlation function and the Euclidean one coincide. Hence,

[TABLE]

Note that

[TABLE]

which leads to

[TABLE]

where $\partial\epsilon/\partial T=c_{v}$ is the specific heat per unit volume.

Appendix B Decomposition of stress-stress correlation function

Consider $\bar{G}^{ij,lm}$ which is the correlation function of $\hat{T}^{ij}$ and $\hat{T}^{lm}$. When the space is isotropic, there are following five independent tensors that respect the symmetry

[TABLE]

which are composed of $\delta_{ij}$ and $\hat{k}_{i}=k_{i}/|\bf{k}|$. Hence, there are altogether five independent functions that can appear in $\bar{G}^{ij,lm}$. Since we are interested in the shear and the bulk responses, it is more convenient to define the transverse metric

[TABLE]

The stress-stress correlation function can be then decomposed as LGY

[TABLE]

This expression must be well defined in the ${\bf k}\to 0$ limit. Hence, the coefficient of, say, $\delta_{ij}\hat{k}_{l}\hat{k}_{m}$ must be $O({\bf k}^{2})$ as ${\bf k}\to 0$. Collecting the coefficients, one sees that

[TABLE]

The coefficient of $\delta_{il}\hat{k}_{j}\hat{k}_{m}$ is

[TABLE]

The coefficient of $\hat{k}_{i}\hat{k}_{j}\hat{k}_{l}\hat{k}_{m}$ must be $O({\bf k}^{4})$,

[TABLE]

Hence in the ${\bf k}\to 0$ limit, only two functions are independent. Choosing those to be $\bar{G}_{1}$ and $\bar{G}_{L}$, we get

[TABLE]

From Eq. (217), one can easily see that $\bar{G}_{R}^{xy,xy}(\omega,k_{y})=-P+\bar{G}_{1}(\omega,k_{y})$ and $\bar{G}_{R}^{xy,xy}(\omega,k_{z})=-P+\bar{G}_{2}(\omega,k_{z})$. From Ref. Baier:2007ix , we know that in the small $k_{z}$ limit,

[TABLE]

However, since $\bar{G}_{1}(\omega,{\bf k})=\bar{G}_{2}(\omega,{\bf k})+O({\bf k}^{2})$, one cannot in general say that $\bar{G}_{R}^{xy,xy}(\omega,k_{z})$ behaves the same as $\bar{G}_{R}^{xy,xy}(\omega,k_{y})$ nor the same as $\bar{G}_{R}^{xy,xy}(\omega,k_{x})$.

To get the Kubo formulas for the bulk viscosity and the bulk relaxation time, the pressure-pressure correlation function is needed. In the limit ${\bf k}\to 0$, Eq. (217) becomes, with $\hat{P}=\delta_{ij}\hat{T}^{ij}/3$,

[TABLE]

Appendix C Closed time path formalism

Here we briefly describe the closed time path or Keldysh-Schwinger formalism, which is studied in more detail, for example, in Chou:1984es . The main object of the formalism is the contour Green function which has four components of the real-time arguments. Here, we define them for the scalar field operators $\phi$ but these definitions may be directly generalized to any composite field operators, such as $\hat{T}^{ij}$ discussed in the main body of the paper. The components of the contour Green function are given as

[TABLE]

where $a_{1},a_{2}\in\{1,2\}$ and the indices 1 and 2 refer to the two branches of the Keldysh contour the field operator $\phi$ is located on. The operator $\mathcal{T}$ represents an ordering of the operators along the contour; $\mathcal{T}=\mathcal{T}_{c}$ chronologically orders the operators on the upper branch and $\mathcal{T}=\mathcal{T}_{a}$ sets anti-chronological ordering on the lower branch. For the angle brackets we use the following notation:

[TABLE]

with $\hat{\rho}(t_{0})$ being a density operator and the trace is understood as a summation over all states of the system at a given initial time $t_{0}$. The averaged products of unordered operators are commonly known as the Wightman functions. All the functions in the (1,2) basis satisfy the relation

[TABLE]

which reflects the fact that only three out of four components are independent of each other.

Going to the $(r,a)$ basis, we define

[TABLE]

and then the four components are defined by

[TABLE]

where $\alpha_{1},\alpha_{2}\in\{r,a\}$ and $n_{r}$ is a number of $r$ indices among $\alpha_{1}$ and $\alpha_{2}$.

For further purposes it is useful to know the relations between the Green functions of $(r,a)$ and $(1,2)$ bases, which read

[TABLE]

where $\theta(x_{0})$ is the step function. A general transformation law which holds for any $n$-point Green function is then

[TABLE]

where the repeated indices are summed over and $Q_{a1}=-Q_{a2}=Q_{r1}=Q_{r2}=\frac{1}{\sqrt{2}}$ are the four elements of the orthogonal Keldysh transformation.

The inverted relations (229)–(232) read

[TABLE]

and the general transformation from $(1,2)$ to $(r,a)$ basis for any $n$-point function is then

[TABLE]

The functions $\Delta_{ra}(x,y)$ and $\Delta_{ar}(x,y)$ are the usual retarded and advanced Green functions. In case of massless theory they are of the following forms in the momentum space

[TABLE]

and they satisfy

[TABLE]

$\Delta_{rr}(k)$ is the correlation function and it is the only function where a distribution function enters. In thermal equilibrium all these three functions are related via the fluctuation-dissipation theorem

[TABLE]

where $n(k^{0})=1/(e^{\beta k^{0}}-1)$ is the Bose distribution function with $\beta$ being the inverse of temperature $T$.

By means of the relation (241), we immediately find real and imaginary parts of the retarded propagator

[TABLE]

Appendix D Four-point Green functions

Here we provide some useful derivations arising during the analysis of the Bethe-Salpeter equation.

Relations among different four-point Green functions

The four-point Green functions which contribute to the real and imaginary parts of $\bar{G}_{R}^{xy,xy}$ are $G_{arrr}$, $G_{rarr}$, $G_{rrar}$, and $G_{rrra}$. The analysis of the Bethe-Salpeter equation gets easier when one realizes that these functions may be expressed by $G_{aarr}$ and $G_{rraa}$. Let us first notice that by amputating two external legs out of the $G_{aarr}$ function from the left hand side one may write

[TABLE]

and by amputating two external legs of $G_{rraa}$ from the right-hand side one gets

[TABLE]

Using the expressions (245) and (246), we can express $G_{rrra}$, $G_{rrar}$, $G_{rarr}$, and $G_{arrr}$ as follows

[TABLE]

where $N_{p}=N(p_{0})$ and we also have used the identity $N(-p_{0})=-N(p_{0})$. Then, we see

[TABLE]

Kernel of the Bethe-Salpeter equation

The general expression of any rung contributing to the kernel is

[TABLE]

and the bare 4-point vertex is given by

[TABLE]

where $n_{a}$ is the number of $a$ indices among $\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}$. With the definitions (253) and (254), the three rungs under interest are given by

[TABLE]

Then, the following relations hold

[TABLE]

which are helpful in investigation of the Bethe-Salpeter equation for $G_{aarr}^{*}$. In the vanishing frequency and momentum limits, the kernel may be denoted as

[TABLE]

The kernel, as explained in Wang:2002nba , may be also expressed in a more convenient form as

[TABLE]

where $\rho(s)=i[\Delta_{ra}(s)-\Delta_{ar}(s)]$. If we denote $\mathcal{K}(p,-p,-l,l)=\mathcal{K}(p,l)$, the following relations for the kernel hold:

[TABLE]

and

[TABLE]

both of which can be shown by changing the integration variable from $s$ to $-s$. One can also show

[TABLE]

by changing the integration variable $s$ to $s-l+p$.