Bulk Viscosity at Extreme Limits: From Kinetic Theory to Strings

In this paper we study bulk viscosity in a thermal QCD model with large number of colors at two extreme limits: the very weak and the very strong 't Hooft couplings. The weak coupling scenario is based on kinetic theory, and one may go to the very strong coupling dynamics via an intermediate coupling regime. Although the former has a clear description in terms of kinetic theory, the intermediate coupling regime, which uses lattice results, suffers from usual technical challenges that render an explicit determination of bulk viscosity somewhat difficult. On the other hand, the very strong 't Hooft coupling dynamics may be studied using string theories at both weak and strong string couplings using gravity duals in type IIB as well as M-theory respectively. In type IIB we provide the precise fluctuation modes of the metric in the gravity dual responsible for bulk viscosity, compute the speed of sound in the medium and analyze the ratio of the bulk to shear viscosities. In M-theory, where we uplift the type IIA mirror dual of the UV complete type IIB model, we study and compare both the bulk viscosity and the sound speed by analyzing the quasi-normal modes in the system at strong IIA string coupling. By deriving the spectral function, we show the consistency of our results both for the actual values of the parameters involved as well for the bound on the ratio of bulk to shear viscosities.


Introduction and summary
The wide and thoroughgoing experimental programs pioneered at the Relativistic Heavy Ion Collider (RHIC) and pursued at the Large Hadron Collider (LHC) offer a unique opportunity to study properties of a most exotic state of matter: the quark-gluon plasma (QGP). Although there is a common agreement that droplets of QGP are produced in heavy ion collisions in the pursued experiments, a unequivocal and quantitative determination of the properties of such a state is still the topic of much research. The time evolution of the plasma, its transport properties, the parameters of the transition to the confined phase, are some of the features that are currently being addressed along with many others. The difficulties in extraction of QGP properties owe much to the fact that the excited nuclear matter produced by colliding heavy ions at currently achievable energy scales is strongly coupled. Accordingly, the applicability of known fundamental methods and approaches to study the system in this regime is very limited, and hence all obtained findings in this limit have to be examined critically. On the other hand, this situation may also provide an opportunity to explore new technical facets of known tools and to explore new directions.
One of the methods that have proven useful to study the properties of QGP in the domain accessible experimentally is viscous hydrodynamics -a low-frequency long-wavelength effective theory. The application of the hydrodynamic framework to heavy ion collisions [1,2,3,4,5,6,7,8,9] and its use in the interpretation of a wide range of experimental observables [10,11,12,13] allowed to conclude that the experimentally produced QGP is a strongly coupled system. In particular, the studies on the hadronic flow and the emergence of other collective phenomena in the hydrodynamic description of QGP were taken as an indication of its fluid-like nature. Moreover, the success of hydrodynamics seemed to necessitate a fast nearthermalization of the QGP. All these arguments led to the conclusion that the created quark gluon plasma must be strongly coupled [14]. For reviews on hydrodynamic applications and formulations see [15,16,17,18,19].
Another powerful tool to study systems in the limit of strong 't Hooft coupling originated with the discovery of the AdS/CFT correspondence [20]. Even though it was hydrodynamic predictions and analyses that provided the empirical evidence that the shear viscosity to entropy density ratio is small [21], it was the AdS/CFT conjecture that established an analytical bound of η/s = 1/4π [22] 1 .
Transport coefficients are valuable elements of the hydrodynamic description as they carry information about the microscopic properties of a medium. In the case of the shear viscosity of strongly interacting matter, numerous phenomenological studies, the AdS/CFT result, kinetic theory calculations in the high-temperature weakly coupled regime of QGP η ∝ 1/(g 4 Y M log(1/g Y M )) [24], and non-perturbative estimates [25] allow a schematic global understanding of the shear viscosity to entropy density ratio. It is known that shear viscosity is large in the perturbative, hightemperature, limit, smaller near the phase transition temperature [26,27], and large again in the confined, pion gas domain [28]. However, the physics of the bulk viscosity is less satisfactorily understood. There are strong indications that bulk viscosity behaviour follows a trend opposite to that of the shear viscosity. In the limit of hightemperature QCD, bulk viscosity was found to have a very small value [29]: this is to be expected, as the coefficient of bulk viscosity can be written as a correlator of the trace anomaly (see Section 2.1), and QCD is known to be approximately conformal at high temperatures. Although it may seem that, in the very large coupling regime, a direct application of AdS/CFT techniques to bulk viscosity exploration is not relevant as the conjecture relies on the N = 4 super Yang-Mills theory (which is perfectly conformal and where the bulk viscosity vanishes identically), this is not quite true. Approaches based on holography in fact have proven useful by providing a lower bound on the ratio of bulk to shear viscosities: ζ/η ≥ 2(1/3 − c 2 s ) [30] 2 . In the vicinity of the transition from QGP to hadronic degrees of freedom, the bulk viscosity should, in principle, be calculated from the equation of state extracted the lattice data [33,34]. It is expected to be proportional to the trace anomaly ( − 3P )/T 4 and hence be notably peaked. Various investigations, both formal and phenomenological [35,36,37,38,39,40], confirm this expectation. Recently, it was demonstrated that the presence of a coefficient of bulk viscosity is important in hydrodynamical simulations, as it has a significant impact on the elliptic flow coefficients [40,41,42] and other heavy-ion observables, strongly interacting and otherwise [43,44,45,46,47]. However, it is fair to write that the precise behavior of the bulk viscosity for systems in extreme conditions of temperature and density is not yet firmly established and therefore needs further studies.
Understanding the behavior of bulk viscosity and knowing how it changes when the coupling strength varies is important for several reasons. First, bulk viscosity is an inherent property of nonconformal systems, and finite-temperature systems obeying QCD are good examples of such environments. The behavior of the bulk viscosity is fixed by parameters that break conformal symmetry. These include, at least at the perturbative region, finite masses of plasma constituents and the Callan-Symanzik β-function which expresses the coupling constant as a function of an energy scale [29,31]. Equivalently, these parameters enter the definition of the speed of sound, and bulk viscosity can be conveniently expressed as some function of 1/3 − c 2 s as well. From the phenomenological point of view, expressing bulk viscosity via the speed of sound is practical as this enables a direct connection with the lattice QCD equation of state. Second, bulk viscosity plays an essential role in the hydrodynamical description and modelling of hot and dense strongly-interacting matter. One could attempt to compute the coefficient within a theory which captures the microscopic interactions, and then insert it into hydrodynamic equations. Alternatively, fluid dynamics may be viewed as an effective theory of the long wavelength behaviour, and its transport coefficients are to be extracted empirically. Either way, viscous hydrodynamics serves as a powerful tool to investigate the strongly coupled nuclear medium produced in RHIC and LHC experiments. It provides information on the dynamics of the plasma, informs how the plasma evolves and also helps to extract, or at least constrain, other plasma characteristics. In addition, bulk viscosity studies have the potential to further the development of new theoretical methods to study the conformal anomaly of QCD. Because of different system dynamics at different coupling regimes one may expect a different dependency of the bulk viscosity on the factor 1/3 − c 2 s . This is what was observed by comparing the bulk to shear viscosity ratios at perturbative and very strong-coupling limits: ζ/η ∝ (1/3 − c 2 s ) 2 [29] and ζ/η ∝ (1/3 − c 2 s ) [30] respectively. Analyzing this difference is one of the main objectives of our studies.
We examine the bulk viscosity of systems governed by the SU(M ) theories with the interaction strength determined by the 't Hooft coupling λ = g 2 Y M M , where g Y M is the gauge coupling and M is the number of colors. The 't Hooft coupling may be thought as an effective coupling of QGP. We distinguish here 3 regions of the 't Hooft coupling: the weak coupling region, the intermediate coupling region (near the phase crossover temperature), and the strong (infinite) coupling one. In each region a different microscopic approach is applicable. The extreme limits are discussed comprehensively while the intermediate coupling part includes a brief summary and discussion on conceptual difficulties preventing one from determining bulk viscosity in this domain.
As already mentioned, the weak-coupling studies on bulk viscosity for QCD were done within kinetic theory in [29]. In our work we intend to adjust the kinetic theory result to the 't Hooft coupling. In this way we provide the form of bulk viscosity which can be directly confronted with its strong-coupling counterpart discussed via string theory methods. In this approach quark contribution is always suppressed by a factor 1/M and may be neglected in the leading order analysis. Kinetic theory [24,48,49] is an effective theory which is commonly and successfully used to compute transport coefficients. Its correspondence to fundamental microscopic theory was directly shown for scalar theory in [31] and then also for QED [50,51]. In this manuscript we undertake the task to justify the validity of kinetic theory for SU(M ) theory by providing power counting of microscopic processes contributing to the collision kernel of the Boltzmann equation. Since a derivation of the transport equation from diagrammatic representation of any non-Abelian theory is highly non-trivial, we intend to present a procedure on how to represent the collision kernel diagrammatically. We discuss how the pinching poles and nearly pinching poles control power counting of elastic and inelastic processes, respectively. The consequences of soft physics on power counting are emphasized. We also show how the integral equations emerge by discussing all topological structures of planar diagrams contributing to them. We believe that this examination may provide solid arguments to prove an equivalence of the Boltzmann equation with the analysis based on the loop expansion.
The intermediate coupling region is considered mostly to summarize the status of studies of bulk viscosity done with microscopic analyses. The bulk viscosity in this regime can be obtained if one can extract the low frequency behavior of the corresponding spectral density [35,36,37,52,53,54,55]. Although these approaches provide some constrains, they do not yet allow definite conclusions on the behaviour of the bulk viscosity to be drawn. We briefly discuss the difficulties.
On the other hand, the strong 't Hooft coupling behavior of bulk viscosity is an interesting playground to study sting theory and gauge theory because of the use of gauge/gravity duality. In fact since the bulk viscosity should truly be studied in a theory with running couplings, the famed AdS/CFT duality is not very useful, as discussed above. Going beyond CFT will require us to find the right gravity dual to answer any questions related to running couplings, and especially questions related to bulk viscosity. The gravity dual that we seek has been first proposed in [56,57] and the full UV completion was given from the type IIB side in [58,59,60] and more recently from the type IIA side in [61].
At this stage one might ask as to how a gravitational background, which has hitherto no connection to gauge theory, could in principle enter the picture to help us solve strongly coupled system like the one that we concentrate on here. There are two ways to answer this question, but none are completely satisfactory. The first one is to relegate this to the magic of duality. However this duality is special because all dualities studied so far have either been between two different gauge theories or between two different supergravity theories. There has never been a duality between a gauge theory and a gravitational theory before AdS/CFT [20].
The second one is to view the gauge theory as to be somehow contained inside a gravitational background. To elaborate this viewpoint, let us consider a four-dimensional Minkowski spacetime. This serves as an arena for gauge theory interactions, and for simplicity we decouple all gravitational interactions by tuning the Newton constant. The gauge theory interactions can happen at various energy scales, and we can assume that a specific slice of four-dimensional Minkowski spacetime is associated with a specific energy configuration. We can stack up all the slices together such that the low energy slices are at the bottom and the high energy slices are kept on top of one another in an increasing fashion. Clearly the topmost slice will be at infinite energy.
The above construction immediate provides a five-dimensional space and if we assume the energy direction to be parametrized by a radial coordinate r, then at r = 0 we have IR physics and at r → ∞ we have UV physics. This is also by construction a five-dimensional gravitational background, and by this simple assumption we seem to have got a five-dimensional gravity theory that captures the dynamics of the fourdimensional gauge theory from IR to UV! Of course this is a very simple construct and does not answer all questions related to gauge/gravity duality but it is instructive to see how two seemingly unrelated physics, one of gauge theory and other of gravity, may be united in a framework like above.
A few quick checks may be easily performed at this stage. If the gauge theory is a CFT, i.e scale independent, then the slicing idea will tell us that we need not worry too much of the physics at any scale r, and instead study the dynamics of the corresponding gauge theory from the boundary at r → ∞. Of course this is what AdS/CFT dictionary tells us, but what is lacking in our simple construct is the justification that the gravitational background is indeed an AdS 5 space. Maybe the idea of scale invariance, combined with decoupling and the supergravity EOMs could uniquely fix that, but this has not been checked.
On the other hand if we are dealing with a gauge theory that is not scale invariant, then every point on the slice matters. At every r we have the corresponding gauge theory dynamics at that scale 3 . Indeed in the Wilsonian sense at this scale all high 3 This argument entails the fact that if we keep r fixed and move along the remaining fourdimensions, nothing should change. However we can envision more generic scenario where the energy scale is mapped to a certain combination of r and the other three directions. In this case the Wilsonian effective action will be sensitive to where we are on a given slice. Of course it should be possible to redefine coordinates in such a way to find a new radial coordinate that will again correspond to the energy scale. For this paper we will however stick to the simplest case where r is mapped to the energy scale, r c to the UV cut-off, and r h , the horizon radius, to the temperature. energy degrees of freedom are integrated out and we are left with a set of relevant, marginal and irrelevant operators. This is of course the premise of our construction in [58], and the UV completion in [59,60] is done by introducing new degrees of freedom from the so-called Region 2 of [59] onwards.
Other checks, that include the exact mapping of the gauge theory operators to supergravity states, are much harder to perform and in fact the dictionary for gauge/gravity duality for the non-conformal case is not yet fully developed compared to the conformal case. Nevertheless one thing is for sure: to have any control on the computations on the supergravity side we need small g s . For a background with a constant dilaton − an example would be the Klebanov-Strassler background [56] − a little bit of numerology can tell us that g s may be made arbitrarily small 4 .
There is also an additional requirement of large number of colors. For a SU(M ) gauge theory, the corresponding supergravity theory will make sense if λ ≡ g s M is very large. In this limit all computations can be restricted to classical supergravity alone, and stringy corrections can be entirely ignored. However if we want to study an actual large M QCD we will have to explore string coupling g s = O(1). How can we ignore stringy corrections now and restrict ourselves to supergravity alone?
A way out of this conundrum was first proposed in [63] by performing a sequence of two stringy dualities: mirror transformation and M-theory uplift. The mirror transformation is a special kind of duality that takes a type IIB background to a type IIA background by simply interchanging the Kähler and the complex structures of the internal manifolds on both sides of the duality. In [63] this was implemented by performing three T-dualities along the isometry directions of the internal manifold in the type IIB side [64]. Being T-dualities they do not change the behavior of the dilaton too much, and therefore takes a weakly coupled background into another weakly coupled one.
The second duality is when we increase the type IIA coupling. At strong coupling a new internal direction opens up and the theory goes to eleven-dimensional M-theory where the dynamics is now miraculously governed by eleven-dimensional supergravity. All the type IIA stringy corrections are now captured succinctly by classical supergravity analysis in M-theory [65], and therefore g s = O(1) can again be studied using supergravity, albeit from eleven-dimensions. Such a dual description was termed as the MQGP limit of thermal QCD with large number of colors in [63].
The above considerations tell us that the strong 't Hooft coupling regime may be studied from the perspectives of both weak and strong string couplings. In the presence of N f flavors, it means we are exploring both g s N f → 0 as well as g s N f = O(1) limits 5 . This in turn boils down to saying that we can have analytic control on the transport coefficients − and here we will concentrate only on bulk viscosities − for pure glue as well as for flavored large M thermal QCD. Section 4 of the paper will therefore be dedicated to studying the bulk viscosity using weak string coupling and with vanishing number of flavors, whereas section 5 will be dedicated to studying the bulk viscosity using the other limit, namely strong string coupling and non-vanishing number of flavors.
There is yet another limit where we can remain at weak string coupling, but explore strong YM coupling. In the type IIB side such a scenario becomes possible once N f flavor degrees of freedom are switched on. That this could happen is a consequence of two conspiracies: one, the dilaton picks up O(g s N f ) corrections forcing it away from being a constant, and two, the NS 2-form field, through the vanishing two-cycle on which we have the M wrapped D5-branes, also picks up O(g s N f ) corrections. These corrections provide additional structure to the already non-constant field, but more importantly they add to the dilaton factor constructively to provide the full structure of the YM coupling.
Interestingly, from either of these limits at strong 't Hooft coupling, the ratio of the bulk to shear viscosities remains proportional to linear power of 1 3 − c 2 s . The difference however lies in the precise coefficients that control the lower bounds at weak and strong string couplings. For example at strong string coupling, the lower bound is almost 9 times bigger than the Buchel bound [30] as we will discuss in section 5. Of course nowhere we see any violation of the Buchel bound, so presumably the violation can only occur once we dimensionally reduce the four-dimensional theory to two-dimensions. This is much like the scenario presented in [66], but we will not discuss it any further here.
What we will discuss however is the appearance of the linear power of the deviation factor, 1 3 − c 2 s , when we study spectral function using the weakly coupled type IIA theory. The spectral function is an important aspect in the study of QGP, and its derivation is rather complicated at weak 't Hooft coupling. At strong 't Hooft coupling there is a way to derive it from the gauge/gravity duality, but the derivation is technical and involves various manipulations of the background. Nevertheless an answer can be found in the present set-up and the final result shows a linear dependence on the deviation factor. In the limit of vanishing frequency, the result matches well with actual QGP, despite the presence of a large number of colors. Such a success points towards some inherent universality, and it will be interesting to explore this further.

Organization of the paper
The paper is organized as follows. In section 2 we study bulk viscosity at weak is because we always want to keep the combination (g s N f ) k gsM 2 N m << 1 even for m = 1 and k ∈ Z. See also footnote 19. 't Hooft coupling. After short introductory remarks on the definition of bulk viscosity and the applied microscopic theory, in section 2.1 we discuss the Kubo formula which provides a general and first-principle method of computation of the coefficient. In section 2.2 we briefly summarize results on the leading order bulk viscosity calculation performed within kinetic theory by solving the Boltzmann equation. Sections 2.3 − 2.6 contain a diagrammatic, but qualitative only, analysis which is to justify validity of the effective kinetic theory formulation for transport coefficients studies. In section 2.3 we consider the one-loop diagram to find a typical size of the bulk viscosity. This step shows also that fermionic contributions are subleading in favor of the gluonic ones. Then, in section 2.4, the power counting of the relevant self-energies is done. Section 2.5 is devoted to an evaluation of typical sizes of multi-loop diagrams which represent scattering processes. Both particle number conserving and particle number changing processes are studied and the role of the soft physics is emphasized in subsections 2.5.1 and 2.5.2. In section 2.6 a schematic form of the relevant integral equations needed for a diagrammatic bulk viscosity computation is presented.
In section 3 an intermediate coupling regime is discussed. The section consists of a brief overview of literature on the approaches aiming at an extraction of bulk viscosity from lattice QCD results by studying mostly QCD sum rules and finding constraints on the spectral density. The difficulties in the quantitative determination of the bulk viscosity are pointed out.
The strong coupling results are discussed in sections 4, 5 and 6. In section 4, the weak string but strong 't Hooft coupling regime is discussed. We start by giving a detailed description as to where the string theory techniques would fit in in the study of bulk viscosity. The various domains of compatibility as well as the UV completion are emphasized and the consistency of the background is shown from both type IIB and its dual type IIA pictures. In section 4.1, a slightly simplified background is taken to quantify various parameters associated with the computation of bulk viscosity. For example, one of the important parameter is the fluctuation associated with the vielbeins. This is elaborated in section 4.2. The fluctuation modes can be divided into positive and negative frequencies, and we show that there are pieces of the fluctuations, called p nk , that are related to certain sources ∆ (n) ab in the gravity dual picture. The analysis of the sources is rather complicated and therefore in section 4.2.1 we first take a toy example to study the equations connecting p nk fluctuations with the ∆ (n) ab sources. The toy example is based on a simplifying constraint, and using this the simplest zero and the non-zero modes of the fluctuations are shown to satisfy equations that relate them to the sources. In section 4.2.2 we go beyond the simple toy example by studying the equations governing the fluctuating modes in a generic setting. As before, the zero and the non-zero modes satisfy equations relating them to certain sources.
Once we have the fluctuations, we can use them to compute the transport coefficients. In section 4.3 we perform two important computations: one, the sound speed, and two, the ratio of the bulk to the shear viscosities. The former is given by an equation which takes into account not only the scale dependence of the temperature, but also the background fluctuations. Needless to say, the ratio of the bulk to the shear viscosities should depend on the sound speed, and we elucidate this by first computing the precise ratio and then showing that the ratio is indeed bounded below by the deviation of the sound speed from its conformal value.
The remaining two sections are devoted to studying bulk viscosity at strong string and strong 't Hooft couplings. The first, i.e section 5, has to do with obtaining a Buchel-like bulk-viscosity-to-shear-viscosity bound by looking at scalar modes of metric perturbations and the associated quasi-normal modes. The second, i.e section 6, has to do with obtaining the same result from spectral functions. Here is a more detailed plan of these two sections.
In section 5, we first briefly review the Strominger-Yau-Zaslow (SYZ) type IIA mirror of [58]'s top-down type IIB holographic dual of large-N thermal QCD, as well as its M-theory uplift as constructed in [63]. This is followed by a discussion on obtaining the EOM for a linear combination of scalar modes of metric perturbations gauge invariant under infinitesimal diffeomorphisms and obtaining the associated quasi-normal modes in section 5.2; it is noted that with a non-zero bare resolution parameter, the horizon turns out to be an irregular singular point, a fact that proves in fact to be quite helpful in obtaining the aforementioned quasi-normal modes. In section 5.3, we show that one cannot avoid non-normalizable modes if one were to turn off the bare resolution parameter resonating well with similar non-normalizable perturbation modes obtained in section 4. A Buchel-like bound for the ratio of the bulk and shear viscosities in terms of the linear power of the deviation of the square of the speed of sound from its conformal value is finally obtained, both for N f = 0 and N f = 0 in section 5.4.
In section 6, we follow a different route − that of spectral function involving correlation function of gauge fluctuations in background value of gauge fields on the world volume of the flavor D6-branes of the aforementioned SYZ type IIA mirror. In section 6.1, we obtain the background value of a D6-brane world volume gauge field A t (r), r being the radial coordinate and set up the EOM for fluctuations about the same. We obtain and explicitly solve the EOM − there turns out to be only one linearly independent EOM − in the zero-momentum limit in section 6.2. From the on-shell action, the gauge-fluctuation correlation function and hence the spectral function per unit frequency in the vanishing frequency limit, is worked out in section 6.3 and it is explicitly seen that the difference between the same at non-zero and zero temperatures is precisely proportional to the deviation of the square of the speed of sound from its conformal value. In section 6.4, we argue that unlike sections 6.1 − 6.3 wherein one had considered weak-string-coupling strong-'t-Hooft-coupling limit, the result of section 6.3 goes through even for the strong string and strong 't Hooft couplings' limit. We argue therein that the g s → 0 limit alongwith non-trivial B-field along the vanishing two-cycle conspires to produce a g 2 Y M in the gauge theory side that is no longer a small number.
Finally in the appendices we discuss three topics. The first one is on a gauge invariant combination of the scalar modes of metric fluctuations. Such a combination is useful to study the quasi-normal modes.The second one is on the derivation of the on-shell action and Green's function required to study the spectral function. The third one is on an estimation of the horizon radius.

Bulk viscosity at weak 't Hooft coupling
When the system exhibits a small deviation from thermal equilibrium, its evolution is well described by the equations of hydrodynamics. These are given in terms of conservation laws of currents accompanied by the equation of state. Here, we focus only on the energy and momentum currents which are encoded in the stress-energy tensor T µν . Its spatial part is: where ζ and η are the bulk and shear viscosities, u i is the fluid flow velocity and the metric is mostly negative. A many-body system can be driven out of its equilibrium state through uniform compression or rarefaction and both processes lead to changes in the energy density , the increase or decrease, respectively. The pressure P also changes but its change is different than that provided by the equation of state P ( ). The trace of the stress tensor carries information on changes in pressure. The deviation from the equilibrium pressure when the system is expanding or contracting is characterized by the bulk viscosity ζ: where ∇ · u is the expansion parameter. Bulk viscosity, as well as other transport coefficients, is determined by microscopic dynamics. Here we discuss how bulk viscosity emerges when the system is governed by the non-Abelian SU(M ) gauge theory with the Lagrangian: Here ψ is the quark field with M × N f degrees of freedom, where M is the number of colors and N f is the number of flavors, D µ = ∂ µ + ig Y M A µ is the covariant derivative with the gluon field A µ , which has M 2 − 1 degrees of freedom, and F µν = 1 ig Y M [D µ , D ν ] is the field strength tensor. The strength of interaction is fixed by the gauge coupling g Y M .
Classically, this theory has conformal symmetry as long as the quarks are massless. Quantum mechanically, renormalization breaks the conformal symmetry since the Callan-Symanzyk β-function is non-zero. Therefore, it is expected that the bulk viscosity of the massless SU(M ) gauge theory is directly related to the β-function. This is manifestly shown within the effective kinetic theory analysis in Ref. [29]. In the rest of our analysis, we mainly consider the large M limit. In this limit, the relevant interaction strength turns out to be the 'tHooft coupling λ = g 2 Y M M and then β ∼ λ 2 /M [67].
In principle, to study bulk viscosity comprehensively one should consider massive fermion fields, since a constant mass is a parameter that breaks conformal symmetry as well. In the light of the forthcoming discussion it is, however, not necessary here as the quark contribution will be M suppressed compared to the gluon contribution.

Kubo formula for bulk viscosity
The first-principles prescription to compute bulk viscosity is given by the Kubo formula [68]: where ρ P P (ω, k) is the spectral function of the pressure-pressure correlator and ω is the frequency of the hydrodynamic mode. In the following discussion, we will often omit the the k dependence from the correlation functions and spectral densities. The common k → 0 limit should be understood in those cases. The spectral function is related to the imaginary part of the pressure-pressure retarded correlation function: where we used G A = G * R . In the rest frame of the fluid cell, the pressure operator is given by the trace of the stress-tensorP = − 1 3T i i . Because of the energy-momentum conservation, one can easily show that the spectral functions ρ P (ω, k) and ρ (ω, k) must vanish in the same limit [19], whereˆ =T 00 is the energy density operator. For theoretical analysis, it is often more advantageous to use the trace of the full stress-energy tensorΘ/3 =T µ µ /3 = 1 3ˆ −P which is Lorentz invariant, or the more kinetic-theory-friendly combinationP =P − c 2 sˆ which reduces to −Θ/3 in the conformal limit. Here c 2 s = ∂P/∂ is the speed of sound squared. Therefore, with the retarded correlation function given in coordinate space as Here the operatorÔ can beP ,Θ/3 orP . This correlation function contains all the essential information about the physics of bulk viscosity and their structures are fixed by the Lagrangian (2.3) and thermal medium effects. Although the Kubo formula (2.6) is general, in this section we focus on the regime of the sufficiently high energy scale, where the expansion in the small 't Hooft coupling λ may be applied. We consider here the limit where the 't Hooft coupling λ = g 2 Y M M remains small and the number of flavours N f is fixed while M → ∞. In this limit one should, in principle, be able to calculate bulk viscosity perturbatively. Due to very complex multi-scale nature of the non-Abelian theory, a comprehensive quantitative computation of bulk viscosity using field theoretical tools is not an easy task. To date, a complete diagrammatic analysis of the bulk viscosity in QCD has not been carried out (for other transport coefficients of QED, see [50,51]). However, an equivalent approach to compute the coefficient is offered by using effective kinetic theory.

Bulk viscosity from kinetic theory
The foundations of the effective kinetic theory of the SU(M ) theory were formulated in Refs. [24,48,49]. The scattering processes governing transport properties of the medium are embedded in the collision kernel of the Boltzmann equation and their sizes in terms of the gauge coupling g Y M , the numbers of degrees of freedom and the Casimir operators are explicitly shown in Ref. [48]. Further this formulation was used in Ref. [29] to calculate bulk viscosity of QCD. Here we briefly summarize these results in the leading order in the 't Hooft coupling λ = g 2 Y M M . In the large-M limit, the bulk viscosity in the leading order in λ is governed only by the pure gluodynamics since quarks are suppressed by at least a factor of M . This can be clearly seen from the following analysis. Bulk viscosity depends on two factors. First, it must be proportional to the nonconformality parameter reflecting the incompressibility of the system. Second, it is controlled by the mean free path carrying the information on the microscopic properties of the medium, in particular, on the nature of interaction, and relevant degrees of freedom. From Ref. [29] one observes that the same dependence of bulk viscosity on the nonconformality parameter is obtained either for quark and for gluon contributions. The mean free path of the quark contribution and of the gluon one is parametrically the same but it is associated with the corresponding numbers of degrees of freedom, which are different. While the number of gluons scales as M 2 , the number of quarks scales as M . This dependence occurs for both the number conserving and number changing processes and can be extracted when analyzing all matrix elements and associated degrees of freedom shown explicitly in [48]. Hence we ignore the quark contribution at every step of the forthcoming analysis.
In kinetic theory one focuses on the evolution of the distribution function of relevant quasiparticles. The evolution of the gluon distribution function f (p, x) is governed by the Boltzmann equation of the form: (2.8) Since f (p, x, t) is slightly out of equilibrium it can be expressed as f = f eq +f 1 , where f eq is of the form f eq (p, f eq is therefore a function of time-space dependent quantities: β(t) being the inverse of temperature T (t) and E p (x) = p 2 + m 2 th (x) -the energy of a gluon where the x dependence appears through the thermally fluctuating mass m th (x). f 1 is the nonequilibrium correction, which includes both the action of hydrodynamic forces and the correction due to thermally fluctuating mass. C[f ] is a collision term, that contains processes involving only gluons, namely, the number conserving gg → gg scatterings and the number changing g → gg splittings. Its explicit form can be found in [48]. The left-hand side of the Boltzmann equation at the linearized order is then: ) and f 0 is the Bose-Einstein distribution function (e βEp −1) −1 . The form of the quantity q(p) is most essential as it establishes the final parametric dependence of bulk viscosity on the nonconformality parameter. It reads: The quantitym 2 is of the form: The formula (2.10) is obtained by taking into account the stress-energy conservation law, thermodynamic relations and space dependence of the quasiparticle energy. Note that as the consequence of the temperature dependence of the quasiparticle mass, given by m 2 th = g 2 Y M (T )M T 2 /6, the beta function of SU(M ) theory β λ = −11λ 2 /(48π 2 M ) arises in the formula (2.11) and, consequently, in Eq. (2.10). The β λ -function is just the parameter that breaks conformal symmetry in the system and the factor 1/3 − c 2 s , with the speed of sound squared c 2 s = ∂P/∂ , is equivalent to it through the relation: (2.12) Due to such a dependence, q(p) can be expressed in a simple form: In all formulas, terms which are suppressed by any powers of M were omitted. The form of left-hand side of the Boltzmann equation, Eq. (2.9), dictates also the form of the correction f 1 which, in turn, fixes the form of the linearized collision kernel. The correction is f 1 = β 2 f 0 (1 + f 0 )χ∇ · u, so that both sides of the Boltzmann equation are proportional to ∇ · u. By dropping this scalar factor, the Boltzmann equation can be expressed in a convenient form S(p) = [Cχ](p). Bulk viscosity may be then found as: where the matrix isC mn = 2M 2 p φ m (p)[Cφ n ](p) and the column vector isS m = 2M 2 p φ m (p)S(p), with the basis functions φ m (p) = p m T K−m−1 /(T + p) K−2 and m = 1, ..., K. The numerical procedure relies on the variational method. Since ζ ∝ S 2 ∝ q 2 , the bulk viscosity is clearly expressed by the nonconformality parameter squared, (1/3 − c 2 s ) 2 or equivalently β 2 λ and the inverted collision kernel introduces the mean free path. The final expression then scales as: where a and b should be obtained by solving the integral equation (2.14). The whole procedure of finding bulk viscosity coefficient of QCD is comprehensively discussed in [29] for different values of the number of flavors N f . One can then reproduce the dependence of bulk viscosity of the SU(M ) theory on the coupling constant λ from Fig. 1 of Ref. [29] by setting all quark masses to 0, taking N f = 0, and rescaling the coupling 4πM α s → λ. Due to the same sizes of the nonconformality parameter and the 't Hooft coupling constant squared, given by the relation (2.12), one can write: where we used the entropy density s = (P + )/T ∝ M 2 T 3 . The formula (2.16) shows that in the very weak coupling regime the leading order bulk viscosity over entropy density ratio is a linear function of the nonconformality parameter 1/3−c 2 s , up to the logarithm. This occurs due to the fact that β λ function is of the same order as the inverse of the mean free path. This behavior is characteristic for the theories when the conformal symmetry is broken only by the β λ function. These are, for example, SU(M ) in the large M limit or massless QCD. Also, the shear viscosity coefficient of QCD with the effective coupling λ was studied in [69] and the result is: with A and B being numerical constants. Combining Eqs. (2.17) and (2.15), one finds that the ratio of ζ/η is characterized by the quadratic dependence of the non-conformality parameter:

One-loop diagram and power counting
So far kinetic theory has been the only utilizable method allowing for a quantitative computation of transport coefficients of non-Abelian weakly coupled gauge theories. However, it is an effective description of quasiparticle dynamics and its equivalence to quantum field theoretical approach has not been fully shown for the SU(M ) theory.
In particular, a complete diagrammatic analysis of the bulk viscosity of that theory has not yet been carried out (for transport coefficients of QED, see [50,51]). As was shown in Refs. [31,70,71,72], the equivalence of the diagrammatic method and the kinetic theory description can be established when the ladder diagram resummation dominates the leading order result. To carry out a qualitative analysis of the weakly coupled large M Yang-Mills theory, it will be therefore enough to confirm that the planar ladder diagrams dominate in the viscosity calculations. The goal of the forthcoming subsections is to do just that for the purpose of establishing the qualitative behavior of the bulk viscosity in the weakly coupled theory. We will not, however, attempt to carry out the full analytical computation necessary for the quantitative analysis as this is beyond the scope of this work.
To perform the qualitative analysis we need to establish the necessary basic ingredients dictated by the Kubo formula (2.6). The full stress-energy tensor of the SU(M ) gauge theory is given by: To have some insight into the parametric form of the bulk viscosity and to establish a starting point for evaluating the size of microscopic processes governing its behavior, it is illuminating to consider only the kinetic terms of the stress-energy tensor, that is, the first term in Eq. (2.19). Since quarks are subleading we focus only on the gluonic contribution to the stress-energy tensor; we briefly comment on this issue later. Power counting of the gluon one-loop diagram is most conveniently accomplished using the (r, a) basis of the thermal field theory. This was shown for the scalar field theory in [73,74] and also for gauge theories in [70]. In this basis, the elementary gluon propagators are the retarded propagator G ra , advanced one G ar and the auto-correlation function, which introduces information on the medium momentum distribution, G rr = (1 + 2n B )(G ra − G ar ), where n B is the Bose-Einstein statistics. These propagators carry indices related to color, spin or the Lorentz structure, but within this analysis we will not explicitly show them.
Since all these propagators describe a propagation of a given particle in a thermal medium they are dressed with self-energies. The retarded self-energy is given by Π = Re Π − iIm Π and the retarded propagator is: , (2.20) where A g (p) carries the necessary color and tensor indices. The advanced propagator then is given by G ar = G * ra . In the weakly coupled limit, the retarded propagator has poles at p 0 ≈ ±E p −iΓ g p where the quasi-particle energy is given by E p = p 2 + m 2 th with the thermal mass m 2 th = ReΠ(p). The thermal width is given by the imaginary part of the self-energy at the on-shell momentum Γ p = Im Π(E p , |p|)/2E p . The resummed propagator can be then expressed as: (2.21) In using the propagators in the (r, a) basis to evaluate the Kubo formula, we encounter two different types of singularities: The pinching pole singularity and the collinear singularity. Both are regulated by the thermal self-energies but they complicate the power counting. In this section, we discuss the pinching pole singularity and its ramification. The effect of collinear singularity is discussed in the later section.
Using the operatorP defined below Eq. (2.5) one finds the gluonic one-loop contribution to the bulk viscosity in the pinching pole approximation as: Note that the propagator part is written in a symbolic way as all internal indices and traces over them are not written explicitly. The retarded propagator has poles at p 0 = ±E p − iΓ g and the advanced one at p 0 = ±E p + iΓ g . Hence the two poles at Re p 0 = E p , for instance, are separated by 2iΓ g p in the imaginary direction on the opposite side of the integration contour. When integrated, these "pinching poles" result in a large 1/Γ g factor leading to the following power counting: This expression requires a few comments. First, the factor in the square bracket has analogous form to the quantity q(p) found within the kinetic theory and given by Eq.
(2.10), up to the thermal mass term. The expression (2.23) is obtained, however, only from the one-loop analysis and it does not include all effects. We expect that when the Lagrangian part and the interaction terms of the stress-energy tensor operator are included in the computation, the term d(m 2 th )/d(logT 2 ) will emerge in Eq. (2.23). This term, when subtracted from m 2 th in (2.23) will be analogous to the expression (2.11) and therefore will result in the β λ function emergence, or equivalently (1/3 − c 2 s ), analogously to what was obtained within the kinetic theory. The inclusion of the temperature dependence of the thermal mass was justified in Ref. [31] and explicitly incorporated to formulate fluid dynamic equations in Ref. [75], but for scalar theories only. We expect that performing full analysis of the spectral function of the SU(M ) theory will result in this dependence of the nonconformality parameter, but we do not intend to derive it. We do focus on discussing the consequences of the presence of 1/Γ g factor in the formula (2.23), which governs the mean free path behavior. Before doing that let us point out that M 2 − 1 in (2.23) reflects the number of degrees of freedom and since the number of colors is large we will be neglecting further the constant "−1". To represent the expression (2.23) diagrammatically, it is convenient to use the 't Hooft notation [76] so that a double line corresponds to a gluon propagator and any fermion propagator is represented by a single line. In this representation power counting relies on the simple formula [76]: where L is the number of closed loops, V 3 is the number of the 3-body interaction vertices and V 4 is the number of 4-body interaction vertices. In case of a fermion occurrence there is an extra factor of N f and L f is the number of fermion loops. Using the 't Hooft notation, the one-loop diagram corresponding to the expression (2.23), together with its typical size, is depicted in Fig. 1a), where the crossed vertices stand for the insertion of the renormalized operator of the trace of the stress-energy tensor. For a comparison, in Fig. 1b) we also present the fermionic one loop with its typical size given in terms of the corresponding degrees of freedom and the fermionic thermal width Γ f being given by the imaginary part of the fermionic self-energy. Therefore, the gauge boson contribution to the correlation function at the leading order scales as M 2 /Γ g , since the diagram is made of two closed loops and the fermionic one scales as M N f /Γ f . As may be implied from Fig. 1 each factor of the thermal width is associated with the presence of a pair of propagators. Thus one observes that the fermion contribution is subleading by a factor of N f /M comparing to the gluonic one as long as the same parametric dependence in the parameter 1/3 − c 2 s holds and Γ f is of the same order as Γ g . From the kinetic theory findings in Refs. [29,48] one finds the parameter (1/3 − c 2 s ) 2 being common for gluons and fermions and we rely on this result. Given that, the estimates of the sizes of the fermionic and gluonic thermal widths are still needed. What is more, to fully estimate the leading order 2-point correlation function for bulk viscosity, one also needs to know typical sizes of the corresponding thermal masses, which are essential for number changing processes.

Self-Energy Power Counting
Both the real part and the imaginary part of the self energy plays a key role in the calculation of transport coefficients. The role of the imaginary part (the thermal width) as a regulator for the pinching pole singularity has already been discussed in the previous section. The role of the real part (the thermal mass) is to regulate the infrared and collinear singularities that occur at finite temperature. Hence, the size of the thermal mass defines the soft scale while the temperature itself defines the hard scale. In QCD, we know that the thermal mass is of O(g Y M T ) while the thermal width is of O(g 4 Y M T ) when the particle momentum is hard (For instance, see [50]). In the large M limit, these need to be re-expressed using the 't Hooft coupling.
The thermal mass is determined by the real part of the self-energy of a particle. In our case, the leading order contribution comes from one-loop diagrams. The corresponding diagrams contributing to the gluon thermal mass in the double-line representation are shown in Fig. 2. For a systematic comparison, the fermion leading contribution to the real part of self-energy is shown in Fig. 3. The coupling dependence comes from counting the interaction vertices and number of degrees of freedom using the formula (2.24). As in case of one-loop diagrams contributing to the spectral density, when M → ∞ the gluon loops, Figs. 2a) and 2b), dominate over the fermion one by a factor of M/N f . Thus, the leading order of the gluon thermal mass as well as the fermion one ( Fig. 3) in the large M limit is: For explicit expressions, see [77,78,79]. The imaginary part of the one-loop self-energy vanishes when bare propagators are used due to kinematic constraints. It does not vanish when the resummed propagators are used, but that is equivalent to the two-loop self-energy which we discuss next. The relevant two-loop self-energy diagrams and their sizes for both gluons and quarks are shown in Fig. 4 and 5, respectively. It is then apparent again that gluon contributions dominate over the fermion ones by a factor of M for both the gluon and quark self-energies. The size of the imaginary parts of the self energies is: which leads to the thermal widths being of the same size, Γ g ∼ Γ f ∝ λ 2 T in the leading order. This is already enough to justify that quark loops do not have to be considered any more since the gluon contribution to a given quantity is always M times bigger than the quark one, up to the factor of N f , which is fixed constant and much smaller than M . This justifies omitting all quark contributions in the forthcoming analysis. One can then observe that the typical size of the propagator part of the correlation function is: up to the logarithm. The parametric estimate of the self-energy is of significant importance since it controls power counting of the scattering processes establishing the bulk viscosity coefficient. Consequently, as we discuss in the next subsection, it is the form of the self-energy at the one and two-loop order that controls the form of the collision kernel of the Boltzmann equation. In particular, by studying the self-energy one is able to find which processes contribute to the collision kernel and what is their sensitivity to different scales. In contrast to the shear viscosity, where the hard scale dominates its typical size, the bulk viscosity is sensitive to the soft scale as well. Since the soft scale is dictated by the size of the thermal mass, whenever we refer to it we mean

Diagrammatic justification of the processes contributing to the Boltzmann equation
Although the parametric estimate of bulk viscosity can be found by considering the one-loop diagram of the correlation function, in the thermal medium infinite number of multiple processes have to be included in the leading order. Equivalently, the infinite number of relevant loops need to be resummed. In the kinetic theory formulation this procedure is supposed to be captured by the collision term of the Boltzmann equation. The equivalence between the two approaches can be established by showing that at a given order of the coupling constant microscopic processes obtained in the diagrammatic representation have their counterparts in the collision kernel of the Boltzmann equation. Here we discuss this issue. There is a twofold source of the need for resummation of infinite number of diagrams. Each of them is related to the presence of a different type of singularity. The first case has already been discussed and it is the pinching pole singularity, regulated by the thermal width Γ g , where no other singularity occurs. The diagrams reflecting this type of singularity correspond to the number conserving, 2 → 2, processes. Then, within the one loop already discussed, any number of gluon exchanges between the side rails is possible. Any possible insertion of the permissible gluon exchanges, meant as rungs, is of the order of λ 2 and it is compensated by an accompanying 1/Γ g factor coming from pinching poles of the pair of retarded and advanced propagators. Infinitely many such combinations are possible. The other type of singularity, characteristic for gauge theories, is the collinear singularity associated with the small angle between the scattering and scattered particles. This type of singularity governs the number changing processes, 1 + N → 2 + N , where N is a number of hard particles taking part in a splitting of one hard gluon into two hard gluons. The collinear processes contribute to the bulk viscosity computation at the same order of the coupling constant as the number conserving processes. The splitting process occurs when a hard gluon traversing the medium interacts with another hard gluon via a soft momentum exchange and then emits an additional hard gluon so that both the emitted and the emitting gluons move almost collinearly within an angle θ ∼ O( √ λ). The collinear region of propagating particles is always associated with the corresponding product of retarded and advanced propagators, which, if not dressed, has a singular behavior. The collinear singularities, similarly to the pinch singularities, are regulated by the thermal width, but in this case, the soft scale fixed by the thermal mass plays an essential role as well. As discussed in detail later, all hard gluons taking part in the process can interact infinitely many times via the soft exchanges with the thermal background and they have to be coherently resummed.
Since each type of singularity involves a resummation of the corresponding set of infinitely many diagrams, there are two integral equations that need to be solved, each of them associated with the corresponding type of singularity. We first focus on the physics of number conserving processes which can be represented by a set of diagrams involving the pinch singularities only. The case of collinear processes is discussed later.

2.5.1
The case without collinear singularities 2 → 2 processes are represented by rungs, which have to be inserted in the one-loop spectral function and then resummed. Finding the structures of rungs is not trivial and to do so one needs to rely on a few constraints: Ward identities, power counting and kinematic boundaries. The most essential constraint is imposed by the Ward identities which provide relations between the effective vertex and dressed propagators and also dictate the way to maintain gauge invariance. Thus, the Ward identities should be used to obtain relations between the full on-shell imaginary part of the selfenergy and possible rung insertions. This is discussed for QED transport coefficients in [50,51] and for any SU(M ) one should expect similar relations. Accordingly, one can reproduce corresponding rungs by cutting the two-loop self energy diagrams in all possible ways and then opening one line in every diagram in all permissible ways. In Fig. 6 we show one schematic example. The two-loop diagram in Fig. 4 a) is cut through the two loops and the cross denotes lines which are open. In this way one gets two possible topologies of a rung. The lines which are cut, but not open, represent particles put on shell. Similarly, the external lines represent thermal on shell excitations. By opening the cut lines one reproduces 2 → 2 scattering processes shown in the right column in Fig. 6. Therefore, the first row shows how to obtain t and u channels of scattering events, they are obtained from the same rung but with a different momentum flow along the rung, which is not shown explicitly. The second row presents how one gets the s channel. The topological structures obtained by using this procedure are depicted in Fig. 7. First, only the diagrams a)-e) in Fig. 4 of the two loop gluon self energy have to be examined when looking for the topological structures of rungs. The diagrams g) and h) are tadpole diagrams with one-loop corrections and they contribute to the real part of self energy. The diagram f) is a one loop diagram with a tadpole correction and it also provides a contribution to the resummed propagator. Therefore, in general, all diagrams containing tadpoles do not have to be investigated any more for this qualitative analysis. For the power counting analysis one should include all rungs which have g 4 Y M factor coming from the interaction vertices and which have one closed loop contributing a factor of M . The other factor of M , expected for the proper 't Hooft coupling order is obtained when the external lines of rungs on the right-hand side are joined with other rungs or with each other. What is more, in Fig. 7 we present all topological structures arising only from the use of the Ward identity. The final relevant contributions to the kernel of the integral equation can be, however, found by using kinematic constraints and power counting arguments. The kinematic constraints are schematically represented by the dashed and the dotted lines. The dashed lines represent the cuts through rungs which are allowed by kinematics. All possible dashed lines that lead to the number conserving processes are shown. The dotted lines, although coming from the Ward identity analysis, reflect forbidden processes since one on-shell massless particle cannot decay into two on-shell masslesss particles. Therefore, the structures j) and k) in Fig. 7 do not contribute to the kernel. Also, any diagram involving crossed lines do not have to be considered here since it is always suppressed by some powers of M . Accordingly, only diagrams a)-i) constitute the kernel of the integral equation determined by the pinching singularity and they all contribute at the order O(λ 2 ). It is also easy to observe that all the contributing rungs with the associated cuts may be converted to reproduce matrix elements in the scattering amplitude defining the collision term of the Boltzmann equation. The rungs a) and b), shown also in Fig. 6, represent a contribution to the scattering amplitude squared given by t, u, and s channels. The diagram c) leads to the respective contribution from the contact interaction. The diagrams d)-i) reflect all possible interference terms.
At hard scale all allowed rungs contribute at the order O(λ 2 ), where it is enough to count the number of interaction vertices and the number of color loops. To see the relevant M -dependence it is more convenient to count closed loops of the spectral function shown in Fig. 1 a) with the rungs inserted.
When the soft scale starts to play a role power counting of the diagrams presented in Fig. 7 can change and not all diagrams are of the same size. The rung c) is not affected by the soft physics since all lines must be hard and on-shell. To do power counting of other diagrams with momenta of the order O( √ λT ), we use the (r, a) basis. It is important to notice that there can be many diagrams of the same topology but with different a and r assignment. In Fig. 8 we show exemplary diagrams from Fig. 7, where the a and r positions and the momentum convention are shown. The expression corresponding to the rung a) is: 4 G ra (l)G ar (l)G rr (k + l)G rr (p + l). (2.28) The size of the rung is estimated as follows. All incoming and outgoing momenta are hard, k ∼ p ∼ O(T ), and on-shell, while the loop momentum l is soft l ∼ O( √ λT ) and off-shell. In this case both G ra (l) and G ar (l) propagators are of the order O(λ −1 T −2 ). Additionally, since both G rr (k + l) and G rr (p + l) are on-shell they contain delta functions to maintain energy-momentum conservation. When the loop momentum integration is performed the phase space d 4 l combined with the delta functions reduces to d 2 l, which is O(λT 2 ) when l is soft. Combining all these factors one gets λ 2 from the explicit interaction vertices, λ from a phase space suppression and λ −2 from two soft propagators which make this rung to be of the order O(λ). The rung has therefore a different size at the soft scale than at the hard one, which is due to the Coulomb divergence characteristic for these scattering processes. This is, however, only a superficial difference since there is an additional mechanism which makes this rung contribute to the integral equation at the expected O(λ 2 ) order. The best way to see it is to refer to the 2 → 2 collision kernel of the Boltzmann equation [49], which is: (2.29) where the functions χ represent a small nonequilibrium deviations from the Bose-Einstein distribution function. When 2 → 2 scattering processes represented by rungs a) and b) in Fig. 7 occur in the medium via the soft momentum exchange, which is when k − k = l, with l ∼ √ λT , then one encounters the following cancellation between the χ functions: (2.30) The prescription dictated by the Kubo formula has similar structure to the Boltzmann equation [31], where the term [χ(k) + χ(p) − χ(k ) − χ(p )] needs to be squared to compute any transport coefficients from the Boltzmann equation. This introduces additional power of λ from the soft momentum and softens the contribution of rung a) so that its final size is O(λ 2 ). An analogous mechanism applies to the diagram c) where the two vertical soft lines cause λ −2 enhancement, the explicit vertices and the phase space introduce λ 3 and the Boltzmann equation structure (2.30) -the factor λ, which altogether give the size O(λ 2 ). We also need to evaluate the interference terms, that is, the rungs d) -i). They are all of O(λ 2 ) order when the off-shell exchange momentum is soft. To see this we first consider the rungs f) and g) (for the notation of the rung f) see Fig. 8). They both contain one propagator with a soft momentum l whose contribution is O(λ −1 ), but this is canceled by an additional phase space suppression. Due to combination of two delta functions in the propagators G rr (k + l) and G rr (p + l) with the phase space d 4 l, the latter one is reduced to d 2 l which leads to d 2 l ∼ O(λ). When assessing the size of the rungs h) and i) the same arguments hold as before. The rung h) is shown in Fig. 8) and the corresponding expression is: The soft propagator G ra (l) introduces O(λ −1 ) and the number of integrals over the loop momentum is reduced as previously so that we are left with d 2 l ∼ O(λ). These two factors cancel each other leaving the rungs O(λ 2 ). In rungs d) and e) the vertical line represents the soft propagator which is λ −1 . Including further the phase space suppression and the couplings from the explicit interaction vertices, one gets this rung of the expected O(λ 2 ) size.

The case with collinear singularities
Number changing processes contribute at the same order as 2 ↔ 2 processes (up to logarithm). They are entangled in the same topological structures as number conserving processes, shown in Fig. 6, but emerge under different kinematic conditions. The mechanism responsible for their occurrence is also more complicated than the one discussed above and it is fully controlled by soft physics. Here we briefly and qualitatively discuss how they emerge and evaluate their sizes.
Collinear processes occur when one hard particle splits into two hard particles with an accompaniment of a soft gluon exchange with the thermal medium [70,71,72]. The topological structures corresponding to these processes can be obtained in the procedure shown in Fig. 9. As presented, the rungs representing collinear processes are reproduced by opening one outer line of the two-loop self-energy. The line which is open is denoted by the black cross in the figure. The internal (shaded) lines of the self-energy represent propagators with soft momenta. They contain the hard thermal loop corrections, which is not shown explicitly in the two first columns of Fig. 9. Thus, whenever the cut is through the soft line it means that the hard thermal loop is cut. Consequently, all the cut lines and the external lines are hard and nearly on-shell. Specifically, in contrast to the number conserving processes, the thermal masses in the respective propagators must be included. To evaluate the size of the processes in Fig. 9 we consider in detail the rung shown in Fig. 10, reproduced with a and r positions and momentum convention. As before there is more than one layout of the a and r assignment and a complete analysis of the kernel of the spectral function has to include all possibilities. The size of this rung can be evaluated similarly to the case where collinear singularities are absent, but the power counting is more subtle. First, it is important to point out that whenever a soft line appears in the rung, it must be G rr propagator since it carries the distribution function to account for the interaction with the medium. G rr propagator, by contrast to other propagators, introduces 1/ √ λ enhancement in the soft momentum region. Moreover, the process under consideration is in the collinear regime when there is a pair of the adjacent retarded and advanced propagators with respect to a given momentum. If these propagators were bare their product would produce a singular behavior as their poles would nearly pinch the real axis in the contour integration. This is, however, cured by the inclusion of the self-energies, which leads to a finite expression. As in case of pinching pole approximation, diagrams containing the products G ra G ra or G ar G ar instead of G ra G ar for the same momentum give much smaller contribution to the whole expression and can be neglected in the leading order analysis. The expression corresponding to the rung shown in Fig. 10 is: where (. . . ) means the contribution from the external propagators, which is not needed to be shown explicitly. As mentioned, the external momenta are hard and nearly on-shell, k ∼ p ∼ T and k 2 ∼ p 2 ∼ O(λT 2 ), while the loop momentum is soft l ∼ √ λT . In this kinematic region the integral over the loop momentum is dominated by dl 0 ∼ O(λT ) in the frequency region, and d 3 l ∼ O(λ 3/2 T 3 ). What is more, G ra (l + k) and G ar (l + k) propagators are both O(λ −1 ) since they are dressed with the self-energies to cure pinch singularities. Additionally, since G rr (l + k − p) is in the collinear regime with G ar (l + k), it is also dressed and is of the order O(λ). The properties of propagators impose that (l + k) 2 and (l + k − p) 2 are O(λT 2 ) and the same holds for (l) 2 , which is soft and dressed with the HTL correction. These conditions are, in turn, equivalent to the fact that the angles between all participating particles are parametrically small so that they all propagate collinearly. The small angles are therefore θ kl ∼ θ pl ∼ O( √ λ). The constraints on the angles impose constraints on the phase spaces, which is, d 3 p ∼ |p| 2 d|p|sinθ pl dθ pl dφ ∼ O(λT 3 ) and d 3 k ∼ |k| 2 d|k|sinθ kl dθ kl dφ ∼ O(λT 3 ). The loop momentum l must be spacelike and since it is soft there is an additional Bose-Einstein enhancement making G rr be of the order O(λ −3/2 ). Combining all these powers of λ and the couplings coming from the explicit interaction vertices combined with the closed color loops one finds this rung to be O(λ 2 ).
The presence of the self-energy in G ar (l + k) and G ra (l + k) propagators signals further interactions, which have not yet been explicitly shown nor discussed. In fact one can attach other lines to the side rails of the rung to reproduce processes involving a larger number of participating excitations. For example, one could add a hard line so that to obtain a double gluon emission. Such a process is however subleading [49]. One could also add many soft lines to reflect the process of a hard excitation interacting many times with the medium via a soft exchange and then ending up with splitting into two hard particles. Attaching any number of soft lines to the side rails is possible and they all contribute O(1) corrections. These processes, however, do not need to be explicitly included inside the diagram in Fig 10 since they are resummed within the integral equation for the bulk viscosity. Also, apart from the pair of propagators with pinching poles, there is also a pair of G ar (l + k) and G ra (l + k − p) propagators, which contain nearly pinching poles, where other insertions are possible. We do not examine them here since they are a part of the forthcoming collinear analysis, which investigates the emergence of an effective vertex in the collinear regime. In Fig. 9 we depicted how one can reproduce the collinear processes when one soft line appears in the two-loop self-energy. The rightmost column in Fig. 9 presents the squares of amplitudes of these processes. For the entire analysis to be completed one also needs the interference terms. The leading order interference terms, which contain number changing processes, are those with only 3-gluon vertices; diagrams with 4-gluon vertices are suppressed. The diagrams in question are shown in Fig.  11. The figure presents the procedure of opening one cut line of the self-energy to reproduce the topological structures representing interference terms between collinear processes. To reproduce number changing processes the cutting line has to go through 3 lines of the 2-loop self-energy including the soft (shaded) line (as shown in Fig.  11). The line which can be open is the cut line which is a part of the side rails. Opening the internal (vertical) line would mean opening all color loops and it would lead to emergence of nonplanar diagrams, which are suppressed in large M limit. One can also cut the self-energy so that the soft line remains uncut. Diagrams obtained in this way could only represent number conserving processes. In Fig. 11 we show only a few typical topologies with respect to the position of the soft line. One can also realize that the same structures, but inverted upside down, are also possible. The latter ones would be the complex conjugates of those shown in Fig.  11. The interference terms are essentially the sums of the rungs shown in Fig. 11 and their complex conjugates. Additionally, all structures shown can have different momentum and a and r assignment and the full computation of the imaginary part of the spectral function requires summation over all possibilities. To show that collinear splittings occur at the same order as 2 → 2 processes we examine one representative rung depicted in r, a basis in Fig. 12. Other rungs with collinear singularities can be considered analogously. The expression corresponding to this rung is where (. . . ) stands for the insertion of propagators corresponding to the incoming and outgoing states and we also included integrals over all momenta since pairs of propagators with respect to all momenta can have nearly pinching poles in the collinear regime. In this particular case there are three such pairs: G rr (l − k + p)G ar (l + p) ∼ G ra (l − k + p)G ar (l + p), G ra (k)G ra (p − k) and G ar (p)G ra (p − k), which have singularities with respect to momentum l, k, and p, respectively. Notice that in case of k momentum integration the propagators which have pinching poles are both denoted as G ra due to the notation and assignment of r and a with respect to p in Fig. 12. However, taking into account that G ra (p − k) = G ar (k − p) one obtains the expected G ra (k)G ar (k−p) responsible for the emergence of the singularity. Due to all these constraints G rr (l), G rr (l−k+p), and G ra (p−k) have to be dressed and therefore they all are of O(λ −1 ). In the leading order analysis l has to spacelike, and thus G rr (l) is O(λ −3/2 ). All these properties of propagators have their equivalence in the kinematic constraints, which reflect collinearity conditions, namely, small scattering . These, in turn, limit the respective phase spaces to . Additionally, the integral over dl 0 is dominated by the narrow frequency width, ∼ O(λ). Collecting all powers of coupling constant one finds this rung to be O(λ 2 ).
When assessing the size of the rung in Fig. 12 one can realize that the enhancement in the rung's size coming from the collinear singularities is always balanced by the suppression coming from the phase space caused by the small scattering angle. Given that, there are more effects that need to be included in the leading order evaluation. One can attach infinitely many soft lines to a given pair of propagators with nearly pinching poles and still get the rung at the same order. This is schematically shown in Fig. 13, where a few exemplary insertions of a soft line are shown (the r, a positions and the momentum convention is the same as in Fig. 12). The insertion of soft lines in the leading order is governed by a few rules. The lines have to be G rr propagators and they cannot cross each other since they must be ordered in time and coherent. Moreover, their insertion must follow the standard a and r assignments so that one has to have an odd number of a in a given vertex. Also, a pair of propagators with the nearly pinching poles must appear, otherwise the rung is suppressed by some powers of λ. If all these rules are kept then attaching G rr soft line to the pair of lines with the nearly pinching poles always introduces λ −3/2 from the very size of the propagator, λ from the two explicit interaction vertices and a closed color loop, a pair of new propagators with the nearly pinching poles with the contribution O(λ −2 ) and a phase space suppression d 4 l ∼ λ 5/2 . Altogether the insertion is O(1) and thus infinitely many soft lines can be added in this way without changing the size of the rung. All possibilities have to be resummed and such a procedure reflects the diagrammatic representation of the LPM effect.
The resummation of all possibilities of adding a soft line to a given rung is most efficiently done by finding an effective vertex. The vertex involves three hard and nearly on-shell lines, where all possible insertions of a soft line are included. One exemplary vertex in (r, a) basis is shown in Fig. 14, where one soft line can be added in 3 possible ways. There are more such combinations since r and a can have a  An inclusion of all possible combinations reflects the need for an integral equation, which needs to be solved to find a form of the effective vertex. The solution should be then inserted in the kernel of the integral equation established by the pinch singularities. This approach is, however, demanding within quantum field theory approach. The only essential point to notice is that insertion of any number of soft lines does not change the size of the vertex nor, consequently, the size of a rung where the effective vertex appears.

Integral equations
The bulk viscosity is controlled by the elementary scattering processes entangled in rungs discussed above and both 2 → 2 and effective 1 → 2 processes between gluons contribute at the same order in the 't Hooft coupling λ. For a quantitative computation of the bulk viscosity coefficient all diagrams representing scattering events have  to be resummed, which leads to the relevant integral equations. For the prescription given by the Kubo formula the integral equation is shown schematically in Fig. 15. The kernel of the equation is presented in Fig. 16. It includes 2 → 2 processes and effective 1 → 2 processes as well. For number changing processes another integral equation needs to be solved. It is the equation for the effective vertex and is shown schematically in Fig. 17. The shaded regions denote resummed parts. For this schematic representation of the integral equations we do not distinguish between the propagators with hard and soft, HTL resummed, momenta, but it can be easily done taking into account the discussion in Sec. 2.5. In general, the leading order analysis requires all propagators to be dressed with the self-energies. In principle all bare vertices could be replaced by the effective ones, but since their contribution in all diagrams apart from the last one in Fig. 16 is subleading we do not show them explicitly. In the last diagram in Fig. 16 effective vertices must be used as arbitrary many interactions with the medium through the soft momentum exchange occur at the same order. The rung is responsible for the interference terms and the coherent resummation of all contributions reflects the LPM effect for the effective 1 → 2 processes.
The analytical computation of the bulk viscosity spectral function in terms of quantum field theory tools is very challenging and only qualitative picture can be sketched. The same physics is, however, embodied in the Boltzmann equation as long as the same elementary processes govern its collision kernel. As has been examined in this section, both 2 → 2 and 1 → 2 processes can be reproduced from the planar diagrams of the spectral function. Both classes of processes occur at the same order and so contribute to the kernel of the Boltzmann equation, discussed in detail in Refs. [24,48,49]. It therefore justifies that the collision kernel of the Boltzmann equation captures the same physics as the kernel of the spectral function, shown in Fig. 16, and serves as a convenient way to compute transport coefficients. In particular, the analysis justifies the employment of the Boltzmann equation to calculate the bulk viscosity coefficient of the SU(M ) theory, as carried out in Ref. [29] and summarized in Sec. 2.2 of this manuscript.

Bulk viscosity at intermediate coupling
In the previous section, we have discussed the behavior of the bulk viscosity in the SU(M ) gauge theory in the weak coupling limit. In the next sections, we will discuss the strong coupling behavior. In these two limits, we have well defined calculational tools, perturbation theory in the case of the weakly coupled limit and the AdS/QCD correspondence in the strongly coupled limit. When the coupling is neither weak nor strong the only reliable QCD results are from Euclidean Lattice QCD (LQCD) calculations. Unfortunately, direct extraction of viscosities from LQCD is very nontrivial since viscosities have to do with dissipation in real time while LQCD calculations are inherently static. Of course, if one can calculate full Euclidean correlation functions, they can be analytically continued to real time correlation functions. But since only discrete and finite number of data points are available from LQCD, this procedure in practice introduces large uncertainties.
In literature, efforts were made to extract information on the bulk viscosity from LQCD results using sum rules [35,36,37,52,53,54,55]. In this section, we summarize the main points and point out why it is difficult to get any information on the bulk viscosity from the sum rules in particular and from LQCD results in general.
Extraction of the bulk viscosity from the sum rule relies on the Kubo formula for the bulk viscosity ζ: where ρ P P is the spectral density for the pressure-pressure correlator. Hence one may expect that the bulk viscosity can be extracted from a sum rule involving ρ P P (ω, 0)/ω. Equivalently, it may be able to extract ζ from a sum rule involving the correlation function of the trace operatorΘ =T µ µ =T 00 − 3P because the energy-momentum conservation laws dictate that the zero wavenumber limit ofT 00 correlators vanish. Note that different forms of the trace operator can be used and they were briefly discussed in the previous section.
In Ref. [36], low energy sum rules for the stress-energy tensor trace, Θ = T µ µ , were derived and the following result was established: where Θ G is the gluon contribution to the stress-energy tensor trace and ρ ΘΘ is the spectral density for the ΘΘ correlation function. 6 As we have argued in the previous section, the quark contribution is negligible in the large M limit and we will not consider it here, either. Consequently we will drop the subscript G. The trace average is: where Θ 0 is the vacuum contribution. In Ref. [52], re-derivation of the results with the direct subtraction of the vacuum spectral density led instead to: is the deviation from the vacuum spectral density at finite temperature T . The difference between the two sum rules was attributed to the non-commutability of the limits lim ω→0 and lim k→0 , see [52]. In Ref. [54], the sum rule Eq. (3.4) is re-cast as: where ρ * is the spectral density for the operatorΘ * =T µ µ −(1−3c 2 s )T 00 . The spectral density of the operatorΘ * satisfies the same Kubo formula but has an added benefit that the limits lim ω→0 and lim k→0 commute.
The right hand side of the sum rule (3.5) can be evaluated using the LQCD results. If one can then show that the left hand side is a well defined function of the bulk viscosity then these sum rules may be used to determine the bulk viscosity in the region of temperature where LQCD calculations can be performed.
A first attempt at relating the sum rule integral (the left hand side of Eq.(3.5) to the bulk viscosity was carried out in [35,36]. In Ref. [36] the following ansatz was introduced: which does satisfy the Kubo formula and makes the sum rule integral in the left hand side of Eq. (3.5) proportional to ζω 0 . However, this form lacks contribution from frequencies higher than the unknown parameter ω 0 , see Ref. [80]. It turned out that the high frequency contribution is actually negative that largely cancels the low frequency contribution. The fact that this ansatz is not adequate has been shown by [37,52,54] both perturbatively and non-perturbatively. The biggest problem is that the right hand side of Eq. (3.5) is negative while Eq.(3.6) makes the left hand side strictly positive. This difference can be attributed to the the presence of the glueball [54]. Hence, the sum rule Eq. (3.5) is not particularly useful since it cannot be definitely established that the dominant contribution to the sum rule integral comes from the low frequencies.
If one cannot rely on the sum rule, then one needs to obtain the spectral density directly from LQCD calculations at least in the k = 0 and |ω| T limit. This is not an easy task as it involves analytic continuation when only a finite number of data points in the Euclidean space is known. First attempts in this direction were made in [53,54], which does show that δρ * (ω, 0)/ω has a peak at ω = 0. Unfortunately, actual values and behavior of ζ/s obtained in this way contains too much uncertainty at this point. All one can conclude from the LQCD studies right now is that

Bulk viscosity at weak string and strong 't Hooft couplings with zero flavors
In the previous section we studied the bulk viscosity at weak 't Hooft coupling, and argued how the ratio of the bulk to shear viscosities should be interpreted at both weak and strong couplings. As mentioned therein, the strong coupling result depends on the existence of a gravity dual of the resulting framework. Our aim here is to analyze SU(M ) gauge theory at various values of the 't Hooft couplings and at high temperatures as depicted in Fig. 18. The three regimes of interest are shown in Fig. 18: the yellow box denotes weak, the green box denotes intermediate and the blue box denotes strong 't Hooft couplings, all at high temperatures. The theories governing each of the three dynamics are also varied as we discussed above. The weak and the intermediate couplings are studied using kinetic and LQCD, whereas the strong coupling will be studied using string theory. The latter however is more elaborate because of its UV properties, and in fact differs quite a bit from what we expect from kinetic theory and LQCD. Let us start with the kinetic theory which was discussed in details earlier. The regular RG flow of such a theory is governed by the black lines in Fig. 18. At low energies the YM coupling becomes very large and the theory confines. However if Figure 18: The three different regimes of interest used here to study bulk viscosity. The yellow box denotes the regime of kinetic theory, the green box denotes the regime of LQCD, and the blue box denotes the regime of string theory. All these regimes are analyzed at high temperatures, i.e above the deconfinement temperatures, and the regular RG flows connecting the yellow and the green boxes are denoted by black curves. On the other hand, the cascading RG flows, that specifically arise from string theory, are shown here in the blue box. All the three different RG flows lead to a consistent picture at low energies connecting the weak, intermediate and strong 't Hooft couplings.
we increase the temperature, the coupling can be made smaller. In fact at high temperature the 't Hooft coupling λ ≡ g 2 YM M can become very small even for large M . This is of course the regime where kinetic theory can be studied (see section 2 for more details), and is denoted by the yellow box in Fig. 18. One can similarly go to the intermediate coupling regime, whose dynamics is governed by LQCD..
What we now require is to understand the regime where the 't Hooft coupling λ can be very large for both weak and strong YM couplings. This is the regime where neither kinetic theory nor LQCD can help us, and therefore the only way we can have any analytic control is to use techniques of string theory. Of course when M , the number of colors, is small even string theory cannot provide a controlled laboratory, so it is the large M limit that can be tackled using stringy techniques. This is then the regime of gauge/gravity dualities, i.e the dynamics at strong 't Hooft in the gauge theory side may now be done using a gravity dual description.
Clearly since string theory provides a UV complete picture, it is natural to ask what UV completion would mean in the present set-up. However before we go about exploring this side of the story, let us first elucidate the IR dynamics of the theory directly from the gravity dual description. The gravity dual description was originally provided in [58] (see also [63] for the mirror set-up which will be useful soon). In simple terms, the gravity dual is given in terms of a resolved warpeddeformed conifold with fluxes with an additional black hole that provides the high temperature physics in the gauge theory side, i.e the physics above the deconfinement temperature. In the absence of a black hole we expect minimal four-dimensional supersymmetry (that may be broken too), whose simplest description appear, on one side from wrapped D5-branes on a non-Kähler resolved conifold [62], and on the other side from fluxes on a resolved warped-deformed conifold alluded to above [61]. The "resolution" parameter in the resolved warped-deformed conifold is responsible for the UV completion, that we will discuss soon (see also [82] for a slightly different realization of the same story). In the following we want to discuss the background as well as the issue of supersymmetry, mostly for the IR part of the gauge theory. For simplicity we will concentrate on the Baryonic branch of the gauge theory where is the issue of supersymmetry is most prominently displayed. Later on, in section 4.1, we will concentrate on a more specific point in the moduli space of the corresponding gauge theory.
In the Baryonic branch, generated by M wrapped D5-branes on a non-Kähler resolved conifold [62], the gravity dual for the IR physics may be given by the following type IIB background with three-and five-form fluxes [61] 7 : where (θ i , φ i , ψ) are the angular coordinates, β is the parameter of the Baryonic branch, h is the warp-factor and J is the fundamental (1, 1) form that is not closed. We have also denoted the dilaton by φ, and the five-form by C 5 . The internal metric ds 2 6 can be expressed as: with H i (r) being the additional warp-factors. Note that the two-spheres, denoted by (θ 1 , φ 1 ) and (θ 2 , φ 2 ) have different curvatures governed by H 3 and H 4 respectively, and their inequality will be responsible for UV completion. The complexified threeform flux G 3 then takes the following form [61]: where M i are certain functions expressed in terms of the vielbeins whose form may be ascertained from eq (2.113) of [61]. The E i are defined with the following choice of the almost complex structure: which is integrable 8 for a constant dilaton, otherwise the three-form flux G 3 is defined as an ISD (Imaginary Self-Dual) form with respect to the almost complex structure (4.4). Note that (4.3) is a (2, 1) form as one would expect from a supersymmetrypreserving background. Additionally, the choice of the Baryonic branch tells us that the gauge group is SU(2M ) × SU (M ), which is in fact one cascading step away from the confining SU(M ) gauge group that we seek! In the blue box of Fig. 18 this may be seen as the second-last stage of the cascading RG flow before permanent confinement sets in.
One can also give a physical meaning to the Baryonic branch directly from the wrapped five-brane picture. The SU (2M ) × SU (M ) gauge group implies that, along with M wrapped D5-branes, we have M D3-branes too. The five-branes wrap the two-sphere parametrized by (θ 2 , φ 2 ). For vanishing size of the two-sphere, the M additional D3-branes preserve the same supersymmetries as the M wrapped D5branes. However if the two-sphere is of finite size, supersymmetry is completely broken, and the only way to preserve supersymmetry in this case would be to dissolve the D3-branes on the D5-branes.
Being on the Baryonic branch does not give a well-defined UV picture. We will still need to find the UV completion of the model. This will be discussed a bit later, but note that being in the Baryonic branch does tell us that if we make one Seiberg duality we will land in the confining SU(M ) gauge theory description. At non-zero temperatures, we will require black holes in the gravity side of our story. Since this is the premise on which our calculations in this paper will be based on, let us elaborate the story a bit more. At zero temperature, the duality sequence that we shall use is laid out in Fig. 19. On the bottom left corner, i.e box (a), is the gauge theory configuration discussed in [58,62] with M D5-branes wrapped on the two-sphere parametrized by (θ 2 , φ 2 ). This is a non-Kähler resolved conifold because at r = 0 there is a resolved two-sphere parametrized by (θ 1 , φ 1 ). The usefulness of such a configuration will be spelled out a little later. The wrapped D5-branes on the non-Kähler resolved conifold give rise to the gravity dual background which is a non-Kähler deformed conifold with three-form fluxes, much in the lines of (4.1), (4.2) and (4.3) and is given by box (b) in Fig. 19. The computations performed in section 4 will be based on this configuration, albeit with a black hole that will signify non-zero temperature, but with no flavors. Figure 19: The two configurations, one in type IIB and the other in the M-theory uplift of type IIA, on which all the computations of sections 4 and 5 respectively will be based upon. On the left is the type IIB picture with the gravity dual given by a resolved warpeddeformed conifold with fluxes. On the right is the M-theory uplift of the type IIA gravity dual. The IIA gravity dual involes a non-Kähler resolved conifold with fluxes, whereas the M-theory uplift is a seven dimensional manifold with a G 2 structure. The type IIB computations will be done at high temperatures, i.e above the deconfinement temperatures, but with zero flavors. The type IIA, and also the M-theory uplift, will take into account both high temperatures as well as non-zero flavors.
A mirror transformation, a la Strominger-Yau-Zaslow [64], on both the type IIB boxes of Fig. 19 will produce the IR type IIA background whose gravity dual configuration involves a non-Kähler resolved conifold with fluxes, as shown in box (d) in the figure. The M-theory uplift of this is given in the top right-hand box of Fig.  19, i.e box (e), which is a seven-dimensional G 2 structure manifold with G-fluxes [63]. Our computations in section 5 will be based on this specific M-theory manifold albeit, again, with non-zero temperatures but now including non-zero flavors. Interestingly for the spectral function computation of section 6, we shall resort back to the type IIA picture.
Let us now come to the UV completion of these models that we alluded to earlier. In the type IIB side, this was first discussed in [58], but a full elaborations on the actual ingredients that constitute the UV degrees of freedom were given in [59] and [60]. We expect the UV theory to be a strongly coupled conformal field theory, as this would be the closest to being asymptotically free. The reason for choosing a CFT − and not an asymptotically free theory − as the UV theory is because we require strong 't Hooft coupling to allow for a gravity dual. In fact, a gravity dual description only exists if the corresponding gauge theory is strongly coupled at all scales, i.e strongly coupled from UV to IR. For large but finite number of colors, this means that the requirement for asymptotic freedom is not quite compatible with the existence of a gravity dual. Therefore the closest we can come to asymptotic freedom is to allow for a CFT in the UV. In the limit of infinite number of colors, the 't Hooft coupling can be very large, yet the YM coupling can be made arbitrarily small.
One specific choice of a UV group that could lead to a CFT is SU(N + M ) × SU (N + M ), where we have introduced an extra parameter of N . In the present context, the choice of N has a special meaning. In the type IIB theory, N signifies the number of D3-branes whereas M is the usual wrapped D5-branes. The two-cycle on which the D5-branes wrap, i.e the two-cycle parametrized by (θ 2 , φ 2 ), should now be of vanishing size to preserve supersymmetry. In the blue box of Fig. 18, we have denoted the UV group SU(N + M ) × SU(N + M ) that is shown to get Higgsed to a smaller group SU(N + M ) × SU(N ) at a certain IR scale. This is followed by a series of cascading RG flows that eventually takes us to the confining gauge group SU(M ) at the far IR.
The complete RG flow that is depicted in the blue box of Fig. 18 can be described rather succinctly from both type IIB as well as the type IIA theories. This will also answer all the questions that we put aside earlier. From the type IIB side, the UV CFT may be easily described by allowing additional M anti-D5-branes distributed on the northern hemisphere of the resolved sphere parametrized by (θ 1 , φ 1 ). These anti-D5-branes are stabilized against collapse by using fluxes, details of which have appeared in [59,84]. The string connecting the branes and the anti-branes are heavy, and they are integrated out at low energies. Thus at low energies we only see the cascading SU(N ) × SU(N + M ) theory. At high energies, the anti-brane degrees of freedom are integrated in and the M D5-brane and the M D5-branes combine to give M D3-brane degrees of freedom. Together with N D3-branes localized at the south pole of the resolved sphere, this leads to the UV CFT described above. Therefore the three stages of operation namely, (1) emergence of CFT, (2) Higgsing and (3) the cascading behavior, are all described neatly from the type IIB configuration of N D3, M D5 and and M D5-branes on a non-Kähler resolved conifold with fluxes.
The correctness of our construction may also be ascertained from a T-dual type IIA configuration as shown in Fig. 20. This T-duality is a single T-duality along the ψ direction and therefore should not be confused with the three T-dualities that we performed earlier to determine the mirror configuration. A single T-duality of a conifold along ψ direction, in the type IIB theory, leads to a configuration of two orthogonal NS5-branes in the type IIA side. In the presence of N + M extra D3branes in the type IIB side, the T-dual configuration is shown on the left of Fig. 20. The M D3-branes have five-brane origins as discussed above, and so the configuration on the left of Fig. 20 gives us a CFT with a gauge group SU(N + M ) × SU(N + M ). The reason why this is a CFT comes from the fact that the NS5-branes are not bent. Clearly, any bendings of the NS5-branes would have lead to running couplings of the gauge theories on the D4-branes [97]. These bendings can be achieved by having an unequal number of D4-branes on both sides of the NS5-branes. Such a feature may be achieved independently but does not seem to come naturally from the configuration on the left of Fig. 20: a consequence due to the absence of Coulomb branches in N = 1 gauge theories.
However, all is not lost as these theories do have other branches, namely Baryonic, Mesonic and possible remnants of the N = 2 Coulomb branches. Without going into too much details, which the readers may find in [61,82], one may easily see that a branch in the moduli space arises by putting an extra NS5-brane along the dotted line in the left configuration of Fig. 20. Happily, this does not break any extra supersymmetries but creates the necessary Higgsing effect that we require to jump-start the cascading process! On the right of Fig. 20 we have shown how one may go from UV conformal to IR cascading behavior. As should be obvious from Fig.  20, moving the M D4-branes along parallel NS5-branes bends the NS5-branes, thus creating RG flows on the remaining D4-branes. The far IR physics is then exactly a confining SU(M ) gauge theory with decoupled U(1)'s that we seek here. Switching on a non-zero temperature we can study the various transport coefficients.
In the gravity dual, the IR story is clear: this is given as in (4.1) and (4.2).
The UV degrees of freedom start appearing from Region 2 onwards as shown in [59], and as we go to large r we are effectively in Region 3 where the three-form fluxes vanish and the background asymptotes to an AdS 5 space. In this section we will use a slightly simplified form of this background and mainly concentrate on Region 1 − to be at low energies − to study the bulk viscosity at strong 't Hooft but weak string coupling in the absence of fundamental flavors. In the next section, we will put in the flavors and study the bulk viscosity as well as the ratio of the bulk to shear viscosites at both strong 't Hooft and strong string couplings, again concentrating on Region 1. For an earlier work on bulk viscosity with bottom-up approach, using two different AdS spaces at UV and IR and for a wide class of models, see [83]. Note however that the study of bulk viscosity in [83] differs from our study here in at least two respects. First, the model considered in [83] has two fixed points: one at UV and the other at IR respectively. This differs from the IR confining model that we consider here. Secondly, the study of bulk viscosity in [83] finds violation of the Buchel bound [30]. Although this is possible in our set-up, by choosing a different lower bound for d 1 in (4.85) and (4.86), we do not analyze such cases here.

The type IIB dual background for large N thermal QCD
In [84] we made some preliminary study of bulk viscosity using the UV complete large N thermal QCD model of [58] with N f = 0. The metric that we took in [84] is of the form: Note that the internal space is a warped resolved conifold and not a resolved warpeddeformed conifold as one would have expected from (4.2). This is a simplifying assumption which helped us to study bulk viscosity without worrying about the far IR regime of the gauge theory. Recall that the far IR regime of the gauge theory is governed by the blown-up three-cycle of the resolved warped-deformed conifold.
However since the small r regime of the geometry is covered by the horizon radius r h , our choice of the metric (4.5) is not too far from the correct answer. The resolution parameter a 2 (r) is not the resolution parameter used in the brane side to control the UV behavior of the theory. In the brane side, i.e in the gauge theory description, the M D5-branes wrap the vanishing two-cycle of the resolved conifold parametrized by (θ 2 , φ 2 ) in a way that the D5-branes (and the N D3-branes) are at the south pole of the resolved 2-cycle, parametrized by (θ 1 , φ 1 ), and the anti-D5-branes are distributed over the upper hemisphere of the 2-cycle.
In the language of the metric (4.5), this means putting the M D5-branes on the (θ 2 , φ 2 ) 2-cycle has caused an asymmetry quantified by the resolution parameter a 2 . From the discussion above this would mean that a(r) 2 = O( ) and should have no terms that are zeroth order in . This can be confirmed by plugging the metric into the equations of motion.The Einstein's equations are: here G 3 , F 5 and τ are the complexified three-form flux, five-form flux and the axiodilaton respectively as defined in (4.1) and (4.3). To figure out how the wrapped D5-branes, inserted in the non-extremal system, affect the warp-factor , we can express the change as: where A 0 ≡ − 1 4 log L 4 r 4 is the conformal value and = 3gsM 2 2πN is our expansion parameter. The resolution parameter a 2 now may be expressed in the following way: where, as we emphasized above, to zeroth order in , the D5-branes wrap vanishing two-cycle. We start seeing non-zero resolution only from the first order in . The two functions, P (r) and Q(r), are related via the following set of equations: where D 1 , D 2 and D 1 are constants that may be fixed from the boundary conditions. This has been discussed in details in [84], and after the dust settles, the functional form for P (r) and Q(r) can be explicitly represented in the following way: where both behave well in the limit r → r h as one would have expected. In fact knowing the functional form for Q(r) immediately tells us what the black-hole factor, e 2B , in the metric (4.5) should be. This may be expressed by the following integral form: which reproduces the conformal result for vanishing . Finally, plugging (4.10) to (4.7) and (4.8) gives us the O( ) corrections to the conformal values for the resolution and the warp-factors: The functional forms for a 2 , e 2B and e −4A are consistent with the general picture developed in [58], [59] and [85]. In particular knowing the O( ) correction to the black-hole factor is consistent with the O( ) corrections to the two black-hole factors g 1 and g 2 in [58]. There are also O(g s N f ) corrections, from the N f flavors, that we do not consider here. This is relegated to section 5.

Details on the bulk viscosity computations from gravity dual
Bulk viscosity appears, in a system with a SO(3) spatial symmetry, from the correlation of T xx at two different points in four-dimensional space-time with one point fixed at the origin. This means, as discussed earlier, in the gravity dual bulk viscosity may be computed from the fluctuations of the vielbeins e k with k = 0, x and r. These fluctuations may be divided into positive and negative frequencies, and are expressed as: where is the same non-conformality factor as before and ω is the frequency. The other parameters appearing in (4.13) may be defined in the following way. The coefficients p nk are in general functions of r as well as |ω| but not constants. With constant p nk , the bulk viscosity would vanish despite the existence of a complex piece in (4.13). Note however that δe k ≡ δe k (r, t) are all real functions of r and t.
The coefficients Γ 0k (r, |ω|) and Γ 1k (r, |ω|) capture the essence of the bulk viscosity computations here. In a system with SO(3) symmetry, Γ 0x takes the following form: where B(r, 0) is given in (4.11) and T , the temperature, is proportional to r h , the horizon radius. We also expect Γ 0y = Γ 0z to be equal to Γ 0x . On the other hand, Γ 00 and Γ 0r take the following form: Although (4.14) and (4.15) are related to conformal theory 9 , we will use them to analyze the non-conformal regime of our model as the imaginary parts of the fluctuations in (4.13) depend on Γ 0k as well as p nk . The latter are associated with extra sources coming from the distribution of anti-D5 branes in Regions 2 and 3. We can quantify these sources in the following way: where we see that the imaginary part involves three infinite series of modes specified by the sources ∆ where A 0 is defined in (4.7) and B 0 ≡ B(r, 0) is given in (4.11). Note that (4.17) involves five fluctuation modes, p n0 , p nx , p nr , p (n−1)x and p (n−1)r ; as well as the three Γ 0k 's defined in (4.14) and (4.15). This means knowing ∆ (n) 2k , we will need at least five equations to solve for the fluctuations p nk . One may also construct the following recursion relations from (4.17): and so on. Note that there are two types of non-derivative terms in the first equation of (4.18): (a) the terms proportional to p 00 , p 0x and p 0r , and (b) terms proportional to p (−1)x and p (−1)r . The latter have no dynamics, so maybe we could use them to cancel the former terms in the following way: This would mean that by knowing p 00 , p 0x and p 0r , one may not only build the next set of fluctuation modes from (4.18) but also determine the functional forms for the non-dynamical modes p (−1)k . Unfortunately such an identification would either overconstrain the dynamics or lead to some apparent contradictions. To avoid this, we will set: In any case, identification like (4.19) can never be used to cancel the non-derivative terms p 1k with p 0k as both set of fluctuations are dynamical now. Thus generically we should assume the existence of p (n−1)k modes along with the p nk modes. The next series of sources appear from ∆ (n) 2x and would follow similar strategy as above. These sources may be expressed in terms of the fluctuation modes p nk in the following way: where this time four, instead of five, modes p nx , p n0 , p nr and p (n−1)x are needed. As before, the zeroth and the first order recursion relations may be written as: where this time an input of p 0x is needed to build the first order fluctuation equation.
In a similar vein, we now construct the third series of sources associated with ∆ (n) 2r in the following way: with the input of four fluctuation modes p nx , p n0 , p nr and p (n−1)r governing the dynamics. The recursion relations for the zeroth and the first order can be easily expressed in terms of the fluctuation modes as: At this stage let us ask whether the above three set of equations, (4.17), (4.21) and (4.23), are enough to determine the five unknown functions 10 , p n0 , p nx , p nr , p (n−1)x and p (n−1)r . It would seem we need at least two more equations. However a careful look tells us that the first equations in each of the three recursion series, (4.18), (4.22) and (4.24), are enough to determine the three functions p 00 , p 0x and p 0r provided the sources ∆ 2k and the boundary conditions are adequately specified. Similar arguments apply for the next three functions, p 10 , p 1x and p 1r : once we specify the sources ∆ (1) 2k and the boundary conditions, this would in principle fix the functional forms for all p 1k . Thus it seems that the above three equations (4.17), (4.21) and (4.23) should suffice.
For the present case we will work out the equation satisfied by p 0x as this is the only component relevant for bulk viscosity. This will be explained soon (see also [84]). In fact what is required is not p 0x , rather p 0x , and therefore we will work out the equation for Y x (r, |ω|) ≡ p 0x . In the process we will also see how to write the equations for p 00 and p 0r . To start, let us define a few variables f i , g i and h i using which the zeroth order equations in (4.18), (4.22) and (4.24) may be re-expressed in the following way: where we have identified p 01 ≡ p 0x and p 02 ≡ p 0r to avoid clutter. One may define similar equations for the first order fluctuation equations, namely for the f 1k using the recursion relations. The various coefficients appearing in (4.25) may be written as: where Γ 0k have been defined in (4.14) and (4.15); and A 0 and B 0 are defined in (4.7) and (4.11) respectively. It is interesting to note that the LHS of the equations in (4.25), i.e the coefficients of p 0k and p 0k , are mostly functions of Γ 0k whereas the RHS of (4.25), i.e the coefficients of p 0k , are all functions of the derivatives of Γ 0k . The set of equations (4.25) are highly non-linear and solving them will in general be a non-trivial exercise. Therefore it might be instructive to first solve a slightly simpler system than (4.25) to gain some familiarity with the solutions and then proceed to address the full set of equations. In the following subsection we analyze a simpler case, and in the next subsection we will study the full system.

A toy example in full details
To study a toy example from (4.25), the first question is: how can we simplify the set of equations in (4.25)? This is where the observation that we made above could become useful, namely, we can assume that the derivatives of Γ 0k are much smaller than Γ 0k at some r >> r h . This would immediately imply: making the RHS of all the equations in (4.25) to only depend on the sources ∆ 2k . Note that (4.27) does not imply absorbing the p 0k terms in the definition of the sources because the sources are independent of the bulk fluctuations. Nor does this imply invoking relations like (4.19), since such a procedure is generically prone to errors. Thus (4.27) would be the only way to simplify (4.25).
With this in mind the next set of procedures may be elaborated in the following way. Using (4.26), let us define another set of functions as: which will help us to avoid cluttering of formulae later when we write the equations for the fluctuations p nk . Note that these functions are all expressed in terms of certain definite integrals (the lower bounds of these integrals could be r h or r = 0, but these details will be irrelevant). There are also four other functions that are not expressed in terms of integrals. They may be expressed as: where ∆ 2k are the zeroth order sources that appear in (4.25). Note that G 2 and G 4 are functions of r as well as |ω| because they depend on the sources ∆ (0) 2k (r, |ω|). Therefore with (4.26), (4.28) and (4.29), we are ready to write the equation governing the fluctuation Y x (r, |ω|) ≡ p 0x as: which is a second order differential equation and therefore would require boundary conditions, both at the cut-off r = r c as well as at the horizon radius r = r h , to determine the functional behavior precisely. The coefficients a I1 appearing in (4.30) are non-trivial functions of F i and G i variables, defined in (4.28) and (4.29), and can be written as: where a 41 = a 41 (r, |ω|) and all other a I1 are functions of r. This implies Y x = Y x (r, |ω|) as expected. Solving (4.30) would provide the fluctuation mode p 0x . Once p 0x is known, we can use it to determine the next fluctuation mode, p 00 . Let us now define p 00 ≡ Y 0 (r, |ω|), instead of p 00 , and write the equation for Y 0 in the following way: where one may use either of the two set of expressions on the RHS of (4.32) to solve for Y 0 . The equality between the two expressions can be argued easily from (4.25). Finally, knowing Y x and Y 0 , one may use any of the three equations in (4.25) to solve for Y r (r, |ω|) ≡ p 0r .
Let us now work out the first order fluctuations for our case invoking (4.27). Again we expect three set of fluctuations of the form p 10 , p 1x and p 1r , similar to the three set of fluctuations p 00 , p 0x and p 0r respectively for the zeroth order case. The equations satisfied by the first order fluctuations are a slight variations of (4.25), namely: where f i , g i and h i are exactly the ones appearing in (4.26); Γ 0k are as in (4.14) and (4.15); and A 0 and B 0 are the zeroth order values in (4.7) and (4.11) respectively. However not everything remain the same: the RHS of the equations (4.33) have two kind of sources, (a) the first order sources ∆ (1) 2k , and (b) sources appearing from the zeroth order in fluctuations, p 0x and p 0r . These changes in sources imply that G 4 (r, |ω|) and G 2 (r, |ω|) in (4.29) may be replaced by: 2k in G 2 and G 4 . Note that there are no additional changes to G 1 (r) and G 3 (r) in (4.29). The above observation immediately tells us that the equation satisfied by p 1x ≡ Y 1x (r, |ω|) should be: where we see that the coefficients appearing in the LHS of (4.35) are the same as the ones appearing in (4.30) with a 11 , a 21 and a 31 as given in (4.31). The only difference from (4.30) is the replacement of a 41 by a 41 , where: Similarly the equation for p 10 ≡ Y 10 (r, |ω|) will be similar to (4.32) with the replacement of Y x by Y 1x and G 2 and G 4 by G 2 and G 4 respectively. Once we know Y 1x and Y 10 , we can use (4.33) to determine the equation for Y 1r . This way the first order fluctuations may be completely determined. The picture is now clear for the generic order fluctuations. If we want to study the n-th order fluctuations Y nx , Y n0 and Y nr , all we need is to rewrite the sources, ∆ 20 − e −4(A 0 +B 0 ) r 3Y (n−1)x (y, |ω|)Γ 0x (r) + Y (n−1)r (y, |ω|)Γ 0r (r) dy.
Once these sources are specified we can construct G 4 using ∆ (n) 2r and ∆ (n) 2x ; and G 2 using ∆ (n) 20 and ∆ (n) 2r using the definitions in (4.34); and finally a 4 using (4.36). The equations for Y nx , Y n0 and Y nr would then follow the steps outlined above.

Towards exact solutions for the fluctuations
To study exact solutions for the system of equations in (4.25), one way would be to eliminate the p 0k pieces on the RHS by rearranging the set of equations there. However a slightly simpler approach is to keep the RHS only as a function of p 0x and eliminate the others. This leads to the following set of equations: where, as mentioned above, we kept the RHS as functions of p 0x only. The set of equations (4.38) are in some sense more symmetrical than the earlier set of equations (4.25). The coefficients are expressed in terms of brackets which may be defined as: This formalism has some distinct advantages that will be clear soon. Note also that, in (4.38), there are no second derivatives of p 0r which in turn will help us to rearrange the set of equations further. But before we do so, let us define the coefficients appearing in (4.38). The k i are defined in the following way: where k 6 and k 7 will be used to describe the sources ∆ 1 and ∆ 2 in (4.38) below. All the k i are in turn constructed out of the (f i , g k ) coefficients defined earlier in (4.26).
In a similar vein, the l i coefficients are defined as: where the (h i , g k ) coefficients, used here to define l i , are given in (4.26). As before, the (l 6 , l 7 ) coefficients will be used below to describe the sources ∆ 1 and ∆ 3 . Finally the m i coefficients may be defined in the following way: where again the (h i , f k ) coefficients are given in (4.26), and m 6 and m 7 will be used to describe the sources ∆ 2 and ∆ 3 . The sources ∆ k may now be expressed as: The new sources are combinations of the original sources ∆ 2k , the coefficients defined in (4.40), (4.41), (4.42) and (4.26); and p 0x . These equation explicitly take us away from the simplifying assumption (4.27), and so are only valid when no approximations are made 11 . Additionally, the dependence of all the sources only on p 0x means that any further rearrangements of the sources will not have new dependences on other fluctuation modes. This means one may eliminate p 0r pieces from (4.38) to simplify them further in the following way: which mix the sources ∆ 1 and ∆ 2 as well as ∆ 2 and ∆ 3 . We could also write another equation, parametrized by γ i coefficients, that mix the sources ∆ 1 and ∆ 3 , but that won't be necessary for us. The new sources may be expressed in the following way: which explicitly show that they are not only linear with respect to the fluctuation mode p 0x but also that no other modes show up in the definition (4.45). The precise coefficients of p 0x appearing in the sources above are respectively: which do not vanish generically, although special cases with vanishing coefficients could appear. Of course in the limit (4.27) everything vanishes, but since we are no longer considering the simplifying condition (4.27), we will assume non-zero coefficients. This consideration also allows us to express the other coefficients in (4.44), namely α i and β i , in the following suggestive way: which again do not generically vanish. At this stage the signs of the various α i and β i coefficients are not important, but they could be worked out by carefully studying the relative terms. The relative terms depend on the (k i , l i , m i ) coefficients defined in (4.40), (4.41) and (4.42) respectively which in turn are expressed in terms of coefficients given in (4.26). We also expect α i = β i as well as α i α j = β i β j for all i = j, which may be inferred from (4.47).
Something interesting happens here. Eliminating p 00 from (4.44) lands us directly to an equation for p 0x ≡ Y x whose form is similar to what we had earlier when we analyzed a toy example. This means, as in (4.28) therein, we can define the following functions: using the integrals of the functions defined in (4.47), assuming neither α i nor β j vanish. If any of the α i or β j vanish, the analysis has to be changed completely to get the requisite equation for Y x . We can also define another set of functions using α i , β j and the sources ∆ [a,b] that do not involve integrals, much like the ones in (4.29). They are: where as before we have − similar to G 2 and G 4 in (4.29) − P 1 and P 2 that are functions of both r and |ω| because of their dependences on the sources ∆ [1,2] and ∆ [2,3] respectively. Thus using (4.48) and (4.49), we can write the equation for Y x in the following way: similar to (4.30). The coefficients a I2 are defined in somewhat similar form to (4.31) in the following way: Note that, although the analysis is similar to what we had for (4.30), there is an important difference now. The RHS of the equation (4.50), defined using a 4 is constructed with P 3 and P 4 which are in turn defined in (4.49). Both P 3 and P 4 are linear in p 0x as may be seen from (4.45) and (4.43). Thus a 4 in (4.51) differs from a 4 in (4.31) by the presence of p 0x , implying (4.50) to be a third order equation in p 0x . We can use the above set of equations to formulate the equation for p 00 , instead of p 00 , as we had in (4.32). Needless to say, the equation for Y 0 (r, |ω|) ≡ p 00 follows similar route as before, and we can write the equation for Y 0 in the following way: where the equality between the two sides is the consequence of (4.50). The way we have constructed the sources P 3 and P 4 in (4.49), Y 0 do not appear on the RHS of (4.52) and therefore knowing Y x we would not only know: respectively. These may be worked out with some effort, but we will not do so here as these fluctuation modes are not important for computing the bulk-viscosity to the order that we want to analyze here. The story however does not end here as there are additional constraints on the p nk modes that appear from the flux EOMs, namely the five-form, the three-forms and the axio-dilaton EOMs. We can also get another equation from the cross-term in the metric, namely the rt component of the metric. All these should further constrain the fluctuation modes, and there is a worry that these additional EOMs may over-constrain the system rendering them inconsistent. The scenario is subtle, so let us proceed carefully. First, and to O( ), we may ignore the three-form EOMs as they start changing the equations only to O( 2 ). Similarly, once we switch off the g s N f corrections we are also effectively switching off the contributions from the axiodilaton EOMs. On the other hand we cannot ignore the five-form and the rt EOMs. They will constrain the p nk modes, and it is easy to see how the rt component of the metric EOM does this: where the summation convention for k follows the same as in (4.25). The other coefficients appearing in (4.54) are defined in the following way: c nk are constants that one may determine from the way the sources arrange themselves in the rt EOM, whereas c [q]nkm are functions of r such that: In fact this is where the above mentioned constraint show up: one can determine the functional forms of the coefficients c [q]nkm and the constants c nk by comparing with the LHS of (4.54). One may also get these coefficients directly from the rt EOM. We expect these two ways of getting these coefficients to match because in the absence of the sources i.e for the conformal case, the extra rt equation did not over-constrain the system [84]. Motivated by the above discussions, one may now give similar arguments for the five-form EOM, where the constraint equation takes the following form: for every choice of n, and with d nk and d [q]nkm being the coefficients similar to c nk and c [q]nkm respectively in (4.54) with d [q]nkm vanishing for m ≥ 3 as (4.55). As before, the RHS of (4.56) may be expressed in terms of the modes p nk and their derivatives which may be compared with the LHS of (4.56). The system will be consistent when all the coefficients on both sides match. There is also a simpler way to see why the coefficients on both sides of the equations in (4.54) and (4.56) would match, once the RHS of these equations have been specified in terms of the sources and the modes. This is because, all the three equations in (4.17), (4.21) and (4.23) may be expressed as: with f (k) [q]nlm being constrained in the same way as in (4.55), implying that the RHS of either of the two equations (4.54) and (4.56) take the following form: where b can be either c or d for (4.54) and (4.56) respectively. In this form (4.58) may easily be made to match with the LHS of the respective equations. Finally, let us give a reason why the RHS of the two equations (4.54) and (4.56) are expressed in terms of the sources ∆ (n) 2k and the modes p nk and p (n−1)k . For (4.54) it is easy to justify since it is the Einstein equation for the rt component and therefore should depend on the sources and the fluctuation modes. To O( ) we expect only a linear combination of the form given as the RHS of (4.54). On the other hand, in the five-form EOM (4.56), the fluxes used to balance the system against any collapse [84] would in turn induce three-brane sources on the anti-D5 branes. The fluctuation modes should also affect these sources, and therefore the RHS of (4.56) is expressed as a linear combination of the sources and the fluctuation modes to O( ), justifying the above analysis.

The speed of sound in the strongly coupled plasma
We are now ready to do the two set of computations related to bulk viscosity: the speed of sound and the bound on the ratio of bulk viscosity to shear viscosity. The latter is again related to the speed of sound [30], so it will suffice to compute the speed of sound in the strongly coupled plasma. However before we go about computing the sound speed, let us present the formula for the ratio of the bulk viscosity ζ to the entropy density s, which was already derived in [84] for an appropriate choice of the quadrant: where Y x (r, ω) ≡ p 0x (r, ω) satisfies the differential equation given in (4.50). The result for the ratio of bulk viscosity to entropy density in a different quadrant can also be written down, and even their equivalence may be shown as in [84], but we will not do so here. Instead we will analyze the sound speed in the medium using all the ingredients we have collected so far.
One of the ingredients that we shall use extensively to compute the sound speed is the entropy density s. This has already appeared in (4.59) above, but the s appearing above is only the conformal result as the ratio (4.59) is already proportional to ≡ 3gsM 2 2πN . What we now need is the non-conformal correction to s. This may be written as: where the non-conformal corrections to s enters through the resolution parameter a 2 given in (4.12). We have chosen zero bare resolution parameter for simplicity and therefore, as evident from (4.12), a non-zero resolution already implies nonconformality in this set-up. One may worry that a zero bare resolution parameter may fail to capture the essential ingredients for a UV completion [58,84]. However that is not much of a concern here as we are not exploring the UV physics. Thus a cut-off r c will prominently feature in our results, as evident from (4.59) already. These all may be easily rectified, and we will elaborate it somewhat in the next section, but since no essential IR physics is lost in this simplified set-up, we will continue with this construction. There is however one issue that we do want to emphasize at this point and it has to do with the sign of the first expression in (4.12). Of course we naively expect a 2 to be positive, but the expression (4.12) involves various functions of log and dilog, so it will be instructive to check the sign of (4.12). Let us therefore start by defining x ≡ r 2 h r 2 c << 1, using which we can express (4.12) by: which is negative definite. This may trigger an alarm because a now becomes imaginary. Note that this problem does not arise if there is a bare resolution parameter a 0 , however small (as one may tune to be smaller than the smallest a 0 ). The way out of this conundrum is to notice that all expressions of fluxes etc involve a 2 and not a. Further, a 2 appears in the metric (4.5) as a combination r 2 + 6a 2 , and since we are only exploring the region r ≥ r h , the sign of a 2 does not create any problem here too. On the other hand, when there is no black-hole, r h vanishes, and so does a 2 (4.61). All this has also appeared in [85] − see discussions around figure 3 therein − for a more generic choice of a 2 given as eq. (2.63) in [85]. We can of course resort to a more conservative approach by writing an expression for |a| instead, and we shall do so in (5.13) in the next section wherein a non-zero bare resolution parameter will also be taken into account. Coming back, the entropy density computed above in (4.60) is proportional to powers of r h , so it vanishes when r h → 0. Additionally when r → r h , the entropy density receives corrections that take us away from the conformal value. These corrections may be easily quantified as powers of but we won't analyze it here 12 . 12 For example, to first order in and for r → r h we can sum up the series (4.61), or use (4.12), to show that the entropy density may be expressed as: where s 0 is the conformal value for the entropy density that can be read up from (4.60). To this order we can see that there is no r h dependence at the horizon.
Instead, at this point it may be instructive to point out the steps that went in the computation of the entropy density s. This would in turn effect the computation of the sound speed c s , since it depends upon s via: The entropy density may be determined directly from supergravity by first computing the energy-momentum tensors and then dividing the result by the temperature T . The energy-momentum tensor, on the other hand, arises from the variation of the action of the the form given by eq (3.120) of [58]. One may add a Gibbon-Hawkings term to it to control the boundary behavior, as evident from equations (3.121) and (3.123) of [58], but that does not alter the required linear term for our case. One may also add counter-terms to holographically renormalize the subsequent action, but since we are using a finite cut-off r c , it is not necessary to add them at this stage. This aspect has already been alluded to earlier, and here we see a more concrete realization of this. Putting everything together, the sound speed for r > r h will be given by: where x is the same parameter used in (4.61) before. Expectedly, the sound speed reduces to c s = 1 √ 3 in the conformal limit, and is smaller than 1 √ 3 when non-conformal corrections are included. One may justify this by looking at either of the two expressions in (4.63): the two terms, that account for the non-conformal corrections, are negative definite 13 when x < 1 (or r h < r c ). In the limit r h << r c , the sound speed (4.63) may be approximated by: (4.64) The above limit is not without its merit as we expect r c to be much bigger than r h , even if we restrict the dynamics completely to Region 1 of [58]. We can now use (4.59) and (4.64), to express the ratio of the bulk viscosity to shear viscosity in the 13 This may be easily seen from the first expression in (4.63) written in terms of the variable x in the following way: which is by construction smaller than c s = 1 √ 3 . Note that, for vanishing we get back the conformal result for the sound speed as one would expect.
following suggestive way: where η = 1 4π is taken at its conformal value to this order in , x = r 2 h r 2 c as before, and α x ≡ Yx(rc,0) Yx(r h ,0) is the ratio if the two fluctuations. We have also defined, without loss of generality, Y x (r h , 0) ≡ 1 r h for x << 1. We expect the ratio (4.65) to be positive definite, as (4.59) is positive definite. The second term is already positive, and the first term can become positive if α x is constrained in the following way: There is something puzzling about (4.65) that we should clarify right now. The way we have expressed (4.65) would seem to put an additional constraint on the ratio α x of the fluctuations as evident from (4.66). However such a constraint does not seem to follow from (4.59). In fact as long as x 2 < 13 16 both (4.59) and (4.65) should be positive definite. Since the expression (4.65) is basically a rewriting of (4.59) using the expression (4.64), it implies that (4.65) should not introduce any additional constraint of the form (4.66) on the ratio α x . Then why is there a new constraint? One way to argue for this would be to observe that the expression (4.65) is generic in the sense that it may be re-expressed as: where a(x) and b(x) are variations of the coefficients appearing in (4.65). This generalization however suffers from the appearance of explicit cut-off dependences of the respective variables. There could also be O( ) corrections to (4.63) and (4.64) that may change the coefficients of (4.65). These corrections appear from O( ) corrections to the temperature T , which we had identified to the horizon radius r h . To see this first let us take the cut-off temperature T c used in [58], which may be expressed as: where e 2B and e 2A are defined in (4.11) and (4.12) respectively. To O( ) the functional forms for the various parameters, i.e g(r) ≡ e 2B(r, ) and h(r) ≡ e −4A(r, ) , appearing in (4.68) may be determined exactly as: where L 4 ≡ 27gsN 4 , and note the appearance of higher powers of x in the black-hole factor g(r) defined at the cut-off r c . This series has similarity with the series (4.61) defined for the resolution parameter a 2 . The connection is of course spelled out earlier in (4.11), and once we plug (4.69) in (4.68), the temperature may be expressed as: ∞ n=1 x 2(n+1 n 2 (n+1)(2n−1) where the second line is in the limit x << 1. The corrections are exactly the ones that one would expect from switching on non-conformalities in the system. However note that even in the limit → 0, our expression for T c seems to have an additional factor of the form: which implies the cut-off dependence of the temperature. Clearly when we make r c → ∞ we recover the conformal result, but the appearance of x in (4.71) as well as in (4.70) means that UV completion is necessary to argue for the physical value of T c here. Naively taking r c → ∞ for non-zero will not give us the correct answer here, which of course resonates well with the UV completion discussed in [58]. Thus there is a way to holographically renormalize the system, following the procedure given in [58], that would take care of the log pieces in the metric and other variables in the problem. Once this is accomplished one may, in some restrictive sense, take r c → ∞. This is a specific UV completion wherein the UV cap gives rise to an asymptotically conformal theory. For such a case the temperature does take a physical value which may be expressed as: where Λ is related to the QCD scale for this model. The above is the so-called boundary temperature of [58] that we define at far UV. We will however need to define the temperature at any given scale, not just the UV, to avoid issues like (4.71) in the absence of any non-conformalities. Let us therefore take the following definition of the temperature: where a 1 and a 2 , which can be functions of x = r 2 h r 2 c , will be determined below. Note that T and T c are similar when a i take specific values extracted from (4.72). In general however, T should be the temperature that would occur naturally in this framework. This means we need to change slightly the formula for entropy in (4.60) by replacing r 3 h in (4.60) by r 4 h /T , with T given by (4.73). The sound speed will also change from (4.63) to the following: where the expected cut-off dependence appears from x as before. Clearly when a 1 = 1 πL 2 and a 2 = 0, we recover the sound speed computed in (4.63). However now, when both a 1 and a 2 are functions of x, the = 0 limit gives us: which takes us away from the conformal value of c 2 s = 1 3 in the conformal limit. This is not what we expect here, so we can use (4.75) to determine the functional form for a 1 (x). There are clearly two possible solutions for a 1 (x), namely: where b is yet an undetermined constant. The second choice is not acceptable in a theory that is holographically renormalizable, as it blows up when the cut-off is taken to infinity. This implies that T in (4.73) can only be: with constant b. What value can a 2 (x) take? To determine this we will need to study the full holographically renormalized temperature. This is in general a tedious exercise, but we can get a hint from the renormalized boundary temperature T c that we determined earlier in (4.72). To the first order in , the renormalized boundary temperature depends on log r h . This tells us that we can make the following ansatze for a 2 (x): where c 1 (x) and c 2 (x) are polynomials in x that do not have either log x or x −n pieces. The two functions c 1 (x) and c 2 (x) contribute to the full sound speed in the following way: x 2n n(2n − 1) where we see that the result is not so different from our earlier value for sound speed (4.63). The difference lies in the additional term proportional to da 2 /dx, which in turn would depend on how c 1 (x) and c 2 (x) depend on x. If c n (x) = −|c n | with constant c n , and x << 1, the sound speed is simple and is given by: where b > 0, and the signs are dictated by the fact that the beta function is negative and so the sound speed is smaller than 1/ √ 3. The additional term in the sound speed (4.80) means that the ratio of bulk to shear viscosities, i.e (4.65), changes to: where expectedly when |c 2 | = 0 we recover (4.65). The RHS now crucially depends on α x i.e on the ratio of fluctuations Y x (r c , 0) and Y x (r h , 0) satisfying (4.50). The equation (4.50) is difficult to solve, partly because of our ignorance of the precise sources a 4 defined in (4.49) using ∆ [1,2] and ∆ [2,3] via (4.45). Nevertheless, using the constraint (4.66), allows us to make the following ansatze for α x : where d 1 is another positive definite quantity. Note that we have not included a term proportional to x 3/2 in (4.82). This is done precisely to bring the ratio of bulk to shear viscosities (4.81) into the following suggestive form: where the cut-off dependence, compared to (4.81), now appears only through the last term. In the absence of the precise knowledge of d 1 and c 2 , this is the best we can do at this stage. However note that the negative terms in (4.83) cannot offset the sign of ζ/η in (4.83) because we have already established the positivity of (4.83) from the original expression (4.59). The concern however is the choice (4.82). How are we justified in the selective choice of the coefficients in (4.82)? How do we even know that such a choice will solve the EOMs? The answer to both the questions lies in the specific UV completion, or more appropriately on the distribution of the anti D5-branes in Regions 2 and 3. Once we plug (4.82) in (4.50), we can in principle determine the form of the sources ∆ [1,2] and ∆ [2,3] in a 4 , given via (4.49). One can then re-arrange the anti D5-distributions to match with the functional forms of ∆ [1,2] and ∆ [2,3] . This way the ansatze (4.82) may be justified.
Once this is settled, note that the negative definite second-term cannot be very large as all the three constants appearing there, namely d 1 , |c 2 | and b, are finite numbers. In fact for large cut-off r c , it is easy to establish the following upper bound on the coefficient d 1 : without loss of generality, where r h is the horizon radius. This also means that the cut-off dependent term in (4.83) will dominate over the negative definite second term. This is still consistent with the overall positivity of the ratio (4.83). However we now need the lower bound on d 1 . To determine this, we first note that the ratio of bulk to shear viscosities (4.83) satisfy the following bound: which would eventually control the behavior of the fluctuation modes studied in section 4.2. In the next section, we will study the sound speed and viscosity bound with non-zero fundamental flavors and with string coupling of order 1. We will rederive some of the above results, but in a different regime of the parameter space.
Such an analysis will hopefully shed light on the underlying universality of the results derived here.

Bulk viscosity at strong string and strong 't Hooft couplings with non-zero flavors
In the previous section we saw how one may study bulk viscosity, sound speed and the bound on the ratio of the bulk to shear viscosities at strong 't Hooft coupling using a gravity dual in type IIB theory. At this stage one may attempt few improvements in the present scenario by including both the flavor degrees of freedom as well as the UV regions. One may even ask the questions in the regime where the string coupling itself is of order 1, which of course still maintains strong 't Hooft coupling in the gauge theory side. The latter is however harder to study because it is the regime where even S-duality does not help. The question then is whether we can say something concrete in this regime of parameter space. One simple answer to the enigma may be to T-dualize the system to type IIA, by including the flavor branes, and then lift the configuration to M-theory. This should in principle accomplish the task, except that the T-dual scenario leads to a configuration of intersecting NS5-branes with the intersection region being blown up to a diamond [86,87]. This is not necessarily bad, and in fact in the past useful results have been drawn out of this configuration [88], but the requirement of keeping track of the NS5 degrees of freedom may thwart a simple analysis of the system. What we are looking for is a configuration with manifold and fluxes that we could use to succinctly address similar set of questions as in the previous section, avoiding the unnecessary requirement of including extra degrees of freedom. This is exactly where the mirror dual of the type IIB framework becomes handy. In fact, latticecompatible results pertaining to glueball spectroscopy were obtained in [101] and P(article)D(ata)G(roup)-compatible results pertaining to meson spectroscopy were obtained in [102], by working with the mirror dual.

The mirror type IIA model and its M-theory uplift
As discussed above, and also alluded to in Fig. 19, the model that we want to use here is the M-theory uplift of the type IIB scenario that we studied earlier. This is the MQGP model of [63,89] where at weak string coupling we have a type IIA description. One of important procedure that goes in the construction of [63,89] is the so-called delocalized mirror symmetry via the Strominger-Yau-Zaslow (SYZ) prescription [64]. This prescription involves a two-step procedure: one, by viewing the Calabi-Yau manifold as a special Lagrangian T 3 fibered over a base that is taken very large, and two, by performing three T-dualities over the T 3 fiber. In this subsection, we will provide some discussions on the details of the procedure.
The first requirement of a large base is important. This has to do with nullifying the contributions from open-string disc instantons with boundaries that appear as non-contractible 1-cycles in the special Lagrangian (sLag) T 3 fibered over the base. To see this more clearly, let us define three delocalized T-dual coordinates (x, y, z) which are basically proportional to (φ 1 , φ 2 , ψ) coordinates respectively that we encountered earlier. These coordinates are valued in the fiber torus T 3 via [63]: where s i are constants whose values may be derived from [91]. Interestingly, the choice of the coordinates (x, y, z) allows us to study the local geometry of the un-derlying manifold. Furthermore, using the results of [90] the following conditions, as shown in [89,103], are satisfied: for the underlying T 2 -invariant special Lagrangian manifold of [90] for resolved and deformed conifold. This immediately implies that, if the underlying resolved warpeddeformed conifold is predominantly either complelely resolved or deformed, the local sLag T 3 of (5.1) is then the required sLag to allow for the SYZ mirror construction via local T-dualities. Let us analyze thus further by taking the type IIB background given in (4.5) but now with e B = 1. The latter requirement is to just simplify the ensuing discussion. As we saw above, to enable use of SYZ-mirror duality via three T dualities, one is required to take a large base. This immediately means taking large complex structures of the aforementioned two two-tori of the sLag T 3 (x, y, z) fibration. One may easily implement this via the following considerations [92]: for appropriately chosen large values of f k (θ k ) with k = 1, 2. This choice does not change the local NS three-form flux, as was shown in [91,92]. Globally the underlying manifold can be a non-Kähler manifold as we discussed earlier. This is the advantage of using the (x, y, z) coordinates. On the other hand, the fact that one may be allowed to choose large values of f k (θ k ), was justified later in [63]. The main idea is basically the requirement that the metric obtained after SYZ-mirror transformation, applied to the non-Kähler resolved warped-deformed conifold, should resemble, at least locally, a non-Kähler warped resolved conifold. This means after incorporating (5.3) to (4.5), the (x, y, z) coordinates discussed in (5.1) will parametrize the local behavior succinctly. The global considerations will follow afterwards as shown in [93,92] 15 . 15 To justify the delocalization method while constructing the type IIA mirror a la SYZ triple-Tduality prescription [64] and its subsequent M-theory uplift one may argue the following. Consider the example of the mirror of a D5-brane wrapping the resolved S 2 with fluxes as studied in the first reference of [91]. The M-theory uplift can be made free of delocalization ensuring that one can construct a permissible G 2 structure manifold for the entire domain of validity of the delocalized coordinates. For example, in the delocalized large-complex structure limit and after a fixed ψ coordinate rotation, one obtains the SYZ mirror to be D6-brane wrapping a non-Kähler deformed conifold. Now, as shown in section 6 of [91], one can define an appropriate set of vielbeins to construct an explicit G 2 structure in terms of which the M-theory uplift of the previously obtained type IIA mirror could be rewritten, and which is valid for all values of ψ. In other words, the mirror In the local geometry, we can now perform three T-dualities 16 , first along coordinate x, then along coordinate y and finally along coordinate z, to get the local mirror manifold. The details of this construction, utilizing the results of [91], [92] was first worked out in [63]. The local mirror captures all the right properties of the expected dual configuration in the type IIA side, and then one may use the coordinates (φ 1 , φ 2 , ψ) to express the global metric as ds 2 IIA (see [91,93,92,63] for details). An additional ingredient that appears naturally from the SYZ procedure from the type IIB three and five-form fluxes as well as the axio-dilaton, is the one-form type IIA potential A. Such a one-form is useful to construct the M-theory uplift of the mirror type IIA as was shown in [63] 17 . The global M-theory metric takes the following form: where ϕ is the type IIA dilaton that appears from the mirror transform of the type IIB dilaton. Once the dilaton is allowed to take a non-trivial value, both in type IIB as well as in the mirror type IIA side, one starts seeing the effects of the flavors. This is simply because, in the type IIB side, non-trivial axio-dilaton shows up only when we switch on N f seven-branes. Of course, not all the N f seven-branes are required to be local D7-branes, but having D7-branes make the mirror picture more transparent as these would eventually contribute to the dilaton ϕ in the type IIA side. Once the dust settles, the g R 3 and g tt components appearing in (5.4) may be defined in the following way 18 : for ψ = ψ 0 coincides with the triple-T-dual-fixed-ψ rotated type IIA mirror obtained assuming delocalization. This essentially states that the type IIA mirror in equation (6.23) of the first reference in [91] obtained by descending to type IIA from arbitrary-ψ M-theory uplift will be the same as the fixed-ψ 0 type IIA mirror of equation (5.64) obtained using delocalization for ψ = ψ 0 . Hence we could just replace ψ 0 by ψ in the type IIA mirror obtained assuming delocalization.This therefore implies that the type IIA mirror is effectively free of the delocalization restriction. 16 Now also switching on e B in (4.5). 17 As is standard in such constructions, the one-form A may not be globally defined, although it's field strength will be. In the type IIB side such one-form will lead to either a RR two-form field or the axion depending on the T-duality direction. 18 Note that, unless mentioned otherwise, we shall always assume log r, in expressions like (5.5), is written as log r rc with r c being the cut-off radius. To avoid clutter, we will also take r c = 1 so that r remains dimensionless.

(5.5)
where r h is the horizon radius, and both g s N f as well as gsM 2 N are expectedly small 19 . Note also that both the metric components are independent of the resolution parameter a 2 . In fact the only metric component that depends on the resolution parameter would be the g rr component, whose explicit value is given by: where the full structure for a 2 will be given later. The functional forms for the coefficients appearing in (5.5) and (5.6) are determined by mapping the local metric to the warped resolved conifold metric 20 with a resolution parameter a 2 . In addition to that, and in the MQGP limit of [63], the α θ k factors for k = 1, 2 are angular coordinates such that around: we can allow the decoupling of the five-dimensional spacetime M 5 (t, x 1,2,3 , u) from the internal six-dimensional space M 6 (θ 1,2 , φ 1,2 , ψ, x 10 ). This decoupling is affected by making the Kaluza-Klein (KK) modes very heavy. The above discussions more or less summarizes the mirror construction as well its M-theory uplift. However it would be instructive to compare this with the type IIA brane construction of Fig. 20, which deals with both the UV and the IR brane configurations. The IR picture is of course the Klebanov-Strassler construction which is got by making a single T-duality along a direction orthogonal to the wrapped D5brane world volume, i.e along z of (5.1). This yields the RHS of Fig. 20, if we ignore the parallel NS5-brane. In other words, we get M D4-branes straddling between a pair of orthogonal N S5-branes whose world-volume directions are parametrized by (θ 1 , x) and (θ 2 , y) [95,96]. The mirror picture discussed here is then got by making two further T-dualities along x and y directions. Each of these T-dualities would yield Taub-NUT spaces from the corresponding NS5-branes [98]. The N f flavor D7-branes would yield N f D6-branes that are then uplifted to M-theory as KK monopoles [99]. These are also Taub-NUT spaces. Combining everything together 19 In section 4 we took g s → 0 with N, M very large and N f vanishing such that gsM 2 N << 1 and g s N f = 0. Here we take g s < 1 and N f = O(1) with N, M still very large. Again gsM 2 N << 1, but g s N f < 1. The latter can be implemented, for example, by choosing g s ∼ 0.4 and N f ∼ 2. Such a choice will guarantee that (g s N f ) m gsM 2 N n << 1 even for n = 1 and m = Z. Note however that g s → 0 does not always imply g 2 Y M → 0. We can have g 2 Y M = O(1) when N f = 0. This will be elaborated in section 6.4. 20 Recall that globally we can only put a non-Kähler metric on the resolved conifold [62].
then leads to a seven-dimensional manifold with a G 2 structure and with G-fluxes. This configuration is precisely equivalent to the uplift of the wrapped D5-branes on a warped resolved conifold of [92,91,62].

Quasi-normal modes, attenuation constant and the sound speed
Let us now discuss the main ingredient of our construction, namely the quasi-normal modes in the dual gravitational background. The procedure involves a few steps that we lay down in the following. Building up on the ideas developed in [100] and [103], and using gauge-invariant combinations of metric perturbations invariant under infinitesimal diffeomorphisms, in other words: as discussed in [100], the gauge-invariant combination of scalar modes of (M-theory) metric 21 perturbations was constructed in [103]. A discussion of the same appears in Appendix A.
Next, we work near the decoupling limit prescribed in (5.7), and choose the other three angular coordinate (ψ, φ 1,2 ) in the mirror metric (5.4) as ψ = 2nπ, with n = 0, 1, 2 and small φ 1,2 . We also choose our radial variable henceforth as u ≡ r h r . Using these we can define: (5.9) in the context of the gravitational dual of large-N thermal QCD with N f = 0, where N f is the number of flavors. This functional form of B(u) appears in the construction of the gauge-invariant Z s (u) in the following way: (5.10) where the H ab functions are given in appendix A, and T is the temperature whose form will be given below. Note that the upshot of appendix A is essentially the construction of the gauge-invariant Z s (u) that will satisfy certain EOM to be elaborated in the following 22 .
In obtaining an EOM for Z s (u), we will make use of q 3 = q πT , w 3 = w πT where T is temperature that appears in (5.10) above. We will express T in terms of all the 21 As discussed above, this corresponds to the local uplift of the delocalized Strominger-Yau-Zaslow [64] type IIA mirror of the holographic type IIB dual of [58] of large-N thermal QCD, having integrated out the six angular directions as in [104], up to NLO in N in the MQGP limit of [63]. 22 Our emphasis here would be to determine the EOM up to NLO in N .
variables that appear in the metric. To proceed, and for later brevity, we start by defining the following quantity: where (k, j) will be integers. Now assuming the resolution to be larger than the deformation in the resolved warped-deformed conifold in the type IIB background of [58] in the MQGP limit, and using the decoupling limit (5.7), the temperature T may be expressed in the following way (see also [103]): where G µν is the M-theory metric (5.4), and C 11 (u) may be extracted from (5.11). We can also go to the limit where α θ i are O(1) numbers. This way the temperature may be written completely in terms of the resolution parameter a 2 and the horizon radius r h . Interestingly when a 2 >> r 2 h , the temperature is expressed in terms of inverse r h . Otherwise, the temperature is proportional to r h . In the limit of vanishing flavors, small bare resolution parameter, and large cut-off, the expression for the temperature becomes identical to what we took on the type IIB side (see (4.77) and (4.78)). The bare resolution parameter in type IIB side, as given in (4.12), was taken to be zero. A natural question then is to ask what happens if we take non-zero bare resolution parameter 23 . A particular choice of a(u) can be: (5.13) this way b may serve as the bare resolution parameter in (5.12) and c 1 (u), c 2 (u) are some slowly varying functions of the u parameter (not to be confused with b and c 1 , c 2 taken in (4.77) and (4.78)). One may compare (5.13) with the type IIB resolution parameter (4.12) in the limit b → 0. The functional forms in the two cases are similar, but not identical. This is intentional because the choice (5.13) allows us to perform computations in the mirror side more efficiently compared to the choice (4.12). This in turn will also effect some of our final results, so comparison with the type IIB side will have to be done more carefully.
With these definitions at hand, we are now ready to write down the equation of motion for Z s (u) appearing in (5.10). This may be expressed in the following way: (5.14) 23 Note that allowing a bare resolution parameter in the type IIB side allows us to perform the SYZ mirror transformation more efficiently [91]. Here however we will use the word bare to denote the part of the resolution parameter that is independent of g s N f and gsM 2 N . Of course the r h independent piece of a(u) vanishes.
which is a second-order differential equation in u whose solutions will tell us the precise gauge-invariant variables that we seek here. This equation depends on two non-trivial functions of u, namely m(u) and l(u), whose functional form will be important. Both these functions may be expressed in terms of C kj (u) in (5.11) as well as certain other functions that we shall elaborate in the following. We start with m(u) which may be written as: where note the appearance of C 21 (u) and C 23 (u) defined in (5.11) as well as A i (u) that form the various coefficients above. The function A 1 (u) may be written as: which expectedly simplifies for vanishing b. On the other hand, at the boundary when u vanishes, A 1 (0) is proportional to (ω 2 3 − q 2 3 ) 2 which is now expressed in terms of T defined at u → 0 instead. Similarly, A 2 (u) may be written in the following way: which is again proportional to (ω 2 3 − q 2 3 ) 2 at the boundary u → 0. When b vanishes, A 2 is unaffected, but the term itself comes multiplied with b in (5.15), so decouples completely. The third term in (5.15) is proportional to b 2 , so we expect it to decouple in the limit of vanishing b. To see what happens at the boundary, i.e when u → 0, we express A 3 (u) as: 18) which is expectedly proportional to (ω 2 3 − q 2 3 ) 2 , but the term itself decouples because it appears together with a factor of u in (5.15), much like the previous term in (5.15). The remaining two coefficients, A 4 (u) and A 5 (u), are in similar vein as (5.16), (5.17) and (5.18) and share much of the same properties as above. They take the following form: and become proportional to (ω 2 3 − q 2 3 ) 2 when u → 0 at the boundary, but do not decouple in a simple way as before. In this sense they share the property of the first term in the definition (5.15). The boundary behavior of m(u) can then be given by the following limiting expression: 12u , (5.20) where ≡ 3gsM 2 2πN is the same expansion parameter that we used in section 4. Both the C 2j factors behave as log u, but are suppressed by as well as g s N f (5.11) (any constant factors get suppressed by from (5.15)). Thus m(0) seems to blow up as 1 u or log u u . This preliminary analysis however is naive because precisely in this regime the UV cap modifies the boundary behavior appropriately to avoid any such pitfalls. Therefore a more relevant question to ask is the behavior of m(u) at the horizon, i.e when u → 1. We will analyze this below, but before that let us discuss the behavior of the other function l(u) appearing in (5.14).
The expression for l(u) turns out to be very large so we shall suffice ourselves by demonstrating that the horizon u = 1 is an irregular singular point whenever N = 0. In the following we give below the expansion of l(u) about u = 1 to see the same: where C 21 (1) may be extracted from (5.11) by putting u = 1 therein. For vanishing bare resolution parameter (5.21) vanishes, so a minimal resolution is necessary to see the behavior at the horizon. The above expression for l(u) near the horizon is what we need, and we could also go to the u → 1 limit for m(u) in (5.15) to determine its behavior at the horizon. However the results are expressed in terms of both q 3 and ω 3 . To elaborate further, we need to first express ω 3 in terms of q 3 and then identify the subsequent behavior of m(u) and l(u) at the horizon. To this effect, we make the following ansatz: 22) and substitute into m(u) and l(u). Here α(u) and β(u) are certain functions whose values will be determined near the horizon. We then first perform a small q 3 expansion, followed by an expansion around u = 1 and lastly a large N expansion. The procedure is straightforward albeit a little tedious. After the dust settles, we come up with the following expansions for m(u) by keeping terms up to O q 3 , gsM 2 N and the most singular term in u near u = 1, namely: with C 21 as in (5.11). For vanishing bare resolution parameter, there is a further simplification: m(u → 1) behaves simply as 1 1−u as may be easily seen from (5.23). On the other hand, the behavior of l(u) at the horizon may be read more directly from (5.21) as: which expectedly vanishes for vanishing bare resolution parameter, and has the required irregular singular point. The functional forms for m(u) and l(u), expressed using the dispersion relation (5.22), and analyzed near the horizon u → 1 is essentially the regime that we want to concentrate here. We can also study the system at the boundary by attaching an appropriate UV cap controlling, in turn, the behavior of m(u) and l(u), but this will not be the emphasis of this section. Our aim would be to explore the near horizon behavior where one sees u = 1 as an irregular singular point of (5.14). To proceed, let us make the following ansatz for Z s (u): where we shall assume [S (u ∼ 1)] 2 > |S (u ∼ 1)|. This derivative requirement essentially converts (5.14) to a simple quadratic equation in S (u) with coefficients m(u) and l(u). The solutions are: At this stage it would be interesting to ask what happens when the derivative condition is not satisfied. Clearly in this case we will get a second order inhomogeneous differential equation which becomes homogenous when the bare resolution parameter vanishes. Generically it is harder to deal with the inhomogeneous case, because of the complicated forms of m(u) and l(u), and the homogenous form is not a suitable choice for the system undergoing SYZ transformations [64,91]. Thus the simplification and calculability attained from the derivative requirement guarantee not only analytic control, but also solutions not far from the regime of interest. With this in mind, the next set of steps are standard 24 . Choosing the minus sign in (5.26), one obtains the following: which has a simple pole structure of − 1 2(u−1) in the limit of vanishing bare resolution parameter b. The other parameters appearing in (5.27) are C 21 (u) defined in (5.11), and the two functions B 1 and B 2 defined in the following way: Let us now assume that q 3 → 0 as N −1−κ with κ > 0. One might worry that imposing this one would obtain, near u = 1 − which is an irregular singular point − a solution of the type e S(u) = (1 − u) γ F (u) implying u = 1 to be a regular singular point. This does not happen, and therefore demanding the vanishing of the residue of S (u) at u = 1 gives the following values for b, α and β: where C 21 (u) is defined in (5.11). In fact this is all we needed to determine both the sound speed c s as well as the attenuation constant Γ because the first term in (5.22), i.e the term proportional to q 3 , gives us the sound speed as 25 : Although we do not undertake here, a more generic analysis away from u = 1 can be performed and from there the limiting form of (5.27) can be ascertained. Needless to say, the results match. 25 Recall that we can express the dispersion relation (5.22) in terms of sound speed c s , shear viscosity η and bulk viscosity ζ as: where one can see that the result differs from the conformal answer of 1 √ 3 expectedly by the gsM 2 N and g s N f factors. Even in the absence of the fundamental flavors N f , the sound speed deviates from the conformal answer. The form of the deviation is consistent with what we had earlier in (4.63) although the precise factors differ. This is understandable in the light of the different choices of the supergravity parameters in the type IIB and the M-theory pictures.
The attenuation constant Γ may now be easily extracted from (5.22) by plugging in the value of β from (5.29). To NLO in N , Γ may be written as: where T is the temperature, and again we see that even in the absence of fundamental flavors, the attenuation constant differs from the conformal value of 1 6πT . The parameter C 21 (1), defined in (5.11), becomes identity when N f = 0, so the deviation from the conformal value is solely governed by gsM 2 N .

The case with a vanishing bare resolution parameter
Let us now discuss what happens if one sets b = 0 in (5.13). We briefly dwelt on this earlier, wherein we saw how (5.15) and (5.21) behave when b vanishes: (5.21) completely decouples but some remnants of (5.15), as seen from (5.16), (5.17), (5.18) and (5.19), survive. Interestingly, this now makes u = 1 a regular singular point of (5.14). To proceed, let us then rewrite Z s (u) using an analytic function F (u) in the following way: where C 23 (1) can be evaluated from (5.11) by putting u = 1 in the required expression and we have defined ≡ 3gsM 2 2πN as the non-conformality factor. With the definition (5.32), the EOM (5.14) becomes: which is a second order homogenous differential equation with coefficients defined by parameters a 1 , a 2 and b 2 . The 1 u and log u u terms are remnants of the equivalent terms in (5.20). The a 1 and a 2 coefficients take the following form: where s is the entropy density. This means α(u) in (5.22) is related to c s and β(u) in (5.22) is related to shear viscosity and bulk viscosity combination, or the attenuation constant Γ. However since we are analyzing the system close to the horizon, i.e u → 1, the relevant parameters for us will be α(1) and β(1).
where C 21 (1) is extracted from C 21 (u) in (5.11). It is interesting to note that the combined expression with a 1 and a 2 may be succinctly expressed as: which simply converts C 21 (1) in (5.34) to C 21 (u). This is expected from the way we represented the EOM for Z s in (5.14). On the other hand, the form for b 2 in (5.33) may be expressed as: where C 23 (u) is given in (5.11) (note the appearance of C 23 (u) instead of C 23 (1), much like what we have in (5.35)) and D is defined in the following way:  (1). Implementing this, F (u) takes the following form: where d 1 and d 2 are constants. We can also motivate the replacement C 21 (u) → C 21 (1) in both (5.35) and (5.36) in the following way. In (5.35), this amounts to dropping the log u term near u = 0 compared to the log N term in the large N limit i.e making a 2 = 0 in (5.33). In (5.36), this amounts to just keeping terms of O gsM 2 N as the log u term in C 21 (u) is already suppressed by g s N f (see (5.11)).
The u → 0 limit may seem a bit puzzling because so far we have analyzed the system near the horizon i.e near u → 1. However in (5.14), for vanishing bare resolution parameter, as we saw earlier, l(u) vanishes and the EOM is solely governed by m(u) (5.15). This may be defined both at the boundary (5.20) and at the horizon (5.23). Thus extrapolating F (u) to the boundary is still well defined, modulo the subtlety of including a UV cap.
Let us now go to the various choices of the parameters d 1 and d 2 in the solution (5.38). If one sets d 2 = 0 then the small-u expansion of F (u) will be given by the small-u expansion of the first part of the solution (5.38), i.e the d 1 part of the solution in (5.38), as: To analyze the boundary conditions, first let us make Im b 2 = 0. There is already a problem at this stage, but let us still carry on. At the boundary u = 0, if Re b 2 < 0 then (5.39) blows up as exp . This behavior will persist even if we include the a 2 log u piece in (5.33). Thus to be able to impose Dirichlet boundary condition on F (u), i.e impose F (u) = 0 at the boundary, one needs to set b 2 = 0. Now, substituting b 2 = 0 in (5.33), one obtains: where d 3 and d 4 are constants. We can fix d 4 in terms of d 3 by demanding Dirichlet boundary condition on F (u). This immediately gives us: This is almost what we need, except for an important caveat. Putting b 2 = 0 (or Im b 2 = 0) in (5.36) gives us: where D is defined in terms of α, β, q 3 as well as gsM 2 N in (5.37). Since the RHS of (5.42) is a c-number, and β defined in (5.22) is a pure imaginary number (at least at the horizon), the equation (5.42) can only be solved if: The above forms of α and β are unfortunately not acceptable as they will not only lead to the wrong sound speed and attenuation constant, but also take us away from the perturbative regime where all our computations were focussed. One might think that this could be rectified if we had started off with a non-zero Im b 2 , but unfortunately the conclusions don't change much as F (u) would still oscillate infinitely fast or blow up. The above failure is a near miss, but it teaches us an important lesson about the choice of the function F (u): the boundary conditions are subtle and important, but one still needs to choose the function carefully, as any arbitrary choice may take us away from the perturbative regime of interest. Therefore at this stage there are two ways to fix the function F (u). One, we may not impose a Dirichlet boundary condition, and allow a non-normalizable functional form for F (u). Two, we again allow for a Dirchlet boundary condition, but choose the functional form for F (u) a little differently from the previous choice (5.39). The latter case is easier to implement, so we start by setting d 1 = 0 in (5.38). This gives us: We see that if Im b 2 = b 2 = 0 then we would encounter similar problem as in (5.43).
On the other hand, if we allow Re b 2 < 0, then we can control the amplitude of oscillation from the Im b 2 piece, provided: The above set of conditions does help us to solve for F (u) as before allowing the required Dirichlet boundary condition at u = 0, although the procedure for getting the exact functional form for F (u) is not as straightforward as in (5.40). However the condition Re b 2 < 0 now leads to the following condition on α and β ≡ iγ: (5.46) where = 3gsM 2 2πN is the non-conformal factor. Although the above condition gets further refined by (5.45), getting α and γ satisfying (5.46) can at least indicate the behavior of α and β with respect to gsM 2 N . A careful look at (5.46) tells us that if both α and γ are proportional to , then q 3 gets constrained. This cannot be right, so it seems the only way to satisfy (5.46) would be to take α and γ to be inversely proportional to , much like what we had in (5.43) before. Such a choice will again take us away from the perturbative regime of interest. Thus it seems the only way to analyze the behavior of α and β from the boundary u = 0 point of view is to allow for a non-normalizable F (u). This resonates well with the analysis of fluctuation modes of the metric in section 4 where p nk and Γ 0k functions were both non-normalizable functions. Note that we did not encounter these issues while studying the b = 0 case because the analysis was performed at the horizon u = 1 where these subtleties were not visible.

Shear viscosity, entropy and the bulk viscosity bound
After this detour, it is time to go back to our analysis of bulk viscosity and the bound on the ratio of the bulk to shear viscosities. To proceed, we will first quantify the functional forms of f 1 (θ 1 ) and f 2 (θ 2 ) in (5.3) in the following way [94]: 47) where α N and the choice (5.7) ensure large base for implementing the SYZ [64] mirror transformation. Recall the necessity of a large base in our set-up to nullify certain disc instanton contributions. The choice (5.47) is essential to compute transport coefficients and entropy in the M-theory uplift of the mirror set-up. We can now combine this with the value of the bare resolution parameter b = √ 6 that we got in (5.29), and using results of [63], [103], we can show that the shear viscosity near the horizon takes the following form: 48) whereĈ 21 (1) is defined as C 21 (u) with u = α θ i = 1 in the definition (5.11). Note also the appearance of α θ 1 and not α θ 2 in (5.48). This is because θ 1 and θ 2 defined in (5.7) approach zero at different rates so the former got selected in the computation 26 . We have also introduced a coefficient Υ in the formula (5.48) for η, whose value will be fixed soon. It is now time to compute the entropy density s. The procedure for computing s remains similar to what we did in section 4, although the choice of the mirror variables differ from the type IIB case. This implies that the entropy density at the horizon may now be expressed as: whereĈ kj (1) is defined for C kj (u) with u = α θ i = 1 in (5.11). One may now compare (5.49) with (4.60) as well as the entropy computed in [58] where we see similar suppressions with respect to gsM 2 N and g s N f . The precise coefficients understandably differ because of the different choices of variables alluded to above. One may get away from this by choosing a uniform definition of the variables in all the models. However this suffers from a reduction in the efficiency of computations of physical quantities in some models and increase in others 27 .
The above discussion however does not spell out a failure to compare the physical quantities in different models; rather one should interpret the validity of different results to be at different range of parameter values. For the choice of parameters in the mirror set-up, and using (5.48) and (5.49), we can now express the ratio between shear viscosity and entropy density as: where note the absence of the parameter α N from (5.47). We also wrote the first term as 1 4π . In the absence of the gsM 2 N correction, this should be the conformal result [106], and therefore we have used this to fix the parameter Υ in (5.48) as: There are a few issues regarding the ratio (5.50) that we should take into account now. First, observe the appearance of an inherent cut-off in (5.50). This appears through the log r h term above as log r h rc , where r c is the cut-off radius. Physical result should not depend on the cut-off so we should interpret (5.50) carefully.
The r c dependence in (5.48) for example should remind us of a similar r c dependence of shear viscosity in the type IIB side as given in eq (3.198) of [58]. The introduction of UV cap to the geometry contributed an additional piece as eq (3.200) in [58]. This eventually led to the ratio of the shear viscosity to the entropy density being given by eq (3.222) therein that depended upon the UV degrees of freedom N uv as e −Nuv . The result for infinite UV degrees of freedom was exactly 1 4π , so we should expect similar result for our case too. However the analysis of η s in (5.50) is done at the horizon with a UV cut-off r c , and one may easily see that the cut-off dependence is log r c which is an expected answer for a QCD like model. This means that, even with a UV cut-off, we expect η s to be at least bigger than 1 4π so that the KSS bound [106] is not violated. In (5.50) it is easy to see that the r h dependent terms are positive definite because log r h rc < 0. Thus if we define c 1 , which is as yet an unfixed function, as: with σ as another undetermined function, then η s > 1 4π . The UV cap can then change the result accordingly, but we will not elaborate on this here anymore. At this stage, it will simply suffice to see that the KSS bound is not violated.
All the ingredients are at hand now to compute both the bulk viscosity ζ as well as the ratio of the bulk to shear viscosities i.e ζ η . As we saw earlier, the shear and the bulk viscosities are connected by the following relation: where s is the entropy density (5.49), T is the temperature (5.12) and Γ is the attenuation constant (5.31). One can therefore use (5.53) to express the ratio ζ s in terms of Γ and η s as: 54) where κ 0 ≡ 364π √ 6 45 , C jk (1) is given by (5.11) for u = 1, andĈ jk (1) is given by (5.11) with u = α θ i = 1. The overall factor of gsM 2 N is interesting and crucial: it tells us that the ratio (5.54) is zero for conformal theories. This is of course consistent with what we had in section 4, and we note that the bulk viscosity may be easily derived, to this order in gsM 2 N , by simply multiplying (5.54) by the conformal entropy density. Any non-conformal corrections to s will change the bulk viscosity only to higher orders in gsM 2 N . Note also that, in the limit of vanishing fundamental flavors i.e N f = 0, the ratio (5.54) takes the following form: where we have inserted back the cut-off radius r c (which was taken to be 1 so far). Looking at (5.55) one might be tempted to compare it with the bulk viscosity that we got in (4.59), which was expressed using the fluctuation mode Y x satisfying (4.50).
Both have somewhat similar r h /r c dependence, but the exact factors differ. This has already been alluded to earlier because of the different choices of the parameters in the two theories. However, as we discussed in section 5.3, the ratio (4.59) is derived for vanishing bare resolution parameter whereas (5.54) is derived with non-zero bare resolution parameter. This of course is not the only difference. The zero bare resolution case, according to section 5.3, involves study of quasi-normal frequencies whereas the result (4.59) is derived from the study of fluctuation modes Y x . The point of comparison between the two results maybe that both involve certain nonnormalizable functions at the boundary u = 0. Plugging in the non-normalizable function F (u) in (5.32) will help us find α and β in (5.22), which in turn may be compared to (5.54).
Finally, the ratio of bulk to shear viscosities may now be determined from (5.54), to first order in gsM 2 N , by taking the conformal limit of (5.50). The result is similar to what we have in (5.54) up to a factor of 4π: 56) where κ 0 is defined earlier and we have inserted back r c , the cut-off radius. To see whether (5.56) does not violate the Buchel bound [30] we will have to determine c 1 and c 2 in (5.56). In (5.52) we expressed c 1 in terms of a negative definite function −|σ| at the horizon u = 1, assuming c 2 to be a positive definite quantity there. However underneath this choice was the assumption that both the bare resolution parameter b and the full resolution parameter a in (5.13) remain positive definite. As long as b > 0, this could still be made true with c 2 > 0. However b can be zero, as we saw in sections 4 and 5.3, and in this case c 2 > 0 will make a < 0 in (5.13) with the choice of c 1 in (5.52), rendering the whole construction meaningless. One might think that c 1 could be changed, but then the KSS bound [106] will be affected. Therefore it seems the only way to avoid any contradictions is to take c 2 = −|c 2 | with: at the horizon u = 1. By construction this preserves the KSS bound [106], and keeps the resolution parameter (5.13) positive definite 28 . With this at hand, it is now time to see if the ratio of bulk to shear viscosities (5.56) preserve the Buchel bound [30]. We will start with the simplest case of vanishing flavor i.e N f = 0. Referring back to sound speed (5.30) and the ratio (5.56), we get: where c 1 and c 2 now satisfy (5.52) and (5.57) respectively. Since log r h rc < 0, all terms in ζ η in (5.58) are positive definite. In the limit where r c > r h , the ratio of bulk to shear viscosities may be related to the sound speed as: which clearly satisfies the Buchel-bound [30]. Interestingly, for r c >> r h exp 11+|c 1 | |c 2 | , one may ignore the second piece in (5.59) and the ratio (5.56) may solely be expressed in terms of 1 3 − c 2 s . In either case, one may easily infer from (5.58) that the Buchelbound is always satisfied, at least for vanishing fundamental flavors N f .
What happens when N f = 0, i.e when we switch on fundamental flavors? Both, the ratio of bulk over shear viscosities and sound speed, have been computed above in (5.56) and (5.30) respectively. So it's time to combine them to see whether the specified combination of them satisfy the Buchel-bound. It is easy to see that the bulk to shear ratio (5.56) may now be expressed as: where the C ij andĈ ij terms may be extracted from (5.11) using the limits u = 1 and u = α θ i = 1 respectively. By switching off the g s N f terms one gets (5.59) from (5.60), so the question now is whether the gsM 2 N terms in (5.60) can again be positive definite.
It turns out, with some algebraic manipulations, one may rewrite the gsM 2 N terms appearing on the RHS of (5.60) in the following suggestive way: where σ 0 ≡ 201 20π log 4 + 603 20 ≈ 100.86 is a positive coefficient. Since N is very large and r c >> r h , every term in (5.61) can be shown to be positive definite, and the negative piece σ 0 g s N f does not change anything as long as log N log rc r h >> 160. The latter is not a constraint as we saw above 29 . Thus generically our model satisfies the Buchel-bound [30], and comparing (5.59) and (5.61), we see that there is in fact a new bound on the ratio of bulk to shear viscosities given by: To see that the terms on the RHS of (5.61) can be positive definite, choose N ≡ exp (n 0 + 60.3) with n 0 being a very large number approaching infinity, and r c ≡ r h exp (n 1 + 2.5125) with n 1 being another large number, not necessarily infinite. The condition for positivity of the RHS of (5.61) is n 0 n 1 > 160. This is easily achieved because going to the gravity dual description forces us to choose both n 0 and n 1 very large.

Type IIA spectral function and the viscosity bound at strong coupling with non-zero flavors
In the above section we found an interesting bound (5.62) for the ratio of bulk to shear viscosities at strong string and strong 't Hooft couplings. In fact the form of the bound seems consistent over the whole strong 't Hooft coupling regime as is obvious from the weak string coupling bound that we got earlier in (4.85): both (5.62) and (4.85) are proportional to 1 3 − c 2 s , although their coefficients differ. On the other hand, the weak 't Hooft coupling bound differs by being proportional to the square of the strong coupling bound as shown in (2.18). The reason for the different results at the two ends of the couplings has been motivated in section 2.2. Loosely, it is the ratio of the shear viscosity over entropy density that creates the difference at the two ends. At weak 'tHooft coupling the ratio is not a constant and is given by (2.17), whereas at strong 't Hooft couping we expect the ratio to be a constant [106]. This is a reasonably strong argument for justifying the difference between the two bounds, despite the fact that we have no control on the dynamics at the intermediate coupling regime as argued in section 3. The very weak coupling results have been justified in great details, and in sections 4 and 5, we provided some justifications for the strong coupling results. However one might be interested in deriving the bound at strong coupling directly from the spectral function, as such an analysis will be in line with the discussions of section 3. Further, we make the following observations: • The ratio of the bulk viscosity ζ to entropy density s is of O gsM 2 N as we saw in (5.54), and the ratio of the shear viscosity η to the entropy density is dominated by the conformal result plus an O gsM 2 N correction term from (5.50). This means up to O 1 N the ratio ζ/η would mimic ζ/s.
• The gauge and the metric perturbations may be required to be considered simultaneously − see subsection 4.2 of [81] and references therein.
• The correlation of gauge fluctuations, A x i A x i for i = 1, 2, 3, along the same direction could hence mimic the spirit behind the correlation of the metric perturbations, h x i x i h x i x i , along the x i axis relevant to the evaluation of bulk viscosity as, for example in [38] or in section 4.
The above three observations provide the necessary motivation for this section. Therefore we would like to evaluate the aforementioned gauge-field correlation function (in the hydrodynamical limit using the prescription of [107]) and see if one obtains a linear bound seen in (4.85) and (5.62). Even if this may not be explicitly tied to ζ/η, we feel the result obtained in this section, in itself, is sufficiently interesting.

Background gauge fluxes and perturbations on the flavor branes
Our starting point is configuration of N f D6-branes in the type IIA mirror set-up. For our purpose, we will however isolate one D6-brane and study world-volume dynamics on it. Alternatively, one can view this as D6-brane probing a non-Kähler warpedresolved conifold with N f flavor D6-branes. The DBI action for a single D6 brane is given as: with 2πα = 1. Here the worldvolume directions of the D6 brane are denoted by the coordinates: (t, x 1 , x 2 , x 3 , Z, θ 2 , ϕ 2 ), with (t, x 1 , x 2 , x 3 ) as the usual Minkowski coordinate, Z as the newly defined radial direction and two angular coordinate (θ 2 , ϕ 2 ); Z is related to the usual radial coordinate r as r = r h e Z and ϕ 2 is the local value for the angle φ 2 (for more details, see sections 3 and 4 of [102]).
In the above, and as before, ϕ denotes the type IIA dilaton which is the triple T-dual version of type IIB dilaton. The pullback metric and the pullback of the NS-NS B field on the worldvolume of the D6 brane are denoted as g and B in (6.1). F is the field strength for a U(1) gauge field A µ , where the only nonzero component of the same is the temporal component A t . In the gauge A Z = 0, the only nonzero component of F is F Zt = −F tZ . Combining together the symmetric g field and the anti-symmetric B field as G ≡ g + B and expanding the DBI action up to quadratic order in A, we get: The second term in (6.2), is coming because of the anti-symmetric B field in G.
As none of the fields in the DBI action depends on the angular coordinates ϕ 2 , the integrand in equation (6.2) is independent of the same. Also we choose to work around the same small values of both θ 2 and θ 1 given by (5.7) earlier. Hence, the upshot is that the integration over θ 2 and ϕ 2 is trivial and we denote Ω 2 as the factor one gets after the integration over θ 2 and ϕ 2 , such that: 3) The equation of motion for the temporal gauge field A t (Z) as obtained from the above lagrangian in (6.3) is given as: At this point we can use the precise functional forms for the background data, namely G tt , G ZZ as well as the dilaton ϕ 2 , to rewrite the above equation in the following form: where a 2 is the resolution parameter, (α θ 1 , α θ 2 ) are the two angular values in (5.7) and C is the integration constant. In the large Z and small a 2 limit, (6.5) yields: which appears from the fact that the second line in (6.5) dominates over the first line. The large Z limit is also the large r limit where one might be concerned about UV issues appearing from AdS cap. This is however not much of a worry at this stage because as long as r h e Z >> a, (6.6) continues to hold. With this in mind, the solution to (6.6) is: where we have used A t to express the background value to avoid confusion. The other parameters appearing in (6.7) are C 1 , which is yet another constant and Ei, which is the exponential integral 30 . In the second line of (6.7) we have shown the first dominant piece in the large Z limit. Higher powers of 1 Z can then be ignored. This background value also prepares us to study the fluctuation of the gauge field components. For example we can express the gauge field appearing in (6.3) as: where the fluctuation A µ only exists along the directions µ = (t, x 1 , x 2 , x 3 ) due to the particular gauge choice and depends only on the radial variable Z. Including the perturbations in the lagrangian of the DBI action (6.1), one gets: with F as the field strength for the gauge field fluctuations. Now defining G ≡ g + B + F and again expanding the above lagrangian up to quadratic order in the gauge field fluctuation one gets: This definition can be used for positive values of x, but the integral has to be understood in terms of the Cauchy principal value due to the singularity of the integrand at zero.
Writing the field strength F in terms of the gauge field fluctuation A µ and after doing some simplifications in terms of the interchange of indices, one can write the above lagrangian as: The second line in equation (6.11) is a total derivative term and equating the last line to zero for any arbitrary A β gives the equation of motion for the gauge field fluctuation: The total derivative term in (6.11) does not necessarily have to vanish at Z → ∞, as there could be non-normalizable modes serving as sources for the dual gauge theory operators. Our EOM in (6.12) however is not affected by this, and in the following section we will discuss possible solutions of (6.12).

Equation of motion for gauge field fluctuations
To derive the equation of motion for the gauge field fluctuation, we first need to write down the fluctuating field as the following Fourier decomposed form: where we assume the fluctuation to have momentum along x 1 direction only, with k 0 = ω, k 1 proportional to q and k 2 , k 3 arbitrary. Now, the equation (6.12) has a free index β and for β = (t, x 1 , x 2 , x 3 , Z), one gets a total of five equations of motion. For example for β = Z, plugging in (6.13) in (6.12) yields: where the RHS vanishes because of the antisymmetry of G [αβ] . The dilaton does not appear because it is independent of the four-dimensional spacetime coordinates. The above equation relates A t with A x 1 . On the other hand if we take β = t, we get the following EOM: which now relates A t with A t and A x 1 . A somewhat similar equation appears when we choose β = x 1 in (6.12), namely: At this stage one can easily verify that pluging in (6.14) in (6.15), we can get (6.16). This shows that the above three equations (6.14), (6.15) and (6.16) are self-consistent. Finally, one may find the equations for β = x 2 and β = x 3 . We expect them to be equivalent, and are given by: To proceed further, we will have to define gauge invariant variables. For our case, there would be two such variables E x 1 and E β with β = x 2 or x 3 , expressed in the following way: With these new variable the three equations in (6.14), (6.15) and (6.16) can be cast into a single second order equation involving E x 1 . Even more obviously, the fourth one for β = x 2 or x 3 , can be rewritten in terms of the new variable E T . Moreover, in the zero momentum limit, i.e in the limit q → 0, it can be shown that the equation involving E x 1 is the same as the one involving E T , given as: 19) implying that in the zero momentum limit all we need is to solve one second order differential equation. This is of course a huge simplification, and one can even rewrite (6.19) in the following suggestive way: where all functions appearing above are only functions of the Z variable. In fact P (Z) and Q(Z) may be easily seen from (6.19) to take the following form: The suggestive way alluded to above is that the above equation (6.20) can be recast in a Schrödinger like form by certain redefinition of the variables involved in the following way: with E T defined as E T (Z) ≡ P (Z)E T (Z), and V E T is the potential term that is expressed as: The Schrödinger like equation is a valid description in the zero momentum limit. Once we go away from that limit, we will have more equations for the fluctuations with different choices of potentials. This is a more complicated scenario and fortunately our present analysis does not call for that. Nevertheless, the potential (6.23) is still highly non-trivial, as both P (Z) and Q(Z) take non-trivial values when expressed in terms of the background metric and dilaton in (6.21). For example, the function P (Z) that may be written as: where C is a constant that appeared first time in (6.5). The other parameters that appear above are G i and P i . All the G i 's depend only on the fixed parameters of the theory, and are defined by: where a is the resolution parameter (5.13) and α θ i are defined in (5.7). The other variables appearing in (6.25) are the P i 's out of which only P 1 is a constant. They are defined in the following way: where we have laid out clearly the g s N f dependences of each of the coefficients. One may see that the g s N f independent terms appear only from P 2 and P 3 . In a similar vein, we can also work out the Q(Z) piece in (6.21). This is given by: where P 2 and G 1 have already been defined in (6.26) and (6.25) respectively, but G 4 and G 5 are new. They can be related to, say, G 3 in the following way: With these set of definitions, the functional forms for P (Z) and Q(Z) are fully determined, although there is one issue that one may want to clarify at this point. This has to do with the presence of terms with relative minus sign inside the square root in (6.24). To avoid (6.24) to develop complex values, we require: where we have used the fact that g s N f P 2 ≥ 4π in the limit g s → 0 (see also (6.26)). This seems to constrain short distances, but since r > r h (6.29) do not put strong constaints. In fact we can take small Z, large N and vanishing momentum limits to re-express the potential (6.23) in the following way: The way we have expressed the above potential, one may clearly see how the various terms in the sum are increasingly suppressed by g s N f and gsM 2 N . The constant β appearing above is related to c i in (5.13) as β = c 1 = c 2 for simplicity 31 ; and m 0 ++ is the mass of the lightest glueball given via: and parametrized by the scale m 0 . This is computed using M-theory metric perturbation, much like the analysis we had in section 5, and is further detailed in [101]. We have also used (6.31) to define d o as: which is a constant. Note that in (6.31) the first term in independent of ω 2 and only depends on Z, b 2 and the glueball mass. The glueball mass also features in the definition of α(Z), that appears in (6.30), in the following way: where b is the bare resolution parameter that is defined in (5.13). Comparing the definition of the glueball mass in (6.31), we see that α(Z) is proportional to g s N , but suppressed by 1 Z 2 . This clearly indicates that the potential (6.30) goes as 1 Z 2 for small Z.
Note that Z = 0 (horizon) is a regular singular point of (6.22) and the exponents of the indicial equation near Z = 0 can then be written as 1 2 ± iI, where I is defined as: The functional form for I shows that it is suppressed by both g s N f and gsM 2 N so to zeroth order there is only a piece that depends on the bare resolution parameter b, the frequency ω, the horizon radius r h and the 't Hooft coupling g s N . We can also express the solution for the Schrödinger type equation (6.22) using I as: (6.35) with F T (Z) being a function that is analytic at the horizon radius r h (or at Z = 0). The above equation tells us the precise behavior of the eigenfunction E T (Z) at Z = 0. This is useful but not exactly relevant for the present case, as what we actually need is the form for E T (Z) when Z 1. The question then is how will the Z = 0 analysis be useful for the large Z domain.
The answer lies in our choice of the ansatze (6.35) that in fact serves as a good ansatze even when Z 1. In other words, the exponent of the indicial equation I that we computed in (6.34) still remains a valid solution for large Z. What changes for large Z is the functional form for F T (Z).
Of course there is yet another change in (6.22): it is the functional form for the potential V T (Z) that we computed earlier in (6.30). Naturally since (6.30) was for small Z, this should change. The change is easy to work out, and may be written in the following way: where we are ignoring higher powers of e −2Z that would appear from the relevant higher powers of e −2Z in P (Z) and Q(Z) in (6.24) and (6.27) respectively. The A and B appearing in (6.36) are not constants, with A defined as: where P 2 is given in (6.26). The function P 2 is defined with N, r h and Z, and one may take appropriate limits in terms of either of these parameters. Before we do this, let us write the expression for B in terms of the background parameters: where the successive suppressions with respect to gsM 2 N as well as g s N f are shown. The term independent of all these is proportional to b/β o where b is the resolution parameter and β o ≡ β log(er h ) with β = c 1 = c 2 in (5.13). The other parameters appearing in (6.38) are defined in the following way: One can now take the form of the potential, given in (6.36), and the wave-function ansatze, given in (6.35), and plug them in the Schrödinger-type equation (6.22) to obtain the following equation for F T (Z): where I is still given by (6.34). The above second order differential equation is rather hard to solve because of the presence of the exponential term e −2Z . However, since we seek the spectral function only in the limit of large Z where e −2Z vanishes, we can easily remove the problematic term from our equation (6.40). Doing this yields the following form for E T (Z): where C + and C − are two integration constants whose values will be determined later. To extract the actual fluctuation E T (Z) from (6.41), we need the functional form for P (Z) in the large Z limit. This is easy to extract from (6.24) and may be written as: which is as expected proportional to g s N f , and P 2 is defined in (6.26). The other coefficient G 6 appearing above can be extracted from some combinations of G i and P i in (6.25) and (6.26) respectively at large N . Here we write this simply as: where n o is a numerical constant given by n o = 4 √ 23 1/3 π 3/2 ≈ 45.43. Combining (6.41), (6.42) and (6.43) together and looking at the definition of E T (Z) given just after (6.22), we can finally determine the form of the fluctuation at large Z as: Few comments are in order related to the form of (6.44). First we see that the suppression factor is (g s N f ) −1/2 . From here it would seem like this does not have a natural zero flavor i.e N f = 0 limit. However when combined with P 2 , g s N f P 2 does have a zero flavor limit, and is given by: (6.45) which one may also verify directly at the level of the Schrodinger equation (6.22). Secondly, the integration constants C ± appearing in (6.44), can in principle be complex valued. So this will require us to investigate few possibilities associated with various choices of C ± satisfying the boundary conditions. Let us start by investigating the form for A given in (6.37). First let us assume that Z goes to infinity as: The above would make sense because N → ∞ and g s N f → 0. In this limit P 2 may be replaced by −3 log r h . In other words A in (6.37) becomes: for large log r h so that the inverse suppression in (6.47) makes sense. Assuming this is possible, plugging (6.47) in (6.44) would imply the following form for fluctuation E T (Z): where we have suppressed inverse log 2 r h dependences. Note that the functional form for E T (Z) is not the only way to express E T (Z) from (6.44). For example if the horizon radius goes as: (6.49) in the limit of very large Z uv and vanishing g s N f , then one may rewrite P 2 simply in terms of log N and not log r h . This means A in (6.37) in-turn will be expressed in terms of log N and not log r h , implying: From the multiple ways of expressing (6.44), for example (6.48) and (6.50), one might worry that the final result would be dependent on our approximation scheme. However we will show in section 6.3 that this will not be the case. Finally, the functional form for E T (Z) in the zero momentum limit matches with the functional form for E x 1 as may be seen from (6.18). This will help us to express the on-shell action completely in terms of known parameters appearing in (6.44), allowing us to compute the spectral function more efficiently. This is the topic that we turn to next in the following section.

On-shell action and the strong coupling spectral function
In the previous section we managed to find the functional form for the gauge field fluctuation E T (Z) in the large Z and in the zero momentum limit. What we now want is the four-dimensional on-shell action. This can be easily extracted from the boundary piece of the Lagrangian (6.11). Earlier we had used (6.11) to determine the EOM for E T (Z) and subsequently for E T (Z). Plugging in the EOM in (6.11) then leaves us only with the boundary term, that we shall label as the on-shell four-dimensional action S 4 . This takes the following form: where x 0 ≡ t and Ω 2 is the same volume of the two-sphere that we had in (6.3). Note that we took Z h to be the lower limit of Z to be consistent with the lower bound (6.29) 32 . However what we seek here is in fact the on-shell action at the boundary Z = Z uv , so the near-horizon geometry is not too relevant for us. At the boundary F tZ = −F Zt = 0, so we must set G tZ = 0 and replace √ −G by √ −G. Incorporating these changes, the boundary value of the on-shell action is now given as: Using the gauge field EOM (6.14), but now resorting to the metric G µν instead of G µν , and the result in Appendix B, the above action can be rewritten in terms of the gauge invariant variables E x 1 , E x 2 and E x 3 as: with k 2 a given in (B.9), and at this point we will be concerned about the x 1 part of the fluctuation. In other words, we only want to study the behavior of E x 1 at zero momentum. At zero momentum, according to (6.18), the fluctuations E x 1 and E T follow the same equation (6.19). Using such an identification, we can define: where one may match the Lorentz indices using (6.18). Plugging (6.54) in (6.53) and using E 0 (k)E 0 (−k) = 1, it is easy to see that the zero momentum limit yields the following action for the x 1 piece of the fluctuation: Before moving ahead, let us make couple of observations. One, E T (Z) is exactly the fluctuation (6.44) that we derived earlier and is therefore subjected to take either of the two possible limits (6.48) and (6.50) that we mentioned above. Two, the coefficient of E T (Z)/E T (Z) looks very similar to P (Z) in (6.21), so one might think that it can take the functional form (6.24). This is unfortunately not the case because P (Z) in (6.21) and (6.24) involves det G ab whereas the coefficient of E T (Z)/E T (Z) in (6.55) involves det G ab . The former differs from the latter by the presence of F ab . The above discussion more or less sets out the tone for the rest of the computations. There are two parts to the computation that we will indulge in the following. One is the coefficient of E T (Z)/E T (Z) in (6.55) and the other is E T (Z)/E T (Z) itself. To condense some of the subsequent formulae, let us define: where P 2 is given in (6.26). The coefficient of E T (Z)/E T (Z) can then be represented in the following way: 57) where we are suppressing higher orders 1/N terms, and n o is a numerical constant that appeared in (6.43). Note that both the denominators are suppressed differently with respect to N, α θ 1 and e 2Z . The numerators are non-trivial functions of e 2Z , and they will govern the behavior of the spectral function. Let us therefore study them carefully by first writing out the form for Σ 11 : where we see that it is proportional to g s N f . This makes sense because in the absence of the flavor D6-branes we won't see this contribution. The forms for P 1 and P 2 are given earlier in (6.26), where P 2 is a function of Z and r h but P 1 is independent of both of them. At this stage we can make (6.58) vanish by choosing: Few questions immediately arise from (6.59). What is the logic behind the choice (6.59), instead of making the other bracketed terms in (6.58), to vanish? What would happen if we make the other bracketed terms in (6.58) vanish? The answers to both the questions lie on the following observation: since b 2 as well as α θ i pieces cannot be large, the first three brackets in (6.58) cannot vanish. Making them zero would lead to contradictions. Therefore from (6.46) we see that Z can be very large, and we can use this to fix the value of r h using (6.59). This gives us: where is a small number that can be derived from above. One may also see that (6.60) cannot be related to (6.49). This is because of our choice between (6.46) and (6.49): we are allowed either of them, but not both. Coincidentally, we can choose either (6.48) or (6.50), but not both. As one may easily verify that the choice (6.46) only, and therefore (6.48), can be consistent with (6.59). The caveat however is that, since log r h is no longer a large number, the expansion in (6.47) cannot be terminated and we shall require the exact form for A in (6.47). We will discuss a way out of this soon. After the dust settles, there will be no Σ 11 term, and so we have to go to the next term given by Σ 22 . The next term incorporates both g s N f as well as gsM 2 N , and takes the following form: where we see that the term is dependent on r 2 h as well as various other factors of log r h . There are also e Z and N dependences that will take large values, so we will have to be careful taking the limits at large Z and large N . The various other quantities appearing above are α b , K(Z) and L(Z) that will be defined in the following. First let us start with α b : where on the right we have shown its behavior at large Z: the resolution parameter b 2 being small, does not contribute anything to α b . In the same vein, K(Z) is defined in the following way: where we have defined P 7 in (6.56) above. Using this definition for P 7 , and (6.46) for Z that we took earlier, one can easily show that: leading to some simplification in (6.63). It also means that for large Z, K(Z) goes as −3e 6Z . This is consistent with the other coefficient for α 4 θ 1 as evident from (6.61). Finally, the last term L(Z) takes the following form: where the large Z behavior is solely governed by the vanishing of P 7 in (6.64). In fact plugging the limiting values of (6.62), (6.63) and (6.65) in (6.61) and then in (6.57), leads us to the following behavior of (6.57) for large Z: (6.66) with κ o being a constant that depends on N as κ o ≡ bβN 0.1 no , where n o remains the same numerical constant that appeared in (6.43).
Before moving ahead, let us pause briefly to examine the situation at hand. The crucial outcome of (6.66) is the dominance of Σ 22 over Σ 11 because of the imposed constraint (6.59). This further lead to the form of the horizon radius given in (6.60) that is of order 1. This in turn gives us the high temperature limit, and so one might ask if there is a way to analyze the spectral function for small r h . Otherwise an expression like A in (6.47) does not have a good expansion in terms of inverse log r h . The situation is subtle because we would still have to impose (6.59) to eliminate Σ 11 piece in (6.57). How can we then avoid the outcome (6.60) for the horizon radius?
A way out of this conundrum is to not impose (6.46) that determines Z from the start, and instead use (6.59) to fix Z. This means (6.47) for A does not hold anymore although the form of A in (6.37) continues to hold. Z then satisfies: which is extracted from (6.59). The RHS has log r h and, as discussed above, we cannot use either (6.60) or (6.49) for r h . Instead we will use a different way, as shown in Appendix C, to determine the horizon radius by demanding the vanishing of the effective number of the three-brane charges in the original type IIB side. Solving (6.67) then gives us the following value for Z: where W n is the analytic continuation of the product log function with integer n. By construction this is a large positive number because N is large whereas r h is a very small number. Plugging (6.68) in (6.37) then gives us the following value for A: which is expectedly different from (6.47). The form of A shows that it is in fact a very large number because in addition to it being inversely proportional to a small number, i.e r h << 1 as mentioned above, it is also exponentially dependent on a large number as g s N f → 0. This will be useful for us because large A can simplify the expression for E T (Z) in (6.44). We will come back to this soon. Let us now compute the coefficient (6.57), which in turn means computing Σ 11 and Σ 22 . As mentioned earlier, Σ 11 vanishes, so we only need to compute Σ 22 at large Z. For this we will need the limiting values for α b , K(Z) and L(Z) in (6.62), (6.63) and (6.65) respectively. The limiting value for α b remains e 6Z as before, but the limiting values for K(Z) and L(Z) change because we can no longer apply (6.64) anymore. They now take the following values: where P 7 is given in (6.56). Plugging (6.70) in (6.61) and using (6.59), gives us the following value for the coefficient in (6.57): 71) where κ 1 = κo g 3/2 s and κ o is the same constant that appeared earlier in (6.66), and in the last line we have used the large Z limit (6.68) to eliminate the e −2Z piece. The above result differs clearly from (6.66), which was computed for r h as in (6.60). Here we expect r h to be small − as shown in Appendix C − and so (6.71) will finally be proportional to r 2 h log r h .
Having done the first part of the computation in (6.55), let us now investigate the second part which is the ratio E T (Z)/E T (Z). The functional form for E T (Z) is given in (6.44) and is expressed in terms of coefficients C ± which could in principle be complex. The ratio then can be written as: where we have introduced three functions α, g and Q that are in general complex.
In fact what we require is that the α, g and Q functions remain complex for large Z and small r h . The precise forms for these functions are: where I and A are defined in (6.34) and (6.37) respectively. Note that in the limit of large N , small g s N f , and small r h , A is large number, implying a small (but non-zero) complex piece in g. On the other hand, I is a large number being proportional to g s N and inversely proportional to the horizon radius r h . However to avoid contradictions, we will not take any limit at this stage and continue with the operations with exact expressions. This gives us: where now there are three distinct sources of imaginary pieces from (6.74): they can come from the C ± coefficients, the exponential term e iZ/ √ A and the bracketed terms in (6.74). The bracketed terms are defined with respect to two new functions P o and Q o , which may be written as: (6.75) The limit that we are looking for now, and as mentioned earlier, is the large Z, large N and small r h limit where Z becomes large as (6.68). Essentially then it is the large N and large |log r h | limit. In this limit P 2 can be expressed using Z as (6.59), which would tell us that it is a small number 33 . Plugging the values of A from (6.69), and Z from (6.68) now implies that Q o may be approximated by: where on the RHS we have shown the behavior of the function as it approaches zero, by ignoring a constant additive factor as the term in the bracket on the LHS of (6.76) will always dominate. We have also defined x and then Z as a function of x in the following way: where the latter should be viewed as an alternative expression for (6.68). For large N , small r h and g s N f → 0, it is easy to see that x vanishes whereas Z becomes very large. However Q o will always go to zero in this limit. What we now want to claim, in this limit, is that: which is easy to justify from the form on Z in (6.77) and the fact that multiplying x with any powers of log x will always approach zero in the above limit. The dominance of I/Z over Q o is a huge simplification for us because this will not only render the expression (6.74) manageable without worrying about contributions from the exponential pieces, but also remove the ambiguity of its dependence on the constants C ± whose values have not been explicitly determined. In fact after plugging in all the values from (6.75) and (6.34) in (6.74) and using the limiting conditions (6.76) and (6.78), it is easy to see that: where m 0 ++ is the mass of the lightest glueball expressed in terms of scale m 0 and is given in (6.31). In fact this is all we need, because the imaginary part of (6.79) can then take the following form: where we have used (6.31) to write it in this form. One may note its linear dependence on ω, the frequency parameter that we encountered earlier. It is also inversely proportional to the horizon radius r h , a fact that will be useful soon. The logic behind the above series of computations should be clear now. What we are looking for is the retarded Green's function in the zero momentum limit. This is now easy to extract from (6.55), and can be written as: which precisely contains the two pieces of computations that we performed above, namely the coefficient of E T /E T in (6.71) and the ratio E T /E T itself in (6.79). One additional input was the imaginary piece in (6.79) that we extracted in (6.80). The reason for this extra bit of work is apparent: the spectral function is exactly the imaginary piece of the retarded Green's function, i.e: where T is the temperature that will be related to the horizon radius r h . Since (6.71) is all real, the imaginary piece in the retarded Green's function can only come from (6.79). Any other contributions to the imaginary piece will be suppressed by higher powers of 1/N so does not concern us here. Putting everything together then gives us the required expression for the spectral function: where expectedly this is proportional to g s N f and g s M 2 /N . It is also proportional to r h (and also log r h ), so at zero temperature ρ(0, ω) = 0. We can use (5.13), or the footnote 31, to express the combination br h in terms of the resolution parameter as . This way, the pre-factor multiplying log r h in (6.83) is not explicitly but only implicitly dependent on r h and it brings out the resolution in the gravity dual rather succinctly. The two other functions appearing in (6.83) may be defined in the following way: where n o is a numerical constant defined after (6.43), and β is defined in footnote 31. Note that if we use the strong string coupling result, as opposed to the weak string coupling analysis presented here (both a strong 't Hooft coupling of course), β can be defined from (5.13) with β = c 1 = c 2 . The coefficient c 1 appears in (5.52) and c 2 is bounded by (5.57). Following this logic, what we now need is the g s N f independent pieces to define β. Thus if we take the negative definite constant piece of c 1 from (5.52) and use this to define both c 2 and β then we can ignore higher order g s N f dependences. Thus essentially, from both strong and weak type IIA couplings, β will be another constant to O(g s N f ), which in turn would make F b to be another constant 34 , that we shall call f b . However the worrisome feature is the other function in (6.84), i.e the function F a that depends on N, g s and Z uv . Both N and Z uv , with Z uv defined in (6.68), go to infinity whereas g s approaches zero. If we define ζ 1 ≡ g s , ζ 2 ≡ 1/N and ζ 3 ≡ 1/Z uv , then we can choose the behavior of each of these parameters such that: 85) 34 Recall that the parameters α θ1 and α θ2 are constants.
with a constant f a . 35 . As T → 0, r h vanishes, and from the expression of the spectral function in (6.83), this also vanishes. Therefore we can finally put everything together and argue that: where we have used (5.58) to express the RHS in terms of the sound speed. Of course, as mentioned above, (5.58) is a strong coupling result, so the comparison has to be done with c 1 and c 2 being proportional to g s N f and not constants (as opposed to the weak string but strong 't Hooft coupling answer). Taking all these into consideration we see a clear linear dependence on ( 1 3 − c 2 s ) at strong 't Hooft coupling, perfectly consistent with the results of sections 4 and 5.
Few comments are in order now. Our analysis is based on small r h as derived in Appendix C, so the natural question is what happens when r h is of order 1, i.e the one given in (6.60). When the horizon radius is of order 1, it means we are at the point where new degrees of freedom are about to enter, i.e we are in Region 2 of [58]. Therefore unless we know the detailed metric configuration of Region 2 and beyond, we cannot perform the analysis as clearly as we have done here because of our definition the radial coordinate as r = r h e Z . When r h is small we are still in Region 1 of [58] and so precise computations may be performed (as shown here).
Secondly r h itself is bounded below by (6.29). This bound is of course to prevent any appearances of unphysical imaginary pieces in the computations. Clearly for the range of Z that we are concerned here, this poses no constraints. Thus happily all the results lead to the following conclusion: 6.4 The strong string coupling limit and pure classical supergravity Most of the analysis section 6 is done with g s → 0 and with large M . This differs a bit from what we did in section 5 where g s = O(1), so that natural question is whether we can work through the analysis of sections 6.1 − 6.3 assuming (g s , N f ) ∼ O(1) and N 1 as part of the MQGP Limit of [63] 36 . This is an unusual large N limit but still warrants the use of pure classical supergravity. To see this, one 35 In the MQGP limit wherein g s < ∼ 1, one can argue that f a will be a finite non-zero constant as follows. As r h < r 0 or |logr h | > |logr 0 | (r 0 being the r where the D3-branes have been entirely cascaded away, and noting min(r) = r h ), hence instead of choosing r h to satisfy (C.9), assume |logr h | = N 1/3 κ 1 f , 0 < f < 1 and κ = n b g s M 2 /3 from (C.9). As Z U V ∼ |logr h | + logN 1/3 ∼ notes that by including terms higher order in g s N f in the RR and NS-NS threeform fluxes than those considered in [63] and the NLO terms in the angular part of the metric, one sees that in the IR in the MQGP limit, there occurs an IR colorflavor enhancement of the length scale as compared to a Planckian length scale in the Klebanov-Strassler (KS) model [56] for large M , thereby showing that stringy corrections will be suppressed. To see this more explicitly, we summarize the main ideas of [94,103] here. Using [58] let us define an effective number of color in the following way: where M eff and N eff f are the effective number of bi-fundamental and fundamental flavors respectively that are defined for our background in the following way: where (m, n) indices are summed from (m, n) = (0, 0) onwards, and henceforth to avoid clutter we will use Einstein summation convention. The coefficients k mn ≡ k mn (r, g s ) and f mn ≡ f mn (r, g s ) and therefore the effective flavors are constructed from the higher orders g s N f and gsM 2 N corrections [58]. Combining these together, it was argued in [94,103] that the length scale in the IR at r = Λ will be dominated by: In the IR, relative to KS geometry, we thus see that (6.90) implies the abovementioned color-flavor enhancement of the length scale. Therefore in the IR, even for g s = 0.45, M = 3 and N f = 2, upon inclusion of of n, m > 1 terms in M eff and N eff f in (6.89), the characteristic length scale in the MQGP limit [63] involving g s ≤ 1 satisfy: where L KS is the characteristic length scale for the Klebanov-Strassler model [56] in the far IR. Because of this enhancement, the stringy corrections are suppressed implying that one can still trust classical supergravity. It is however interesting to note that in the IR, one can obtain g 2 Y M = O(1) even for g s → 0, provided N f = 0). To see this let us first consider vanishing N f . The NSVZ RG flow equation for the SU(M ) gauge group that survives at the end of the Seiberg duality cascade, gives us: (6.92) where the RHS appears from the integral of the NS two-form field over a vanishing two cycle S 2 in the type IIB side. This is of course the same two-cycle discussed at the beginning of section 4, parametrized by (θ 2 , φ 2 ), on which we have M wrapped D5-branes. The question is whether (6.92) can allow g 2 Y M = O(1). Solving the equation (6.92) gives the inverse YM coupling in terms of M and log r, for r = Λ. It is easy to see that, with M = O(1) this is only possible if Λ is proportional to the UV cutoff itself. Clearly since we want to concentrate on far IR physics, such a choice is not feasible. Additionally, since near the UV cutoff we expect the theory to become scale invariant, M automatically vanishes there.
On the other hand, when N f = 0, the above conclusion can change because the dilaton on the gravity side is no longer a constant. Recall that, with N f flavors, the dilaton takes the following form [57,58]: where a 2 is the resolution parameter that we encountered earlier. Using the fact that we have an almost vanishing resolution parameter, and the angular coordinates (θ 1 , θ 2 ) are parametrized by (5.7), the inverse of the YM coupling now satisfy: at the scale r = Λ measured with respect to the cutoff scale Λ ∞ . What we are looking for now is a Λ in the IR whereat g 2 SU (M ) = O(1). The scenario is more subtle now because of the additional O(g s N f ) pieces appearing in (6.94). These pieces come from carefully looking at the NS B-field threading the vanishing two-sphere on which we have the wrapped D5-branes. The B-field is more non-trivial than what we had above, and is given by: where the first term is precisely what we had on the RHS of (6.92) for the case with vanishing N f , and the second term involves the g s N f corrections. These correction terms have been worked out in [58], and may be expressed as: where we have removed any dependence on the resolution parameter when writing (6.96) from [58] 37 . In fact the O(g s N f ) term alluded to in (6.94) comes precisely from Q in (6.96). There is however one subtlety associated with the angular variables θ i and φ 2 . Since the integral of the B 2 field over the two-sphere parametrized by (θ 2 , φ 2 ) contributes to the YM coupling g 2 Y M , one needs to be careful while imposing (5.7). One way would be to impose (5.7) to θ 1 in (6.95) and then integrate over θ 2 . In that case an additional N dependence would appear from the second term in (6.96). Alternatively we could also insert the value of θ 2 from (5.7) after integration over the two-sphere. The latter would imply that the integration of B 2 field over the two-sphere is concentrated mostly near the regime defined in (5.7). After the dust settles, the equation that we need to solve to determine Λ can be derived from (6.94) as: which is a quadratic equation to the first order in g s N f . To higher orders in g s N f the equation starts becoming more complicated. The various coefficients of (6.97) are defined as 38 : Let us pause a bit to see what are the dominating terms in the above set of coefficients. We want g s → 0, and small N f , but we also want very large N . Let us therefore take the following limiting values for g s , N f and N : There is one subtlety that we are putting under the rug. A part of the B-field in (6.95) goes as 9 4π (g s M ) (g s N f ) log r log |a|, where a is the resolution parameter. This blows up in the limit a → 0, so one might be worried that Q given in (6.96) is not well defined in this limit. This is however not the case because the derivation of the B-field in [58] was done with non-zero resolution parameter, and for zero resolution parameter we have to do the analysis separately. The result then is of course independent of the log |a| piece, and is as given in (6.96). 38 We have used the following values of the integrals governing the B 2 field using the θ i values given in (5.7): where α N could be a large number and 1 < b < 2. This clearly shows that the g s N f log N term in B dominates and B 2 4AC. Using this criteria, and solving (6.97) immediately gives us: implying that Λ can be in the deep IR. Hence, one can obtain an O(1)g Y M in the IR without requiring an O(1) g s , in the presence of flavors but not in their absence. In the IR, of course N f = 0. Before ending this section, let us make a few observations. First, if we also take M to be very large, then the first term of C in (6.98) will be suppressed by 1/M . This of course will not change the conclusion of (6.100). Secondly, in Section 6, (6.29) will be replaced by the observation that for large Z the argument of the square root in (6.24) is obviously positive and for small Z: where P 2 is given in (6.26); as long as logN, |logr h | 1. This is obviously true from our earlier considerations. Therefore, the argument of the square root in P (Z) in (6.24) is always positive.

Conclusions and discussions
In this work we have studied bulk viscosity in the whole range of the 't Hooft coupling constant λ. One of the main goal of our studies was to express the bulk viscosity as a function of the speed of sound within well-established first-principle theories. Our efforts were put into clarifying possible differences in the parametric form of the ratio ζ/η, obtained at different coupling constant. Apart from the discussion on the final forms of the bulk viscosity, we have also elaborated on the relevant employed analytic methods. This was to adapt them to the examined regimes and to justify their relevance. We focused on the extreme limits of the coupling where analytical methods are applicable. At weak coupling, kinetic theory was used, which is currently the most common effective approach to calculate the bulk viscosity and other transport coefficients. To confirm the validity of kinetic theory we provided its justification from a more fundamental diagrammatic approach. At strong ('t Hooft) coupling, the UV complete type IIB holographic dual (and its M theory uplift when addressing also the strong string coupling limit) of large-N thermal QCD was employed. The intermediate coupling behavior, most relevant for the quark gluon plasma produced experimentally in the heavy ion collisions, was also briefly discussed. We mainly summarize known challenges related to the first-principles extraction of bulk viscosity.
To discuss the weak coupling limit we summarized and matched the analysis of bulk viscosity of QCD done extensively within the effective kinetic theory in Ref. [29] to the case when the interaction is governed by the 't Hooft coupling λ = g 2 Y M M . In such case, the bulk viscosity behavior is controlled by gluons only as the quark contributions are suppressed by a factor 1/M , where M → ∞ is the number of colors. 39 The parametric form of the bulk viscosity as a function of the speed of sound is ζ/s ∝ 1/3 − c 2 s , while the ratio ζ/η ∝ (1/3 − c 2 s ) 2 . Then, starting from the Kubo formula, we performed a multi-loop analysis which enabled us to determine which scattering processes contribute to the collision kernel of the Boltzmann equation and provided a power counting in the weak 't Hooft coupling and high temperature. Collecting all evaluated diagrams we have shown a schematic procedure how to derive an integral equation which may be thought of as a diagrammatic representation of the Boltzmann equation. The integral equation is formed by infinite number of planar diagrams with propagators and vertices being dressed. Both number conserving and number changing processes have to be included in the complete bulk viscosity examination. For the vertices a separate integral equation, governed mainly by the soft physics and capturing the LPM effect, has to be solved. The diagrammatic analysis presented in this work stands for the first explicit justification of the validity of the Boltzmann equation, whose solution is needed for transport coefficients investigation governed by SU(M ) theories.
Within the intermediate coupling region, we have summarized a state of knowledge on the bulk viscosity studies. Although the prescription of calculation of bulk viscosity is given by the Kubo formula, it is difficult to reliably establish the hydrodynamic limit of the spectral function and determine which physical phenomena may be responsible for its shape. Therefore we can only conclude that all compiled findings do not allow one for quantitative determination of the bulk viscosity behavior in this region starting from first principles and new methods and/or perspectives are needed.
After analyzing the weak and the intermediate 't Hooft coupling regimes, we go to the next stage, i.e the strong 't Hooft coupling regime. Clearly neither pQCD, nor Lattice results can help us here. A new paradigm is needed and is given by the so-called gauge/gravity duality. This is a refined form of the famed AdS/CFT duality, constructed precisely to tackle strongly coupled gauge theories that are nonconformal. In section 4 we study a SU(M ) gauge theory in the IR at high temperature (i.e the temperature above the deconfinement temperature) and at strong 't Hooft coupling. We take large M , but keep the string coupling g s very small, such that λ = g s M is still very large. To avoid additional complications, we take no flavor degress of freedom. 39 The number of colors is N + M in the UV and M in the IR; both are kept very large in sections 2.2 and 4 and N f (along with string coupling) could be taken to be O(1) in sections 5 and 6 keeping N to be very large as part of the "MQGP" limit.
In such a setup, the computation of bulk viscosity boils down to the computation of metric fluctuations in the corresponding gravity dual. In section 4.2 we study the equations governing the fluctuations using two steps: one, in section 4.2.1, we relax some of the constraints and study a toy example which in turn provides a nice solvable system; and two, in section 4.2.2, we do a more precise and careful computations of the fluctuation equations. Knowing the precise fluctuations help us to compute both the sound speed as well as the ratio of the bulk to the shear viscosities. In section 4.3 we perform the aforementioned computations and show that the ratio of the bulk to shear viscosities is indeed bounded below by the deviation of the square of the sound speed from its conformal value.
It is believed that QGP is an example of a strongly coupled system at finite temperature wherein unlike as considered in most gravity duals, the gauge coupling and hence the string coupling, is of O(1). Motivated by the same and with the idea of also including the flavor degrees of freedom as well as the UV region, in section 5, we calculate holographically at finite string coupling, the deviation of the square of the speed of sound from its conformal value, the attenuation constant and the ratio of the bulk and shear viscosities and find a Buchel-like bound for the latter. Finite string coupling necessitates addressing these issues from the M-theory uplift of the type IIB construct of [58] which was obtained by the M-theory uplift of the SYZ type IIA dual in [63]. This also enjoys the additional benefit of not having to keep track of the NS5-degrees of freedom that one needs to while working with a single T-dual of the type IIB configuration of [58]. Based on [100,103], an equation of motion (EOM) for a combination of scalar modes of metric perturbations invariant under infinitesimal diffeomorphisms, is constructed. Upon investigating this EOM near the horizon, it is realized that for a non-zero bare resolution parameter, the horizon turns out to be an irregular singular point. Demanding the same of an ansatz for the solution to the same, in section 5.2, the dispersion relation for the quasinormal modes obtains not only the conformal values of the speed of sound and the attenuation constant but also their respective non-conformal corrections. Interestingly, for the case of a vanishing bare resolution parameter, by looking at the solution to the EOM near the asymptotic boundary, in section 5.3 one realizes that one can not consistently impose Dirichlet boundary condtion (at the asymptotic boundary); like section 4, non-normalizable modes are required to propagate. In section 5.4, with a non-zero bare resolution parameter, we first show that the KSS bound on the shearviscosity-to-entropy-density is not violated having incorporated the non-conformal corrections. We then obtain the bulk-viscosity-to-entropy-density ratio and the deviation of the square of the speed of sound from its conformal value, and confirm that the conformal value of both vanish and they are both hence determined entirely by the non-conformality of the theory. One of the main results of this section is a crisp bound: ζ η ≥ 91 5 1 3 − c 2 s with(out) the flavor degrees of freedom. In section 6, we approach the issue of obtaining the deviation of the square of the speed of sound from its conformal value, from two-point correlators involving gauge field fluctuations on the world-volume of flavor D6-branes using the prescription of [107]. To begin with one considers the weak-string-coupling strong-'tHooft-coupling limit. The fluctuations are considered over a background value of the gauge field − worked out in section 6.1 − assumed to be having only a temporal component and radial dependence. In the zero-momentum limit, interestingly and as shown in section 6.2, there is only a single second order equation in a gauge-invariant perturbation field − the 'electric field' − which needs to be and is solved for (in section 6.3). Finally the subtracted (zero temperature from the non-zero temperature) spectral function per unit frequency in the vanishing-frequency limit yields that the same is proportional to the linear power of the deviation of the square of the speed of sound from its conformal value, thereby validating the same as obtained in the previous sections 4 and 5. We conclude section 6 with some remarks (in section 6.4) arguing that this result remains unchanged even in the strong-string-coupling stong-'tHooft-coupling, or the true MQGP limit of [58].
Let us briefly discuss some future directions. It would be rather interesting to probe better the regime of intermediate 't Hooft coupling whereat the number of colors is large, the gauge coupling is small but the 't Hooft coupling is finite, i.e., neither small (weak coupling regime) nor large (strong coupling regime). As discussed, techniques based on QCD do not offer currently a reliable way to explore this region. The ansatz proposed for the spectral function parametrization does not capture properly a high frequency tail and the QCD sum rule cannot be directly applied to constrain bulk viscosity. Since it is not clear how to handle the issues with the QCD tools, the region can be alternatively explored within the supergravity framework. One could invoke higher derivative corrections in the supergravity action which would hence back-react on the background. The same in the context of N = 4 SYM has been studied recently in [108]. For the present case there are two ways to go about it. One, we could start from the type IIB background of [58] and consider corrections to the metric and fluxes in powers of α 3 and solve the modified equations of motion up to O(α 3 ). Two, we could use the MQGP limit (with g s ∼ O(1), large N but finite g s M ) and start with D = 11 supergravity action up to sextic power in the eleven-dimensional Planck length [109], and construct solutions to the EOMs as Planck-length perturbation of the M-theory uplift of [63]. Clearly the latter is a bit more practical because of the reduced number of fields in M-theory. Following this, one can then include metric perturbations and solve for their EOMs and hence see the effects of the inclusion of the aforementioned higher derivative terms on some spectral functions. It would be interesting to evaluate the non-zero frequency contribution to the spectral function per unit frequency and compare with previous studies on this topic as in [110] (which had excluded higher derivative corrections) in N = 4 SYM.
Another possible future direction could be to look at simultaneously turning on gauge and vector modes of metric perturbations [81] and then see the modification in the spectral function of gauge fluctuations considered in Section 6. The same in the context of type IIB for evaluating electrical and thermal conductivities, was considered in [103].

A.1 The equation of motion for the fluctuation mode H tt
To elaborate the implications of the above discussion, let us discuss the EOM for H tt . This can be expressed in terms of the other fluctuation modes in the following way: , (A.7) which at first glance seems to be well defined in the regimes r h > 0 and u ≥ 1.32. The precise regimes of interest however is not important for the kind of details that we are aiming for here. This will be illustrated later. Note also that A i are not constants but certain nested functions whose forms may be given in the following way: where the denominator of the form (a, b) is to be understood as being identified with the subscript bracket A (a,b) so that individual relations for A a and A b may be constructed. The nested function D 1 takes the following form: which is also expressed in terms of B(u) in a somewhat similar form as in (A.9) above. Together they would determine the coefficient A 1 in (A.7). The next coefficient A 2 is now determined in terms of D 3 and D 4 . The former is simple: 11) and expressed in terms of B(u), whereas the latter is more involved and may be expressed in the following way: where the B(u) independent terms appear, in our notation, as B (0) (u) and one may verify the uniqueness of the proposed form (A.13). This is also evident from the next coefficient, namely D 6 which may be determined from c 6nm as: We see that the structure is somewhat similar to (A.14), i.e the coefficient D 5 in the sense that we have B (0) , B (1) , B (2) and B (3) terms distributed in an identical way (although the precise c knm coefficients differ) as the derivations of these terms involve similar manipulations of the Einstein's equations. This is evident from the form of the next coefficients, namely D 8 which may be expressed as: which is structurally similar to D 6 in (A.15). On the other hand, the last two coefficients require a slightly different analysis and therefore we expect them to differ from the above D k coefficients. This becomes clear from the expression of D 9 which is written as: D 9 ≡ 6g s π u 4 − 15 u 4 − 1 q 2 + w 2 B (u)u 5 + 8g s π u 4 − 3 u 8 + 8u 4 − 3 q 2 + 4u 4 + 3 w 2 B(u) + u 4 − 3 − 2r h 2 u 8 − 12u 4 + 3 u 2 + 9g s π u 4 − 1 u 4 − 1 q 2 + w 2 B (u)u 2 + 8g s N π u 8 + 8u 4 − 3 q 2 + 4u 4 + 3 w 2 , (A. 18) that takes the form, although similar to (A.13), different from the other D k coefficients. The final coefficient, D 10 , may be presented in the following to illustrate the same point: This completes our analysis of the EOM (A.7) for H tt (u). Our next step is to analyze the EOM for H s (u) defined above in (A.6).

A.2 The equation of motion for the combined mode H s
The functional form for H s (u), as evident from (A.6), can be expressed as certain linear combination of H xx and H yy . As in (A.7), we can express the EOM for H s (u) in the following way: , (A.20) where we see that both the denominator and the numerator have the same set of factors as in the denominator and the numerator of (A.7). The only thing that would differ are the actual values of B k . The functional forms for B k may be expressed in terms of certain nested functions in the following way: where the F k functions may be compared with the D k functions in (A.8). In fact one may even express the functional forms for F k as a power series in u and derivatives of B(u), much like (A.13), but now with coefficients g knm instead of c knm . The coefficients g knm are independent of u, and may be determined easily as before by analyzing the corresponding Einstein's equations. For example finding g 1nm and g 2nm immediately reproduces the functional forms for F 1 and F 2 in the following way: which as expected takes the form (A.13). One may also easily see the pattern repeating for the next two coefficients, namely F 3 and F 4 , in the following way: We can now go to the other set of coefficients where we can see how we could relate to the D k coefficients studied above. A priori there shouldn't be any apparent connections, but the functional forms for F 5 and F 6 are similar to what we had earlier. For example: which should be somewhat reminiscent of (A.14). Similarly the functional form for F 6 , expressed here as: , we see that they are related via the following relation: Thus knowing F 1 would determine the functional forms for F 2 as well as F 6 . In fact one can show that F 1 or F 6 can also fix the functional forms for two other coefficients, namely F 8 and F 10 , in the following way: . (A.27) The remaining two coefficients, namely F 7 and F 9 , are however more complicated and are not anyway related to F 1 in a simple way. For example the functional form for F 7 may be expressed as: which of course follows the pattern similar to (A.13), but cannot be decomposed in terms of any of the above F k coefficients. A similar thing may also be said for the coefficient F 9 , written as: take values in (t, x, y), (A.30) may also be expressed in terms of H ab (u) and H tt (u). This pattern follows for the next component H yy (u) as: H yy (u) = − q (u 4 − 1) H tt (u) + 2qu 3 H tt (u) + wH tx (u) 2q (u 4 − 1) , (A. 31) implying that solutions may be found once we know the background values. Finally, combining the above set of equations with the defining equation for H s (u), namely (A.6), gives us a way to formulate the EOM for H xx (u) as: Basically this is all we need to construct gauge invariant perturbation modes. For us, following [100], a specific combination of the above set of perturbations is useful to quantify the required perturbation as: ) is a precise reproduction of this fact. Once m(u) is determined, l(u) can also be obtained by equating the coefficient of H tt from Z s (u) to the sum of coefficients of H tt from Z s (u) and Z s (u). In (5.21) we quoted the functional form for l(u) for u → 1. The generic form for l(u) is straightforward but technically challenging, and is therefore left as an exercise for the reader.
After the dust settles, one may verify that the EOM (A.34) is satisfied by the gauge-invariant choice of the perturbation Z s (u) in (A.33).

B. A derivation of the on-shell action and the Green's function
The four-dimensional action that we considered in (6.51) uses the pull-back metric G µν constructed out of the type IIA metric, the NS B-field and the world-volume gauge field background. When the gauge field fluctuation, whose Fourier component is written as A µ in (6.13), is also taken into account, the four-dimensional action takes the following form: where x 0 ≡ t, T D6 is the tension of the probe D6-brane, and Ω 2 is the volume of the two-sphere that we had in (6.3). The presence of Z derivative in the integrand, despite integrating out the Z variable, is from a total derivative term as may be inferred from (6.51). This also explains the two limits of Z in (B.1) Note that we took Z h to be the lower limit of Z to be consistent with the lower bound (6.29). However what we seek here is in fact the on-shell action and the Green's function at the boundary Z = Z uv , so the near-horizon geometry is not too relevant for us. At the boundary F tZ = −F Zt = 0, so we must set G tZ = 0 and replace √ −G by √ −G. Incorporating these changes, the boundary value of the on-shell action may now be re-written from (6.52) as: where we have suppressed the ω dependence, and will have to resort back to G µν component if we want to take Z h , i.e the lower limit of Z. Recall also that we have used EOM to get to the boundary action (B.2), so it makes sense to use the EOM further to simplify the above action. For example we can use (6.14) to rewrite G tt in the following way: for non-vanishing ω. Plugging (B.3) in (B.2) then gives us the following action: which is almost similar to the action (B.2) except with three major differences: one, the appearance of 1 ω as an overall factor; two, the sum over a now being from 1 to 3; and three, the appearance of three new variables E xa for a = 1, 2, 3. The new variables are defined in the following way: which are clearly borne out from (B.3) and explains the appearance of the 1 ω suppression of the full action. We could also use (B.5) to express A y and A z in terms of E x 2 and E x 3 respectively, but we won't do this right way. Instead let us use the first equation in (B.5) to write: where to get the second equality in the above we have used equation (B.3). To complete the picture we will need the ratio of the two metric components. Using the fact that r ≡ r h e Z , we can easily argue that G xx G tt = − 1 − e −4Z . Plugging this in (B.6) gives us: This is all we need, because the derivatives on the other components are straightforward replacements of E xa with a = 2, 3. Therefore combining (B.7) with (B.5) and plugging this in (B.4) gives us the final action: which is the action given earlier in (6.53). The k 2 a appearing in (6.53) are the poles in (B.8) and may be identified with the variables of (B.8) as: Since we are only interested in the x 1 part of the fluctuation, the values of k 2 2 and k 2 3 are not very useful for us. Of course one may perform a more generic study, but we will not do so here. For the simplest case, the next step would be to define (6.54) and then re-write the action as in (6.55). From here the story follows as depicted in section 6.3.

C. Effective number of three-brane charges with background three-forms and the horizon radius
The horizon radius that we computed in (6.60) was typically a O(1) number that was written as r h = 1 − 2 by demanding the vanishing of (6.59). The small parameter is defined as: with both b, the resolution parameter, and g s N f small. This choice of the horizon radius is not very useful for us because this would imply that r h is placed right at the point where new degrees of freedom would appear to UV complete the system. With the definition of the radial coordinate r as r = r h e Z , this means Z only measures geometry beyond r h , i.e the geometry of Regions 2 and 3. Question then is how to place r h deep inside Region 1 where the background is well known. However we cannot make r h arbitrarily small, as there exists a lower bound on r h given in (6.29).
If r

(o)
h denotes the lower bound, then: with C being an integration constant that appeared in (6.5), and we expect the horizon radius to satisfy r h > r h . Such a lower bound is necessary otherwise an expression like P (Z) in (6.24) will develop unphysical imaginary piece.
To find an appropriate r h it would be easier to do the analysis in the type IIB side instead of the mirror type IIA side. Such an analysis won't change the expression for r h as the mirror transformation a la SYZ [64] keeps the radial coordinate unchanged. To proceed then, let us define an effective number of three-brane charge as: where B 2 , F 3 and F 5 are given, for N f = 0 and in the Baryonic branch, in (4.1). The five-dimensional internal space M 5 , with coordinates (θ i , φ i , ψ), is basically the resolved warped-defomed conifold of (4.2), or its simplified avatar given in (4.5).
If we now collectively denote the lower limits of all the angular variables, namely (θ i , φ i , ψ) as R − and the upper limits of all the angular variables as R + ; and also use the fact that at fixed r, dr = 0, then the effective number of three brane charges take the following form: thus simplifying the expression (C.3) tremendously. Here N denotes the integral over F 5 , and is therefore related to the integer D3-branes in the dual gauge theory side at the Higgsing scale. The second term combined with N then denotes the effective number of cascading D3-brane charges at the scale r = r 0 . The functional forms for (a 0 , e 0 , f 0 , b 2 , c 2 , d 2 ) can be extracted from (C.5). Combined together leads to the following expression for N eff : 18a 2 (−1) k log r + r 2 108a 2 log r 2k + 1 + r (C.8) + 5 3a 2 (g s − 1) + r 2 (3g s N f log r + 2π)(9g s N f log r + 4π) 9a 2 g s N f log e 2 r 3 + 4πr 2 , = N 1 + 6πlog r (3g s N f log r + 2π) (9g s N f log r + 4π) where we have only kept terms linear in gsM 2 N , linear and quadratic in g s N f , and ignored higher order terms. Of course one may question the logic of suppressing a term linear in log N . Such a term typically comes with (g s N f ) 2 and with either a 2 or with higher powers of r = r 0 . Since we will be assuming r 0 << 1, we can safely ignore the log N piece. Note that the assumption of small r 0 is crucial here. This implies the domination of g s N f |log r 0 | over other constant pieces in (C.8). Implementing this 40 , and putting N eff = 0, gives the following estimate for r 0 that we will identify with the horizon radius r h as: where n b ≡ 3 (6π) 1/3 . Since both g s N f and gsM 2 N are very small quantities, the horizon radius is indeed deep inside Region 1. Note that this estimate has to be bigger than the lower bound r (o) h which in turn has a range (C.2). 40 Otherwise one will have to solve a cubic equation in log r 0 from (C.8). This will have one real solution that we can identify with the horizon radius r h .