text:
# Nonlinear Response from Transport Theory and Quantum Field Theory at Finite Temperature

## I Introduction

Fluctuations occur in a system perturbed slightly away from equilibrium. The responses to these fluctuations are described by transport coefficients which characterize the dynamics of long wavelength, low frequency fluctuations in the medium . The investigation of transport coefficients in high temperature gauge theories is important in cosmological applications such as electroweak baryogenesis  and in the context of heavy ion collisions .

There are two basic methods to calculate transport coefficients: transport theory and response theory . Using the transport theory method one starts from a local equilibrium form for the distribution function and performs an expansion in the gradient of the four-velocity field. The coefficients of this expansion are determined from the classical Boltzmann equation . In the response theory approach one divides the Hamiltonian into a bare piece and a perturbative piece that is linear in the gradient of the four-velocity field. One uses standard perturbation theory to obtain the Kubo formula for the viscosity in terms of retarded green functions . These green functions are then evaluated using equilibrium quantum field theory. As is typical in finite temperature field theory, it is not sufficient to calculate perturbatively in the coupling constant: there are certain infinite sets of diagrams that contribute at the same order in perturbation theory and have to be resummed .

In this paper, we want to compare these two methods. The response theory approach allows us to calculate transport coefficients from first principles using the well understood methods of quantum field theory. On the other hand, the transport theory approach involves the use of the Boltzmann equation which is itself derived from some more fundamental theory using, among other things, the quasiparticle approximation. In this sense, the response theory approach is more fundamental than the transport theory method. However, the response theory approach can be quite difficult to implement, even for a high temperature weakly coupled scalar theory, because of the need to resum infinite sets of diagrams . These considerations motivate us to understand more precisely the connection between the more practical transport theory method, and the more fundamental response theory approach.

Some progress has already been made in this direction. It has been shown that keeping only terms which are linear in the gradient expansion in the transport theory calculation is equivalent to using the linear response approximation to obtain the usual Kubo formula for the shear viscosity in terms of a retarded two-point function, and calculating that two-point function using standard equilibrium quantum field theory techniques to resum an infinite set of ladder diagrams . This result is not surprising since it has been known for some time that ladder diagrams give large contributions to $`n`$-point functions with ultra soft external lines .

To date, calculations of transport coefficients have been limited to linear response. In some physical situations however nonlinear response can be important, especially for relativistic gauge theories . It is therefore of interest to study nonlinear response. In this paper we study nonlinear response to fluctuations in weakly coupled high temperature scalar $`\varphi ^4`$ theory using both the transport method and the response theory approach. We study the relationship between these two approaches at the level of quadratic response.

This paper will be organized as follows. In section II we define shear viscosity and quadratic shear viscous response using a hydrodynamic expansion of the energy-momentum tensor which includes up to quadratic terms in the gradient of the four-velocity. In section III we calculate viscosity using the transport theory method by performing a gradient expansion on the Boltzmann equation and obtaining a hierarchy of equations. In section IV we derive the quadratic response Kubo formula by using Zubrave’s techniques in nonequilibrium statistical mechanics. We obtain an expression which relates the quadratic shear viscous response coefficient to the three-point retarded green function of the viscous shear stress tensor. Starting from this generalized Kubo formula we calculate the quadratic shear viscous response using standard techniques of finite temperature quantum field theory. We show that the quadratic shear viscous response can be obtained as an integral over a four-point vertex. This four-point vertex satisfies an integral equation involving terms which are quadratic in a retarded three-point function, which itself satisfies a linear integral equation. These two integral equations represent the resummation of an infinite series of ladder and extended-ladder diagrams. We show that these integral equations, which represent ladder and extended-ladder resummations, have the same form as the first two equations in the hierarchy obtained from expanding the Boltzmann equation. We discuss our results and present our conclusions section V.

## II Viscosity

In a system that is out of equilibrium, the existence of gradients in thermodynamic parameters like the temperature and the four dimensional velocity field give rise to thermodynamic forces. These thermodynamic forces lead to deviations from the equilibrium expectation value of the energy momentum tensor which are characterized by transport coefficients like the thermal conductivity and the shear and bulk viscosities. In order to separate these different physical processes we decompose the energy-momentum tensor as,

$$T^{\mu \nu }=ϵu^\mu u^\nu p\mathrm{\Delta }^{\mu \nu }+P^\mu u^\nu +P^\nu u^\mu +\pi ^{\mu \nu };\mathrm{\Delta }_{\mu \nu }=g_{\mu \nu }u_\mu u_\nu .$$
(1)
The quantities $`ϵ`$, $`p`$, $`P_\mu `$ and $`\pi _{\mu \nu }`$ have the physical meanings of internal energy density, pressure, heat current and viscous shear stress, respectively. The four vector $`u_\mu (x)`$ is the four dimensional four-velocity field which satisfies $`u^\mu (x)u_\mu (x)=1`$. The expansion coefficients are given by

$`ϵ`$ $`=`$ $`u_\alpha u_\beta T^{\alpha \beta }`$ (2)
$`p`$ $`=`$ $`{\displaystyle \frac{1}{3}}\mathrm{\Delta }_{\alpha \beta }T^{\alpha \beta }`$ (3)
$`P_\mu `$ $`=`$ $`\mathrm{\Delta }_{\mu \alpha }u_\beta T^{\alpha \beta }`$ (4)
$`\pi _{\mu \nu }`$ $`=`$ $`(\mathrm{\Delta }_{\mu \alpha }\mathrm{\Delta }_{\nu \beta }{\displaystyle \frac{1}{3}}\mathrm{\Delta }_{\nu \mu }\mathrm{\Delta }_{\alpha \beta })T^{\alpha \beta }.`$ (5)
The viscosity is obtained from the expectation value of the viscous shear stress part of the energy momentum tensor. We expand in gradients of the four-velocity field and write,

$`\delta \pi _{\mu \nu }=\eta ^{(1)}H_{\mu \nu }+\eta ^{(2)}H_{\mu \nu }^{T2}+\mathrm{}`$ (6)
$`H_{\mu \nu }=_\mu u_\nu +_\nu u_\nu {\displaystyle \frac{2}{3}}\mathrm{\Delta }_{\mu \nu }\mathrm{\Delta }_{\rho \sigma }^\rho u^\sigma `$ (7)
$`H_{\mu \nu }^{T2}:=H_{\mu \rho }H_\nu ^\rho {\displaystyle \frac{1}{3}}\mathrm{\Delta }_{\mu \nu }H_{\rho \sigma }H^{\rho \sigma }`$ (8)
where $`\eta ^{(1)}`$ and $`\eta ^{(2)}`$ are the coefficients of the terms that are linear and quadratic respectively in the gradient of the four-velocity. The first coefficient is the usual shear viscosity. The second has has not been widely discussed in the literature – we shall call it the quadratic shear viscous response.

Throughout this paper we work with $`\varphi ^4`$ theory. The Lagrangian for this theory is

$`={\displaystyle \frac{1}{2}}[(_\mu \varphi )^2m^2\varphi ^2]{\displaystyle \frac{\lambda }{4!}}\varphi ^4`$ (9)
and the coupling constant is assumed to be small: $`\lambda 1`$.

## III Viscosity from Transport Theory

Kinetic theory and the Boltzmann equation can be used to calculate transport properties of dilute many-body systems. One assumes that, except during brief collisions, the system can be considered as being composed of classical particles with well defined position, energy and momentum. This picture is valid when the mean free path is large compared with the Compton wavelength of the particles. At high temperature the typical mean free path of thermal excitations is $`𝒪(1/\lambda ^2T)`$ and is always larger than the typical Compton wavelength of effective thermal oscillations which is $`𝒪(1/\sqrt{\lambda }T)`$ . We introduce a phase space distribution function $`f(x,\underset{¯}{k})`$ which describes the evolution of the phase space probability density for the fundamental particles comprising a fluid. In this expression and in the following equations the underlined momenta are on shell, since we are describing a system of particles. The form for $`f(x,\underset{¯}{k})`$ in local equilibrium is,

$$f^{(0)}=\frac{1}{e^{\beta (x)u_\mu (x)\underset{¯}{k}^\mu }1}:=n_k;N_k:=1+2n_k.$$
(10)
We study the Boltzmann equation in the hydrodynamic regime where we consider times which are long compared to the mean free time and describe the relaxation of the system in terms of long wavelength fluctuations in locally conserved quantities. For a simple fluid without any additional conserved charges, the only locally conserved quantities are energy and momentum. To solve the Boltzmann equation in this near equilibrium hydrodynamic regime, we expand the distribution function around the local equilibrium form using a gradient expansion. We go to a local rest frame in which we can write $`\stackrel{}{u}(x)=0`$. Note that this does not imply that gradients of the form $`_iu_j`$ must be zero. In the local rest frame (6) becomes,

$`\delta \pi _{ij}=\eta ^{(1)}H_{ij}+\eta ^{(2)}H_{ij}^{T2}+\mathrm{}`$ (11)
$`H_{ij}=_iu_j+_ju_i{\displaystyle \frac{2}{3}}\delta _{ij}(\stackrel{}{}\stackrel{}{u})`$ (12)
$`H_{ij}^{T2}:=H_{ik}H_j^k{\displaystyle \frac{1}{3}}\delta _{ij}H_{lm}H_{lm}`$ (13)
In all of the following expressions we keep only linear terms that contain one power of $`H_{ij}`$ and quadratic terms that contain two powers of $`H_{ij}`$, since these are the only terms that contribute to the viscosity coefficients we are trying to calculate.

We write,

$`f=f^{(0)}+f^{(1)}+f^{(2)}+\mathrm{}`$ (14)
with,

$`f^{(1)}\underset{¯}{k}_\mu ^\mu f^{(0)};f^{(2)}\underset{¯}{k}_\mu ^\mu f^{(1)}.`$ (15)
Using (10) we obtain,

$`\underset{¯}{k}_\mu ^\mu f^{(0)}=\beta n_k(1+n_k){\displaystyle \frac{1}{2}}I_{ij}(k)H_{ij}`$ (16)
$`f^{(1)}:=n_k(1+n_k)\varphi _k;\varphi _k=\beta {\displaystyle \frac{1}{2}}B_{ij}(\underset{¯}{k})H_{ij}`$ (17)
$`\underset{¯}{k}_\mu ^\mu f^{(1)}=\beta ^2n_k(1+n_k)N_k{\displaystyle \frac{1}{2}}I_{ij}(k)H_{ij}{\displaystyle \frac{1}{2}}B_{lm}(\underset{¯}{k})H_{lm}`$ (18)
$`f^{(2)}:=n_k(1+n_k)N_k\theta _k;\theta _k:=\beta ^2{\displaystyle \frac{1}{4}}C_{ijlm}(\underset{¯}{k})H_{ij}H_{lm}`$ (19)
where we define

$`\widehat{I}_{ij}(k)=(\widehat{k}_i\widehat{k}_j{\displaystyle \frac{1}{3}}\delta _{ij});I_{ij}(k)=k^2\widehat{I}_{ij}(k)`$ (20)
and write,

$`B(\underset{¯}{k})_{ij}=\widehat{I}_{lm}(k)B(\underset{¯}{k});C_{ijlm}(\underset{¯}{k})=\widehat{I}_{ij}(k)\widehat{I}_{lm}(k)C(\underset{¯}{k}).`$ (21)
The viscous shear stress part of the energy momentum tensor is given by

$`\pi _{ij}={\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}f(k_ik_j\frac{1}{3}\delta _{ij}k^2)}.`$ (22)
Using the expansion (14) we get,

$`\pi _{ij}={\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}[f^{(0)}+f^{(1)}+f^{(2)}](k_ik_j\frac{1}{3}\delta _{ij}k^2)}.`$ (23)
The lowest order term is identically zero. We obtain the linear and quadratic contributions by substituting in (16) and (21). We use $`H_{ii}:=0`$ and the following results which are obtained from rotational invariance:

$`k_ik_jB(\underset{¯}{k}){\displaystyle \frac{1}{3}}\delta _{ij}k^2B(\underset{¯}{k})`$ (24)
$`k_ik_j\widehat{k}_l\widehat{k}_mB(\underset{¯}{k}){\displaystyle \frac{1}{15}}(\delta _{ij}\delta _{lm}+\delta _{il}\delta _{jm}+\delta _{im}\delta _{jl})k^2B(\underset{¯}{k})`$ (25)
$`k_ik_j\widehat{k}_l\widehat{k}_m\widehat{k}_a\widehat{k}_bC(\underset{¯}{k}){\displaystyle \frac{1}{105}}[\delta _{ab}(\delta _{lm}\delta _{ij}+\delta _{lj}\delta _{mi}+\delta _{li}\delta _{mj})+\delta _{al}(\delta _{bm}\delta _{ij}+\delta _{bi}\delta _{mj}+\delta _{bj}\delta _{mi})`$ (26)
$`+\delta _{am}(\delta _{bl}\delta _{ij}+\delta _{bj}\delta _{li}+\delta _{bi}\delta _{lj})+\delta _{ai}(\delta _{bl}\delta _{mj}+\delta _{bm}\delta _{lj}+\delta _{bj}\delta _{ml})`$ (27)
$`+\delta _{aj}(\delta _{bl}\delta _{mi}+\delta _{bm}\delta _{li}+\delta _{bi}\delta _{ml})]k^2C(\underset{¯}{k})`$ (28)
We obtain,

$`\delta \pi _{ij}={\displaystyle \frac{\beta }{15}}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}n_k(1+n_k)k^2B(\underset{¯}{k})H_{ij}}+{\displaystyle \frac{2\beta ^2}{105}}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}[n_k(1+n_k)N_k]k^2C(\underset{¯}{k})H_{ij}^{T2}}`$ (29)
Comparing with (11) we have,

$`\eta ^{(1)}={\displaystyle \frac{\beta }{15}}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}n_k(1+n_k)k^2B(\underset{¯}{k})}`$ (30)
$`\eta ^{(2)}={\displaystyle \frac{2\beta ^2}{105}}{\displaystyle \frac{d^3k}{(2\pi )^32\omega _k}[n_k(1+n_k)N_k]k^2C(\underset{¯}{k})}`$ (31)
Thus we have shown that the shear viscosity and the quadratic shear viscous response can be obtained from the functions $`B(\underset{¯}{k})`$ and $`C(\underset{¯}{k})`$ respectively. These two functions are the coefficients of the linear and quadratic terms in the gradient expansion of the distribution function. In the next section we will show that these functions can be obtained from the first two equations in the hierarchy of equations obtained from the gradient expansion of the Boltzmann equation.

### A Expansion of the Boltzmann Equation

The Boltzmann equation describes the evolution of the distribution function $`f(x,\underset{¯}{k})`$ and can be used to obtain integral equations for the functions $`B(\underset{¯}{k})`$ and $`C(\underset{¯}{k})`$ defined in (16) and (21). The Boltzmann equation has the form:

$`\underset{¯}{k}_\mu ^\mu f(x,\underset{¯}{k})=𝒞[f]`$ (32)
where $`𝒞[f]`$ is the collision term:

$`𝒞[f]={\displaystyle \frac{1}{2}}{\displaystyle _{123}}𝑑\mathrm{\Gamma }_{123k}[f_1f_2(1+f_3)(1+f_k)(1+f_1)(1+f_2)f_3f_k]`$ (33)
with $`f_i:=f(x,\underset{¯}{p}_i)`$, $`f_k:=f(x,\underset{¯}{k})`$. The symbol $`d\mathrm{\Gamma }_{123k}`$ represents the differential transition rate for particles of momentum $`P_1`$ and $`P_2`$ to scatter into momenta $`P_3`$ and $`K`$ and is given by

$`d\mathrm{\Gamma }_{123k}:={\displaystyle \frac{1}{2\omega _k}}|𝒯(\underset{¯}{k},\underset{¯}{p}_3,\underset{¯}{p}_2,\underset{¯}{p}_1|^2\mathrm{\Pi }_{i=1}^3{\displaystyle \frac{d^3p_i}{(2\pi )^32\omega _{p_i}}}(2\pi )^4\delta (\underset{¯}{p}_1+\underset{¯}{p}_2\underset{¯}{p}_3\underset{¯}{k})`$ (34)
where $`𝒯`$ is the multiparticle scattering amplitude. Using the expansion (14) and (16) produces a hierarchy of equations.

#### 1 First Order Boltzmann Equation

The first order equation is:

$$\underset{¯}{k}^\mu _\mu f^0(x,\underset{¯}{k})=𝒞[f^{(0)};f^{(1)}]$$
(35)
where we keep terms linear in $`f^{(1)}`$ on the right hand side. The left hand side gives,

$`\underset{¯}{k}^\mu _\mu f^0(x,\underset{¯}{k})=n_k(1+n_k)\beta {\displaystyle \frac{1}{2}}I_{ij}(k)H_{ij}`$ (36)
The right hand side is

$$𝒞[f^{(1)}]=\frac{1}{2}_{123}𝑑\mathrm{\Gamma }_{123k}(1+n_1)(1+n_2)n_3n_k(\varphi _k+\varphi _3\varphi _1\varphi _2)$$
(37)
Using the definition of $`\varphi `$ given in (16) and comparing the coefficients of $`H_{ij}`$ on both sides of (35) we obtain,

$`I_{ij}(k)={\displaystyle \frac{1}{2}}{\displaystyle _{123}}𝑑\mathrm{\Gamma }_{123k}{\displaystyle \frac{(1+n_1)(1+n_2)n_3}{1+n_k}}[B_{ij}(\underset{¯}{p}_1)+B_{ij}(\underset{¯}{p}_2)B_{ij}(\underset{¯}{k})B_{ij}(\underset{¯}{p}_3)]`$ (38)
This result is an inhomogeneous linear integral equation which can be solved self-consistently to obtain the function $`B_{ij}(\underset{¯}{k})`$.

#### 2 Second Order Boltzmann Equation

The second order contribution to (32) is,

$$\underset{¯}{k}^\mu _\mu f^{(1)}(x,\underset{¯}{k})=𝒞[f^{(0)};f^{(1)};f^{(2)}]$$
(39)
where we keep terms linear in $`f^{(2)}`$ and quadratic in $`f^{(1)}`$ on the right hand side. Using (16) the left hand side becomes,

$`\underset{¯}{k}^\mu _\mu f^{(1)}=`$ $`n_k(1+n_k)N_k\beta ^2{\displaystyle \frac{1}{4}}I_{ij}(k)B_{lm}(\underset{¯}{k})H_{ij}H_{lm}`$ (40)
and the right hand side gives,

$`C^{(2)}[f^{(0)};f^{(1)};f^{(2)}]=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{123}}𝑑\mathrm{\Gamma }_{123k}(1+n_1)(1+n_2)n_3n_k`$ (44)
$`[(N_1\theta _1+N_2\theta _2N_3\theta _3N_k\theta _k)`$
$`+{\displaystyle \frac{1}{2}}((N_1+N_2)\varphi _1\varphi _2(N_k+N_3)\varphi _3\varphi _k+(N_3N_1)\varphi _1\varphi _3`$
$`+(N_kN_1)\varphi _1\varphi _k+(N_3N_2)\varphi _3\varphi _2+(N_kN_2)\varphi _k\varphi _2)]`$
Using the definitions of $`\varphi `$ and $`\theta `$ given in (16) and comparing the coefficients of $`H_{ij}H_{lm}`$ on both sides we obtain,

$`n_k`$ $`(1+n_k)N_kI_{ij}(k)B_{lm}(\underset{¯}{k})={\displaystyle \frac{1}{2}}{\displaystyle _{123}}𝑑\mathrm{\Gamma }_{123k}(1+n_1)(1+n_2)n_3n_k`$ (48)
$`\{[N_1C_{ijlm}(\underset{¯}{p}_1)+N_2C_{ijlm}(\underset{¯}{p}_2)N_kC_{ijlm}(\underset{¯}{k})N_3C_{ijlm}(\underset{¯}{p}_3)]`$
$`+{\displaystyle \frac{1}{2}}[((N_1+N_2)B_{ij}(\underset{¯}{p}_1)B_{lm}(\underset{¯}{p}_2)(N_k+N_3)B_{ij}(\underset{¯}{p}_3)B_{lm}(\underset{¯}{k})+(N_3N_1)B_{ij}(\underset{¯}{p}_1)B_{lm}(\underset{¯}{p}_3)`$
$`+(N_kN_1)B_{ij}(\underset{¯}{p}_1)B_{lm}(\underset{¯}{k})+(N_3N_2)B_{ij}(\underset{¯}{p}_3)B_{lm}(\underset{¯}{p}_2)+(N_kN_2)B_{ij}(\underset{¯}{k})B_{lm}(\underset{¯}{p}_2)]\}`$
This integral equation can be solved self consistently for the quantity $`C_{ijlm}(\underset{¯}{k})`$ using the result for $`B_{ij}(\underset{¯}{k})`$ from (38).

## IV Viscosity from Field Theory

The Kubo formulae allow us to use quantum field theory to calculate nonequilibrium transport coefficients. The results should be the same as those obtained in the previous section using transport theory. To do calculations in a system that is out of equilibrium, the Hamiltonian is separated into an equilibrium piece $`H_0`$ and a nonequilibrium piece $`H_{ext}`$ which depends on the gradients of the thermodynamic parameters: the four-velocity field and the inverse temperature field. For systems close to equilibrium the nonequilibrium piece of the Hamiltonian can be treated as a perturbation and the deviations of physical quantities from their equilibrium values can be calculated perturbatively. The linear response calculation includes only the first order contribution to this expansion and gives transport coefficients that can be expressed as integrals of retarded two-point Green functions over space and time. One of the results we obtain in this paper is that the quadratic shear viscous response coefficient can be written in terms of a retarded three-point function. Throughout this section we use capital letters to denote four-vectors and small letters for three-vectors. We also define $`d^4p/(2\pi )^4:=𝑑P`$.

To calculate the expectation value of an operator we take the trace over the density matrix. We follow the presentation of . We work with the density matrix in the Heisenberg representation which satisfies,

$$\frac{\rho }{t}=0$$
(49)
and can be written as,

$$\rho =\frac{e^{A+B}}{\mathrm{Tr}e^{A+B}}$$
(50)
where

$`A`$ $`=`$ $`{\displaystyle d^3xF^\nu T_{0\nu }},`$ (51)
$`B`$ $`=`$ $`{\displaystyle d^3x_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}T_{\mu \nu }(x,t^{})^\mu F^\nu (x,t^{})}`$ (52)
with $`F^\mu =\beta u^\mu `$ and $`ϵ`$ to be taken to zero at the end. In this expression $`A`$ is the equilibrium part of the Hamiltonian and $`B`$ is a perturbative contribution that is linear in the gradient of the four-velocity field. Note that in the local rest frame $`u^\mu =(1,0,0,0)`$ and $`[A]^{LRF}=d^3x\beta T_{00}`$. Using the identity for the exponential function of two operators

$$e^{\beta (\widehat{a}+\widehat{b})}=e^{\beta \widehat{a}}[1+_0^\beta 𝑑\lambda e^{\lambda \widehat{a}}\widehat{b}e^{\lambda (\widehat{a}+\widehat{b})}]$$
(53)
we can expand

$`e^{A+B}=`$ $`e^A[1+{\displaystyle _0^1}d\lambda e^{\lambda A}Be^{\lambda A}`$ (55)
$`+{\displaystyle _0^1}d\lambda {\displaystyle _0^\lambda }d\tau e^{\lambda A}Be^{\lambda A}e^{\tau A}Be^{\tau A}+𝒪(B^3)]`$
Thus the density matrix can be written as

$`\rho =`$ $`\rho _0[1+{\displaystyle _0^1}d\lambda (B(\lambda )B(\lambda ))+{\displaystyle _0^1}d\lambda {\displaystyle _0^\lambda }d\tau (B(\lambda )B(\tau )B(\lambda )B(\tau ))`$ (57)
$`{\displaystyle _0^1}d\lambda {\displaystyle _0^1}d\lambda ^{}(B(\lambda )B(\lambda ^{})B(\lambda )B(\lambda ^{})]+O(B^3)`$
where

$$\rho _0=\frac{e^A}{\mathrm{Tr}e^A}$$
(58)
is the local equilibrium density matrix.

### A Viscosity as an Expansion in Green Functions of Composite Operators

#### 1 Linear Response

Using the first three terms of (57) produces the linear response approximation for the deviation of the expectation value of the viscous shear stress part of the energy momentum tensor from the equilibrium value:

$`\delta \pi _{\mu \nu }(x,t)^l`$ $`=\pi _{\mu \nu }(x,t)^l\pi _{\mu \nu }(x,t)_0`$ (60)
$`={\displaystyle d^3x^{}_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}(\pi _{\alpha \beta }(x^{},t^{}),\pi _{\mu \nu }(x,t))^\alpha F^\beta (x^{},t^{\prime \prime })}`$
where the correlation function $`(\pi _{\alpha \beta }(x^{},t^{}),\pi _{\mu \nu }(x,t))`$ is defined as,

$$(\pi _{\mu \nu }(x,t),\pi _{\alpha \beta }(x^{},t^{}))=\frac{1}{\beta }_0^\beta d\tau (\pi _{\alpha \beta }(x^{},t^{}+i\tau )\pi _{\mu \nu }(x,t)\pi _{\alpha \beta }(x^{},t^{}+i\tau )_0\pi _{\mu \nu }(x,t)_0$$
(61)
and $`\mathrm{}_0:=\mathrm{Tr}(\rho _0\mathrm{})`$. From now on we drop the subscript $`0`$ on the correlation functions. We take the initial time $`t_0`$ to minus infinity and we assume the system is in equilibrium at $`t=t_0`$ and that the external forces are switched on adiabatically. We assume that changes in the thermodynamic forces are small enough over the correlation length of the correlation functions that factors $`F_\mu =\beta u_\mu `$ can be taken out of the integral. Using (4) and rotational invariance we have,

$`\pi _{\mu \nu }\pi _{\alpha \beta }={\displaystyle \frac{1}{10}}(\mathrm{\Delta }_{\mu \alpha }\mathrm{\Delta }_{\nu \beta }+\mathrm{\Delta }_{\mu \beta }\mathrm{\Delta }_{\nu \alpha }{\displaystyle \frac{2}{3}}\mathrm{\Delta }_{\mu \nu }\mathrm{\Delta }_{\alpha \beta })\pi _{\sigma \tau }\pi ^{\sigma \tau }`$ (62)
Using this result we obtain,

$`\delta \pi _{\mu \nu }^l`$ $`={\displaystyle \frac{H_{\mu \nu }}{10}}{\displaystyle d^3x^{}_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}}`$ (64)
$`{\displaystyle _0^\beta }((\pi (x^{},t+i\tau )\pi (x,t)\pi (x^{},t+i\tau )\pi (x,t))`$
where we have written $`\pi _{\sigma \tau }\pi ^{\sigma \tau }:=\pi \pi `$. If we assume that correlations vanish as $`t^{}\mathrm{}`$ this expression can be rewritten as,

$`\delta \pi _{\mu \nu }^l`$ $`={\displaystyle \frac{H_{\mu \nu }}{10}}{\displaystyle d^3x^{}_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}}`$ (66)
$`{\displaystyle _0^\beta }{\displaystyle _{\mathrm{}}^t^{}}dt^{\prime \prime }{\displaystyle \frac{d}{dt^{\prime \prime }}}((\pi (x^{},t^{\prime \prime }+i\tau )\pi (x,t)\pi (x^{},t^{\prime \prime }+i\tau )\pi (x,t))`$
Using $`\frac{}{dt^{\prime \prime }}f(t^{\prime \prime }+i\tau )=i\frac{}{\tau }f(t^{\prime \prime }+i\tau )`$ we can perform the integration over $`\tau `$. We obtain,

$`\delta \pi _{\mu \nu }^l`$ (67)
$`={\displaystyle \frac{i}{10}}H_{\mu \nu }{\displaystyle d^3x^{}_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t^{}𝑑t^{\prime \prime }(\pi (x^{},t^{\prime \prime }+i\beta )\pi (x,t)\pi (x^{},t^{\prime \prime })\pi (x,t))}.`$ (68)
Using the KMS condition

$`\pi (x^{},t^{\prime \prime }+i\beta )\pi (x,t)`$ $`=\mathrm{Tr}[e^{\beta H}\pi (x^{},t^{\prime \prime }+i\beta )\pi (x,t)]`$ (70)
$`=\mathrm{Tr}[e^{\beta H}\pi (x,t)\pi (x^{},t^{\prime \prime })]=\pi (x,t)\pi (x^{},t^{\prime \prime })`$
we obtain,

$$\delta \pi _{\mu \nu }^l=\frac{H_{\mu \nu }}{10}d^3x^{}_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t^{}𝑑t^{\prime \prime }D_R(x,t;x^{},t^{\prime \prime })$$
(71)
where

$$D_R(x,t;x^{},t^{\prime \prime })=i\theta (tt^{\prime \prime })[\pi (x,t),\pi (x^{},t^{\prime \prime })]$$
(72)
We extract the shear viscosity using the definition (6). We obtain,

$`\eta ^{(1)}={\displaystyle \frac{1}{10}}{\displaystyle d^3x^{}_{\mathrm{}}^0𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t^{}𝑑t^{\prime \prime }D_R(0;x^{},t^{\prime \prime })}`$ (73)
We can rewrite this result in a more useful form. We rewrite the integrals inserting theta functions: $`_{\mathrm{}}^0𝑑t^{}=\theta (t^{})_{\mathrm{}}^{\mathrm{}}𝑑t^{}`$ and use the integral representation for the theta function:

$`\theta (x)={\displaystyle \frac{d\omega }{2\pi i}\frac{e^{i\omega x}}{\omega iϵ}}.`$ (74)
Inserting $`D_R(0;x^{},t^{\prime \prime })=𝑑Qe^{iQX^{}}D_R(Q)`$ with $`X^{}=(x^{},t^{\prime \prime })`$ we obtain,

$`\eta ^{(1)}={\displaystyle \frac{i}{10}}{\displaystyle \frac{d}{dq_0}}[\underset{\stackrel{}{q}0}{lim}D_R(Q)]|_{q_0=0}`$ (75)
Since it is the imaginary part of the two-point function that is odd in $`q_0`$ we have,

$`\eta ^{(1)}={\displaystyle \frac{1}{10}}{\displaystyle \frac{d}{dq_0}}\mathrm{Im}[\underset{\stackrel{}{q}0}{lim}D_R(Q)]|_{q_0=0}`$ (76)
This is the well known Kubo formula which expresses the shear viscosity in terms of a retarded two-point green function that can be calculated using equilibrium quantum field theory.

#### 2 Quadratic Response

Now we consider corrections to the linear response approximation. We calculate the quadratic shear viscous response from the terms in (57) that are quadratic in the interaction. We show that the result can be written as a retarded three-point correlator. We obtain,

$`\delta \pi _{\mu \nu }(x,t)^q=`$ $`{\displaystyle d^3x^{}d^3x^{\prime \prime }_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t𝑑t^{\prime \prime }e^{ϵ(t^{\prime \prime }t)}}`$ (78)
$`(\pi _{\mu \nu }(x,t),\pi _{\alpha \beta }(x^{},t^{}),\pi _{\rho \sigma }(x^{\prime \prime },t^{\prime \prime }))^\alpha F^\beta (x^{},t^{})^\rho F^\sigma (x^{\prime \prime },t^{\prime \prime }).`$
We use (4) and rotational invariance to write,

$`\pi _{\mu \nu }(x,t)\pi _{\alpha \beta }(x^{},t^{})\pi _{\rho \sigma }(x^{\prime \prime },t^{\prime \prime })^\alpha u^\beta ^\lambda u^\tau `$ $`={\displaystyle \frac{3}{35}}H_{\mu \nu }^{T2}(\pi _{\alpha \beta }(x,t)\pi ^{\beta \lambda }(x^{},t^{})\pi _\lambda ^\alpha (x^{\prime \prime },t^{\prime \prime }))`$ (80)
$`:={\displaystyle \frac{3}{35}}H_{\mu \nu }^{T2}(\pi (x,t)\pi (x^{},t^{})\pi (x^{\prime \prime },t^{\prime \prime }))`$
Writing out the correlation function we obtain,

$`\delta \pi _{\mu \nu }(x,t)^q=`$ $`{\displaystyle \frac{3}{35}}H_{\mu \nu }^{T2}{\displaystyle d^3x^{}d^3x^{\prime \prime }_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t𝑑t^{\prime \prime }e^{ϵ(t^{\prime \prime }t)}}`$ (85)
$`[{\displaystyle _0^\beta }d\tau {\displaystyle _0^\tau }d\lambda \pi (x^{},t^{}+i\tau )\pi (x^{\prime \prime },t^{\prime \prime }+i\lambda )\pi (x,t)`$
$`{\displaystyle _0^\beta }𝑑\tau {\displaystyle _0^\beta }𝑑\tau ^{}\pi (x^{},t^{}+i\tau )\pi (x^{\prime \prime },t^{\prime \prime }+i\tau ^{})\pi (x,t)`$
$`{\displaystyle _0^\beta }𝑑\tau {\displaystyle _0^\tau }𝑑\lambda \pi (x^{},t^{}+i\tau )\pi (x^{\prime \prime },t^{\prime \prime }+i\lambda )\pi (x,t)`$
$`+{\displaystyle _0^\beta }d\tau {\displaystyle _0^\beta }d\tau ^{}\pi (x^{},t^{}+i\tau )\pi (x^{\prime \prime },t^{\prime \prime }+i\tau ^{})\pi (x,t)]`$
Assuming that correlations vanish when the time approaches minus infinity we can rewrite (85) as

$`\delta \pi _{\mu \nu }(x,t)^q=`$ $`{\displaystyle \frac{3}{35}}H_{\mu \nu }^{T2}{\displaystyle d^3x^{}d^3x^{\prime \prime }_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t𝑑t^{\prime \prime }e^{ϵ(t^{\prime \prime }t)}_{\mathrm{}}^t^{}𝑑s^{}_{\mathrm{}}^{t^{\prime \prime }}𝑑s^{\prime \prime }}`$ (90)
$`{\displaystyle \frac{d}{ds^{}}}{\displaystyle \frac{d}{ds^{\prime \prime }}}[{\displaystyle _0^\beta }d\tau {\displaystyle _0^\tau }d\lambda \pi (x^{},s^{}+i\tau )\pi (x^{\prime \prime },s^{\prime \prime }+i\lambda )\pi (x,t)`$
$`{\displaystyle _0^\beta }𝑑\tau {\displaystyle _0^\beta }𝑑\tau ^{}\pi (x^{},s^{}+i\tau )\pi (x^{\prime \prime },s^{\prime \prime }+i\tau ^{})\pi (x,t)`$
$`{\displaystyle _0^\beta }𝑑\tau {\displaystyle _0^\tau }𝑑\lambda \pi (x^{},s^{}+i\tau )\pi (x^{\prime \prime },s^{\prime \prime }+i\lambda )\pi (x,t)`$
$`+{\displaystyle _0^\beta }d\tau {\displaystyle _0^\beta }d\tau ^{}\pi (x^{},s^{}+i\tau )\pi (x^{\prime \prime },s^{\prime \prime }+i\tau ^{})\pi (x,t)]`$
Carrying out the integration over $`\lambda `$, $`\tau `$ and $`\tau ^{}`$ and using the fact that we have symmetry under interchange of $`(x^{},t^{})`$ and $`(x^{\prime \prime },t^{\prime \prime })`$ we obtain,

$`\delta \pi _{\mu \nu }(x,t)^q=`$ $`{\displaystyle \frac{3}{70}}H_{\mu \nu }^{T2}{\displaystyle d^3x^{}d^3x^{\prime \prime }_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t𝑑t^{\prime \prime }e^{ϵ(t^{\prime \prime }t)}_{\mathrm{}}^t^{}𝑑s^{}_{\mathrm{}}^{t^{\prime \prime }}𝑑s^{\prime \prime }}`$ (92)
$`{\displaystyle \frac{1}{2}}([[\pi (x,t),\pi (x^{},s^{})],\pi (x^{\prime \prime },s^{\prime \prime })]+[[\pi (x,t),\pi (x^{\prime \prime },s^{\prime \prime })],\pi (x^{},s^{})])`$
We define

$`G(x,t;x^{},s^{};x^{\prime \prime },s^{\prime \prime })={\displaystyle \frac{1}{2}}([[\pi (x,t),\pi (x^{},s^{})],\pi (x^{\prime \prime },s^{\prime \prime })]+[[\pi (x,t),\pi (x^{\prime \prime },s^{\prime \prime })],\pi (x^{},s^{})])`$ (93)
and write

$`\delta \pi _{\mu \nu }(x,t)^q`$ (94)
$`={\displaystyle \frac{3}{70}}H_{\mu \nu }^{T2}{\displaystyle d^3x^{}d^3x^{\prime \prime }_{\mathrm{}}^t𝑑t^{}e^{ϵ(t^{}t)}_{\mathrm{}}^t𝑑t^{\prime \prime }e^{ϵ(t^{\prime \prime }t)}_{\mathrm{}}^t^{}𝑑s^{}_{\mathrm{}}^{t^{\prime \prime }}𝑑s^{\prime \prime }G(x,t;x^{},s^{};x^{\prime \prime },s^{\prime \prime })}`$ (95)
Using (6) we extract,

$`\eta ^{(2)}={\displaystyle \frac{3}{70}}{\displaystyle d^3x^{}_{\mathrm{}}^0𝑑t^{}e^{ϵt^{}}_{\mathrm{}}^0𝑑t^{\prime \prime }e^{ϵt^{\prime \prime }}_{\mathrm{}}^t^{}𝑑s^{}_{\mathrm{}}^{t^{\prime \prime }}𝑑s^{\prime \prime }G(0;x^{},s^{};x^{\prime \prime },s^{\prime \prime })}`$ (96)
As we did previously, we can rewrite this result in a neater form. We rewrite the integrals inserting theta functions: $`_{\mathrm{}}^0𝑑t^{}=\theta (t^{})_{\mathrm{}}^{\mathrm{}}𝑑t^{}`$ and use the integral representation for the theta function (74). We use $`X^{}=(x^{},s^{})`$ and $`X^{\prime \prime }=(x^{\prime \prime },s^{\prime \prime })`$ and rewrite,

$`G(0;x^{},s^{};x^{\prime \prime },s^{\prime \prime })={\displaystyle 𝑑Qe^{iX^{}Q}𝑑Q^{}e^{iX^{\prime \prime }Q^{}}G(QQ^{},Q,Q^{})}.`$ (97)
We obtain,

$`\eta ^{(2)}={\displaystyle \frac{3}{70}}{\displaystyle \frac{d}{dq_0}}{\displaystyle \frac{d}{dq_0^{}}}[\underset{\stackrel{}{q}0}{lim}G(QQ^{},Q,Q^{})]|_{q_0=q_0^{}=0}`$ (98)
Using the techniques of  to write the three-point function in the spectral representation, it is tedious but straightforward to show that the only contribution to the result above comes from the real part of the three-point function which is retarded with respect to the first leg. In coordinate space this retarded three-point function is written: $`G_{R1}(x,y,z)=\theta (t_xt_y)\theta (t_yt_z)[[\pi (x),\pi (y)],\pi (z)]+\theta (t_xt_z)\theta (t_zt_y)[[\pi (x),\pi (z)],\pi (y)]`$. We obtain,

$`\eta ^{(2)}={\displaystyle \frac{3}{70}}{\displaystyle \frac{d}{dq_0}}{\displaystyle \frac{d}{dq_0^{}}}\mathrm{Re}[\underset{\stackrel{}{q}0}{lim}G_{R1}(QQ^{},Q,Q^{})]|_{q_0=q_0^{}=0}.`$ (99)
This is an interesting new result. We have obtained a type of nonlinear Kubo formula that allows us to obtain the quadratic shear viscous response from a retarded three-point function using equilibrium quantum field theory.

### B Diagrammatic Expansion

We obtain a diagrammatic expansion for the viscosity coefficients given in (76) and (99). We use the closed time path formulation of finite temperature field theory, and work in the Keldysh representation. Several reviews of this technique are available in the literature . The closed time path integration contour involves two branches, one running from minus infinity to positive infinity just above the real axis, and one running back from positive infinity to negative infinity just below the real axis. All fields can take values on either branch of the contour and thus there is a doubling of the number of degrees of freedom. It is straightforward to show that this doubling of degrees of freedom is necessary to obtain finite green functions. We discuss below the structure of correlation functions of field operators.

The two-point function or the propagator can be written as a $`2\times 2`$ matrix of the form

$$D=\left(\begin{array}{cc}D_{11}& D_{12}\\ D_{21}& D_{22}\end{array}\right),$$
(100)
where $`D_{11}`$ is the propagator for fields moving along $`C_1`$, $`D_{12}`$ is the propagator for fields moving from $`C_1`$ to $`C_2`$, etc. The four components are given by

$`D_{11}(XY)`$ $`=`$ $`iT(\varphi (X)\varphi (Y))`$ (101)
$`D_{12}(XY)`$ $`=`$ $`i\varphi (Y)\varphi (X)`$ (102)
$`D_{21}(XY)`$ $`=`$ $`i\varphi (X)\varphi (Y)`$ (103)
$`D_{22}(XY)`$ $`=`$ $`i\stackrel{~}{T}(\varphi (X)\varphi (Y))`$ (104)
where $`T`$ is the usual time ordering operator and $`\stackrel{~}{T}`$ is the anti-chronological time ordering operator. Physical functions are obtained by taking appropriate combinations of the components of the propagator matrix. It is straightforward to show that the usual retarded and advanced propagators: $`D_R=i\theta (x_0y_0)[\varphi (X),\varphi (Y)]`$ and $`D_A=i\theta (y_0x_0)[\varphi (X),\varphi (Y)]`$ are given by the combinations,

$`D_R`$ $`=`$ $`D_{11}D_{12}`$ (105)
$`D_A`$ $`=`$ $`D_{11}D_{21}.`$ (106)
The 1PI part of the two-point function, or the polarization insertion, is obtained by truncating legs. The retarded and advanced parts are given by,

$`\mathrm{\Pi }_R`$ $`=`$ $`\mathrm{\Pi }_{11}+\mathrm{\Pi }_{12}`$ (107)
$`\mathrm{\Pi }_A`$ $`=`$ $`\mathrm{\Pi }_{11}+\mathrm{\Pi }_{21}.`$ (108)
The situation is similar for higher $`n`$-point functions. For example, the three-point function which is retarded with respect to the first leg is given by

$`\mathrm{\Gamma }_{R1}=\mathrm{\Gamma }_{111}+\mathrm{\Gamma }_{112}+\mathrm{\Gamma }_{121}+\mathrm{\Gamma }_{122}.`$ (109)
The other three-point vertices that we will need are:

$`\mathrm{\Gamma }_{R2}=\mathrm{\Gamma }_{111}+\mathrm{\Gamma }_{112}+\mathrm{\Gamma }_{211}+\mathrm{\Gamma }_{212}`$ (110)
$`\mathrm{\Gamma }_{R3}=\mathrm{\Gamma }_{111}+\mathrm{\Gamma }_{121}+\mathrm{\Gamma }_{211}+\mathrm{\Gamma }_{221}`$ (111)
$`\mathrm{\Gamma }_F=\mathrm{\Gamma }_{111}+\mathrm{\Gamma }_{121}+\mathrm{\Gamma }_{212}+\mathrm{\Gamma }_{222}`$ (112)
The four-point function which is retarded with respect to the first leg is given by,

$`M_{R1}=M_{1111}+M_{1112}+M_{1121}+M_{1211}+M_{1122}+M_{1212}+M_{1221}+M_{1222}.`$ (113)
The other four-point vertices that we will need are:

$`M_{R4}=M_{1111}+M_{1121}+M_{1211}+M_{2111}+M_{1221}+M_{2121}+M_{2211}+M_{2221}`$ (114)
$`M_F=M_{1111}+M_{1121}+M_{1211}+M_{1221}+M_{2112}+M_{2212}+M_{2122}+M_{2222}.`$ (115)
There are two relations that we can use to simplify expressions involving these vertices. The first is a consequence of the fact that real time green functions are related to each other through the KMS conditions. These identities are a consequence of symmetries generated through cyclic permutations of field operators. The identity that we will need is ,

$`\mathrm{\Gamma }_F+N_1\mathrm{\Gamma }_{R3}+N_3\mathrm{\Gamma }_{R1}=(N_1+N_3)\mathrm{\Gamma }_{R2}^{}`$ (116)
where each vertex carries arguments $`\mathrm{\Gamma }(P_1,P_2,P_3)`$ and we write $`N(P_1):=N_1`$ etc. The second simplification occurs because of the fact that the equations we obtain contain a specific combination of four-point vertices which we define as,

$`\overline{M}_F=M_F+N_1M_{R4}+N_4M_{R1}`$ (117)
where each vertex carries arguments $`M(P_1,P_2,P_3,P_4)`$. It has been shown that this particular combination of four-point vertices has special properties in a number of other contexts .

In addition, bare vertices that carry a Keldysh index ‘2’ have an extra factor of minus one. This factor is accounted for by inserting a factor of $`\tau `$ for each vertex where $`\tau `$ is the two component vector given by,

$`\tau =\left({\displaystyle \genfrac{}{}{0pt}{}{1}{1}}\right).`$ (118)
We want to obtain a perturbative expansion for the correlation functions of composite operators $`D_R(X,Y)`$ and $`G_{R1}(X,Y,Z)`$ which appear in (76) and (99). We obtain expressions for these correlation functions in terms of the vertices $`\mathrm{\Gamma }_{ij}`$ and $`M_{ijlm}`$ which are defined as the vertices obtained by truncating external legs from the following connected vertices:

$`\mathrm{\Gamma }_{ij}^C=T_c\pi _{ij}(X)\varphi (Y)\varphi (Z)`$ (119)
$`M_{ijlm}^C=T_c\pi _{ij}(X)\pi _{lm}(Y)\varphi (Z)\varphi (W)`$ (120)
where

$`\pi _{ij}(X)=_i\varphi (X)_j\varphi (X){\displaystyle \frac{1}{3}}\delta _{ij}(_m\varphi (X))(_m\varphi (X))`$ (121)
and $`T_c`$ is the operator that time orders along the closed time path contour. These definitions allow us to write the two- and three-point correlation functions as integrals of the form,

(122)
$`D_{ab}(Q)=2i{\displaystyle 𝑑P\mathrm{\Gamma }_{cad}^{ij}(P,Q,PQ)iD_{bc}(P)iD_{db}(P+Q)I_{ji}(p)\tau _b}`$ (123)
$`G_{abc}(QQ^{},Q,Q^{})=8{\displaystyle 𝑑KM_{ebcd}^{ijjk}\tau _aiD_{ae}(K)iD_{da}(K+Q+Q^{})I_{ki}(k)}`$ (124)
where the indices $`\{a,b,c,d,e\}`$ are Keldysh indices and take values $`\{1,2\}`$. These expressions are shown diagrammatically in Fig. . We perform the sum over Keldysh indices using the Mathematica program described in . We obtain,

$`D_R(Q)=i{\displaystyle 𝑑K(N_{k+q}N_k)\mathrm{\Gamma }_{R2}^{ij}(K,Q,KQ)D_A(K)D_R(K+Q)I_{ji}(k)}`$ (125)
$`G_{R1}(QQ^{},Q,Q^{})`$ (126)
$`=4{\displaystyle 𝑑K(\overline{M}_F)_{ikkj}(K,Q,Q^{},KQQ^{})D_R(K)D_A(K+Q+Q^{})I_{ji}(k)}`$ (127)
To rewrite these expressions in a simpler form we use rotational invariance to write,

$`M_{ijlm}:=\widehat{I}_{ij}\widehat{I}_{lm}M;\mathrm{\Gamma }_{ij}:=\widehat{I}_{ij}\mathrm{\Gamma }`$ (128)
Using an obvious notation we write the pairs of propagators $`D_A(P)D_R(P+Q)`$ and $`D_R(K)D_A(K+Q+Q^{})`$ as $`a_pr_{p+q}`$ and $`r_ka_{k+q+q^{}}`$. When $`\{q_0,q_0^{}\}0`$ the dominant contribution to the integral from these pairs of propagators is produced by what is called the pinch effect: the contour is “pinched” between the poles of the two propagators which gives rise to a factor in the denominator that is proportional to the imaginary part of the propagators. We regulate the pinching singularity with the imaginary part of the hard thermal loop self energy and obtain ,

$`r_ka_{k+q}{\displaystyle \frac{\rho _k}{2\mathrm{I}\mathrm{m}\mathrm{\Sigma }_\mathrm{k}}};\rho _k=i(r_ka_k)`$ (129)
where $`\mathrm{\Sigma }`$ is the retarded part of the hard thermal loop self energy.

Now we expand in $`q_0`$ and $`q_0^{}`$. In (125) we keep terms proportional to $`q_0`$ since these terms are the only ones that contribute to (76); in (126) we keep terms proportional to $`q_0q_0^{}`$ since these terms are the only ones that contribute to (99). In each term there are products of thermal factors of the form $`N_xN_{x+q}`$ and $`N_xN_{x+q^{}}`$ where $`x`$ is some combination of $`\{k,p,r\}`$. The expansion of these thermal factors is straightforward:

$`N_xN_{x+q}=2q_0\beta n_x(1+n_x)+\mathrm{}`$ (130)
Consider the behaviour of the vertices when $`\{q_0,q_0^{}\}0`$. Using (125) and (129) in (76) it is easy to see that only the real part of $`\mathrm{\Gamma }_{R2}`$ contributes. Similarly, using (126) and (129) in (99) it is clear that we need only the real part of $`\overline{M}_F`$. First we consider (125). Because of the explicit factor of $`(N_{k+q}N_k)`$ we can set $`Q`$ to zero in the vertex. The expansion of (126) in $`Q`$ and $`Q^{}`$ is more difficult. One can show that the terms in the expansion that contain gradients acting on the vertices $`M_F,M_{R1}`$ and $`M_{R4}`$ do not contribute to the result and that we can make the replacement:

$`\mathrm{Re}\overline{M}_F(K,Q,Q^{},KQQ^{})(N_kN_{k+q+q^{}})\mathrm{Re}M_{R1}(K,0,0,K)`$ (131)
From now on, to simplify the notation, we define $`\mathrm{\Gamma }(K,0,K):=\mathrm{\Gamma }(K)`$ and $`M(K,0,0,K):=M(K)`$.

Using these results to simplify (125) and (126) and substituting into (76) and (99) we obtain,

$`\eta ^{(1)}={\displaystyle \frac{\beta }{15}}{\displaystyle 𝑑Kk^2\rho _kn_k(1+n_k)\left[\frac{\mathrm{Re}\mathrm{\Gamma }_{\mathrm{R2}}(\mathrm{K})}{\mathrm{Im}\mathrm{\Sigma }_\mathrm{k}}\right]}`$ (132)
$`\eta ^{(2)}={\displaystyle \frac{2\beta ^2}{105}}{\displaystyle 𝑑Kk^2\rho _kn_k(1+n_k)N_k\left[\frac{\mathrm{ReM}_{\mathrm{R1}}(\mathrm{K})}{\mathrm{Im}\mathrm{\Sigma }_\mathrm{k}}\right]}`$ (133)
Comparing with (30) and (31) we see that the results are identical if we identify

$`B(\underset{¯}{k})={\displaystyle \frac{\mathrm{Re}\mathrm{\Gamma }_{\mathrm{R2}}(\underset{¯}{\mathrm{k}})}{\mathrm{Im}\mathrm{\Sigma }_\mathrm{k}}}`$ (134)
$`C(\underset{¯}{k})={\displaystyle \frac{\mathrm{ReM}_{\mathrm{R1}}(\underset{¯}{\mathrm{k}})}{\mathrm{Im}\mathrm{\Sigma }_\mathrm{k}}}`$ (135)
with the momentum $`K`$ on the shifted mass shell: $`\delta (K^2m_{th}^2)`$ where $`m_{th}^2=m^2+\mathrm{Re}\mathrm{\Sigma }_K`$.

#### 1 The Ladder Resummation

It has been known for some time that the set of diagrams which give the dominate contributions to the vertex $`\mathrm{\Gamma }_{ij}(P,Q,PQ)`$ are the ladder diagrams. These diagrams contribute to the viscosity to the same order in perturbation theory as the bare one loop graph and thus need to be resummed . This effect occurs for the following reason. It appears that the ladder graphs are suppressed relative to the one loop graph by extra powers of the coupling which come from the extra vertex factors that one obtains when one adds rungs (vertical lines). However, these extra factors of the coupling are compensated for by a kinematical factor. This factor arises through the pinch effect which is described above. The addition of an additional rung in a ladder graph always produces an extra pair of propagators of the form $`a_kr_{k+q}`$ or $`r_ka_{k+q}`$. Products of this form contribute a factor which produces an enhancement. This factor occurs when the contour is “pinched” between the poles of the two propagators, which gives rise to a contribution in the denominator that is proportional to the imaginary part of the inverse propagators.

In order to include ladder diagrams we obtain the vertex $`\stackrel{~}{\mathrm{\Gamma }}_{ij}(P,Q,PQ)`$ as the solution to the integral equation shown in Fig. .

(136)
$`\stackrel{~}{\mathrm{\Gamma }}_{abc}^{lm}(K,Q,KQ)`$ (137)
$`={\displaystyle 𝑑P𝑑R\tau _3^b\tau _3^a\tau _3^cI_{lm}(p)D_{ca}(R)D_{ac}(K+RP)D_{bc}(P+Q)D_{ab}(P)}`$ (138)
$`+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑P𝑑R\stackrel{~}{\mathrm{\Gamma }}_{dbe}^{lm}(P,Q,PQ)D_{ad}(P)D_{ec}(P+Q)D_{ac}(R+KP)\tau _3^cD_{ca}(R)\tau _3^a}`$ (139)
We perform the sums over the Keldysh indices using the Mathematica program in  and simplify the result by taking $`Q`$ to zero, keeping only the pinching terms, and using (129). The vertex function which includes the tree vertex is obtained by shifting: $`\mathrm{\Gamma }_{ij}=I_{ij}+\frac{1}{2}\stackrel{~}{\mathrm{\Gamma }}_{ij}`$. (The factor of 1/2 is a symmetry factor). We obtain,

$`\mathrm{\Gamma }_{R2}^{lm}(K)=`$ $`k_mk_l{\displaystyle \frac{1}{3}}\delta _{ml}k^2{\displaystyle \frac{\lambda ^2}{4}}{\displaystyle 𝑑P𝑑R𝑑P^{}}`$ (141)
$`(2\pi )^4\delta ^4(P+P^{}RK)\rho _p\rho _r\rho _p^{}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P)}{\mathrm{Im}\mathrm{\Sigma }_R(P)}}(1+n_p)(1+n_p^{})n_r/(1+n_k)`$
Note that this integral equation is decoupled: the only three-point vertex that appears is $`\mathrm{\Gamma }_{R2}^{lm}`$. To simplify this expression further we use (128), (134) and the fact that $`\mathrm{\Gamma }_{R2}(P)`$ is pure real, and symmetrize the integral on the right hand side over the integration variables $`\{P,P^{},R\}`$. We multiply and divide the left hand side by $`\mathrm{Im}\mathrm{\Sigma }_\mathrm{k}`$ and replace this expression in the numerator by the HTL result ,

$`\mathrm{Im}\mathrm{\Sigma }_k=`$ $`{\displaystyle \frac{\lambda ^2}{12}}\left({\displaystyle \frac{1}{1+n_k}}\right){\displaystyle 𝑑P𝑑R𝑑P^{}(2\pi )^4\delta ^4(P+P^{}RK)}`$ (143)
$`\rho _p\rho _r\rho _p^{}(1+n_p)(1+n_p^{})n_r`$
Rearranging we obtain ,

$`I(k,k)_{lm}={\displaystyle \frac{\lambda ^2}{12(1+n_k)}}{\displaystyle 𝑑P𝑑R𝑑P^{}(2\pi )^4\delta (P+P^{}RK)\rho _p\rho _r\rho _p^{}}`$ (144)
$`[B_{lm}(P)+B_{lm}(P^{})B_{lm}(K)B_{lm}(R)](1+n_p)(1+n_p^{})n_r`$ (145)
where we have used (21), (128) and (134). When the delta functions are used to do the frequency integrals, this equation has exactly the same form as the equation obtained from the linearized Boltzmann equation (38) with a shifted mass shell describing effective thermal excitations. Comparing (30) and (38) with (132) and (144) we conclude that calculating shear viscosity using effective transport theory by keeping only first order terms in the gradient expansion is equivalent to using the Kubo formula obtained from linear response theory, with a three-point vertex obtained by resumming ladder graphs. The contributions to the viscosity are shown in Fig. .

(146)
#### 2 The Extended-Ladder Resummation

In this section we consider an integral equation for the vertex $`M_{ijlm}`$ which resums an infinite set of graphs that includes ladder graphs and some other contributions which we will call extended ladder graphs. We will show that this integral equation has exactly the same form as the integral equation (48) obtained by expanding the Boltzmann equation to second order. We consider the integral equation shown in Fig. .

(147)
Following (117) and (126) we calculate contributions to $`M_F`$, $`M_{R1}`$, $`M_{R4}`$ for each diagram. We introduce some notation to simplify the equations. For four-point vertices we list the first three momenta only and for three-point vertices we list the first two only. In both cases the last momentum is the one that is determined by the conservation of energy. For example: $`M(K,Q,Q^{},KQQ^{}):=M(K,Q,Q^{})`$ and $`\mathrm{\Gamma }(K,Q,KQ):=\mathrm{\Gamma }(K,Q)`$. In addition we write (as before) $`M(K,0,0,K):=M(K)`$ and $`\mathrm{\Gamma }(K,0,K):=\mathrm{\Gamma }(K)`$. From Fig. \[4a\] we have,

$`(M_F^{ijlm})^{4a}(K,Q,Q^{})=2i[r_{k+q}\mathrm{\Gamma }_{R1}^{ij}(K,Q)\mathrm{\Gamma }_{F2}^{lm}(K+Q,Q^{})`$ (148)
$`+a_{k+q}\mathrm{\Gamma }_{F2}^{ij}(K,Q)\mathrm{\Gamma }_{R3}^{lm}(K+Q,Q^{})+f_{k+q}\mathrm{\Gamma }_{R1}^{ij}(K,Q)\mathrm{\Gamma }_{R3}^{lm}(K+Q,Q^{})]`$ (149)
$`(M_{R1}^{ijlm})^{4a}(K,Q,Q^{})=2ir_{k+q}\mathrm{\Gamma }_{R1}^{ij}(K,Q)\mathrm{\Gamma }_{R1}^{lm}(K+Q,Q^{})`$ (150)
$`(M_{R4}^{ijlm})^{4a}(K,Q,Q^{})=2ia_{k+q}\mathrm{\Gamma }_{R3}^{ij}(K,Q)\mathrm{\Gamma }_{R3}^{lm}(K+Q,Q^{})`$ (151)
We can rewrite this result by substituting in the expanded form of the vertex $`\mathrm{\Gamma }_{ij}=I_{ij}+\frac{1}{2}\stackrel{~}{\mathrm{\Gamma }}_{ij}`$ with $`\stackrel{~}{\mathrm{\Gamma }}_{ij}`$ obtained from the integral equation that corresponds to Fig. . We use,

$`\mathrm{\Gamma }_{R1}^{ij}(K,Q)`$ (152)
$`=I^{ij}(k){\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Rr_pa_{p+q}(N^{}N_r)a_rr^{}[\mathrm{\Gamma }_{F2}^{ij}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{ij}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{ij}(P,Q)]}`$ (153)
$`\mathrm{\Gamma }_{R3}^{ij}(K,Q)`$ (154)
$`=I^{ij}(k){\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Rr_pa_{p+q}(N_rN^{})r_ra^{}[\mathrm{\Gamma }_{F2}^{ij}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{ij}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{ij}(P,Q)]}`$ (155)
$`\mathrm{\Gamma }_{F2}^{ij}(K,Q)`$ (156)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Rr_pa_{p+q}(1N^{}N_r)\rho _r\rho ^{}[\mathrm{\Gamma }_{F2}^{ij}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{ij}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{ij}(P,Q)]}`$ (157)
It is easy to show that the contribution from Fig. \[4a\] is equivalent to the contributions from the four diagrams in Fig. .

(158)
From Figs. \[5a,5b\] we obtain,

$`(M_F^{ijlm})^{5a}(K,Q,Q^{})=iI_{ij}(k)[a_{k+q}\mathrm{\Gamma }_{F2}^{lm}(K,Q)+f_{k+q}\mathrm{\Gamma }_{R1}^{lm}(K,Q)]`$ (159)
$`(M_{R1}^{ijlm})^{5a}(K,Q,Q^{})=ir_{k+q}I_{ij}(k)\mathrm{\Gamma }_{R1}^{lm}(K,Q)`$ (160)
$`(M_{R4}^{ijlm})^{5a}(K,Q,Q^{})=ia_{k+q}I_{ij}(k)\mathrm{\Gamma }_{R3}^{lm}(K,Q)`$ (161)
$`(M_F^{ijlm})^{5b}(K,Q,Q^{})=iI_{ij}(k)[r_{k+q}\mathrm{\Gamma }_{F2}^{lm}(K+Q,Q^{})+f_{k+q}\mathrm{\Gamma }_{R3}^{lm}(K,Q)]`$ (162)
$`(M_{R1}^{ijlm})^{5b}(K,Q,Q^{})=ir_{k+q}I_{ij}(k)\mathrm{\Gamma }_{R1}^{lm}(K+Q,Q^{})`$ (163)
$`(M_{R4}^{ijlm})^{5b}(K,Q,Q^{})=ia_{k+q}I_{ij}(k)\mathrm{\Gamma }_{R3}^{lm}(K+Q,Q^{})`$ (164)
From Figs. \[5c,5d\] we obtain,

$`(M_F^{ijlm})^{5c}(K,Q,Q^{})`$ (165)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p}`$ (166)
$`[\mathrm{\Gamma }_{F2}^{lm}(P+Q,Q^{})N_{p+q+q^{}}\mathrm{\Gamma }_{R1}^{lm}(P+Q,Q^{})+N_{p+q}\mathrm{\Gamma }_{R3}^{lm}(P+Q,Q^{})]`$ (167)
$`[\mathrm{\Gamma }_{F2}^{ij}(K,Q)(f_ra_p^{}+r_rf_p^{})a_{k+q}`$ (168)
$`+\mathrm{\Gamma }_{R1}^{ij}(K,Q)(f_{k+q}(f_ra_p^{}+r_rf_p^{})+r_{k+q}(f_p^{}f_r+a_p^{}r_r+r_p^{}a_r))]`$ (169)
(170)
$`(M_{R1}^{ijlm})^{5c}(K,Q,Q^{})`$ (171)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p}`$ (172)
$`[\mathrm{\Gamma }_{F2}^{lm}(P+Q,Q^{})N_{p+q+q^{}}\mathrm{\Gamma }_{R1}^{lm}(P+Q,Q^{})+N_{p+q}\mathrm{\Gamma }_{R3}^{lm}(P+Q,Q^{})]`$ (173)
$`\mathrm{\Gamma }_{R1}^{ij}(K,Q)r_{k+q}(f_p^{}a_r+r_p^{}f_r)`$ (174)
(175)
$`(M_{R4}^{ijlm})^{5c}(K,Q,Q^{})`$ (176)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p}`$ (177)
$`[\mathrm{\Gamma }_{F2}^{lm}(P+Q,Q^{})N_{p+q+q^{}}\mathrm{\Gamma }_{R1}^{lm}(P+Q,Q^{})+N_{p+q}\mathrm{\Gamma }_{R3}^{lm}(P+Q,Q^{})]`$ (178)
$`\mathrm{\Gamma }_{R3}^{ij}(K,Q)a_{k+q}(f_p^{}r_r+a_p^{}f_r).`$ (179)
(180)
$`(M_F^{ijlm})^{5d}(K,Q,Q^{})`$ (181)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q}r_p[\mathrm{\Gamma }_{F2}^{lm}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{lm}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{lm}(P,Q)]}`$ (182)
$`[\mathrm{\Gamma }_{F2}^{ij}(K+Q,Q^{})r_{k+q}(a_rf^{}+f_rr^{})`$ (183)
$`+\mathrm{\Gamma }_{R3}^{ij}(K+Q,Q^{})(f_{k+q}(a_rf^{}+f_rr^{})+a_{k+q}(f_rf^{}+r_ra^{}+a_rr^{}))]`$ (184)
(185)
$`(M_{R1}^{ijlm})^{5d}(K,Q,Q^{})`$ (186)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Rr_pa_{p+q}[\mathrm{\Gamma }_{F2}^{lm}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{lm}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{lm}(P,Q)]}`$ (187)
$`\mathrm{\Gamma }_{R1}^{ij}(K+Q,Q^{})r_{k+q}(a_rf^{}+f_rr^{})`$ (188)
(189)
$`(M_{R4}^{ijlm})^{5d}(K,Q,Q^{})`$ (190)
$`={\displaystyle \frac{i\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Rr_pa_{p+q}[\mathrm{\Gamma }_{F2}^{lm}(P,Q)N_{p+q}\mathrm{\Gamma }_{R1}^{lm}(P,Q)+N_p\mathrm{\Gamma }_{R3}^{lm}(P,Q)]}`$ (191)
$`\mathrm{\Gamma }_{R3}^{ij}(K+Q,Q^{})a_{k+q}(a^{}f_r+f^{}r_r)`$ (192)
The contribution from Fig. \[4b\] is

$`(M_F^{ijlm})^{4b}(K,Q,Q^{})={\displaystyle \frac{\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p\overline{M}_F^{ijlm}(P,Q,Q^{})[r_ra_p^{}+a_rr_p^{}+f_rf_p^{}]}`$ (193)
$`(M_{R1}^{ijlm})^{4b}(K,Q,Q^{})={\displaystyle \frac{\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p\overline{M}_F^{ijlm}(P,Q,Q^{})[a_rf_p^{}+f_rr_p^{}]}`$ (194)
$`(M_{R4}^{ijlm})^{4b}(K,Q,Q^{})={\displaystyle \frac{\lambda ^2}{8}}{\displaystyle 𝑑P𝑑Ra_{p+q+q^{}}r_p\overline{M}_F^{ijlm}(P,Q,Q^{})[r_rf_p^{}+f_ra_p^{}]}`$ (195)
where $`P^{}=R+KP`$. Fig. \[4c\] gives,

$`(M_F^{ijlm})^{4c}(K,Q,Q^{})=`$ $`{\displaystyle \frac{i\lambda ^2}{4}}{\displaystyle 𝑑P𝑑R(N_pN_{p+q})(N_{rq^{}}N_r)}`$ (197)
$`\mathrm{\Gamma }_{R2}^{ij}(P,Q)\mathrm{\Gamma }_{R2}^{lm}(RQ^{},Q^{})f_p^{}r_pa_{r+q}a_rr_{rq^{}}`$
$`(M_{R1}^{ijlm})^{4c}(K,Q,Q^{})=`$ $`{\displaystyle \frac{i\lambda ^2}{4}}{\displaystyle 𝑑P𝑑R(N_pN_{p+q})(N_{rq^{}}N_r)\mathrm{\Gamma }_{R2}^{ij}(P,Q)}`$ (199)
$`\mathrm{\Gamma }_{R2}^{lm}(RQ^{},Q^{})r_r^{}r_pa_{r+q}a_rr_{rq^{}}`$
$`(M_{R4}^{ijlm})^{4c}(K,Q,Q^{})=`$ $`{\displaystyle \frac{i\lambda ^2}{4}}{\displaystyle 𝑑P𝑑R(N_pN_{p+q})(N_{rq^{}}N_r)\mathrm{\Gamma }_{R2}^{ij}(P,Q)}`$ (201)
$`\mathrm{\Gamma }_{R2}^{lm}(RQ^{},Q^{})a_p^{}r_pa_{r+q}a_rr_{rq^{}}`$
To combine these expressions we use (116) and (117). We keep only the pinching contributions and use (129) to regulate. As discussed previously, we expand in $`q_0`$ and $`q_0^{}`$ and keep the term proportional to $`q_0q_0^{}`$, since that is the only term that will contribute to the quadratic shear viscous response coefficient. As before, one can show that the terms in the expansion that contain gradients acting on the vertices $`M_F,M_{R1}`$ and $`M_{R4}`$ do not contribute to the result and that we can use (131). We simplify further by using (130). We also make repeated use of the set of identities below which hold for momenta which satisfy $`P_1+P_2=P_3+P_4`$:

$`n_1n_2(1+n_3)(1+n_4)=(1+n_1)(1+n_2)n_3n_4`$ (202)
$`n_3n_2+n_3+n_1(n_3n_2)=(1+n_1)(1+n_2)n_3/(1+n_4)`$ (203)
$`n_1(1+n_1)n_3(1+n_3)(N_2N_4)=(1+n_1)(1+n_2)n_3n_4(N_3N_1).`$ (204)
Finally, since $`K`$ is the momentum for an external leg we take the on shell piece:

$$r_k\frac{i}{\mathrm{Im}\mathrm{\Sigma }_k}$$
(205)
We obtain:

$`N_kn_k(1+n_k)M_{R1}^{ijlm}(K)`$ (206)
$`=n_k(1+n_k)N_kI_{ij}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}+{\displaystyle \frac{\lambda ^2}{4}}{\displaystyle 𝑑P𝑑R(1+n_p)(1+n_p^{})n_rn_k\rho _p\rho _p^{}\rho _r}`$ (207)
$`\left[{\displaystyle \frac{N_pM_{R1}^{ijlm}(P)}{\mathrm{Im}\mathrm{\Sigma }_p}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P)}{\mathrm{Im}\mathrm{\Sigma }_p}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(R)}{\mathrm{Im}\mathrm{\Sigma }_r}}(N_rN_p)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P)}{\mathrm{Im}\mathrm{\Sigma }_p}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}(N_pN_k)\right]`$ (208)
We introduce the symmetric notation: $`P:=P_1;P^{}:=P_2;R:=P_3`$ and rewrite the equation above after symmetrizing on the integration variables. Rearranging we obtain,

$`n_k(1+n_k)N_kI_{ij}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}`$ (209)
$`={\displaystyle \frac{\lambda ^2}{12}}{\displaystyle 𝑑P_1𝑑P_2𝑑P_3(2\pi )^4\delta ^4(P_1+P_2P_3K)(1+n_1)(1+n_2)n_3n_k\rho _1\rho _2\rho _3}`$ (210)
$`[{\displaystyle \frac{N_{p_3}M_{R1}^{ijlm}(P_3)}{\mathrm{Im}\mathrm{\Sigma }_{p_3}}}+{\displaystyle \frac{N_kM_{R1}^{ijlm}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}{\displaystyle \frac{N_{p_1}M_{R1}^{ijlm}(P_1)}{\mathrm{Im}\mathrm{\Sigma }_{p_1}}}{\displaystyle \frac{N_{p_2}M_{R1}^{ijlm}(P_2)}{\mathrm{Im}\mathrm{\Sigma }_{p_2}}}`$ (211)
$`+{\displaystyle \frac{1}{2}}\{(N_1+N_2){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P_1)}{\mathrm{Im}\mathrm{\Sigma }_{p_1}}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P_2)}{\mathrm{Im}\mathrm{\Sigma }_{p_2}}}(N_k+N_3){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P_3)}{\mathrm{Im}\mathrm{\Sigma }_{p_3}}}`$ (212)
$`+(N_3N_1){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P_1)}{\mathrm{Im}\mathrm{\Sigma }_{p_1}}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P_3)}{\mathrm{Im}\mathrm{\Sigma }_{p_3}}}+(N_kN_1){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P_1)}{\mathrm{Im}\mathrm{\Sigma }_{p_1}}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}`$ (213)
$`+(N_3N_2){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(P_3)}{\mathrm{Im}\mathrm{\Sigma }_{p_3}}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P_2)}{\mathrm{Im}\mathrm{\Sigma }_{p_2}}}+(N_kN_2){\displaystyle \frac{\mathrm{\Gamma }_{R2}^{ij}(K)}{\mathrm{Im}\mathrm{\Sigma }_k}}{\displaystyle \frac{\mathrm{\Gamma }_{R2}^{lm}(P_2)}{\mathrm{Im}\mathrm{\Sigma }_{p_2}}}\}]`$ (214)
Note that once again we have obtained an integral equation that is decoupled: it only involves $`M_{R1}`$ and $`\mathrm{\Gamma }_{R2}`$. With $`\mathrm{\Gamma }_{R2}`$ determined by the integral equation (144), Equation (214) can be solved to obtain $`M_{R2}`$. From the diagrams in Figs.  we see that the solutions to the integral equation will contain contributions of the form shown in Fig. . Finally, by using (21), (128), (134) and (135) and comparing (133) and (214) with (31) and (48) we see that calculating the quadratic shear viscous response using transport theory describing effective thermal excitations and keeping terms that are quadratic in the gradient of the four-velocity field in the expansion of the Boltzmann equation, is equivalent to calculating the same response coefficient from quantum field theory at finite temperature using the next-to-linear response Kubo formula with a vertex given by a specific integral equation. This integral equation shows that the complete set of diagrams that need to be resummed includes the standard ladder graphs, and an additional set of extended ladder graphs. Some of the diagrams that contribute to the viscosity are shown in Fig. .

(215)
(216)
## V Conclusions

We have studied nonlinear response using two different methods. The first method uses standard transport theory. We start from a local equilibrium form for the distribution function and perform a gradient expansion. We calculate the quadratic shear viscous response coefficient by expanding the Boltzmann equation and obtaining a hierarchy of equations that can be solved consistently. The second technique uses response theory. We work with a perturbative Hamiltonian that is linear in the gradient of the four-velocity field and study quadratic response. We generalize the Kubo formula for linear response and obtain an expression that allows us to calculate quadratic shear viscous response from the retarded three-point green function of the viscous shear stress tensor. The transport theory calculation involves the use of the Boltzmann equation which is itself obtained from some more fundamental theory. The response theory calculation uses the well known methods of equilibrium finite temperature quantum field theory and is, in this sense, more fundamental. However, the response theory calculation is complicated by the need to resum infinite sets of diagrams at finite temperature.

At leading order, it is well known that a correct calculation of the linear response coefficient involves the resummation of ladder graphs. Beyond leading order in response theory it is difficult even to identify which diagrams need to be resummed. We have identified precisely which diagrams need to be resummed by studying the connection between the transport theory calculation and the response theory calculation. We have shown that calculating the quadratic shear viscous response coefficient using transport theory by keeping terms that are quadratic in the gradient of the four-velocity field in the expansion of the Boltzmann equation, is equivalent to calculating the same response coefficient from quantum field theory at finite temperature using the next-to-linear response Kubo formula with a vertex given by a specific integral equation. This integral equation shows that the complete set of diagrams that need to be resummed includes the standard ladder graphs, and an additional set of extended ladder graphs.

There are several directions for future work. It has been shown that the Boltzmann equation can be derived from the Kadanoff-Baym equations by using a gradient expansion and keeping only linear terms . The connection between this result and the work discussed in this paper can probably be understood by studying the dual roles of the gradient expansion and the quasiparticle approximation. In addition, it would be interesting to generalize this work to gauge field theories.

###### Acknowledgements.

This work was supported by the Natural Sciences and Engineering and Research Council of Canada (NSERC), and the National Natural Science Foundation of China (NSFC).


format: text

text:
# 1 Hadron mass corrections.

## 1 Hadron mass corrections.

The leading twist approximation to hard processes in QCD is affected by a number of multiplicative and additive corrections. To the first class one obviously attributes radiative corrections in the strong coupling, while the latter encodes higher twist contributions, which provide power suppressed corrections. In the case when a given large scale $`Q`$, that controls the factorization of a reaction, becomes rather low, one has to take care of power suppressed effects since they can significantly modify the scaling behaviour and the magnitude of the leading twist prediction for the corresponding cross section. Power corrections can be divided into two classes according to their origin: dynamical and kinematical ones. The first one comprises of multiparton correlations inside hadrons  and gives rise to new non-perturbative functions. The second one arises from a separation of composite operators into parts that have definite transformation properties w.r.t. the Lorentz group and thus possess a well defined geometrical twist, i.e. dimension minus spin. This decomposition provides the so-called Wandzura-Wilczek type contributions  arising from components having definite symmetry properties as well as target mass corrections  stemming from the subtraction of trace terms in the afore mentioned operators.

In the present paper we develop a formalism for the resummation of the target mass corrections for off-forward two-photon processes whose amplitudes are given by a chronological product of electromagnetic currents. This method is indispensable for a study of hadron mass effects in the cross section of e.g. deeply virtual Compton scattering , — the processes sensitive to the generalized parton distributions (GPDs). The present day facilities can offer the $`Q^2`$ of the order of a few (2-4) $`\text{GeV}^2`$ (for the reactions in question) and one has to account for power suppressed contributions to the twist-two and -three observables arising from the ratio of the nucleon mass to the hard momentum transfer. Unfortunately, the formalism for resummation of target mass corrections invented by Georgi and Politzer  for the deep inelastic scattering (see also  for a first discussion of the topic), and accepted in all consequent generalizations and applications to forward electroweak scattering amplitudes, is not directly applicable to off-forward processes in the context of GPDs. The reason for that is the appearance of towers of new Lorentz structures in the matrix elements of local operators, arisen due to the non-zero $`t`$-channel momentum transfer, which cannot be handled in the fashion proposed in . Note that a discussion of mass effects in exclusive processes with a much simpler kinematics can be found in Ref. .

In the present context, instead of dealing with single variable functions, e.g. GPDs, we use spectral functions, the so-called double distributions (DDs) , which depend on two variables. This allows to resum the mass corrections in a straightforward manner and to derive all-order results in a compact form. After this is accomplished one can in principle use an inverse integral representation to express DDs by means of GPDs, however, this transformation requires the knowledge of the support of GPDs including an unphysical region. Therefore, this transformation can be hardly used for a successful phenomenological application.

The consequent presentation runs as follows. In the next section we give a short introduction to a procedure of the trace subtraction by means of harmonic polynomials. In Section 3 we resum the target mass effects resulting from the trace subtraction in matrix elements of twist-two operators for zero and spin-$`\frac{1}{2}`$ targets. Then it is shown that the expansion of these results in powers of $`\left(M^2/q^2\right)^j`$ can be cast in a conventional form of GPDs and be used for numerical predictions of target mass effects in physical cross sections. Finally, we conclude.

## 2 Harmonic polynomials and twist decomposition.

The goal of our study is an evaluation of the off-forward Compton amplitude keeping target mass effects stemming from the twist-two contributions in the light-cone expansion of the off-forward matrix element of the chronological product $`\stackrel{~}{T}_{\mu \nu }(\frac{x}{2},\frac{x}{2})P_2|T\left\{j_\mu (\frac{x}{2})j_\nu (\frac{x}{2})\right\}|P_1`$ of electromagnetic currents $`j_\mu =\overline{\psi }\gamma _\mu \psi `$. Since these are entirely kinematical effects, we consider dynamical higher twist effects given by multiparticle operators as independent and, therefore, we safely neglect them. Note, however, that the trace subtraction in twist-three contributions in the Wandzura-Wilczek approximation will generate kinematical mass effects accompanied by leading twist operators. We will comment on this issue in the concluding section. In the approximation, we are considering, the Fourier transformed two-photon amplitude reads

$`T_{\mu \nu }`$ $`=`$ $`{\displaystyle d^4x\mathrm{e}^{ixq}\stackrel{~}{T}_{\mu \nu }(\frac{x}{2},\frac{x}{2})}={\displaystyle \frac{1}{\pi ^2}}{\displaystyle d^4x\frac{x_\sigma }{[x^2+i0]^2}}`$ (1)
$`\times `$ $`P_2|\mathrm{cos}\left(xq\right)S_{\mu \nu ;\rho \sigma }{}_{}{}^{V}𝒪_{\rho }^{}(\frac{x}{2},\frac{x}{2})+\mathrm{sin}\left(xq\right)ϵ_{\mu \nu \rho \sigma }{}_{}{}^{A}𝒪_{\rho }^{}(\frac{x}{2},\frac{x}{2})|P_1,`$
where $`S_{\mu \nu ;\rho \sigma }=g_{\mu \rho }g_{\nu \sigma }+g_{\mu \sigma }g_{\nu \rho }g_{\mu \nu }g_{\rho \sigma }`$ and for the totally antisymmetric tensor we use the normalization $`ϵ^{0123}=1`$. The vector $`q`$ is a semi-sum of incoming $`q_{1\nu }`$ and outgoing $`q_{2\mu }`$ photon momenta, $`q=\frac{1}{2}(q_1+q_2)`$. The non-local operators in Eq. (1) are $`{}_{}{}^{\mathrm{\Gamma }}𝒪_{\rho }^{}(x,x)=\overline{\psi }(x)[x,x]\mathrm{\Gamma }_\rho \psi (x)`$ with a Dirac matrix $`\mathrm{\Gamma }_\rho `$ which equals $`\gamma _\rho `$ or $`\gamma _\rho \gamma _5`$ for vector ($`V`$) or axial-vector ($`A`$) sectors, respectively. The path ordered exponential $`[x,x]`$ will be dropped everywhere in the following formulae.

We are aiming now in a decomposition of $`𝒪_\rho `$ into traceless operators $`_\rho ^2`$ and $`_\rho ^3`$ with definite geometrical twist-two and -three, respectively:

$$𝒪_\rho (x,x)=_\rho ^2(x,x)+_\rho ^3(x,x)+_\rho ^r(x,x).$$
(2)
To this end, we expand first the non-local operator $`𝒪_\rho (x,x)`$ in a Taylor series

$`𝒪_\rho (x,x)={\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i)^j}{j!}}x_{\mu _1}\mathrm{}x_{\mu _j}𝒪_{\rho ;\mu _1\mathrm{}\mu _j},\text{with}𝒪_{\rho ;\mu _1\mathrm{}\mu _j}=\overline{\psi }\mathrm{\Gamma }_\rho i\underset{\mu _1}{\overset{}{𝒟}}\mathrm{}i\underset{\mu _j}{\overset{}{𝒟}}\psi ,`$ (3)
where $`\underset{\mu }{\overset{}{𝒟}}=\underset{\mu }{\overset{}{𝒟}}\underset{\mu }{\overset{}{𝒟}}`$ with the covariant derivative $`𝒟_\mu =_\mu igB_\mu `$, and then we perform the twist decomposition for the local operators $`𝒪_{\rho ;j}=x_{\mu _1}\mathrm{}x_{\mu _j}𝒪_{\rho ;\mu _1\mathrm{}\mu _j}`$. An effective way of trace subtraction is to solve the sufficient condition $`^2H^j(x^2,^2)T_j(x)=0`$ , where $`T_j`$ results from the contraction of a rank-$`j`$ tensor with the product of four-vectors $`x_{\mu _1}\mathrm{}x_{\mu _j}`$, and $`H^j`$ is a subtraction operator. The solution is given by the so-called harmonic polynomials

$`H^j(x^2,^2)={\displaystyle \underset{k=0}{\overset{[j/2]}{}}}{\displaystyle \frac{\mathrm{\Gamma }\left(jk+1\right)}{k!\mathrm{\Gamma }\left(j+1\right)}}\left({\displaystyle \frac{x^2}{4}}\right)^k\left(^2\right)^k.`$ (4)
For operators having a free vector index one has an extra condition so that in total the tracelessness implies $`^2_{\rho ;j}^i=0`$ and $`_\rho _{\rho ;j}^i=0`$, for $`i=2,3`$, where the twist-two operators are completely symmetric while the twist-three operators are antisymmetric in a pair of their Lorentz indices. Employing the equation of motion for the massless Dirac field, one can show that $`^r`$ is proportional to total derivative, $``$, squared, apart from multiparticle operators, $`_\rho ^r^2_\rho ^2+\mathrm{multiparticle}\mathrm{operators}`$. When it is sandwiched between hadronic states, it generates contributions proportional either to the $`t`$-channel momentum transfer squared or given by multi-parton correlations functions. Both of them will be neglected in the following. Thus, our considerations will be limited by the condition $`M^2\mathrm{\Delta }^2`$, with $`M^2=P_1^2=P_2^2`$ and $`\mathrm{\Delta }^2`$, being the hadron mass and $`t`$-channel momentum transfer, respectively.

A straightforward algebra leads to the desired decomposition in operators with definite geometrical twist:

$$𝒪_{\rho ;j}=_{\rho ;j}^2+\frac{2j}{j+1}_{\rho ;j}^3+_{\rho ;j}^r,$$
(5)
where the traceless twist-two and -three operators read (cf. )

$`_{\rho ;j}^2`$ $`=`$ $`{\displaystyle \frac{1}{j+1}}_\rho H^{j+1}(x^2,^2)x_\sigma 𝒪_{\sigma ;j},j0,`$ (6)
$`_{\rho ;j}^3`$ $`=`$ $`{\displaystyle \frac{1}{2j}}\left[g_{\rho \sigma }xx_\sigma _\rho \right]\left[g_{\sigma \tau }{\displaystyle \frac{1}{j+1}}x_\sigma _\tau \right]H^j(x^2,^2)𝒪_{\tau ;j},j1,`$ (7)
and the remainder is

$`_{\rho ;j}^r`$ $`=`$ $`{\displaystyle \frac{1}{j+1}}\{_\rho [1H^{j+1}(x^2,^2)]x_\sigma +[g_{\rho \sigma }xx_\sigma _\rho ][1H^j(x^2,^2)]`$
$`+{\displaystyle \frac{1}{j+1}}[x_\rho xx^2_\rho ]_\sigma H^j(x^2,^2)\}𝒪_{\sigma ;j}.`$
It is obvious that the twist-two operators (6) are traceless. Using $`^2H^j(x^2,^2)𝒪_{\mu ;j}=0`$ and Euler theorem $`x_\mu H^j(x^2,^2)𝒪_{\mu ;j}=(j1)_\mu H^j(x^2,^2)𝒪_{\mu ;j}`$, it is an easy task to check that the operators (7) are traceless too.

## 3 Resummation within DD formalism.

In this section we will give a general framework for the resummation of kinematical mass corrections. For simplicity we consider here only the twist-two sector. The major steps consist of: a parametrization of the off-forward matrix elements of symmetric local operators by means of moments of DDs, the trace subtraction with the harmonic projectors $`H^j`$, then a Fourier transformation to the momentum space and, finally, a resummation of infinite series.

In order to demonstrate the main features of our formalism, let us consider first matrix elements of twist-two operators sandwiched between states of spinless target $`P_2|_{\rho ;j}^2|P_1`$. Here only the matrix elements of the parity even sector are relevant and read

$$P_2|{}_{}{}^{V}𝒪_{\rho ;\mu _1\mathrm{}\mu _j}^{\mathrm{sym}}|P_1=P_{\{\rho }P_{\mu _1}\mathrm{}P_{\mu _j\}}B_{j+1,j+1}+\mathrm{}+\mathrm{\Delta }_{\{\rho }\mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j\}}B_{j+1,0}+\mathrm{}.$$
(8)
They are build from the two vectors $`P=P_1+P_2`$ and $`\mathrm{\Delta }=P_2P_1`$, $`\{\mathrm{}\}`$ denotes symmetrization of corresponding indices and the ellipsis stands for possible terms containing the metric tensor. Since we act on this expression with $`H^j`$ we can neglect the latter. After projection with harmonic polynomials via Eq. (6) one gets

$`P_2|{}_{}{}^{V}_{\rho ;j}^{2}|P_1`$ $`=`$ $`{\displaystyle \frac{1}{j+1}}_\rho H^{j+1}(x^2,^2){\displaystyle \underset{k=0}{\overset{j+1}{}}}\left(xP\right)^{j+1k}\left(x\mathrm{\Delta }\right)^kB_{j+1,j+1k}`$
$`=`$ $`x_{\mu _1}\mathrm{}x_{\mu _j}\underset{\rho \mu _1\mathrm{}\mu _j}{\text{S}}\left\{P_\rho P_{\mu _1}\mathrm{}P_{\mu _j}B_{j+1,j+1}+\mathrm{}+\mathrm{\Delta }_\rho \mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j}B_{j+1,0}\right\},`$
with the operation $`\text{S}`$ standing for symmetrization and trace subtraction. The first line gives a compact expression for subtracted operator which is extremely convenient for consequent considerations and is the basis of our formalism.

On the other hand we can deduce the twist-two operator from Eq. (8) by a contraction of all Lorentz indices with a light-like vector $`n`$ so that the trace terms vanish identically. This tells us that the coefficients $`B_{j+1,j+1k}`$ in front of the Lorentz tensors are related to the reduced matrix elements of light-cone operators and can be represented in terms of moments of conventional leading twist GPDs or DDs via

$$B_{j,jk}=\frac{1}{k!}\frac{^k}{\eta ^k}|_{\eta =0}_1^1𝑑xx^{j1}B(x,\eta )=\left(\genfrac{}{}{0pt}{}{j}{k}\right)_\mathrm{\Omega }𝑑y𝑑zy^{jk}z^kf(y,z)$$
(10)
where $`0kj`$, $`1j`$ and the integration domain is $`\mathrm{\Omega }=\{1y1,1+|y|z1|y|\}`$. Consequently, the parametrization111This is not the definition given in Refs.  which suffered from inconsistencies due to omitted Lorentz structures as noted in . Our solution of this problem differs also from  where a two-component parametrization has been suggested. of the matrix element for the non-local symmetric operator and the light-ray operators coincide up to terms proportional to $`x_\rho `$ or $`x^2`$:

$$P_2|{}_{}{}^{V}𝒪_{\rho }^{\mathrm{sym}}(x,x)|P_1=_\mathrm{\Omega }𝑑y𝑑zf(y,z)𝒫_\rho \mathrm{e}^{ix𝒫}+\mathrm{},$$
(11)
where we introduced a new shorthand notation $`𝒫_\mu yP_\mu +z\mathrm{\Delta }_\mu `$.

The resummation of matrix elements of traceless local operators can be performed in the formalism of DD $`f(y,z)`$. By means of Eq. (10) it is easy to see that the finite sum w.r.t. $`k`$ in Eq. (3) provides

$$P_2|{}_{}{}^{V}_{\rho ;j}^{2}|P_1=_\mathrm{\Omega }𝑑y𝑑zf(y,z)\frac{1}{j+1}_\rho H^{j+1}(x^2,^2)\left(x𝒫\right)^{j+1}.$$
(12)
The evaluation of projection of harmonic polynomials results into conventional Chebyshev polynomials related to the irreducible representation of the orthogonal group $`SO(4)`$

$$H^j(x^2,^2)\left(x𝒫\right)^j=\left(\frac{x^2𝒫^2}{4}\right)^{j/2}U_j\left(\frac{x𝒫}{\sqrt{x^2𝒫^2}}\right).$$
(13)
This result when plugged into Eq. (12) gives our final expression for the twist-two traceless local operators. We will use it after we perform the Fourier transformation in Eq. (1). We accomplish it making use of the following general formula

$$d^4x\mathrm{e}^{ixq}\frac{x_{\mu _1}\mathrm{}x_{\mu _j}}{[x^2]^k}=(i)^{j+1}2^{42k+j}\pi ^2\frac{\mathrm{\Gamma }(2k+j)}{\mathrm{\Gamma }(k)}\frac{q_{\mu _1}\mathrm{}q_{\mu _j}}{[q^2]^{2k+j}},$$
(14)
where on the r.h.s. we dropped all terms involving the metric tensor $`g_{\mu _i\mu _k}`$, since they vanish when contracted with a traceless tensor. Thus, we get

$`{\displaystyle d^4x\mathrm{e}^{ixq}\frac{x_\sigma }{[x^2]^2}\frac{x_{\mu _1}}{2}\mathrm{}\frac{x_{\mu _j}}{2}P_2|{}_{}{}^{V}_{\rho ;\mu _1\mathrm{}\mu _j}^{2}|P_1}`$
$`=i^{j+2}\pi ^2\mathrm{\Gamma }(j+1){\displaystyle \frac{1}{q^2}}\mathrm{\Pi }_{\sigma \mu _1}{\displaystyle \frac{q_{\mu _2}}{q^2}}\mathrm{}{\displaystyle \frac{q_{\mu _j}}{q^2}}P_2|{}_{}{}^{V}_{\rho ;\mu _1\mathrm{}\mu _j}^{2}|P_1,`$ (15)
with projector

$$\mathrm{\Pi }_{\mu \nu }g_{\mu \nu }\frac{2}{q^2}q_\mu q_\nu ,\mathrm{\Pi }_{\mu \rho }\mathrm{\Pi }_{\rho \nu }=g_{\mu \nu }.$$
(16)
Then the Fourier transformation of the twist-two contribution to the Compton amplitude gives

$`{\displaystyle d^4x\mathrm{e}^{ixq}\frac{x_\sigma }{[x^2]^2}P_2|{}_{}{}^{V}_{\rho }^{2}(\frac{x}{2},\frac{x}{2})|P_1}`$
$`=2\pi ^2\left\{{\displaystyle \frac{q_\sigma }{q^2}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}P_2|{}_{}{}^{V}\stackrel{~}{}_{\rho ;j}^{2}|P_1{\displaystyle \frac{1}{2}}\mathrm{\Pi }_{\sigma \tau }{\displaystyle \frac{}{q_\tau }}{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{j}}P_2|{}_{}{}^{V}\stackrel{~}{}_{\rho ;j}^{2}|P_1\right\}.`$ (17)
Here we used a new convention

$$\stackrel{~}{}_{\rho ;j}^2=\frac{q_{\mu _1}}{q^2}\mathrm{}\frac{q_{\mu _j}}{q^2}_{\rho ;\mu _1\mathrm{}\mu _j}^2,$$
(18)
and exploited an identity

$`{\displaystyle \frac{}{\frac{q_\sigma }{q^2}}}=q^2\mathrm{\Pi }_{\sigma \tau }{\displaystyle \frac{}{q_\tau }}.`$
To evaluate the sums in (3) we substitute (13) in (12) and replace $`x_\mu `$ by $`q_\mu /q^2`$. This gives the local traceless matrix elements $`P_2|{}_{}{}^{V}\stackrel{~}{}_{\rho ;j}^{2}|P_1`$ in terms of Chebyshev polynomials. The summation in Eq. (3) can be now done making use of the generating function for the Chebyshev polynomials $`G_U(a,b)_{j=0}^{\mathrm{}}a^jU_j(b)=(12ab+a^2)^1`$. To get a required $`j`$-dependent coefficient in the series we do integration(s) of both sides with appropriate weight. The sums required in the calculation of (3) read

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{a^{j+1}}{(j+1)}}U_{j+1}(b)={\displaystyle _0^a}{\displaystyle \frac{da^{}}{a^{}}}\left(G_U(a^{},b)1\right),`$ (19)
$`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{a^{j+1}}{j(j+1)}}U_{j+1}(b)={\displaystyle _0^a}𝑑a^{}{\displaystyle _0^a^{}}{\displaystyle \frac{da^{\prime \prime }}{\left(a^{\prime \prime }\right)^2}}\left(G_U(a^{\prime \prime },b)12a^{\prime \prime }b\right).`$
Here $`a=\frac{1}{2}\sqrt{𝒫^2/q^2}`$ and $`b=q𝒫/\sqrt{q^2𝒫^2}`$. This immediately leads to

$`{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}P_2|{}_{}{}^{V}\stackrel{~}{}_{\rho ;j}^{2}|P_1=q^2\mathrm{\Pi }_{\rho \tau }{\displaystyle \frac{}{q_\tau }}{\displaystyle _\mathrm{\Omega }}𝑑y𝑑zf(y,z)`$ (20)
$`\times \left\{{\displaystyle \frac{1\sqrt{1+^2}}{2\sqrt{1+^2}}}\mathrm{ln}\left(1+{\displaystyle \frac{1\sqrt{1+^2}}{2\mathrm{\Xi }}}\right){\displaystyle \frac{1+\sqrt{1+^2}}{2\sqrt{1+^2}}}\mathrm{ln}\left(1+{\displaystyle \frac{1+\sqrt{1+^2}}{2\mathrm{\Xi }}}\right)\right\},`$
$`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{j}}P_2|{}_{}{}^{V}\stackrel{~}{}_{\rho ;j}^{2}|P_1=q^2\mathrm{\Pi }_{\rho \tau }{\displaystyle \frac{}{q_\tau }}{\displaystyle _\mathrm{\Omega }}𝑑y𝑑zf(y,z)`$
$`\times \{{\displaystyle \frac{1}{\mathrm{\Xi }}}+{\displaystyle \frac{1\sqrt{1+^2}}{2\sqrt{1+^2}}}(1+{\displaystyle \frac{1\sqrt{1+^2}}{2\mathrm{\Xi }}})\mathrm{ln}(1+{\displaystyle \frac{1\sqrt{1+^2}}{2\mathrm{\Xi }}})`$
$`{\displaystyle \frac{1+\sqrt{1+^2}}{2\sqrt{1+^2}}}(1+{\displaystyle \frac{1+\sqrt{1+^2}}{2\mathrm{\Xi }}})\mathrm{ln}(1+{\displaystyle \frac{1+\sqrt{1+^2}}{2\mathrm{\Xi }}})\},`$ (21)
where we introduced the conventions

$$\mathrm{\Xi }\frac{q^2}{q𝒫},^2\frac{q^2𝒫^2}{\left(q𝒫\right)^2}.$$
(22)
Inserting our findings (20) and (3) into Eq. (3) the next step is a mere differentiation w.r.t. $`q_\mu `$, which is done by the formula

$$q^2\mathrm{\Pi }_{\mu \nu }\frac{}{q_\nu }\tau (^2,\mathrm{\Xi })=\left\{2^2\left(q_\mu +\mathrm{\Xi }𝒫_\mu \right)\frac{}{^2}𝒫_\mu \frac{}{\mathrm{\Xi }^1}\right\}\tau (^2,\mathrm{\Xi })$$
(23)
for a test function $`\tau `$. Finally, we contract the resulting equation with the tensor $`S_{\mu \nu ;\rho \sigma }`$ which leads to the hadronic tensor

$$T_{\mu \nu }^2=\frac{1}{q^2}_\mathrm{\Omega }𝑑y𝑑zf(y,z)\left\{\left(q^2g_{\mu \nu }q_\mu q_\nu \right)_1+\left(𝒫_\mu +\frac{q_\mu }{\mathrm{\Xi }}\right)\left(𝒫_\nu +\frac{q_\nu }{\mathrm{\Xi }}\right)_2\right\},$$
(24)
with mass-dependent coefficient functions

$`_1={\displaystyle \frac{4\mathrm{\Xi }^2\left(1\mathrm{\Xi }\right)}{\mathrm{\Xi }\left[4\mathrm{\Xi }(1+\mathrm{\Xi })^2\right]}}+{\displaystyle \frac{^2\left(2\mathrm{\Xi }^2\right)}{4\mathrm{\Xi }\left(1+^2\right)^{3/2}}}\mathrm{ln}\left({\displaystyle \frac{1\sqrt{1+^2}+2\mathrm{\Xi }}{1+\sqrt{1+^2}+2\mathrm{\Xi }}}\right)+\left(\mathrm{\Xi }\mathrm{\Xi }\right),`$
$`_2={\displaystyle \frac{\mathrm{\Xi }\left[4\mathrm{\Xi }^2\left(1\mathrm{\Xi }\right)\right]}{(1+^2)\left[4\mathrm{\Xi }(1+\mathrm{\Xi })^2\right]}}+{\displaystyle \frac{3\mathrm{\Xi }^2\left(2\mathrm{\Xi }^2\right)}{4\left(1+^2\right)^{5/2}}}\mathrm{ln}\left({\displaystyle \frac{1\sqrt{1+^2}+2\mathrm{\Xi }}{1+\sqrt{1+^2}+2\mathrm{\Xi }}}\right)`$
$`+\left(\mathrm{\Xi }\mathrm{\Xi }\right).`$
In order to get the right $`i0`$-prescription, so as to pick up a correct sheet of the Riemann surface for the logarithm, we notice that in (3) we have to restore the suppressed Feynman prescription as follows $`q^2q^2+i0`$.

Let us comment on our result. As we observe, the leading order massless (generalized) Callan–Gross-type relation (see ) $`_2=\mathrm{\Xi }^2_1`$ is violated by target mass corrections. We checked that our result coincides in the forward limit with the well-known ones . Obviously, we have $`q_\mu T_{\mu \nu }^2=0`$ and therefore current conservation is fulfilled in the forward case. Note that this is only the case in the sum of the direct and the crossed amplitudes. However, the gauge invariance in the off-forward kinematics is violated, i.e. $`q_{2\mu }T_{\mu \nu }^20`$ and, as it was previously studied in the case of massless amplitudes, it is restored once higher twist corrections are accounted for . This must also persist to the case of traceless operators which are not projected on the light cone.

Now let us generalize the formalism to a spin-$`\frac{1}{2}`$ target. The matrix elements $`P_2|_{\rho ;j}^2|P_1`$ are build now from the vectors $`P`$, $`\mathrm{\Delta }`$ and a set of (independent) Dirac bilinears, characterizing the spin content of the target. For our purposes it is convenient to use the following structures

$$(h_\mu ,\stackrel{~}{h}_\mu )=\overline{U}(P_2)\gamma _\mu (1,\gamma _5)U(P_1),(b,\stackrel{~}{b})=\overline{U}(P_2)(1,\gamma _5)U(P_1),$$
(26)
for (vector, axial) sectors, respectively. Thus, we parametrize the matrix element of e.g. vector operator according to

$`P_2|{}_{}{}^{V}𝒪_{\rho ;\mu _1\mathrm{}\mu _j}^{\mathrm{sym}}|P_1`$ $`=`$ $`h_{\{\rho }P_{\mu _1}\mathrm{}P_{\mu _j\}}A_{j+1,j+1}+\mathrm{}+h_{\{\rho }\mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j\}}A_{j+1,1}`$
$`+`$ $`{\displaystyle \frac{b}{2M}}\left\{P_{\{\rho }P_{\mu _1}\mathrm{}P_{\mu _j\}}B_{j+1,j+1}+\mathrm{}+\mathrm{\Delta }_{\{\rho }\mathrm{\Delta }_{\mu _1}\mathrm{}\mathrm{\Delta }_{\mu _j\}}B_{j+1,0}\right\}+\mathrm{},`$
where again terms proportional to the metric tensor are not needed. The non-local twist-two vector operator reads

$`P_2|{}_{}{}^{V}_{\rho }^{2}(x,x)|P_1`$ $`=`$ $`{\displaystyle _\mathrm{\Omega }}𝑑y𝑑z\left\{f_A(y,z)h^𝒫+f_B(y,z){\displaystyle \frac{b}{2M}}𝒫^𝒫\right\}`$
$`\times _\rho {\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i)^j}{j!(j+1)^2}}\left({\displaystyle \frac{x^2𝒫^2}{4}}\right)^{(j+1)/2}U_{j+1}\left({\displaystyle \frac{x𝒫}{\sqrt{x^2𝒫^2}}}\right).`$
Here we have generated the factor of $`x`$ in $`hx`$ by a differentiation w.r.t. $`𝒫`$. The same we did with the $`b`$ form factor for the purpose of a uniform representation, as well as to have a cross check on our previous calculations. Obviously, this step is not required at all and we can readily borrow the results already found for a (pseudo) scalar target. The same equation holds for axial operator with a trivial dressing of symbols with tildes.

The reduced matrix elements $`B_{j+1,j+1k}`$ are represented in terms of GPDs or DDs as in Eq. (10), — replace $`f(y,z)`$ by $`f_B(y,z)`$. For the reduced matrix elements $`A_{j+1,j+1k}`$ we write

$`A_{j,jk}`$ $`=`$ $`{\displaystyle \frac{1}{k!}}{\displaystyle \frac{^k}{\eta ^k}}|_{\eta =0}{\displaystyle _1^1}𝑑xx^{j1}A(x,\eta )=\left({\displaystyle \genfrac{}{}{0pt}{}{j1}{k}}\right){\displaystyle _\mathrm{\Omega }}𝑑y𝑑zy^{jk1}z^kf_A(y,z),`$ (29)
where $`j1`$ and the index $`k`$ now varies in the intervals $`0kj1`$. There is an immediate consequence of the fact that $`j`$-th moment of the function $`A(x,\eta )`$ is a polynomial of order $`\eta ^{j1}`$ only. Namely, translating it to the conventional Ji’s distributions, $`H=A+B`$, $`E=B`$, we immediately find that

$`{\displaystyle \frac{}{\eta ^j}}{\displaystyle _1^1}𝑑xx^{j1}\left\{H(x,\eta )+E(x,\eta )\right\}=0.`$
As a particular example we recall that for the $`j=2`$ moment, which corresponds to Ji’s sum rule, the $`\eta `$-dependence drops off as it was noted before . For the axial channel the GPDs $`\stackrel{~}{A}`$ and $`\stackrel{~}{B}`$ are identical to the conventional $`\stackrel{~}{H}`$ and $`\stackrel{~}{E}`$ and this implies that the $`j`$-th moment of $`\stackrel{~}{H}`$ ($`\stackrel{~}{E}`$) has the expansion to order $`j1`$ ($`j`$) in $`\eta `$.

The Fourier transform and the resummation are done in the same way as before and result in

$`T_{\mu \nu }^2`$ $`=`$ $`{\displaystyle \frac{1}{q^2}}{\displaystyle _\mathrm{\Omega }}dydz\{(\stackrel{~}{f}_A(y,z)\stackrel{~}{h}^𝒫+\stackrel{~}{f}_B(y,z){\displaystyle \frac{\stackrel{~}{b}}{2M}}𝒫^𝒫)iϵ_{\mu \nu \rho \sigma }q_\rho 𝒫_\sigma 𝒢_1`$
$`+\left(f_A(y,z)h^𝒫+f_B(y,z){\displaystyle \frac{b}{2M}}𝒫^𝒫\right)`$
$`\times ((q^2g_{\mu \nu }q_\mu q_\nu )_1+(𝒫_\mu +{\displaystyle \frac{q_\mu }{\mathrm{\Xi }}})(𝒫_\nu +{\displaystyle \frac{q_\nu }{\mathrm{\Xi }}})_2)\},`$
with massive coefficient functions

$`_1`$ $`=`$ $`{\displaystyle \frac{[2\mathrm{\Xi }\left(1+2^2\right)^4]L_{}}{4\mathrm{\Xi }\left(1+^2\right)^{3/2}}}{\displaystyle \frac{L_+}{2\left(1+^2\right)}}{\displaystyle \frac{^2L}{2\left(1+^2\right)^{3/2}}}+\left(\mathrm{\Xi }\mathrm{\Xi }\right),`$
$`_2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Xi }[2\mathrm{\Xi }\left(1+4^2\right)3^4]L_{}}{4\left(1+^2\right)^{5/2}}}{\displaystyle \frac{\mathrm{\Xi }^2(12^2)L_+}{2\left(1+^2\right)^2}}{\displaystyle \frac{3\mathrm{\Xi }^2^2L}{2\left(1+^2\right)^{5/2}}}+\left(\mathrm{\Xi }\mathrm{\Xi }\right),`$
$`𝒢_1`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Xi }L_{}}{2\left(1+^2\right)^{1/2}}}{\displaystyle \frac{\mathrm{\Xi }L_+}{2\left(1+^2\right)}}{\displaystyle \frac{\mathrm{\Xi }^2L}{2\left(1+^2\right)^{3/2}}}\left(\mathrm{\Xi }\mathrm{\Xi }\right),`$ (31)
where

$`L_\pm `$ $``$ $`\mathrm{ln}{\displaystyle \frac{1\sqrt{1+^2}+2\mathrm{\Xi }}{2\mathrm{\Xi }}}\pm \mathrm{ln}{\displaystyle \frac{1+\sqrt{1+^2}+2\mathrm{\Xi }}{2\mathrm{\Xi }}},`$
$`L`$ $``$ $`\mathrm{Li}_2\left({\displaystyle \frac{1\sqrt{1+^2}}{2\mathrm{\Xi }}}\right)\mathrm{Li}_2\left({\displaystyle \frac{1+\sqrt{1+^2}}{2\mathrm{\Xi }}}\right).`$
Here $`\mathrm{Li}_2`$ is the Euler dilogarithm $`\mathrm{Li}_2(x)=_0^x\frac{dy}{y}\mathrm{ln}(1y)`$.

As we mentioned before, the resummed mass corrections can not be expressed in terms of GPDs (at least we did not succeed in doing this). However, every term in the mass expansion, $`M^2/q^2`$, can be converted into the conventional GPDs representation. Let us demonstrate it for the expanded Compton form factor $`F_1`$ for a (pseudo) scalar target given in Eq. (3). To the first non-trivial order, i.e. to $`𝒪\left(M^4/q^4\right)`$ accuracy, the latter reads in DD form

$`F_1`$ $`=`$ $`{\displaystyle _\mathrm{\Omega }}dydzf(y,z)\{(C_1^{(0)}(\mathrm{\Xi }^1i0)+C_1^{(0)}(\mathrm{\Xi }^1i0))`$
$`{\displaystyle \frac{M^2}{q^2}}y^2\mathrm{\Xi }^2(C_1^{(1)}(\mathrm{\Xi }^1i0)+C_1^{(1)}(\mathrm{\Xi }^1i0))\},`$
where in our approximation, $`M^2\mathrm{\Delta }^2`$, we set $`^2=4\frac{M^2}{q^2}y^2\mathrm{\Xi }^2`$ and introduced the coefficient functions $`C_1^{(0)}\left(\mathrm{\Xi }^1\right)=\left(1+\mathrm{\Xi }\right)^1`$ and $`C_1^{(1)}\left(\mathrm{\Xi }^1\right)=\left(1+\mathrm{\Xi }\right)^2+2\mathrm{ln}\left(\frac{\mathrm{\Xi }}{1+\mathrm{\Xi }}\right)`$. To cast the DD into the usual GPD representations, which are related by the equation $`x_\mathrm{\Omega }𝑑y𝑑zf(y,z)\delta (xy\eta z)B(x,\eta )`$, we need the following general result

$$_\mathrm{\Omega }𝑑y𝑑z\{\begin{array}{c}z^n\\ y^n\end{array}\}\left(\mathrm{\Xi }^1\right)f(y,z)=𝑑x\left(\frac{x}{\xi }\right)𝑑x^{}V_1^{(n)}(x,x^{})\{\begin{array}{c}\left(\frac{\stackrel{}{}}{\eta }\right)^n\\ _{k=0}^{n1}\left(\widehat{d}(x^{},\eta )k\right)\end{array}\}B(x^{},\eta ),$$
(32)
where we used a shorthand notation for the differential operator $`\widehat{d}(x,\eta )=x\frac{\stackrel{}{}}{x}\eta \frac{\stackrel{}{}}{\eta }`$ and the kernels

$`V_1^{(n)}(x,x^{})={\displaystyle \frac{(x^{}x)^{n1}}{(n1)!}}\mathrm{\Theta }_{11}^0(x,xx^{}),\text{with}\mathrm{\Theta }_{11}^0(x,y)={\displaystyle \frac{\theta (x)\theta (y)}{xy}}.`$ (33)
Substitution of these expressions back into Eq. (3) results into

$`F_1`$ $`=`$ $`{\displaystyle _1^1}dx{\displaystyle _1^1}dx^{}\{(C_1^{(0)}(\frac{x}{\xi }i0)+C_1^{(0)}(\frac{x}{\xi }i0))\delta (xx^{})`$
$`{\displaystyle \frac{M^2}{q^2}}(C_1^{(1)}(\frac{x}{\xi }i0)+C_1^{(1)}(\frac{x}{\xi }i0))V_1^{(2)}(x,x^{})(\widehat{d}(x^{},\eta )1)\widehat{d}(x^{},\eta )\}B(x^{},\eta ),`$
which has now the desired form. Completely analogous manipulations produce the GPD representation for all other form factors.

## 4 Conclusions.

In the present note, we have presented a formalism for the resummation of target mass corrections stemming from the twist-two approximation to the off-forward Compton scattering amplitude. The machinery is based on the group-theoretical content of the procedure of trace subtraction from local operators with definite symmetry property – the operators with a well defined twist are given by a harmonic projection of the latter. This allows to resum target mass corrections to all orders in a compact form in terms of DDs. However, once the amplitudes are expanded in $`M^2/q^2`$ series, they can be cast in the form involving GPDs.

An important feature of the result is that the mass corrections enter with additional powers of the scaling variable $`\xi `$. Thus, we expect that they are suppressed in the small $`\xi `$ ($`x_{\mathrm{Bj}}`$) region. Since there still might be a numerical enhancement due to coefficient functions, a numerical model dependent study of the target mass effects to the deeply virtual Compton scattering process for the kinematics of HERMES and Jefferson Lab experiments is necessary.

In the twist-two approximation the current conservation is violated and must be restored by taking into account contributions with higher geometrical twists. To make the former manifest, one has to decompose those higher twist operators in a linear independent set by means of the QCD equation of motion . This procedure gives us a part which is expressed in terms of total derivatives acting on operators with a lower twist. Of course, the latter when combined together with the leading prediction must render a manifestly gauge invariant result. However, the trace subtraction and the application of the QCD equation of motion do not commute and one can easily fail to reproduce an amplitude which respects current conservation. This problem goes beyond the scope of the present study and requires further investigations.

We would like to thank A.V. Radyushkin for a clarifying discussion of Ref. .