Second order equilibrium transport in strongly coupled $\mathcal{N} = 4$ supersymmetric $SU(N_c)$ Yang-Mills plasma via holography

A relativistic fluid in 3+1 dimensions with a global $U(1)$ symmetry admits nine independent static susceptibilities at the second order in the hydrodynamic derivative expansion, which capture the response of the fluid in thermal equilibrium to the presence of external time-independent sources. Of these, seven are time-reversal $\mathbb{T}$ invariant and can be obtained from Kubo formulas involving equilibrium two-point functions of the energy-momentum tensor and the $U(1)$ current. Making use of the gauge/gravity duality along with the aforementioned Kubo formulas, we compute all seven $\mathbb{T}$ invariant second order susceptibilities for the $\mathcal{N} = 4$ supersymmetric $SU(N_c)$ Yang-Mills plasma in the limit of large $N_c$ and at strong 't-Hooft coupling $\lambda$. In particular, we consider the plasma to be charged under a $U(1)$ subgroup of the global $SU(4)$ R-symmetry of the theory. We present analytic expressions for three of the seven $\mathbb{T}$ invariant susceptibilities, while the remaining four are computed numerically. The dual gravitational description for the charged plasma in thermal equilibrium in the absence of background electric and magnetic fields is provided by the asymptotically AdS$_5$ Reissner-Nordstr\"{o}m black brane geometry. The susceptibilities are extracted by studying perturbations to the bulk geometry as well as to the bulk gauge field. We also present an estimate of the second order transport coefficient $\kappa$, which determines the response of the fluid to the presence of background curvature, for QCD, and compare it with previous determinations made using different techniques.


INTRODUCTION AND BASIC SET-UP
The holographic principle and its realization via the Anti de Sitter/Conformal Field Theory (AdS/CFT) correspondence [1][2][3] has been a powerful tool in exploring the properties of strongly coupled quantum field theories which are not amenable to a perturbative analysis. The correspondence allows one to study strongly interacting quantum systems via a dual gravitational description in an asymptotically AdS spacetime, and has been used extensively to explore several phenomena in high energy physics as well as condensed matter systems, including transport phenomena at finite temperature (see [4][5][6][7][8][9][10] for reviews, including some applications).
In this article, we employ the AdS/CFT correspondence to study the equilibrium thermodynamic properties of a strongly coupled charged plasma at a finite temperature and chemical potential. The particular holographic quantum system we have in mind is the N = 4 supersymmetric SU (N c ) Yang-Mills plasma at finite temperature and chemical potential in 3+1 dimensions in thermal equilibrium. The theory enjoys a global SU (4) R-charge symmetry, and we consider it to be charged under a U (1) subgroup of the full R-symmetry group. 1 In fact, this U (1) is the diagonal U (1) of the maximal Abelian subgroup U (1) 3 of SU (4) [12]. An important point to keep in mind is that this U (1) symmetry is axial in nature, implying that the associated chemical potential µ and current J µ are axial too. 2 Due to the chiral anomaly of the theory [3,13], the current J µ is not conserved, and its non-zero divergence is proportional to the product of external electric and magnetic fields applied to the system: ∇ µ J µ ∝ E · B. However, in the present article, we will study the transport properties of the plasma in the absence of any external electric and magnetic fields, thereby restoring current conservation.
In the limit of large N c and with the 't-Hooft parameter λ ≡ g 2 YM N c 1 i.e. at strong coupling, the super-Yang-Mills plasma in thermal equilibrium, in the absence of background electric and magnetic fields, is holographically dual to a charged Reissner-Nordström black brane solution of the five dimensional Einstein-Maxwell-Chern-Simons theory in asymptotically AdS 5 spacetime [12,14]. 3 The transport properties of the plasma can be extracted by studying perturbations to this black brane geometry. Studying transport in the super-Yang-Mills plasma provides an indirect route for a better understanding of transport in strongly coupled systems found in nature, such as the quark-gluon plasma (QGP). Produced in relativistic heavy ion collisions, the QGP is strongly coupled and can be considered conformal at the very high energy scales involved, since the quark masses can be neglected at such high energies to a very good approximation. As the strongly coupled N = 4 super-Yang-Mills plasma can be probed via holography, wherein its properties can be studied using a dual gravitational description, it provides a simpler and tractable model for the much more complicated dynamics of the QGP which is hard to probe via direct computations [17]. 4 It has been observed experimentally that QGP behaves like a low viscosity relativistic fluid [19], with the shear viscosity to entropy density ratio close to the Kovtun-Son-Starinets bound η/s ≥ 1/4π [20]. Its dynamics can thus be modeled using the formalism of relativistic viscous hydrodynamics 1 A comprehensive overview of R-symmetry charges and chemical potentials in N = 4 SYM theory appears in [11]. 2 An axial chemical potential is parity odd but charge conjugation even. 3 The dual bulk solution is different, and in general is to be constructed numerically, once the plasma is kept in external electromagnetic fields. See for instance [15,16] for recent discussions where charged magnetic black brane solutions have been constructed, dual to the charged N = 4 super-Yang-Mills plasma kept in a background magnetic field, and the resulting transport properties of the plasma have been studied. 4 A simple holographic model for QCD wherein the breaking of conformality is captured by a non-trivial dilaton profile in the bulk, dual to a marginally relevant operator on the boundary, appears in [18]. [21,22]. Hydrodynamics provides an effective description for the dynamics of many-body systems at or near thermal equilibrium, and is a vast subject in itself (see [19,[23][24][25] for pedagogical reviews). It owes its successes to the idea that the effective dynamics of such near-equilibrium systems can be captured in terms of a small number of variables associated to the system. These hydrodynamic variables are typically the local fluid four-velocity u µ and the temperature T . In case the system also has a conserved U (1) charge, then the associated chemical potential µ also becomes a hydrodynamic variable. The dynamics of the near-equilibrium system is governed by the conservation equation for the energy-momentum tensor T µν , along with the conservation equation for the current J µ . The hydrodynamic approximation asserts that T µν and J µ each admits a derivative expansion, with the transport parameters of the system, such as viscosities and conductivities, appearing as coefficients in this expansion. Note that the transport parameters are functions of the state of the system, parametrized by the temperature and chemical potential, and act as inputs to the macroscopic hydrodynamic description i.e. they have to be computed using the underlying microscopic theory describing the system.
In the limit of thermal equilibrium, as was pointed out in [26,27], the physical properties of the fluid can further be encoded in a static generating functional W[g µν , A µ ]. This generating functional is a functional of time-independent external sources, which source the conserved currents of the system. 5 For instance, the energy-momentum tensor T µν is sourced by the background metric g µν , and the U (1) current J µ is sourced by the corresponding background gauge field A µ . The hydrodynamic equations, corresponding to energy-momentum and current conservation, follow as a consequence of the diffeomorphism and gauge invariance of the generating functional. Here F µν = ∂ µ A ν − ∂ ν A µ is the field strength tensor associated with the background field A µ . The gauge invariance of the generating functional and the ensuing conservation of the current J µ can get violated in case the underlying microscopic system has quantum anomalies, which need to be accounted for in the macroscopic hydrodynamic description appropriately [36,37]. When the background metric and the gauge field vary slowly over length scales much larger than the typical microscopic length scales associated to the system in thermal equilibrium, such as the mean free path or the static correlation length, the generating functional W admits a derivative expansion in the hydrodynamic variables and the sources. The various coefficients that enter this derivative expansion are called thermodynamic susceptibilities, which quantify the response of the equilibrium system to the presence of external time-independent sources. One can obtain the equilibrium one-point functions for the energy-momentum tensor and the current by varying the static generating functional with respect to the sources: (1.2) 5 A more recent approach based on the Schwinger-Keldysh closed time path formalism for non-equilibrium systems has been developed to study fluctuating hydrodynamics [28][29][30][31][32][33] (see [34] for a review). The approach allows one to systematically construct the full hydrodynamic effective action including fluctuations about the state of thermal equilibrium, taking into account the effects of non-linear interactions between hydrodynamic modes, purely based on the symmetries of the system. Of particular importance is the dynamical Kubo-Martin-Schwinger (KMS) condition [31] used in constructing the effective action, which ensures that the retarded correlation functions of hydrodynamic variables that follow from the action have the correct analyticity properties in the static ω → 0 limit [35], which is of interest for the present work. TABLE I. Parity P, charge conjugation C, and time reversal T eigenvalues for the nine second order equilibrium invariants for a U (1) charged relativistic fluid in 3+1D. Since the U (1) charge we have for the super-Yang-Mills plasma is axial, the electric and magnetic fields that appear here are taken to be axial as well. If the U (1) symmetry is vector in nature, then the seven T-even susceptibilities f 1 , . . . f 7 are P-even as well.
The resulting constitutive relations for T µν and J µ automatically inherit a derivative expansion, and the coefficients that appear are the transport parameters of the system in equilibrium. Thus, in thermal equilibrium, the susceptibilities appearing in the generating functional are the fundamental objects associated to the system, with the transport parameters entering the constitutive relations being linear combinations of the susceptibilities and their derivatives [26,27]. 6 For a relativistic fluid in thermal equilibrium, there are no first order terms in the generating functional. However, for a U (1) charged fluid in 3+1 dimensions, there are nine independent thermodynamic susceptibilities at the second order in derivatives [26]. The generating functional has the form [39] where p(T, µ) is the equilibrium pressure of the fluid, the only zeroth order contribution to the generating functional. The S n are nine second order diffeomorphism and gauge invariant quantities constructed from the sources. These are (1.4) Here R is the Ricci scalar for the background metric g µν , a µ ≡ u ν ∇ ν u µ is the acceleration and Ω µ ≡ µναβ u ν ∇ α u β is the vorticity vector. The electric and magnetic fields are given by E µ ≡ F µν u ν and B µ = 1 2 µναβ u ν F αβ . Note that a µ , Ω µ , E µ , B µ are all orthogonal to the velocity u µ . The parity, charge conjugation and time reversal properties of the second order invariants are presented in table I. The f n (T, µ) that appear as coefficients of various second order invariants in the equilibrium generating functional eq. (1.3) are the thermodynamic susceptibilities of the system. These susceptibilities, along with the equilibrium pressure p, characterize and capture the equilibrium response of the fluid to external sources, up to the second order in derivatives. Of the nine second order susceptibilities f 1 , . . . f 9 , seven are time-reversal T invariant, for both the cases of a vector or an axial U (1) charge. These are the susceptibilities f 1 , . . . f 7 . Kubo formulas for these seven susceptibilities in terms of the As mentioned above, the holographic description of the charged strongly coupled N = 4 super-Yang-Mills plasma in thermal equilibrium in the absence of background electric and magnetic fields is provided by the Reissner-Nordström black brane solution in asymptotically AdS 5 spacetime [12,14]. By perturbing this geometry and the gauge field in the bulk, one can compute the retarded two-point functions of the energy-momentum tensor and the U (1) current for the boundary theory purely from the bulk perspective. By computing these retarded correlation functions and making use of the Kubo formulas, it is straight forward to extract the seven T-invariant thermodynamic susceptibilities f 1 , . . . f 7 for the strongly coupled N = 4 super-Yang-Mills plasma. In fact, the behaviour of these second order susceptibilities as functions of the parameters µ and T , extracted using holography, is one of the main results of this paper. These results for the equilibrium transport properties of the super-Yang-Mills plasma serve as a proxy for the behaviour of more realistic strongly interacting systems found in nature, such as strongly coupled QCD matter at a non-zero baryon number chemical potential. 7 Let us briefly comment upon some of the known aspects of second order transport for the strongly coupled N = 4 super-Yang-Mills plasma, which is a conformal theory. For the uncharged state of the plasma, there are five transport parameters at the second order in the derivative expansion. These are usually denoted by τ π , κ, λ 1 , λ 2 and λ 3 , and in the limit of infinite 't-Hooft coupling are given by [21,22] Here η is the shear viscosity of the plasma, which is the only transport parameter at the first order for the uncharged state, and is given by η = π 8 N 2 c T 3 [41]. 8 Note that κ and λ 3 are thermodynamic in nature, and are related to our f 1 , f 3 via κ = −2f 1 and λ 3 = 2 T ∂f 1 ∂T + µ ∂f 1 ∂µ − 4f 3 [39]. Interestingly, the parameters λ 1 , λ 2 and τ π satisfy the Haack-Yarom relation [53] which is satisfied by a large class of conformal gauge theory plasmas. Leading corrections to the results in eq. (1.6) ensuing from a finite but large value of the 't-Hooft coupling λ have been found in [54][55][56][57]. In particular, [57] argues that the Haack-Yarom relation holds true even when leading coupling constant corrections are taken into account. Finally, for the U (1) charged super-Yang-Mills plasma, there are eight additional transport parameters at the second order. This entire set of thirteen second order parameters was computed in [58,59] using the fluid/gravity correspondence [22]. As mentioned earlier, in thermal equilibrium, the generating functional and the susceptibilities provide the fundamental description of the macroscopic state of the system, with the transport parameters being derived quantities. We therefore obtain the second order susceptibilities directly in the present article. Other important references for transport in N = 4 super-Yang-Mills plasma include [60,61].
The paper is organized as follows. In section 2, we provide details of the holographic geometry dual to the charged strongly coupled N = 4 super-Yang-Mills plasma in thermal equilibrium, i.e. the Reissner-Nordström black brane solution in asymptotically AdS 5 spacetime. In section 3 we provide details of the procedure used to extract the seven T-even second order susceptibilities f 1 , . . . f 7 for the charged super-Yang-Mills plasma. We also compare our results to some of the things previously known in literature, and also make an estimate of the transport parameter κ for QCD using our holographic results, comparing it with previous determinations made using different techniques. We end with some concluding comments in section 4. Appendices provide supplementary material relevant for the discussion in the main body of the paper, including details of the numerical procedure in appendix B.

THE HOLOGRAPHIC MODEL
As per the AdS/CFT correspondence [1][2][3], the four-dimensional N = 4 supersymmetric SU (N c ) Yang-Mills theory is holographically dual to type IIB superstring theory on AdS 5 × S 5 . Further, in the limit where N c → ∞ and the gauge theory is strongly coupled g 2 YM N c → ∞, the dual bulk description simplifies considerably and is provided by classical type IIB supergravity. By appropriately compactifying the ten-dimensional supergravity theory on S 5 , one obtains a consistent truncation to the five-dimensional Einstein-Maxwell-Chern-Simons theory [12,62,63], given by the action 9 where R is the Ricci scalar corresponding to the five-dimensional bulk metric G M N , and is the field strength tensor for the bulk Abelian gauge field A M . Also, κ 2 5 = 8πG N , and the cosmological constant is related to the AdS radius L via Λ = −6/L 2 . The coefficient of the Chern-Simons term, α, has the numerical value 2/ √ 3, and is a measure of the strength of the chiral anomaly in the boundary super-Yang-Mills theory. 10 The action eq. (2.1) is supplemented by boundary-and counter-terms [69][70][71][72][73][74], where u denotes the radial direction in Poincaré slicing. The metric γ µν is induced by the bulk metric G M N on the conformal boundary of AdS 5 , and R is the associated four-dimensional boundary Ricci scalar. Furthermore, K M N is the extrinsic curvature, given by with ∇ being the covariant derivative and n M being the outward pointing normal vector to the boundary ∂M. The equations of motion for the metric G M N and the gauge field A M read Here N M OP Q is the totally anti-symmetric Levi-Civita symbol in 4 + 1 dimensions, with txyzu = 1. The four-dimensional N = 4 supersymmetric SU (N c ) Yang-Mills theory carries a global SU (4) R-symmetry. From the AdS/CFT perspective, this is same as the SO(6) symmetry of the dual type IIB supergravity theory on AdS 5 × S 5 , corresponding to the isometries of S 5 . The effective theory eq. (2.1) can be obtained from supergravity by simply providing rotations or twists to the angular directions on S 5 and compactifying, as illustrated beautifully in [12]. In fact, their construction corresponds to 9 The consistent truncation of IIB supergravity on AdS5 × S 5 gives the bulk gauge field coupling to be g 2 = 2κ 2 5 L 2 . See for instance [12]. We have absorbed a factor of L in defining the bulk gauge field to make it dimensionless. 10 The presence of the Chern-Simons term is also required by supersymmetry in N = 2 gauged supergravity theory on AdS5, which has precisely the same action for its bosonic sector as the one in eq. (2.1) [64][65][66][67][68].
introducing a rotation in the diagonal U (1) of the maximal Abelian subgroup U (1) 3 of the R-symmetry group. For the super-Yang-Mills theory, this corresponds to considering states with a non-vanishing expectation value for the dual U (1) current J µ . This precisely is the charged super-Yang-Mills plasma whose transport properties are of our interest here. In particular, the equilibrium state of the plasma at a finite temperature and chemical potential, in the absence of background electric and magnetic fields, is holographically dual to the Reissner-Nordström black brane solution for the action eq. (2.1) [12,14], which in infalling Eddington-Finkelstein coordinates takes the form [9,58,59] Here, for brevity, we denote the null ingoing Eddington-Finkelstein time coordinate by t, while the radial coordinate is denoted by u. The horizon of the black brane is located at u = u h , and the temperature of the dual field theory is given by the black brane temperature, i.e. T = |f (u h )|/4π. The boundary is located at u = 0, where we impose asymptotically AdS 5 boundary conditions, i.e. f (0) = 1. Also, the bulk gauge field A M takes the form with other components vanishing. We fix our gauge by requiring that A t (u h ) = 0 and by imposing A u (u) = 0. In this gauge, the chemical potential of the boundary theory is simply the boundary value of with the gauge field given by From eqs. (2.6) and (2.8), it is straightforward to compute the temperature T and entropy density s, which are given by where A BH is the area of the black brane horizon.
As discussed in section 1, the thermodynamic susceptibilities f 1 , . . . f 7 are given by Kubo formulas involving equilibrium two-point functions of the energy-momentum tensor and the U (1) current (see table II). Now the energy-momentum tensor of the boundary theory is sourced by the bulk metric, while the U (1) current is sourced by the bulk gauge field. One can thus extract the two-point functions for the energy-momentum tensor and the U (1) current by studying perturbations to the background metric and gauge field configurations discussed above. For our purposes, it is useful to consider time-independent spatial fluctuations on top of the background solution (G M N , A M ), i.e. (2.12) to first order in a small parameter . To fix our gauge at the level of fluctuations, we work in the radial gauge by setting h M u , a u to zero. Due to the isotropy of the system, without loss of generality we have chosen the momentum to point along the z-direction, and will drop the subscript z for the remainder of this work i.e. k ≡ k z .

COMPUTING THE THERMODYNAMIC SUSCEPTIBILITIES
In this section, we compute the seven T-even second order thermodynamic susceptibilities f 1 , . . . f 7 for the N = 4 super-Yang-Mills plasma making use of the holographic perturbed Reissner-Nordström black brane model described in the previous section.
We may compute holographic two-point functions of the form G O a O c (k) in terms of the expectation value of the one-point function O a by exploiting the following relation to the first order in the fluctuations [9,75], where the perturbation δφ c is the source dual to the operator O c on the boundary and δ O a is the response for turning on this perturbation. For the two-point functions we would like to compute, O a ≡ T µν , J µ , and the sources correspond to perturbations in the bulk metric and the gauge field. The expressions for the renormalized one-point functions for T µν , J µ follow from the standard holographic renormalization procedure. The renormalized energy-momentum tensor may be extracted from the covariant expression [70] 2κ 2 5 T µν = lim Similarly, the one-point function of the renormalized current is given by [16,63,76] 11 Up to contact terms and terms from the holographic renormalization procedure, the two-point function reads where ∆ is the conformal dimension of O a . We obtain the value of the two-point function by computing the expectation value of O a (k) in the presence of a source term δφ c (k) for the operator O c and computing their ratio. Note that the asymptotic boundary expansion of the gauge field fluctuations read formally as a µ (u; k) ∼ a L µ + u 2 /L 2 a S µ (u; k) + a L µ k 2 log(u/L) , where the superscripts L and S denote the leading and subleading modes. According to the holographic dictionary, the expectation value of the dual operator is encoded in the subleading mode of the asymptotic boundary expansion. However, at the order u 2 , we observe logarithmic divergences which are removed by boundary counterterms. These inherently break conformal invariance and we have to specify the renormalization scale when regularizing the action (see [75,77] for more details). In our context of the second order thermodynamic transport coefficients, this leads to additional finite contributions to the renormalized expectation value at order k 2 and affects the current-current two-point functions G J t J t and G J x J x . The two-point functions involving the energy-momentum tensor are unaffected at order k 2 , since for them the logarithmic terms due to the Weyl anomaly are of order k 4 .
We compute the susceptibility f 1 analytically in appendix A -see eqs. (A.8) and (A.9). Working in units where we set g 2 = 2κ 2 5 = L = 1 for simplicity, we get the analytic result (3.5) Now, for a conformal theory, such as N = 4 super-Yang-Mills, the susceptibilities f 2 and f 7 are not independent of f 1 [39]. In fact, f 2 = 6f 1 and f 7 = −6∂f 1 /∂µ. These relations between the susceptibilities are required to maintain the Weyl invariance of the equilibrium generating functional W[g µν , A µ ] for the conformal theory. Using eq. (3.5), we thus get the following results for f 2 and f 7 , (3.7) The rest of the susceptibilities are computed numerically by means of pseudo-spectral methods as explained in appendix B. To compute the second order in k Green's function, we compute the Green's function for several small values of the momentum k and extract the k 2 dependence by fitting the numerical results to a second order polynomial (in k). Our final results for the seven T-even thermodynamic susceptibilities for the charged strongly-coupled N = 4 supersymmetric SU (N c ) Yang-Mills plasma in the limit of large N c , extracted via the holographic procedure described above, are depicted in figure 1.
Finally, we discuss the behaviour of the susceptibilities under a parity transformation, as well as for small chemical potentials (or high temperatures) |µ|/T 1, summarized in table III. The behaviour under a parity transformation matches exactly the P eigenvalues of the equilibrium invariants given in eq.  0 μ/T  μ/T  reaches its minima f 3 /T 2 ≈ −9.870 at around µ/T ≈ ±10.883, and successively changes sign for large |µ|/T . The transport coefficient f 5 reaches its symmetrically located extremas of f 5 /T ≈ ±1.296 at µ/T = ∓4.60. At around µ/T ≈ ±10.883, f 5 changes sign and starts to behave linearly in µ/T with f 5 /T ∼ ±2π + 0.4837µ/T . The slope is very similar to the behavior at small |µ|/T , even though its sign has now changed. It is interesting to note that f 5 /T changes sign for the exact same values of µ/T where the transport coefficient f 3 exhibits its minima. Note that at large µ/T the transport coefficient f 7 /T changes slope and behaves like f 7 /T ∼ 3.788 sign(µ/T ) + µ/T . Note that in strong magnetic fields, where B µ ∼ O(1) in the derivative counting scheme, the Kubo formula for f 5 is no longer second order but already has first order in derivative contributions. The first order Kubo formula for f 5 in the presence of a strong background magnetic field B 0 reads (compare with M 4 in [78] or M 5 in [16]) 8) and was computed in [16]. The authors found that where c 1 is a negative constant. This is in perfect agreement with our results (to zeroth order in the magnetic field).

Comments on vanishing anomaly α = 0
As discussed in section 2, the consistent truncation of type IIB supergravity leads to the Einstein-Maxwell-Chern-Simons action eq. (2.1), dual to the strongly coupled regime of the N = 4 super-Yang-Mills theory with a U (1) charge [12]. The parameter α, which is the coefficient of the Chern-Simons term in the bulk action eq. (2.1), takes on the value 2 √ 3 . This is related to the strength of the chiral anomaly C for the boundary super-Yang-Mills theory, There are no background electric or magnetic fields applied to the plasma in our setup, i.e. E, B = 0, and thus one would not see the anomaly in the equation for current conservation. However, it will be manifest in other physical quantities, such as triangle diagrams involving the anomalous current, as well as the equations of motion eqs. (2.4), (2.5). In the results we have presented in this section, figure  1, the effects of this anomaly enter via the equations of motion satisfied by the bulk perturbations.
It is still instructive to consider the case when one sets the anomaly to vanish, α = 0. The resulting bulk description in terms of the Einstein-Maxwell action is no more dual to the N = 4 super-Yang-Mills plasma with a U (1) charge on the AdS 5 boundary. It can rather be considered dual to the strongly coupled regime of some non-anomalous U (1) charged conformal theory living on the boundary of AdS 5 . It is important to note that the U (1) symmetry in this case is vector in nature, and hence the charge conjugation and parity transformation properties will be different from the ones in table I. When α = 0, the Chern-Simons term couples the sectors with the fluctuations {a x , h tx , h xz } and {a y , h ty , h yz }, while all other sectors stay unaffected by the chiral anomaly. Thus, from the Kubo formulas in table II, it is clear that the transport coefficients f 1 , f 2 , f 6 and f 7 do not get influenced by the chiral anomaly, and remain the same for the non-anomalous case α = 0. However, the transport coefficients f 3 , f 4 , f 5 get affected, and we present their behaviour when α = 0 in figure 2 below. Note that the transport coefficient f 4 vanishes identically without the Chern-Simons term.   fig. 1.  Similar to the case with the chiral anomaly, the second order transport coefficient f 5 may be computed in the presence of strong magnetic fields by the first order Kubo formula eq. (3.8) when α = 0. The authors of [16,79] found that where c 2 is a positive constant. We again find perfect agreement with our results (to zeroth order in the magnetic field).

Comments on the transport parameter κ at µ = 0
In the uncharged case µ = 0, the susceptibility f 1 for the N = 4 super-Yang-Mills plasma was computed in [21,22]. In fact their result is quoted in terms of the second-order transport parameter κ which is defined in the Landau-Lifshitz frame, 12 and is related to our f 1 via a standard frame transformation (see appendix A of [39] perfectly matches our analytical calculation of f 1 in eq. (A.8).
The transport coefficient κ was also computed using lattice gauge theory techniques for the SU (N c ) Yang-Mills plasma, and it was found that for N c = 3 [80] κ = 0.36(15)T 2 . (3.12) We can also make a numerical estimate of κ for QCD using our results. We first match the entropy density of our holographic model to that of lattice QCD as explained for e.g. in [81]. The Stefan-Boltzmann value of the entropy density of three flavor QCD is Note that this value is only reached at asymptotically high temperatures. We can match this to the black hole entropy of our gravitational model at vanishing chemical potential, where we have used eq. (2.10) in the last equality to express u h in terms of the temperature T . If we want to compare eq. (3.13) with N = 4 super-Yang-Mills plasma in the limit of infinite 't-Hooft coupling λ, eq. (3.14), we have to take into account a factor of 3/4 between the entropy densities [9,82], i.e. 3s SB /4 = s BH , yielding

15)
12 Note that κ here should not be confused with the coefficient κ5 appearing in the action eq. (2.1) for the holographic model, i.e. κ 2 5 = 8πGN .
Using this result for κ 5 in eq. (3.11) gives at µ = 0 an estimate of This is within the error range of the lattice result eq. (3.12). The difference is expected as the result of [80] is for pure Yang-Mills plasma, whereas in QCD one also has to take the quarks and their flavours into account. 13 Additionally, the authors of [80] quote a much smaller value for κ as the AdS/CFT result since they use a different entropy density too. Our computations thus provide an interesting comparison between holography and lattice gauge theory results. 14

DISCUSSION
In this short article, we computed the seven T-invariant second order thermodynamic susceptibilities f 1 , . . . f 7 for the charged strongly coupled N = 4 super-Yang-Mills plasma and illustrated their dependence on the dimensionless parameter µ/T (see fig. 1). To compute the susceptibilities, we employed Kubo formulas that gave these susceptibilities in terms of equilibrium two-point functions of the energy-momentum tensor and the U (1) current. We computed the two-point functions holographically, by studying fluctuations to the asymptotically AdS 5 Reissner-Nordström black brane geometry. As also mentioned in the Introduction, studying the transport properties of the N = 4 super-Yang-Mills plasma via holography provides insights into the expected behaviour of the quark-gluon plasma produced in heavy ion collisions.
Interestingly, we also computed the second order transport coefficient κ analytically, which is related to the susceptibility f 1 , and measures the response of the fluid to the presence of background curvature. Using the analytic result, an estimate could be made for the value of κ for QCD. Our calculations give a result in proximity of the lattice gauge theory computations for Yang-Mills plasma, giving further credibility to the idea that transport properties of strongly coupled QCD matter can be captured via holographic techniques modeling the super-Yang-Mills plasma.
A few comments are in order. The super-Yang-Mills theory is a conformal field theory, and has a trace anomaly in the presence of external electromagnetic fields, given by 15 Here we have used the notation f n = T ∂f n /∂T + µ∂f n /∂µ. For us, f 4 , f 6 = 0. However, since we apply no external electromagnetic fields to the super-Yang-Mills plasma i.e. E = B = 0, the trace anomaly must vanish. That this indeed happens is very easy to see by computing the stress tensor components and evaluating the trace, since we know the background solution analytically.
We computed the behaviour of the second order susceptibilities as a function of µ/T at strong coupling via holographic techniques. At the other extreme, when the plasma can be considered 13 For N f = 0, we find 2κ 2 5 = 15 L 3 π 2 /2 and thus κ = 2 T 2 /15 ≈ 0.1333T 2 . 14 Another estimate may be found in [83,84] for an SU (Nc) gauge theory, κ = (N 2 c − 1) T 2 /18 = 0.44 T 2 , where we used Nc = 3 in the last step. Furthermore, for 3 flavour QCD the authors in [85] found κ = 0.5694 T 2 . 15 There are further contributions to the trace anomaly if the conformal theory is put on a curved background, proportional to R 2 . However, these are fourth order in derivatives and irrelevant for the present discussion.
weakly coupled, it is much simpler to compute the same susceptibilities. The matter content of the The quark-gluon plasma is produced at heavy ion colliders such as the LHC and RHIC in an outof-equilibrium state due to thermal fluctuations. Needless to say, there is a very rich spectrum of physical effects that arise as a consequence of these out of equilibrium fluctuations, and the second order equilibrium transport properties discussed here will receive important corrections due to the fluctuation effects. Large magnetic fields are also produced during heavy ion collisions [86][87][88] making the background anisotropic, and their presence can also significantly affect the transport and relaxation phenomena occurring in the quark-gluon plasma [16,81,89,90]. Also, the remaining two second order susceptibilitiesf 8 , f 9 -which are T-odd and were not the subject of study in the current paper, can in principle be determined using equilibrium three-point functions of the energy-momentum tensor and the current [39], although closed form expressions of Kubo formulas for them are not yet known. 16 All these possible directions to extend the present work are extremely interesting -we leave them out here for near future investigations.
Note that the second order Green's function is unbothered by the logarithmic term scaling with k 4 . We can compute the Green's function G T xy T xy analytically by expanding the fluctuation h xy to second order in k, i.e.
and solving the equation of motion order by order in k, with the boundary condition h (0) xy (0) = 0. With this, the only solution regular at the horizon is Expanding this solution at the boundary u = 0 yields which perfectly matches the asymptotic expansion in eq. (A.3). For c 1 = 1 we find for the expectation value The explicit dependence of f 1 on µ and T can now easily be obtained by eliminating u h from the expression above using eq. (2.10). This leads to the following result for f 1 , This result is presented in eq. (3.5) of the main text, in the simplified set of units g 2 = 2κ 2 5 = L = 1.

Appendix B: Numerical methods
We determine the two-point functions required to compute the second order thermodynamic transport coefficients numerically by means of a pseudo-spectral method. In this appendix, we briefly explain the numerical techniques used; for a detailed introduction see for example [92]. Our discussion follows [91,[93][94][95]. The basic idea of spectral methods is to account for the radial dependence of the unknown functions in the ordinary differential equations governing the fluctuations by expanding them in terms of a (truncated) series in the basis functions. The basis functions are analytically known and so are their radial derivatives. To set up the spectral method, we choose a discrete, finite collocation grid, a set of basis functions {T k } and a finite number N at which we truncate the series in the basis functions, Spectral methods are global methods, which means that we take the function values on the whole domain into account when we compute numerical derivatives. This is in strong contrast to finite difference methods, which require only the values from the neighboring grid-points. We discretize the equations of motion by a Chebyshev-Lobatto grid with N grid-points, In this work, we choose the Chebyshev polynomials T k = cos(k arccos(2u − 1)), where u ∈ [0, 1] as basis functions. The derivatives acting on the functions simply translate into derivatives of the basis functions, With the derivative matrixD, we can compute the discrete derivatives on our grid. Since spectral methods are global methods (in contrast to finite differences), the information on each grid-point is needed for the derivative (instead of just the neighborhood). In the context of holography, this is very convenient since it allows us to impose boundary conditions on both ends of the domain simultaneously. We can determine the coefficients of the derivatives simply from the original coefficients by multiplying the original coefficients with a derivative matrix. In this work, we do not work directly with the coefficients but with the values of the functions on the discretized grid (this is referred to as pseudo-spectral method). We can translate from coefficients to the values of the functions with whereK is the discretized matrix (see step 1 above), X a M N dimensional vector with the discretized functions and φ the right hand side with the values independent of the functions.

Check the convergence of the numerical solution.
To confirm the accuracy of our numerics, we perform several checks. First of all, we check that the Chebyshev coefficients of the numerical solution drop sufficiently fast and below our threshold, as presented in the left plot of figure 3. Furthermore, we check that the equations of motion and constraint equations are satisfied on an equidistant grid. Finally, we show that the functions converge to an accurate solution (right side of figure 3). For the accurate solution we choose N = 130 since the change of the solution for taking even more grid-points is smaller than 10 −40 .
To show that the numerical solutions converge to the accurate solution, we compute numerical solutions with less grid-points. We compare these "less accurate" solutions to the solution with N = 130 on an equidistant grid with z = i/129 for i = 0, . . . , 129 and compute the biggest deviation. We can see that solutions converge exponentially towards the N = 130 solution.