Shear Viscosity at Infinite Coupling: A Field Theory Calculation

I derive an exact integral expression for the ratio of shear viscosity over entropy density $\frac{\eta}{s}$ for the massless (critical) O(N) model at large N with quartic interactions. The calculation is set up and performed entirely from the field theory side using a non-perturbative resummation scheme that captures all contributions to leading order in large N. In 2+1d, $\frac{\eta}{s}$ is evaluated numerically at all values of the coupling. For infinite coupling, I find $\frac{\eta}{s}\simeq 0.42(1)\times N$. I show that this strong coupling result for the viscosity is universal for a large class of interacting bosonic O(N) models.

I derive an exact integral expression for the ratio of shear viscosity over entropy density η s for the massless (critical) O(N) model at large N with quartic interactions. The calculation is set up and performed entirely from the field theory side using a non-perturbative resummation scheme that captures all contributions to leading order in large N. In 2+1d, η s is evaluated numerically at all values of the coupling. For infinite coupling, I find η s 0.42(1)×N . I show that this strong coupling result for the viscosity is universal for a large class of interacting bosonic O(N) models.

PACS numbers:
The shear viscosity is a transport coefficient that encodes how efficiently spatial anisotropies are transmuted into momentum anisotropies (velocity gradients). For relativistic systems, a key dimensionless ratio that quantifies this efficiency is found by the ratio of shear viscosity η to entropy density s. Experimental measurements of momentum anisotropies in heavy-ion collisions together with hydrodynamic modeling constrain the value of shear viscosity in QCD to η s < ∼ 0.2 [1][2][3][4][5]. This numerical value happens to be not too dissimilar from the result η s = 1 4π 0.08 found for the conjectured strong-coupling limit of another gauge theory, N = 4 Super-Yang-Mills (SYM) theory, in the large N limit [6][7][8]. By contrast, in QCD where calculations are limited to weak coupling only, one typically encounters values η s 1 4π [9][10][11]. At intermediate couplings, results for the shear viscosity exist for QED with many fermion flavors [12] and for SU(3) gauge theory from lattice Monte-Carlo simulations [13,14]. At strong coupling, results for transport coefficients have been limited to theories with known holographic duals, with the exception of so-called "thermodynamic transport coefficients" such as κ [15]. In particular, there is no known example of a theory where η s has been determined for all coupling strengths.
The present work is meant to fill this gap, and provide the complete coupling-dependence for the shear viscosity over entropy ratio by directly calculating the relevant low-frequency limit of the energy-momentum tensor correlator from the field theory. The theory to be studied will be the O(N) vector model with quartic interactions in 2+1 dimensions in the large N limit, which exists for all values of the coupling. The choice of this theory is motivated by the fact that the entropy density s is known for all couplings [16] and that an efficient resummation scheme that captures all relevant contributions to the shear viscosity at large N is known [17,18]. It also helps that the most difficult part of the calculation, namely the evaluation of the shear viscosity coefficient for the O(N) vector model, has already been set up for the 3+1d theory using a variant of the 2PI formalism in Ref. [19]. Thus strictly speaking, the only new result in the present work will be the determination of η s at infinite coupling, which can be done in the 3d large N O(N) model (but does not make sense because of the Landau pole in 4d).
A major objection to calculating the shear viscosity in a theory with only two space dimensions is the presence of so-called long-time tails [20], which normally lead to a divergent two-dimensional shear viscosity when naively taking the low frequency limit. However, for the O(N) model it so happens that η s ∝ N , so that in the large N limit long-time tails are suppressed by three powers of N, cf. the discussion in Ref. [21]. For this reason, the calculation of the shear viscosity as the zero-frequency limit of the relevant energy-momentum correlator reported in this work is well defined.

Preliminaries
In the interest of brevity, I will not review the set up of finite temperature correlators in quantum field theory (the interested reader is referred to the Supplemental Material for this topic). Denoting Minkowski momenta as K = (ω, k) and space-time coordinates as X = (t, x), I will use the relation of retarded real-time Greens functions G R and the spectral density ρ given by This relation can be used to derive the analytic continuation of the corresponding Euclidean correlator to real time (cf. Refs. [22,23]) In order to connect the retarded correlator with transport coefficients, one employs the low-momentum expansion of G R provided by fluid dynamics. Fluid dynamics is the effective field theory of conserved quantities for small moment. For a theory that conserves energy and momentum, good EFT variables are the energy density and fluid 4-velocity , u µ , fulfilling u µ u µ = −1 (see e.g. Ref. [22][23][24] for reviews of modern fluid dynamics). Using fluid dynamics, it is straightforward to derive the relation G R (ω, k = 0) = P − iωη + O(ω 2 ) for the T xy correlator in d ≥ 3 spacetime dimensions, from which the so-called Kubo relation follows: Including thermal fluctuations in the fluid dynamic calculations leads to a long-time tail contribution for d = 3 of the form G R (ω, k = 0) ∝ iωT 2 η/s ln ω (see Supplemental Material). This term in general invalidates the Kubo formula (3), except in the large N limit where it is suppressed by η s ∝ N , thereby allowing the use of (3) to extract the shear-viscosity for d = 3 to O(N 0 ).
The O(N) model The field theory I consider in this work is that of a massless ("critical") N-component scalar field φ a , a = 1, 2, . . . , N with Euclidean action where at finite temperature T the X 0 = τ direction is compactified on a circle with radius β. The partition function for this theory Z = Dφe −S may be rewritten by inserting 1 = Dσδ(σ − φ a φ a ) = DσDζe i ζ(σ−φaφa) . Integrating out the σ field gives Splitting the field ζ = 1 2 ζ 0 + ζ into a zero mode and fluctuations, the action in (5) becomes where S R0,I = i X ζ φ a φ a . At leading order in large N, only the zero-mode from the ζ-field contributes to the partition function, so S R0,I may be ignored. This corresponds to a particular resummation of finite-temperature Feynman diagrams (certain tadpoles), dubbed the R0-level resummation in Ref. [17]. One finds for R0 where V is the d − 1 dimensional volume of space and In the large N limit, the remaining ordinary integral in (7) can be evaluated exactly from the saddle point located at iζ 0 = m 2 , so that Z R0 = e −N V β(m 4 −J(m)) , where [25,Eq. (2.44)]. Using dimensional regularization for d = 3 − 0 + , the saddle point condition becomes π λm 2 +m + 2 ln 1 − e −m = 0 , wherem ≡ βm,λ ≡ βλ. Note that in the strong coupling limitλ → ∞ (9) has a universal solutionm = 2 ln 1+ √ 5 2 [16,26]. Basic thermodynamic relations give the exact large N entropy density as s = ∂ ∂T ln Z R0 βV . Using (9), integration by parts, and a variable change to E = E k , for d=3 this relation leads to For some correlation functions, additional interactions not included in the R0 resummation may contribute at leading order in large N, making it necessary to go beyond the R0-level. To this end, slightly changing notation from Ref. [17] rewrite the action (6) as S R0,0 +S R0,I = S R2,0 + S R2,I with .
(12) The R2-level resummation then consists of calculating Σ, Π self-consistently up to (including) one-loop level using S R2,I , finding At finite temperature, more subtle ways generating contributions at leading order in large N arise, making it necessary to go beyond the R2 resummation. To this end, rewrite the action (6) again using S R2,0 , S R2,I but express where Γ 3 (P, K) = 1 + δΓ 3 (P, K) is the resummed threevertex function and P, K are the incoming momenta for the ∆-propagators. The R3 level resummation then consists of calculating Γ 3 self-consistently to one-loop level, and Σ, Π self-consistently up to (including) twoloop level. Diagrammatically, one finds [17] δΓ 3 = −4 , Π = 2 , Σ = 4 (15) where the full vertex is denoted by wiggly lines denote D(K) and straight lines denote ∆(K). While the R0 resummation contains all leading order large N contributions for the zero-point function, R2 and R3 contain all large N contributions for the two-and three-point function, respectively. The energymomentum tensor is a four-point function, so including all large N contributions requires going to R4. The R4 resummation is found by adding and subtracting a term Calculating the one-loop 4-point vertex self-consistently, one finds diagrammatically whereas the large N three-point vertex in R4 becomes and the self-energies Π, Σ are still diagrammatically given by (15). Note that since D ∝ N −1 , the 4-point vertex and the triangle diagram contribute at the same order at large N.

Energy-Momentum Tensor Correlators
For the action (4), the Euclidean operator for the energy-momentum tensor is given by T xy (X) = ∂ x φ a (X)∂ y φ a (X) so that the Euclidean energy-momentum tensor correlator is defined by In the R0-level resummation, following (6) and neglecting S R0,I , the R0-action is quadratic in φ and hence G xy,xy E can be calculated using Wick's theorem. At finite temperature, φ a (X) = T k E k e iK·X φ a (K) where k E = 2πnT are the bosonic Matsubara frequencies with n ∈ Z. In Fourier space with p = px, one finds [15,16] G xy,xy E,R0 (P ) = 2N where the mass m is determined by (9). The result may be analytically continued to real frequencies ω using Eq.
(2). The R0-level result G xy,xy R,R0 (ω, p = 0) constitutes the correct large N result for the retarded correlator except for the region βω To appreciate this statement, let us reconsider G xy,xy E in the R2-level resummation. Using (11) and neglecting S R2,I , the R2-action is once again quadratic in the fields, and hence G xy,xy E can be calculated using Wick's theorem. Setting again p → 0 gives with Σ = Σ R2 in (12) and Σ R2 specified by (13). Using the spectral representation of the propagator the thermal sum and angular integral is straightforward, Note that since Σ ∝ 1 N , in the naive large N limit (21) becomes the spectral function for a free massive particle, One finds that the product ρ(k 0 )ρ(k 0 − ω) has contributions where poles "pinch" the real k 0 axis from above and below in the large N limit. Focusing on the limit βω 1 N , this contribution becomes which is proportional to N in the large N limit. Note that the other contributions (as well as the whole product ρ(k 0 )ρ(k 0 − ω) for βω 1 N ) give contributions of order O(N 0 ), bringing us back to the R0 result (19). The non-commutative nature of the limits βω → 0, 1 N → 0 imply that for the calculation of transport coefficients, contributions that are naively subleading at large N become important. This enhancement process, known as "pinching poles", has a long history in nuclear physics literature, cf. Refs. [19,27].
In the small frequency limit, the R2-level expression for the shear viscosity from Eq. (3) is The R2-level expression for the shear viscosity is welldefined for all values of the coupling, but it does not contain all contributions to leading order in large N. To this end, let us reconsider the correlator G xy,xy E in R4-level resummation, specified by (11) with (17). Since the R4level action contains a non-trivial vertex, Wick's theorem can no longer be used to evaluate the correlator; instead, interactions must be included (See the Supplemental Material for an example of how a generic 4-point function is evaluated beyond R2.) For the energy-momentum correlator G xy,xy E,R4 , the situation is slightly simpler than for a generic 4-point function, because the spatial derivatives ∂ x φ, ∂ y φ in the definition imply that some contributions do not contribute after angular integration (see Supplemental Material). However, instead of the regular 3-vertex Γ 3 , the corresponding contribution to G xy,xy E,R4 contains momenta k x , k y inside the vertex. One thus finds where Γ 3,xy = k x k y + δΓ 3,xy denotes the resummed 3vertex, and the propagator ∆(K) contains Σ in the R4level resummation, cf. Eq. (15). The resummed vertex is the only modification with respect to the R2 result (20), so one needs to check if δΓ 3,xy , which naively is order 1 N gets enhanced in the limit P → 0. One finds (27) where the integral kernel is W (P, Q, K) = D(P −K−Q)+ 2N Γ 4 (Q, P − Q, K, P − K). The structure of (27) indeed suggests a "pinching pole" similar to (20) in the limit P → 0, whereas for other kinematic regions δΓ 3,xy ∝ 1 N . Two of the internal vertices in (27) therefore do not receive corrections in the limit P → 0. One may verify that in this limit, δΓ 3,xy (K, −K) = kxky k 2Γ3 (K), such that after doing the angular average one obtains where to leading order in large N The analytic structure of the vertexΓ 3 (K, P − K) is as follows: expressing the propagators in terms of their spectral functions ρ(µ), ρ D (µ), one can perform the analytic continuation to real frequencies. In a first step, set-tingΓ 3 (Q, P − Q) = q 2 , the thermal sum over q E in (27) can be done explicitly. Since the integrations over the arguments of the spectral functions range over the whole real axis, one finds thatΓ 3 (ik 0 , ip 0 − ik 0 ) has branch cuts along the whole real line for Im p 0 = 0, Im k 0 = 0, Im (p 0 − k 0 ) = 0. This structure is unchanged when iterating the vertex.
With the analytic structure of the vertex known, one proceeds to evaluate (28). First, write the thermal sum as with C encircling the poles at the imaginary Matsubara frequencies k 0 = ik E in a counter-clockwise fashion. Next, the propagators ∆(K)∆(P − K) have branch cuts at Im k 0 = 0, Im(k 0 ) = −p E . The additional branch cut from the vertex is independent from k 0 , and hence not relevant for the thermal sum (30). Deforming the contour C such that it runs along both sides of the two branch cuts, one encounters terms such as The leading large N contributions are then given by combinations such as ∆ R (k 0 )∆ A (k 0 −ω) which have singularities on either side of the real k 0 axis ("pinching poles"), whereas the others can be neglected. Letting k 0 → k 0 − p 0 in the second branch cut contribution, only a particular analytic continuation of the vertex function contributes, which for p 0 → 0 and k 0 = ±E k becomes (32) Clearly, this corresponds to use standard analytic continuation ofΓ 3 (K, P − K) first for k E and then for p E . Assuming that F (E k ) is real (which will be shown below), one can then follow the same procedure that led to (25), with the only modification arising from the resummed vertex function. Using (10), one finds where again k = √ E 2 − m 2 . This relation is exact in the large N limit, as it contains all leading order large N contributions to the shear viscosity and entropy density, cf. Ref [19]. Note that because Σ ∝ 1 N , one finds that η s ∝ N in the large N limit. This N-scaling is generic for vector or fermionic theories [12,19], but is qualitatively different from large N gauge theories where η/s ∝ O N 0 [6,9] In order to get a result for η s , one needs to know the functions Im Σ(E) and F (E). These follow from finitetemperature field theory calculations, see Supplemental Material and Refs. [19,29]. For the evaluation of η s , it is convenient to use quadrature to recast the integrals in for n = 0, 1, . . . , K. Expanding s because k 2 in the numerator of (33) only involves polynomials P n (E) up to degree two.
Results and Universality at Infinite Coupling Numerical evaluation of (33) for all values of the dimensionless couplingλ = βλ is shown in Fig. 1. For weak couplingsλ, η sN is large because the thermal width is small. As the coupling is increased, I find that the ratio of shear viscosity over entropy density is dropping monotonically, but is finite in the limit ofλ → ∞. In this strong coupling limit, the numerically calculated result becomes One may ask about the universality of the strong coupling result. To this end, consider a modification of the action (4) to where now λ is dimensionless. This action has the property that it is a CFT for all values of λ at large N [16].
Using the same replacement σ = φ a φ a as before, and integrating out σ one finds (37) At large N, the asymptotic properties of the Airy function then give the form of the partition function as (38) While it may be challenging to construct R4 for this Lagrangian for general values of λ, the strong coupling limit λ → ∞ of this theory is exactly equal to the strong coupling limitλ → ∞ of Eq. (5). For this reason, one can explicitly write down R4 for this theory in the strong coupling limit where the only difference is the form of the propagator cf. Eq. 11. Hence in the infinite coupling limit the result for the shear viscosity for (36) is identical to that of (4). It is not hard to generalize this proof to theories with other interactions φ a φ a → U (φ a φ a ) for a large class of potentials U (x), which demonstrates that (35) is the universal strong coupling shear viscosity over entropy ratio for a large class of bosonic QFTs. The same universal behavior was found [16] for the weak-strong ratio s λ=∞ s λ=0 = 4 5 and for the boson in-medium mass lim λ→∞m = 2 ln 1+ √ 5

. Summary and Discussion
In this work, I derived an exact large N expression for the ratio of shear viscosity over entropy density for the interacting O(N) model, using well-established field theory techniques. Evaluating the expression numerically in the case of 2+1 dimensions, I found the strong coupling result (35).
Regardless of the numerical value, the present work demonstrates that it is possible to calculate transport properties at infinite coupling directly from quantum field theory, without invoking dualities or conjectures of any kind. While the field theory studied here may not be of interest to most readers, it can nevertheless serve as a test bed for strong coupling transport which would otherwise be inaccessible or very hard by any other means. For instance, having access to exact full energy-momentum tensor correlation functions for all values of the coupling allows to study the onset/breakdown of hydrodynamics from first principles, cf. Refs. [30][31][32]; exact real-time correlators also can be used to test analytic continua-tion techniques employed in lattice Monte-Carlo studies [13,14]; exact results for transport coefficients can be used as a rigorous test case for approximation schemes that are used for e.g. QCD [9][10][11].
In addition to serving as a well-defined test bed for general-purpose tools, the present calculation may be generalized in several ways. By, for instance, calculating other transport coefficients such as the bulk viscosity ζ as well as relaxation time τ π for the O(N) model; calculating transport coefficients for other large N theories [33][34][35]; calculating exact far-from equilibrium real-time dynamics at large N for a quantum field theory.
For these reasons, I am optimistic that the present result can become useful in the future.
Acknowledgments I am indebted to Gert Aarts for clarifying some questions I had concerning Ref. [19], as well as providing numerical data for η s in the 4d O(N) model from this reference. Also, I thank Scott Lawrence and Max Weiner for fruitful discussions, and Marcus Pinto for pointing out a typo in (38). This work was supported by the Department of Energy, DOE award No DE-SC0017905.

Finite Temperature Correlators
Given the operatorφ(X ), one can define the Wightman functions G > (X ) = φ (X )φ(0) and G < (X ) = φ (0)φ(X ) , which in turn are related to the retarded real-time Greens function [22, Eq. (2.9)] . As shown in [25, Eq. (8.10)], for thermal equilibrium, the Wightman functions fulfill the KMS relation G > (t − iβ) = G < (t), where β = 1 T is the inverse temperature. Using mostly plus metric convention, the Fourier transform G(K) = X e −iX ·K G(X ) of the KMS relation becomes G > (K)e −βω = G < (K). The spectral function is defined as In thermal equilibrium, all correlators can be related to each other, and can be expressed through the spectral function.
In particular, using the representation of the step function θ(t) = i dν 2π e −iνt ν+i0 + one proves the relation which is called the spectral representation of the retarded correlator. Using the definition of the advanced correlator G A (K) = G * R (K), this relation can be used to show One can relate these real-time Greens function to the imaginary-time Greens function obtained in the Euclidean formulation. Euclidean momenta and coordinates will be denoted by K = (k E , k, k 4 ) and X = (τ, x), respectively. Defining the Euclidean correlator for imaginary time τ as this matches G > (t, x) if formally identifying τ → it. Therefore, Using these relations, and using (41), one may then prove which is called the spectral representation of the Euclidean correlator. Comparing (46), (42), the retarded correlator for real frequencies ω is obtained by analytic continuation of the Euclidean frequencies (2), such that ρ(P) = Im G E (ω E → iω − 0 + , p).

Fluid Dynamics in two space dimensions
In fluid dynamics, using the EFT variables , u µ as building blocks, one constructs the energy-momentum tensor in a gradient expansion as where g µν is the metric tensor, p = p( ) is the pressure related to the energy density via the equilibrium equation of state, η, ζ are the shear and bulk viscosity coefficients, respectively, ∆ µν = g µν + u µ u ν and σ µν = ∇ µ ⊥ u ν + ∇ ν ⊥ u µ − 2 d−1 ∆ µν ∇ ⊥ λ u λ , ∇ µ ⊥ = ∆ µν ∂ µ and d denotes the number of space-time dimensions, cf. Ref. [23,Eq. (2.30)]. Higherorder versions of fluid dynamics exist, but will not be considered here. In order to calculate the real time correlator of the energy-momentum tensor, one considers the linear response of T µν with respect to metric fluctuations δg µν around Minkowski space, G 0i,0j where and c 2 s = dP d is the speed of sound squared. Using ∂ µ T µν = 0, neglecting contact terms and letting k = kx, this leads to In particular, find for the spectral functions for d = 3 and k = 0 one finds ρ xx,xx (ω) = (η + ζ) ω , ρ xy,xy (ω) = ηω .
The poles of the retarded correlator (49) show that the fluid dynamics EFT contains collective modes (sound modes, shear modes). These modes will interact, thereby modifying the form of (47). It is possible to estimate the effect of these interactions for low momenta by calculating the contribution to G R arising from hydrodynamic self-interactions. Since this is an effective theory calculation, the resulting momentum integrals have to be cut-off at a scale p max ∝ γ −1 η below which the EFT can be trusted. Specifically, focusing on the shear mode contribution, one finds for d = 3 [22] G xy,xy R,1−loop,shear−shear (ω, which needs to be added to G xy,xy R in Eq. (50). Additional contributions to G xy,xy R result from the sound-sound mode contribution (same order as (52)), as well as higher-order loops (suppressed by powers of the temperature).
Example: 4-point correlator in R2, R3 As an example for the use of the R-level resummations, consider the zero-temperature 4-point function In R2, recognizing φ 2 i (X) = σ(X), integrating out the σ field as with the partition function (5) leads to At weak couplingλ 1, one may expand D(X) in a power series, finding C(X) = 2N ∆ 2 (X). For d = 3, the zero-temperature propagator is proportional to ∆(X) ∝ 1 X , so that C(X) ∝ N X 2 for λ 1. Similarly, Π R2 (X) ∝ 1 X 2 so that for strong coupling limλ →∞ D(X) ∝ 1 N X 4 . Neglecting the contact term, this leads to C(X) ∝ N λ 2 X 4 forλ 1. which explicitly proves the statement that the dimension of the operator changes from the free theory limit ("UV fix point") to the strong coupling limit ("IR fix point"), cf. [36].
It is instructive to reconsider this correlator at higher R-level resummation. At R3, one has as well as many other contributions that mutually cancel in R3, such as where boxes indicate insertions of self-energies or vertex functions in the R3-level resummation.
In addition, one finds that the first line in (55) can be summed-up in closed form using the resummed vertex Γ 3 , finding with the only difference to the R2 result (54) being that the Π is to be evaluated in the R3-resummation, cf. (15).

Thermal width
From (15), the self-energy Σ is given by The vertex function is O 1 N for generic momenta, with a pinching-pole enhancement only when K = P . Since the corresponding integration region has measure 1 N , the resummed vertex does not contribute to Σ(P ) to leading order in large N.
Using the spectral representation of the propagators, performing the thermal sum is straightforward. Upon analytic continuation p E → ip 0 − i0 + , one has This expression can be further simplified by using the large N free particle result for ρ(K) and coth βx 2 = 1 + 2n(x) to find Since only the on-shell limit will be needed, I denote Im Σ(E p ) ≡ Im Σ(p 0 = E p , p), and define the angular averages so that Note that while n(E k − E p ) has a pole at E k = E p ,ρ D (E p − E k ) has a zero there, so that the pole is integrable.

Polarization Tensor
To evaluate (60), the spectral function for the auxiliary propagator ρ D (K) is needed. This in turn depends on the polarization tensor (15) This expression has pinching poles and an order O(N 0 ) vertex correction for P → 0. However, the corresponding contribution to (58) has vanishing integral weight in the large N limit, so to leading order in N one may neglect both the resummed vertex as well as the pinching poles in (61). As a consequence, the required contribution for the polarization tensor is given by the simple sum-integral Π(P ) = 2 with m given by the solution from (9). I write the integral as two parts Π(P ) = Π V (P ) + Π T (P ) where Π T (P ) contains the Bose-Einstein distribution factor after evaluating the thermal sum, while Π V does not. For Π V = 2 K ∆(P − K)∆(K) one finds This may be analytically continued to give where E 1 = E k , E 2 = E p−k and I have used the symmetry k → p − k to simplify the integrand. Further simplification yields Taking apart the integrand, one can use [37, Eq. (B6)] to integrate over angles and after analytically continuing to real frequencies I find The remaining integrals can be done numerically, where it is useful to employ integration by parts. Using Π(P) = Π V (P) + Π T (P) the spectral function for the auxiliary correlator is Since Im(P) = 0 for P 2 ∈ [−4m 2 , 0], the spectral function is also vanishing in this region.

Vertex Resummation -R3
The relevant part of Eq. (27) is the pinching-pole contribution to the three-vertex. For pedagogical reasons, let me first evaluate this contribution in the R3-level resummation, where the integral kernel is W (P, Q, K) = D(P − K − Q). One findsΓ with the angular weight For (28), only the case p = 0 is needed. Using the spectral representation for D( , the relevant thermal sum can be written as (30) where all the spatial arguments of the propagators and vertex are q. The contour is encircling the poles of coth βq 0 2 on the imaginary axis. Apart from these poles, the integrand possesses a pole at q 0 = −µ − i(p E − k E ), a branch cut at Im q 0 = 0 from ∆(iq 0 ) and Γ 3 (iq 0 , p E + iq 0 ), and a branch cut at Im q 0 = −p E from ∆(p E − iq 0 ) and Γ 3 (iq 0 , p E − iq 0 ). The pole contribution to (68) is According to (32), first the analytic continuation for k E and then for p E must be done, leading to The branch cut contributions can be calculated in complete analogy to the steps leading up to (33), finding [29] lim p 0 →0 Both contributions are manifestly real, proving the above statement that F (E q ) is real. Using (24), (59), after integrating this gives Note the typo in [19,Eq. (103)]. Hereρ D denotes the angular average with weight (67), cf. Eqns. (59), (60). Using the orthogonal polynomials defined for the quadrature (34), it is convenient to parametrize so that (72) has the solution b n = A −1 nm a m where The values of A nm are calculated numerically.
In the R4 resummation, one needs the 4-vertex (16) in the resummation of the three vertex (17). Writing the three vertex (17) using a similar momentum convention as in (66), one has where W (P, Q, K) = D(P − K − Q) + 2N Γ 4 (Q, P − Q, K, P − K). In the large N limit, the four-vertex is given from (16) as There are no pinching poles in this expression, so the thermal width in the propagators ∆ can be neglected to leading order in large N. Using (30) to calculate thermal sum, one gets contributions from each of the propagators in the loop. These can be calculated by using (46) sequentially for each propagator, finding where r E = ir 0 . This expression is to be evaluated inside the resummation for the three vertex (76). The propagators ∆(Q)∆(P − Q) and the vertex Γ 3 (Q, P − Q) have branch cuts for Imq 0 = 0, Imq 0 = −p E . In addition to these, the thermal sum over q E also picks up the singularities of Γ 4 , which from (78). These are singularities on the q 0 real axis (for ∆(R − Q), D(Q + R)), on the line Imq 0 = −p E (for ∆(R + P − Q), D(Q + R − P )), as well as on the line ). For the branch cut at Imq 0 = 0, there is a pinching pole contribution to Γ 3 (K, P − K) which arises when analytically continuing q E → iq 0 − 0 + . After performing the q 0 integration, the analytic continuation for k E , p E are fixed by (32). So one finds that for the branch cut at Imq 0 = 0, Γ 4 must be analytically continued by taking in this sequence. For the branch cut at Imq 0 = −p E , the pinching pole contribution to Γ 3 (K, P − K) arises when analytically continuing q E → p E + iq 0 + 0 + , and contributes with a minus sign. Finally, the contribution for Imq 0 = −k E are simple poles which also give rise to a pinching pole contribution to Γ 3 (K, P − K). Adding the two analytic continuations of Γ 4 for the branch cuts, one finds that the terms proportional to ρ D in (78) cancel (matching the observation in Ref. [19]). Needing only the p 0 → 0 limit for the viscosity, one obtains branch−contributions where D R (R) is the retarded D correlator. This is to be contrasted with the branch-cut contributions for D(P −K −Q) from (76), which is −2iρ D (Q − K). In addition to the branch-cut contribution, there is the pole contribution arising from ∆(R + K − Q), ∆(R + K − Q). After shifting q 0 → q 0 − k 0 , the pole contribution gives whereas the corresponding pole contribution from D(P − K − Q) is 2iρ D (Q − K). All told, the three vertex contribution (76) thus becomes One can now do the angular integrations that are part of Q , R , respectively. Writing f q ·k = f (q ·r) f k ·r + sin(2q ·r) sin(2k ·r) , the only non-vanishing angular integrals are of the form Therefore, the angular average of Λ(Q, K) becomes The corresponding expression may be evaluated numerically by performing the same expansion as in (74).

Quadrature
This section provides some detail on how to construct the orthogonal polynomials P n (x) needed for Eq. (34). Finding the roots x n from P K+1 (x n ) = 0 allows to recast any integral with f (x) a polynomial of degree 2K + 1 or lower, where w n are the weights. In practice, to find P n (x), first define Constructing polynomials as P n (x) = x n + a n−1 x n−1 + . . . + a 1 x + a 0 , the orthogonality (34)  which leads to the coefficients a n after simple matrix inversion. From the coefficients a n , one obtains the roots x n of the polynomial P K+1 (x) as the eigenvalues of the companion matrix Once the roots have been obtained, the weights w n are found from requiring (86) with f (x) = 1, x, x 2 , . . . , x K , requiring one more matrix inversion. This can be achieved efficiently using standard numerical matrix inversion packages.

Shear viscosity in the 4d O(N) model
As a sanity-check on the methodology and numerics, it is prudent to compare the results found as part of this work to the previous work in Ref. [19] for the O(N) model in 3+1 dimensions. There are several modifications in the formulae presented in this work when considering three spatial dimensions, which will be pointed out in this section.
First, the O(N) model in 4 − 2 dimensions at large N requires non-perturbative renormalization of the coupling constant which in MS is given by As a consequence, the theory possesses a Landau pole at large N, effectively rendering the theory ill defined at strong coupling. In terms of the renormalized coupling λ, the saddle point condition (9) and the entropy density for d = 4 reads cf. Eq. (10) and Ref. [15]. (Note that (89) possesses two solutions, one of which can be recognized as unphysical as long as the coupling λ 1.) The Kubo formula for the T xy correlator is the same for d=3 and d=4, cf. Eq. (50), so (3) holds for both cases. Also, expressions such as the thermal width (58) are formally the same, except that k ≡ d 3 k (2π) 3 for d=4. The vacuum contribution of the polarization tensor (61) is divergent in d=4 in dimensional regularization, but Π only appears in the auxiliary field propagator (12) in conjunction with 1 λ0 , so that (88) renders this combination finite. Effectively, this means there is a additional contribution 1 8π 2 lnμ Horizontal axis is compactified in order to fit values λ ∈ [0, ∞). For comparison, results from Aarts and Resco in Ref. [19] are shown. One finds reasonable agreement between the numerical results from Ref. [19] and this work.
The quadrature rules (34) are modified in d=4 to read ∞ m dx x 2 − m 2 n (x)P i (x)P j (x) = δ ij , and are otherwise constructed analogously to the d=3 case. The vertex resummation is formally identical to (66) except for the angular weight (67) which for d=4 becomes also known as the second Legendre polynomial. As a consequence, the shear-viscosity over entropy density ratio is modified from (33) to read in d=4 Evaluating (93) using the same numerical techniques as for d=3, one finds the results shown in Fig. 2. As can be seen from this figure, there is reasonable agreement between the present work and the published results from Ref. [19], serving as a sanity-check of both the algorithm and the numerical technique used here.