The Chern-Simons Diffusion Rate in Improved Holographic QCD

In (3+1)-dimensional SU(Nc) Yang-Mills (YM) theory, the Chern-Simons diffusion rate, Gamma_{CS}, is determined by the zero-momentum, zero-frequency limit of the retarded two-point function of the CP-odd operator tr[F ^ F], with F the YM field strength. The Chern-Simons diffusion rate is a crucial ingredient for many CP-odd phenomena, including the chiral magnetic effect in the quark-gluon plasma. We compute Gamma_{CS} in the high-temperature, deconfined phase of Improved Holographic QCD, a refined holographic model for large-Nc YM theory. Our result for Gamma_{CS}/(sT), where s is entropy density and T is temperature, varies slowly at high T and increases monotonically as T approaches the transition temperature from above. We also study the retarded two-point function of tr[F ^ F] with non-zero frequency and momentum. Our results suggest that the CP-odd phenomena that may potentially occur in heavy ion collisions could be controlled by an excitation with energy on the order of the lightest axial glueball mass.


Introduction
(3+1)-dimensional SU(N c ) Yang-Mills (YM) theory has an infinite number of degenerate classical vacua distinguished by a topological invariant, the Chern-Simons number, N CS . Normalizing the YM kinetic term as − 1 4g 2 tr[F µν F µν ], N CS is where i, j, k = 1, 2, 3 and the trace is over gauge indices. A change in the Chern-Simons number is thus where x µ = (t, x). In a state invariant under translations in space and time, the rate of change of N CS per unit volume V per unit time t is called the Chern-Simons diffusion rate, denoted Γ CS , where the subscript W denotes the Wightman function. In an equilibrium state with non-zero temperature T , let G R (ω, k) denote the retarded Green's function of q(x µ ) in Fourier space, with frequency ω and spatial momentum k. In such states, eq. (1. 3) can be rewritten as 2T ω Im G R (ω, k = 0). (1.4) Gauge field configurations for which d 4 x q(x µ ) is non-zero produce a non-zero ∆N CS . At zero temperature, such gauge field configurations, called instantons, represent quantum tunneling events between vacua. At both zero and non-zero T , the contribution of instantons to Γ CS is exponentially suppressed [1,2]. When T is nonzero, however, classical thermal fluctuations can also produce a non-zero ∆N CS , for example by exciting unstable gauge field configurations called sphalerons [3,4] which generate non-zero ∆N CS upon decay. Such classical thermal processes are not exponentially suppressed [5][6][7]: in YM perturbation theory Γ CS ∝ λ 5 t log(λ t ) T 4 , where λ t ≡ g 2 N c is the 't Hooft coupling [8][9][10][11].
In YM coupled to fundamental-representation fermions, d 4 x q(x µ ) also contributes to chiral anomalies in global symmetries. In the electroweak theory, gauge field configurations with non-zero d 4 x q(x µ ) play a role in electroweak baryogenesis [12,13], while in Quantum Chromodynamics (QCD), for sufficiently high T they may play a role in generating bubbles of net chirality (more left-handed than righthanded quarks, for example), in which parity, P, and charge conjugation times parity, CP, are broken [14]. 1 Such CP-odd domains in hot QCD may have observable consequences in heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). These collisions produce a hot soup of QCD matter with T on the order of two to four times the QCD crossover temperature. The resulting state behaves as a nearly-ideal fluid of strongly-interacting quarks and gluons, the quark-gluon plasma (QGP) [15][16][17]. A non-central collision may produce a QGP with non-zero angular momentum and hence a magnetic field, both pointing perpendicular to the reaction plane (spanned by the beam axis and impact parameter). In the presence of a magnetic field, a net chirality will produce an electric current parallel to the magnetic field, due to the axial anomaly. This is the Chiral Magnetic Effect (CME) [18,19]. A detection of the CME in heavy ion collisions would thus be a detection of CP-odd processes in QCD.
One observable consequence of the CME in a heavy ion collision is charge separation: positive charges will move to one side of the reaction plane, negative charges to the other. We know from experiment that the strong interactions preserve P and CP, however, so any charge separation from CP-odd sources will, over many events, average to zero. An observable that could serve as a "smoking gun" for the CME is thus hard to find. For heavy ion collsions at RHIC and LHC the focus so far has been on three-particle correlations [20][21][22], which indeed indicate that charge separa-tion occurs in heavy ion collisions. These correlations are sensitive to event-by-event charge separation from both CP-odd and CP-even processes, however, making a positive identification of a signal from the CME difficult [23]. In short, to date the experimental evidence for the detection of the CME in heavy ion collisions at RHIC and LHC is inconclusive.
The experimental situation raises a number of urgent questions for theorists. Can we compute the size of the signal from the CME, relative to backgrounds? How will that signal depend on temperature, magnetic field, centrality, etc.? Clearly an auxiliary question is: how big is the rate of chirality production, which is ∝ Γ CS , in a heavy ion collision?
Unfortunately, Γ CS is difficult to calculate for the QGP, for the same reasons that the shear viscosity, η, is difficult to calculate. The quarks and gluons are stronglyinteracting, so perturbation theory is a priori unreliable. Calculating transport coefficients, such as Γ CS and η, from lattice QCD requires a problematic analytic continuation from Euclidean signature. 2 Currently no reliable technique exists to compute Γ CS or η for QCD at the temperatures reached in the QGP.
An alternative approach is holography [26][27][28], which equates certain stronglycoupled gauge theories in the large-N c limit with weakly-coupled theories of gravity in spacetimes of one higher spatial dimension. A deconfined thermal state of the gauge theory is dual to a black hole spacetime [29], and transport coefficients are relatively straightforward to calculate [30][31][32]. Remarkably, the ratio of η to entropy density, s, for any theory dual to higher-dimensional Einstein gravity is η/s = 1/(4π) [33], which is close to the estimate for η/s for the QGP extracted from data [34]. Such universality serves as encouragement for computing other transport coefficients, like Γ CS , from holography.
Previous calculations of Γ CS in holography employed "top-down" models, i.e. models descending from a known string theory or supergravity construction. The best-understood example is N = 4 supersymmetric YM (SYM) with large N c and large λ t , dual to supergravity in an Anti-de Sitter (AdS) space, for which Γ CS ∝ λ 2 t T 4 [30]. Other holographic calculations included the effects on Γ CS due to a magnetic field [35] or confinement [36]. To our knowledge, in all previous cases the holographic results for Γ CS were ultimately fixed by some underlying (perhaps "hidden" [36]) conformal symmetry.
In this paper we compute Γ CS in Improved Holographic QCD (IHQCD) [37][38][39][40][41][42][43][44], a holographic model of large-N c YM theory. The model is "bottom-up," i.e. does not descend from a known string theory or supergravity construction, but is tailored to model string theory systems very closely, unlike other bottom-up models. The bulk theory is Einstein-dilaton gravity, where the dilaton Φ is dual to tr F µν F µν . A nontrivial dilaton solution will describe non-trivial running of the YM coupling, hence the choice of dilaton potential is crucial. The simplest choice involves only two free parameters, which can be adjusted such that the model reproduces both the T = 0 glueball spectrum and the thermodynamics of large-N c YM, including a first-order deconfinement transition at a critical temperature T c . In particular, the model has no (hidden) conformal symmetry. We briefly review IHQCD in section 2.
In IHQCD the operator q(x µ ) is dual holographically to an axion field in the bulk [37,38,40,42,44]. Defining for convenience a holographic 't Hooft coupling λ ≡ e Φ , the normalization of the axion's kinetic term includes a dilaton-dependent factor, Z(λ). In principle, Z(λ) could be fixed by matching to lattice results for the Euclidean correlator of q(x µ ), as we explain in section 2. We work instead with several simple choices for Z(λ), in part to study the generic behavior of Γ CS in holographic models. Specifically, we consider a Z(λ) with two free parameters, which we fix by demanding that the model match large-N c YM lattice results for the topological susceptibility and for axial glueball mass ratios to within one sigma.
In section 3 we compute Γ CS in the high-temperature, deconfined phase of IHQCD. Letting s denote the entropy density and λ h the value of λ at the black hole horizon, our result for Γ CS is of the form Figs. 4, 5, and 6 show our numerical results for Γ CS /(sT /N 2 c ) = Z(λ h )/(2π). For our choices of Z(λ), the value of Z(λ h ) is bounded from below as a function of T by its value in the T → ∞ limit, and increases monotonically as T approaches T c from above, with most of the increase occurring between 2T c and T c . We will argue that such behavior is generic in a large class of confining theories with classical gravity duals. In a scan through various choices of Z(λ), each of which reproduces the first two axial glueball mass ratios to within one sigma, we find that the increase can be as large as 60%. For our optimal choice of Z(λ), which provides the best fit to the lattice results for the first two axial glueball mass ratios, the increase is only 0.01%.
To obtain Γ CS , we compute the low-frequency limit of G R (ω, k = 0) holographically. In section 4 we initiate the study of G R (ω, k) at non-zero ω and | k|, in the T ≥ T c regime. We focus in particular on Im G R (ω, k), which is proportional to the spectral function of q(x µ ). After suitably subtracting the high-frequency asymptotics, by computing the difference in the value of the correlator at two temperatures, our results suggest the presence of a reasonably long-lived excitation with energy on order of the lightest axial glueball mass at T = 0. That is sufficiently light to prompt the speculation that perhaps such an excitation could dominate many CP-odd phenomena in the QGP created in heavy ion collisions.
In section 5 we summarize our results and discuss directions for future research.

Improved Holographic QCD
The holographic model that we consider as the dual to pure large-N c YM is (4+1)dimensional Einstein-dilaton gravity with a well-chosen dilaton potential [37][38][39][40][41][42][43][44]. In terms of the holographic 't Hooft coupling λ ≡ e Φ , the bulk action is where M p is the Planck Mass, related to the (4+1)-dimensional Newton's constant G 5 as M 3 p = 1/(16πG 5 N 2 c ), g and R are the determinant and Ricci scalar of the bulk metric, V (λ) is the dilaton potential, and S bdry represents all boundary terms, including the Gibbons-Hawking term as well as the counterterms needed for holographic renormalization [45].
If V (λ) = 12/ℓ 2 with a constant length scale ℓ, then the equations of motion arising from eq. (2.1) admit a solution with constant λ and an AdS metric with radius of curvature ℓ, Here r is the holographic radial coordinate, dual to the field theory energy scale: the region near the AdS boundary at r → 0 is dual to the ultra-violet (UV) of the field theory, while the region near the Poincaré horizon at r → ∞ is dual to the infra-red (IR). Such a solution describes a conformal field theory. For non-trivial V (λ), the equations of motion admit vacuum solutions in which λ depends only on r and the metric takes the form with warp factor b 0 (r). In IHQCD we demand that as r → 0 the metric approach that of AdS, b 0 (r) → ℓ/r (up to corrections logarithmic in r) and that λ vanish logarithmically, λ → −1/ log r, to mimic the running of the large-N c YM coupling. Large-N c YM approaches a free theory in the UV, so we expect the holographic dual in the r → 0 region to be a string theory, not just a classical gravity theory like IHQCD. On the other hand, in large-N c YM, λ t diverges in the IR, so a classical gravity theory may be a reliable description in the r → ∞ region. IHQCD is intended to be such a low-energy effective description of large-N c YM, reliable in the r → ∞ region. In practice, the role of the r → 0 region in IHQCD is simply to provide boundary conditions for the fields in the r → ∞ region. We impose those boundary conditions at a cutoff, i.e. at some small but finite r = ǫ. We then compute low-energy quantities that are insensitive to the cutoff, some of which we match to large-N c YM, while the rest are predictions of the model. A more detailed discussion of these issues appears for example in ref. [41].
Using classical gravity in the r → 0 region has an important consequence, however: IHQCD will actually be dual to a theory that flows to a non-trivial UV fixed point. Generically, the UV physics of large-N c YM and IHQCD will thus be different. For example, in IHQCD, η/s = 1/(4π) [33], which is much smaller than the high-T perturbative result for η/s in large-N c YM [46]. Nevertheless, in order to match IHQCD to known results for IR quantities in large-N c YM, we must match some quantities in the UV. For example, to reproduce lattice results for the free energy of large-N c YM for T T c with the correct normalization, we must demand that at high T the free energy of IHQCD obey a Stefan-Boltzmann law. That requirement fixes the value of ℓ in the asymptotic AdS region in units of the Planck mass: (M p ℓ) −3 = 45π 2 [39,40].
By matching to another UV quantity, the perturbative large-N c YM β-function, we can also constrain V (λ). In the r → 0 region, where λ is small, V (λ) has a regular series expansion Committing to an identification of the field theory renormalization scale E ≡ E 0 b 0 (r), where E 0 can be fixed by matching to the lowest glueball mass or to the result for T c from lattice large-N c YM, we can fix the coefficients v 0 and v 1 in terms of the coefficients of the perturbative large-N c YM β-function [37,38,40,44,45]: In the vacuum solutions, generically λ diverges as r → ∞. The large-λ expansion of V (λ) must take the form V (λ) ∝ λ 4 3 √ log λ in order for the glueball spectrum to be gapped and discrete with asymptotically linear trajectories [37,38,41,44]. With this asymptotic form for V (λ), as r → ∞ the warp factor and λ(r) take the form where L is a length scale determined by the value of λ at r = ǫ. The form of b 0 (r) in eq. (2.6) is sufficient to guarantee that the dual field theory is confining [37,38]. The metric actually has a mild singularity 3 at r = ∞ that can be cloaked by a regular horizon and hence is a "good" singularity [47]. Moreover, the singularity is repulsive [38,40], which guarantees that the low-energy spectrum and other observables are insensitive to the details of the resolution of the singularity.
The black hole solutions of the model defined by the action in eq. (2.1) have non-trivial λ(r) and a metric of the form [39,40] The surface r = r h is the horizon, where f (r h ) = 0, and the corresponding Hawking temperature is T = 4πf ′ (r h ). Black hole solutions only exist for temperatures above a value T min , and in fact two branches of solutions exist, the large and small black holes (comparing r h to ℓ). Fig. 1 depicts the typical form of T as a function of r h , including the two branches of black hole solutions. For both large and small black holes, as r → r h , the warp factor b(r) asymptotes to a constant whose value determines the entropy density, s = b(r h ) 3 /(4G 5 ), and as r → 0, b(r) → r/ℓ, up to O(r 4 ) (times logarithmic) corrections, indicating that in the field theory the thermal energy density and pressure are both of order N 2 c .
Big black holes Small black Holes 0 r min r h T min T Figure 1: Schematic plot for the typical form of the black hole Hawking temperature T as a function of the horizon position r h , for a generic choice of V (λ) (with the correct smalland large-λ asymptotics). The temperature exhibits a minimum, T min , at r min , which separates the large black hole (r h < r min ) from the small black hole (r h > r min ) branches.
In large black hole solutions, λ(r) decreases monotonically as T increases, so that λ → 0 as T → ∞. In small black hole solutions, λ(r) increases as T increases. In particular, as discussed in refs. [40,42], the value of λ(r) at the horizon, λ h ≡ λ(r h ), is a monotonically increasing function of r h , so a plot of T versus λ h is qualitatively similar to fig. 1: in the T → ∞ limit, λ h → 0 on the large black hole branch and λ h → ∞ on the small black hole branch (see for example fig. 2 (a) of ref. [40]).
If we Wick-rotate to a compact Euclidean time direction of length 1/T , then for T ≥ T min three bulk solutions exist: the Wick-rotated version of eq. (2.3), which describes a thermal gas of gravitons and is dual to a confined state, and the Wickrotated large and small black holes, which are dual to deconfined states. To determine which solution is thermodynamically preferred at any given T , we must determine which has the smallest on-shell Euclidean action, dual to the field theory's free energy (times 1/T ). As shown in refs. [39,40,42], the small black hole solutions are never thermodynamically preferred, but at some T c > T min the large black hole solutions become thermodynamically preferred. Indeed, the system exhibits a first-order Hawking-Page type transition at T c , dual to a confinement-deconfinement transition.
In general, for a given potential V (λ) we cannot solve the equations of motion arising from eq. (2.1) exactly, so we resort to numerics. Here we will only sketch our numerical procedure, which is described in detail for example in ref. [42]. At the cutoff r = ǫ we impose a Dirichlet condition on each field, and in particular we demand that the metric take the AdS form. We then fix the remaining integration constants, including λ h , by a shooting algorithm. Given a choice of V (λ) and the Dirichlet conditions at r = ǫ, we obtain a one-parameter family of solutions labeled by T , or equivalently by λ h . Following refs. [42,43], in our numerics we use a simple form for V (λ) with the correct small-and large-λ asymptotics, Expanding eq. (2.8) about λ = 0 and matching to eq. (2.4), we The coefficients V 0 and V 2 can be determined in terms of V 1 by imposing the conditions in eq. (2.5b). The potential thus has two free parameters, V 1 and V 3 . We fit these two parameters by matching to lattice results for two thermodynamic quantities in large-N c YM: the latent heat of the deconfinement transition, which is proportional to the entropy density at the transition, [48], and the pressure at T = 2T c [48][49][50]. Upon fixing V 1 and V 3 in this fashion, IHQCD describes very well both the T = 0 glueball spectra (0 ++ and 2 ++ ) as well as the finite T thermodynamics of large-N c YM [42,51].
The instanton number density operator, q(x µ ) in eq. (1.2b), is dual to a bulk pseudoscalar, the axion α (as in many top-down models). In the field theory, the source for q(x µ ) is an angular variable, the θ-angle. As a result, the action of the bulk axion must be invariant under shifts of α, and hence must depend only on derivatives 4 ∂α. General arguments in string theory and in YM theory, including the argument that θ dependence should appear in the YM vacuum energy only at order one in the large-N c limit rather than at order N 2 c [52], imply that the axion action is suppressed by O(1/N 2 c ) compared to the action S in eq. (2.1) [37,38,41,44]. We thus add to the model an axion with an action S α of the form [37,38,41,44] where, following the rules of effective field theory, we have included a dimensionless, λ-dependent normalization function, Z(λ), consistent with the symmetries. Being a massless pseudo-scalar, in an expansion of α(r) about r = 0, the leading, non-normalizable term is a constant, which is proportional to the YM θ-angle defined in the UV, where in top-down models the proportionality constant κ will be fixed, but not in bottom-up models. In other words, in our model the normalization of the operator dual to α is ambiguous: α is dual to q(x µ )/κ. Nevertheless, by fixing the normalization of the topological susceptibility we will be able to compute two-point functions of q(x µ ) unambiguously, as we explain below.
To specify S α completely we must specify Z(λ). In principle, Z(λ) can be fixed as follows. First, perform a lattice calculation of the Euclidean two-point function of q(x µ ) with non-zero T for some set of frequencies. Second, compute the same Euclidean two-point function holographically for all frequencies for some choice of Z(λ). A least squares fit of the holographic results to the lattice results should then determine Z(λ). To study generic behavior of holographic models, we will instead proceed by using simple forms for Z(λ) that we constrain by matching to lattice results for the topological susceptibility and axial glueball mass spectrum. Notice that matching to any lattice data will always have room for improvement: lattice definitions of q(x µ ) generically suffer from power-law divergences that dominate in the continuum limit, making lattice calculations of correlators of q(x µ ) noisy [24]. Accurate calculations may be possible in the near future 5 .
We can constrain Z(λ) as follows. Since Z(λ) is the coefficient of a kinetic term, we demand that Z(λ) ≥ 0. We can also constrain Z(λ)'s small-and large-λ asymptotics [37,38,41,44]: where Z 0 is a dimensionless constant. The small-λ form follows from the rules of effective field theory: a constant is the most general allowed term. The large λ behavior is fixed by glueball universality [38]. Various towers of glueballs have linear asymptotic trajectories: for large excitation number n, their squared masses go as (m i n ) 2 = c i n + · · · , with constants c i , where the integer i labels different towers. Glueball universality is the statement that all the slopes c i are similar, i.e. do not depend on i. That is automatic for the 0 ++ and 2 ++ glueballs. Requiring the same for the 0 −+ glueballs forces Z(λ) to go as λ 4 at large λ [38].
We will use the simplest form of Z(λ), also used for example in ref. [42], where c 4 is a dimensionless constant. To fix Z 0 we match to the large-N c YM lattice result for the Euclidean topological susceptibility, χ, defined in terms of the T = 0 vacuum energy density E(θ) as where the subscript E denotes the Euclidean correlator. The holographic result for χ is 6 [38], . (2.14) Clearly χ will be proportional to κ 2 Z 0 . Thus, for any These values can then be used to predict the masses in the full tower of 0 −+ glueballs. As shown in refs. [37,38,41,44], the holographic result for the first excited 0 −+ glueball mass, m 0 * −+ , agrees very well with the lattice result [54], Crucially, notice that by fixing the normalization of χ we have fixed the normalization of any two-point function of q(x µ ), and thus have eliminated the normalization 6 The holographic calculation of χ in ref. [38] assumed κ = 1. Here we allow for arbitrary κ. 7 For a recent lattice study of the glueball spectrum at large N , see ref. [55]. We prefer to use the older results of ref. [54] because in the latter work an excited state of the 0 −+ tower is given. 8 The result for κ 2 Z 0 in ref. [42] was too large by a factor of four, producing an erroneous result, κ 2 Z 0 = 133. In eq. (2.16) we present the correct value, κ 2 Z 0 = 133/4 = 33.25. ambiguity mentioned below eq. (2.10). In other words, the holographic calculation of the two-point functions of q(x µ ) will only depend on the combination κ 2 Z 0 (as we will see explcitly in section 3), which we have fixed to the value in eq. (2.16).
Solutions for α as a function of T were studied in refs. [40,42]. When T = 0, a non-trivial UV θ-angle forces α(r) to be non-trivial. The resulting normalizable solution then indicates that the non-zero UV θ-angle flows to zero in the IR, as shown in fig. 2, and additionally triggers a non-zero q(x) /κ. The T = 0 solution for α(r) is unchanged when T < T c : Wick-rotating the metric in eq. (2.3) to a compact Euclidean time does not affect the static solution α(r). Such behavior is expected in a confined phase at leading order in N c , due to large-N c volume independence. When T > T c , however, the only non-singular solution for the axion is a constant, α(r) = κ θ, indicating that q(x) = 0, in agreement with evidence from lattice data for large-N c YM [24]. What is the behavior of the topological susceptibility as a function of temperature, χ(T )? When T < T c , χ(T ) is independent of T , i.e. takes the same value as at T = 0, eq. (2.14), again due to large-N c volume independence. When T > T c , the holographic result for the topological susceptibility is .
(T > T c ) (2.18) The denominator on the right-hand-side of eq. (2.18) diverges at the black hole horizon, so in fact χ(T ) = 0 when T > T c , up to O(e −Nc ) corrections [40,42]. We will also consider a form for Z(λ) more general than that of eq. (2.12). On the large black hole branch, if T is large then λ is small, in which case we expect the largest polynomial correction to the Z(λ) in eq. (2.12) to be a term linear in λ, hence we consider where c 1 is a dimensionless constant, which we choose to be positive. If we continue to fit only to the lattice result for the lowest 0 −+ glueball mass, we find a substantial degeneracy (which is not surprising, given that we have introduced an additional parameter, c 1 ). Specifically, for any positive value of c 1 , a value of c 4 exists such that, upon matching to the lowest axial glueball mass in eq. (2.15), the value of the first excited axial glueball mass is in rough agreement with the value in eq. (2.17), exhibiting at most a 3% discrepancy, as shown in fig. 3. To constrain c 1 we will demand that our holographic results for the axial glueball masses fall within one sigma of the lattice values in eqs. As we have seen, the function Z(λ) must be non-negative and is constrained in the λ → 0 and λ → ∞ limits. For intermediate values of λ, the most natural assumption is that Z(λ) is monotonic. At least, we are not aware of any compelling evidence for the existence of maxima or minima in Z(λ). Our choices for Z(λ) were thus monotonic functions of λ, namely polynomials in λ with strictly positive coefficients. To test the effect of maxima and minima in Z(λ), we considered two changes to the Z(λ) in eq. (2.19). First, we allowed slightly negative c 1 , while maintaining Z(λ) ≥ 0. Second, we introduced a maximum by a adding a Gaussian peak to Z(λ). In each case we computed the axial glueball mass spectrum. After matching to the lattice result for the lowest axial glueball mass, we found that the fit to the first excited axial glueball mass was worse, deviating from the lattice result by about 10%. We consider that a preliminary indication that monotonic Z(λ) may indeed be the best choice. We leave more thorough tests for future research.
Our assumption that Z(λ) is monotonic in λ determines the qualitative behavior of Z(λ) as a function of T . On the large black hole branch, as T → ∞, λ → 0, and as T decreases towards T c , λ increases monotonically. As a result, for our choices of Z(λ)-simple polynomials in λ with positive coefficients-when T → ∞, Z(λ) → Z 0 , and when T → T c , Z(λ) will increase monotonically. As functions of T , our Z(λ) are thus bounded from below by their value in the T → ∞ limit: Z(λ) ≥ Z 0 . The behavior of Z(λ) as a function of T will translate directly into the behavior of Γ CS as a function of T , as we will show in the next section. In particular, the dimensionless combination Γ CS /(sT ) will be bounded from below by its value in the T → ∞ limit, and will increase as T → T c . In the next section we will also present a more general argument that Γ CS /(sT ) must increase as T approaches T c from above.  [54]. Only the mass of the n = 2 state is appreciably sensitive to changes of c 1 and c 4 , differing from the lattice result by 3% at most.

The Chern-Simons Diffusion Rate
We will compute Γ CS using eq. (1.4), rewritten as whereĜ R (ω, k) is the retarded two-point function of q(x µ )/κ, the operator dual to our axion α.
In holography, the on-shell bulk action is the generating functional for field theory correlation functions [27,28]. To compute the two-point functionĜ R (ω, k) in the high-temperature, deconfined phase of IHQCD, we must solve the linearized equation of motion of the axion in the black hole spacetime with metric in eq. (2.7), with T ≥ T c . We thus introduce a fluctuation of the axion, δα(r, x µ ), where x µ = (t, x). When T ≥ T c , the background solution for the axion is trivial, hence the linearized equation of motion for δα(r, x µ ) is simply where the metric is that of eq. (2.7). Notice in particular that δα will not couple to the fluctuations of any other fields because the background solution preserves CP and the axion is the only CP-odd field in the bulk. We must solve eq. (3.2) with Dirichlet boundary condition at the asymptotically AdS boundary and with in-going wave boundary condition at the horizon [30]. The solution takes the form where k µ = (ω, k) and where a(k µ ) is fixed by the Dirichlet boundary condition, while δα(r, k µ ) obeys the equation 1 Z(λ(r)) √ −g ∂ r Z(λ(r)) √ −g g rr ∂ r δα(r, k µ ) − g µν k µ k ν δα(r, k µ ) = 0, (3.5) with unit normalization at the asymptotically AdS boundary, lim r→0 δα(r, k µ ) = 1, and in-going wave boundary condition at the horizon. The on-shell axion action is then The retarded Green's function is then [30] To compute Γ CS , we need to solve eq. (3.5) with k = 0 and with small ω. We will do so in two ways, first using near-horizon matching and second using the membrane paradigm, following ref. [32]. In each case we can determine Γ CS analytically, essentially because δα is a massless fluctuation.
In the near-horizon matching technique, we first solve eq. (3.5) with ω = 0 and then expand the solution near the horizon. We then reverse the order of operations, solving the equation in the near-horizon region and then expanding the solution in ω. Finally, we match the two solutions to obtain F (r, k µ ).
When k = 0 and ω = 0 the solution of eq. (3.5) is with constant coefficients C 1 and C 2 . The second term on the right-hand side of eq. (3.9) diverges as r → r h . As a result, when ω = 0 a normalizable solution must have C 2 = 0. When ω is small but non-zero, a normalizable solution may have C 2 ∝ ω. Plugging eq. (3.9) into eq. (3.7), we find (ω ≪ T, k = 0) (3.10) We will choose C 1 = 1 so that our δα has unit normalization at the asymptotically AdS boundary. Our task is thus to determine C 2 . Expanding the solution in eq. (3.9) around the horizon, we find where f ′ (r h ) = 4πT . Now we reverse the order of operations. Expanding eq. (3.5) in (r h − r), we find the solution in the near-horizon region, with coefficients C ± that depend on ω but not on r. We set C + = 0 so that the nearhorizon solution is an in-going wave [30]. Now we expand the solution in eq. (3.12) for small ω: By matching the constant and logarithmic terms in eqs. (3.11) and (3.13), we find (3.14) Setting C 1 = 1, we obtain lim r→0 F (r, k µ ) via eq. (3.10) and thenĜ R (ω, k) via eq. (3.8),Ĝ We thus obtain our main result for Γ CS , where we have used M 3 p = 1/(16πG 5 N 2 c ) and where s = b 3 (r h ) 4G 5 is the entropy density. Notice that the normalization of this result is fixed by the product κ 2 Z 0 , which we fixed in section 2 by matching to the topological susceptibility at T = 0.
The second equivalent, but more efficient, method that we will use to obtain G R (ω, k) is the membrane paradigm [32]. Kubo's formula for the retarded Green's function is Π(ω, k) =Ĝ R (ω, k)δα(ω, k), (3.17) where Π(ω, k) is the one-point function of q(x µ )/κ in Fourier space. Following ref. [32], we extend eq. (3.17) into the bulk by defining an r-dependent response function, , (3.18) where Π(r, ω, k) is the canonical momentum of δα(r, ω, k) with respect to the rfoliation of the bulk space-time, The retarded Green's function is then proportional to the boundary value of ζ: An equation of motion for ζ is straightforward to derive using eq. (3.19) and δα's equation of motion, eq. (3.5), ∂ r ζ = ω Z(λ(r)) √ −g g rr ζ 2 + Z(λ(r)) 2 g g rr g tt 1 + g xx g tt To obtain the retarded Green's functionĜ R (ω, k), we must impose regularity at the horizon, meaning ∂ r ζ is finite there [32], hence the term in brackets in eq. (3.21) must vanish 9 at r = r h : We can now easily derive Γ CS . In eq. (3.21) we take k = 0 and observe that if ω → 0 then ζ becomes independent of r. The value of ζ for all r is then the same as the value at the horizon, eq. (3.22), and via eq. (3.20) we trivially obtainĜ R (ω, k = 0), which is identical to eq. (3.15). We thus find again Our result suggests a natural dimensionless quantity to study, which has implicit dependence on T through Z(λ h ), and is constant in T if and only if Z(λ) is a constant in λ, that is, if the axion does not couple to the dilaton. Indeed, as we mentioned at the end of section 2, the behavior of Z(λ) as a function of T determines the behavior of Γ CS /(sT /N 2 c ) as a function of T . In particular, on the large black hole branch, Γ CS /(sT /N 2 c ) is bounded from below by its value in the T → ∞ limit, If we use the preferred value κ 2 Z 0 = 33.25 [42] then κ 2 Z 0 /(2π) ≃ 5.29. Moreover, Γ CS /(sT /N 2 c ) will increase monotonically as T approahces T c from above. For the simplest choice of Z(λ), given in eq. (2.12), Γ CS /(sT /N 2 c ) has extremely mild dependence on T : as T approaches T c from above, Γ CS /(sT /N 2 c ) is nearly constant, experiencing an increase of only about 0.01%, mostly between 2T c and T c , as shown in fig. 4. In bulk terms, the reason for this mild T dependence is that between T → ∞ and T = T c , λ h increases from zero up to only λ h ≈ 0.14, which for the Z(λ) in eq. (2.12) translates into a very small change in Γ CS /(sT /N 2 c ). increases near T c by more than a factor of six. For all values of c 1 and c 4 that we considered, most of the increase occurs between 2T c and T c . Fig. 6 shows Γ CS /(sT /N 2 c ), normalized to the T → ∞ value κ 2 Z 0 /(2π), as a function of T /T c for values of c 1 and c 4 that reproduce the lattice results for axial glueball mass ratios to within one sigma, eq. (2.20). At the upper limits of the allowed (c 1 , c 4 ) values, namely (c 1 , c 4 ) = (5, 50), we find that as T approaches T c from above, Γ CS /(sT /N 2 c ) increases by about 60%, with most of the increase occuring between 2T c and T c . 2). In all of these cases, as T approaches T c from above Γ CS /(sT /N 2 c ) increases by anywhere from 0.01% up to a factor greater than six. The increase occurs mostly between 2T c and T c .
In heavy ion collisions at RHIC and LHC, T reaches two to four times the QCD crossover temperature. We would thus like to know the value of Γ CS in QCD near the crossover temperature, which is a key ingredient determining the magnitude of any current produced via the CME [18]. 10 No controlled calculation of Γ CS from QCD at these temperatures exists, hence we turn to holography. Suppose we use N = 4 SYM as a holographic proxy for QCD near the crossover temperature. The result for Γ CS in large-N c , strongly-coupled N = 4 SYM is [30], Being a conformal field theory, N = 4 SYM has no phase transitions at non-zero T , so to obtain a sensible result we should consider the dimensionless quantity Γ CS /T 4 . As a crude estimate we take α s ≡ g 2 /(4π) = 0.5 and we use N c = 3, so that λ t = 6π, in which case we find Γ N =4 CS /T 4 ≈ 0.045. (λ t = 6π) (3.27) For a better estimate, let us consider Γ CS (T c )/T 4 c in IHQCD. As discussed above, if Z(λ) is monotonic in λ, then Γ CS /(sT /N 2 c ) is bounded from below by its value in the T → ∞ limit, eq. (3.25). We can obtain a lower bound on Γ CS (T c )/T 4 c by using the large-N c YM lattice result for the entropy density at T c [48], s(T c ) = 0.31N 2 c T 3 c . Letting λ c denote the value of λ h at T c , we find Finally, we have also calculated Γ CS using the small black hole solutions [39,40,42]. Our results for those cases appear in the appendix. Although the small black hole branch is always thermodynamically disfavored, we can actually use the results for Γ CS /(sT /N 2 c ) on the small black hole branch to argue quite generally that on the large black hole branch Γ CS /(sT /N 2 c ) should increase as T approaches T c from above. A similar argument also applies for the bulk viscosity, as discussed in ref. [43]. The key result, shown in fig. 11 in the appendix, is that for T > T min Γ CS /(sT /N 2 c ) is larger on the small black hole branch than on the large black hole branch, but the two branches meet at T min . On the large black hole branch, then, Γ CS /(sT /N 2 c ) must increase as T → T min from above, in order to meet Γ CS /(sT /N 2 c ) from the small black hole branch. In fact, we can show in full generality that on the large black hole branch Γ CS /(sT /N 2 c ) must increase as T → T min : we simply take (d/dT )(Γ CS /(sT /N 2 c )) = (dλ h /dT )(d/dλ h )(κ 2 Z(λ h )/2π) and observe that by definition (dλ h /dT ) diverges when T → T min , while (d/dλ h )(κ 2 Z(λ h )/2π) remains finite. Notice also that Γ CS /(sT /N 2 c ) itself remains finite when T → T min . Given that T min is generally very close to T c , we are then guaranteed that Γ CS /(sT /N 2 c ) will be increasing as T → T c from above, if we assume that Z(λ) is monotonic as a function of T between T min and T c . In principle, Z(λ) could exhibit maxima or minima for T ∈ (T min , T c ), although such behavior seems un-natural. On the large black hole branch an increase of Γ CS /(sT /N 2 c ) as T → T c from above seems to be the generic behavior. We thus learn that the increase in Γ CS /(sT /N 2 c ) in the vicinity of T c on the large black hole branch is tied to the existence of T min , and hence to the existence of small black hole solutions. As argued in ref. [40], the existence of small black hole solutions follows from the fact that the zero-temperature theory is confining. These arguments suggest that perhaps any confining, strongly-interacting, large-N c gauge theory with a (4+1)-dimensional holographic dual 11 may exhibit an increase in Γ CS /(sT /N 2 c ) in the vicinity of T c .

The Spectral Function
We now turn our attention to G R (ω, k) with non-zero ω and k. Generically G R (ω, k) is a complex-valued function of the real variables ω and k. A pole in G R (ω, k) indicates a large response to an infinitesimal source for q(x µ ), and is thus associated with a resonant excitation of the system. Being complex-valued, G R (ω, k) is not directly observable. To study the excitations of our system, we thus turn to the spectral function, −2 Im G R (ω, k), which is real and hence observable in principle. 12 Typically, a pole in G R (ω, k) produces a peak in the spectral function. In this section we initiate the study of these peaks in our system.
To be precise, we will compute Im G R (ω, k). To do so, we will compute G R (ω, k) using the membrane paradigm [32], as explained in section 3. In particular, we must solve eq. (3.21), which we reproduce here for convenience with the boundary condition in eq. (3.22), and then obtain G R (ω, k) via eq. (3.20), We have not been able to solve eq. (4.1) exactly for all values of ω and k, hence we turn to numerical solutions. In this section we exclusively use the Z(λ) in eq. (2.12), with c 4 = 0.26. We consider first the case k = 0. Fig. 7 shows our numerical result for Im G R (ω, k = 0)/(T c M 3 p ) at T c as a function of ω/T c . As we saw in section 3, for ω sufficiently small, Im G R (ω, k = 0) ∝ ω. On the other hand, at asymptotically large ω we expect Im G R (ω, k = 0) ∝ ω 4 because in the UV the theory is conformally invariant and q(x µ ) is dimension four. Our results are consistent with that expectation: fig. 7 shows that the function (1.6 × 10 −7 ) × (ω/T c ) 4.051 provides an excellent fit to our data.
The ω 4 scaling of Im G R (ω, k = 0), and hence of G R (ω, k = 0), at asymptotically large ω is a divergence in the coincidence limit of the two-point function that prevents the correlator from obeying the sum rules and dispersion relations typically used to give physical meaning to the poles of G R (ω, k) in the complex ω plane, which require G R (ω, k) to vanish at large frequency. Such a divergence may overwhelm peaks in Im G R (ω, k), rendering them practically invisible.
One way to improve the large-ω behavior of Im G R (ω, k) is to consider subtracted correlators. For example, one possible option is to determine the form of Im G R (ω, k) at large ω exactly by solving eq. (4.1) in a WKB approximation, and then subtracting that large-ω form from all subsequent calculations of Im G R (ω, k). That approach encounters ambiguities in sub-leading divergences in ω, as discussed for example in ref. [56]. We will instead eliminate the large-ω divergence by computing G R (ω, k) at two different temperatures, T 1 and T 2 , and then taking the difference, (4.4) We could also imagine subtracting the T = 0 result for G R (ω, k), that is, by taking T 1 = 0, but that is difficult to do numerically. When T = 0, G R (ω, k) is a sequence of delta-functions whose locations and amplitudes correspond to the masses and wave-function normalizations of axial glueballs. We would need to subtract the enveloping function of this sequence of delta-functions, which is difficult to implement numerically. We will thus always consider T 1 , T 2 ≥ T c . Fig 8 shows our numerical results for Im G R (ω, k = 0) at two different temperatures, T c and 2T c , while fig. 9 shows our numerical results for ∆Im G R (ω, k = 0; T c , 2T c ). In each figure we observe that the difference in Im G R (ω, k = 0) between T c and 2T c approaches zero as ω/T c → ∞, at least within our numerical precision. Our numerical subtraction thus appears to be reliable, so we may interpret peaks in Im G R (ω, k) as physical excitations.
From figs. 8 and 9, we see that as T increases from T c to 2T c , ImG R (ω, k = 0) changes by at most 10%. Fig. 9 also clearly reveals a minimum in ∆Im G R (ω, k = 0; T c , 2T c ) near ω/T c ≈ 10 and a maximum near ω/T c ≈ 22, indicating a shift in spectral weight towards higher ω as T increases. Indeed, fig. 9 strongly suggests that a peak in the spectral function is moving to higher ω as T increases. The location of the peak, at ω on the order of twenty times T c , is roughly the same as the scale of the lightest 0 −+ glueball mass at T = 0, around 2600 MeV [54]. In other words, fig. 9 provides evidence that the plasma supports an excitation with roughly the same energy as the lightest 0 −+ glueball at T = 0. The width of the peak in fig. 9 is about 10T c ≈ 1300 MeV, so the excitation is reasonably long-lived.  using the same two temperatures as above. We observe that as | k| increases up to | k|/T c ≈ 10, the largest peak shifts from ω/T c ≈ 22 up to ω/T c ≈ 30. Although this change in the position of the peak is roughly order one, the change in the shape of the peak is very mild. In particular, the width of the peak changes very little, indicating that the lifetime of the excitation stays nearly constant as | k| increases.
The typical time scale for dynamical processes in the QGP created in heavy ion collisions is about 1 fm/c ≈ (200 MeV) −1 . Our results suggest the existence of a relatively long-lived excitation with energy on the order of 2600 MeV, corresponding to a time scale of about 0.1 fm/c. We cannot resist speculating that perhaps such an excitation, if present in the QGP, could dominate correlators of q(x µ ) and hence many dynamical CP-odd phenomena. Regrettably, we will leave a detailed analysis of this excitation, and its effect on CP-odd physics, for the future.

Discussion and Outlook
IHQCD is a state-of-the-art bottom-up holographic model for the low-energy physics of (3+1)-dimensional large-N c YM theory. In this paper we computed the retarded Green's function of the instanton density operator q(x µ ) in the high-temperature, deconfined phase of IHQCD. Our primary motivation was to compute the Chern-  Simons diffusion rate, Γ CS , with the result in eq. (1.5). In particular, our result for Γ CS is proportional to Z(λ h ), where Z(λ) is the normalization factor of the bulk axion action, and λ h is the value of the holographic 't Hooft coupling at the black hole horizon. A combination of available data for the topological susceptibility and axial glueball spectrum of large-N c YM, and glueball universality, are sufficient to determine the small and large λ limits of Z(λ) [37,38,41,44]. We considered several forms for Z(λ). Assuming that Z(λ) is a monotonic function of λ, we found quite generally that Γ CS /(sT /N 2 c ) is bounded from below by its value in the T → ∞ limit and increases monotonically as T → T c from above. Indeed, we presented an argument that the same will be true in many (3+1)-dimensional, confining, strongly-coupled, large-N c theories with holographic duals. For the Z(λ) producing our optimal fit to the lattice results for the axial glueball spectrum, we found that the increase was only 0.01%. Fixing Z(λ) completely by a least-squares fit to lattice results for the Euclidean two-point function of q(x µ ), as explained in section 2, is an important task for the future. We also presented evidence for a relatively long-lived excitation in the system with energy roughly on the order of the mass at T = 0 of the lightest 0 −+ glueball, which prompted us to speculate that perhaps such an excitation could dominate CP-odd phenomena in the QGP created in heavy ion collisions.
IHQCD is dual to pure large-N c YM, so an important goal for the future is to include the effects of quarks in the holographic calculation of Γ CS . Some key questions are how the quark mass and chiral symmetry breaking affect Γ CS . The axial and vector flavor U(1) currents are dual to two U(1) Maxwell fields in the bulk, and the quark mass operator is dual to a complex scalar field, a tachyon, that is bi-fundamental under these two gauge fields. In the bulk, the axion couples to the axial U(1) gauge field and the to phase of the tachyon, as explained in refs. [57,58]. A solution for the tachyon describing either a non-zero quark mass or chiral symmetry breaking can thus influence the axion and affect Γ CS .
Introducing flavors fields would also enable us to compute holographically the current produced via the CME. A preliminary requirement is a bulk solution describing a magnetic field and a net chirality.
We plan to study these and other related issues in the future.

Appendix: The Small Black Hole Branch
As discussed in section 2, when T > T min IHQCD admits two branches of black hole solutions, large black holes and small black holes [40]. In this appendix we compute Γ CS using the small black hole solutions.
(a) (b) Figure 11: (a.) Our numerical result for Γ CS /(sT /N 2 c ), divided by Z 0 κ 2 /(2π), as a function of T /T c , for the Z(λ) in eq. (2.12) with c 4 = 0.26. The upper dot-dashed blue curve is our result obtained from small black hole solutions, while the lower solid blue curve is the result obtained from large black hole solutions. Both curves begin at T min , indicated by the vertical dashed black line, which is slightly below T c . The result on the small black hole branch increases as T increases, and in the T → ∞ limit approaches the form in eq. (1). (b.) Close-up of (a.) near T min .
We can determine the dependence of Γ CS on T in the large-T limit of the small black hole solutions as follows. For generality, we will consider a dilaton potential V (λ) whose large-λ asymptotic form is V (λ) ∝ λ 4/3 (log λ) P , with P a non-negative real number. In the body of the paper we used P = 1/2. From fig. 1, we observe that for the small black hole solutions, when T is large, r h is also large. When r h is large, λ h is also large, in which case we can approximate Z(λ h ) ≈ Z 0 c 4 λ 4 h and hence Γ CS /(sT /N 2 c ) ≈ κ 2 Z 0 c 4 λ 4 h /(2π). As shown in refs. [37,38,41,44], in the r → ∞ limit, λ(r) ∝ exp(r 1/(1−P ) ) r 3 4 P 1−P . Evaluating at r h gives us λ h in terms of r h . From ref. [40] we know r h in terms of T on the small black hole branch in the r h → ∞ limit, r h ∝ T (1−P )/P . We thus find where C is a dimensionless positive constant that depends on the choice of V (λ).
To compute Γ CS in the entire range T min < T < ∞, we resorted to numerics. For the Z(λ) in eq. (2.12) with c 4 = 0.26, our results appear in fig. 11, where we see clearly that the result grows as T increases, and in the T → ∞ limit approaches the form in eq. (1). Fig. 11 also shows that the result for Γ CS /(sT /N 2 c ) is always larger on the small black hole branch than on the large black hole branch, except at T min where the two are equal. This result is important for our argument at the end of section 3 that Γ CS /(sT /N 2 c ) computed on the large black hole branch will increase as T approaches T c from above.