Hydrodynamic transport coefficients for the non-conformal quark-gluon plasma from holography

In this paper we obtain holographic formulas for the transport coefficients $\kappa$ and $\tau_\pi$ present in the second-order derivative expansion of relativistic hydrodynamics in curved spacetime associated with a non-conformal strongly coupled plasma described holographically by an Einstein+Scalar action in the bulk. We compute these coefficients as functions of the temperature in a bottom-up non-conformal model that is tuned to reproduce lattice QCD thermodynamics at zero baryon chemical potential. We directly compute, besides the speed of sound, 6 other transport coefficients that appear at second-order in the derivative expansion. We also give an estimate for the temperature dependence of 11 other transport coefficients taking into account the simplest contribution from non-conformal effects that appear near the QCD crossover phase transition. Using these results, we construct an Israel-Stewart-like theory in flat spacetime containing 13 of these 17 transport coefficients that should be suitable for phenomenological applications in the context of numerical hydrodynamic simulations of the strongly-coupled, non-conformal quark-gluon plasma. Using several different approximations, we give parametrizations for the temperature dependence of all the second-order transport coefficients that appear in this theory in a format that can be easily implemented in existing numerical hydrodynamic codes.


Introduction
After the discovery of the quark-gluon plasma (QGP) in ultra-relativistic heavy ion collisions [1], a lot of effort has been put towards understanding how the spatial anisotropies present in the initial state are converted into the final flow of hadrons. Relativistic dissipative hydrodynamics has played an important role in our current view of the complicated spacetime evolution of the QGP formed in heavy ion collisions (for a recent review see [2]). The overall picture that is consistent with experimental data is that before hadronization the QGP evolves in time and space as a relativistic fluid with minimal dissipative effects. Indeed, current estimates [2] for the shear viscosity to entropy density ratio, η/s, of the QGP obtained by comparison to data are in the ballpark of the very small value η/s = 1/(4π) [3] found in a broad class of strongly-coupled non-Abelian plasmas described by the gauge/gravity duality [4][5][6]. This suggests that the gauge/gravity duality may be useful for the study of the non-equilibrium properties of strongly interacting plasmas that are similar (if not quantitatively at least qualitatively so) to the QGP and, thus, several applications have been studied over the last years (see, for instance, the review [7]).
In fact, we shall show in this paper that a simple bottom-up holographic model that is able to describe (some of) the thermodynamic properties of the QGP 1 near the crossover phase transition [8] can be instrumental in providing estimates for the temperature dependence of a large number of second order transport coefficients that appear in consistent theories of (non-conformal) dissipative relativistic hydrodynamics. This paper is organized as follows. In the next section we shall review the dissipative hydrodynamic theory obtained at second-order in gradients in the case of a non-conformal, relativistic plasma (in the absence of conserved charges such as the baryon number) in a curved spacetime [9]. At this order in the gradient expansion, there are 17 coefficients (besides the speed of sound) that may possess some nontrivial temperature dependence (especially near the phase transition) 2 . In Section 3 we present our method to holographically compute the second order transport coefficients κ and τ π in strongly-coupled plasmas that are described by a bulk action including the metric and a dynamical scalar field (see also Appendix A). In Section 4 we give the details of the bottom-up holographic model we use and fix its parameters through a comparison to lattice QCD thermodynamics. In Section 5 we compute the temperature dependence of several transport coefficients for this holographic bottom-up model. In Section 6 we use the 2nd-order gradient theory defined in Section 2 to write an Israel-Stewart-like hydrodynamic theory in flat spacetime with 13 transport coefficients that could be implemented in numerical hydrodynamics. Also, a guide to the temperature dependence of the several 2nd-order transport coefficients considered in this paper (given in terms of fitting functions that could be easily used in numerical hydrodynamics) can be found in Appendix B. Furthermore, in Appendix C we perform a linear stability analysis around the static equilibrium for the non-conformal, 2nd order gradient expansion theory discussed in Section 2. Our conclusions and outlook can be found in Section 7.
The reader that is mostly interested in the hydrodynamic discussions and the specific temperature dependence of the transport coefficients (shown in Figs. 3 to 12) may want to focus on Sections 2, 6, and Appendices B and C. The other sections are devoted to more detailed calculations involving the gauge/gravity duality.
Throughout this paper we use natural units c = = k B = 1 and a mostly plus metric signature. Also, we use capital Latin indices to denote the bulk coordinates x M = (t, x, y, z, u) while Greek indices x µ = (t, x, y, z) denote 4-dimensional coordinates. 1 As it will be clear later in the paper, this model cannot describe aspects of QCD that are directly related to chiral symmetry. 2 We note that many more coefficients would be needed in the case where spatial isotropy in the equilibrium state is lost, as it occurs in anisotropic hydrodynamics [10][11][12][13][14][15][16] and in fluids in the presence of strong magnetic fields (see, for instance, [17]). However, these interesting generalizations will not be pursued here.

Second-order non-conformal hydrodynamics via the gradient expansion
Relativistic dissipative hydrodynamics can be viewed as a type of effective theory for the long wavelength, low frequency behavior of an interacting system at finite temperature and/or chemical potential [18,19]. Such an effective theory may be constructed at weak coupling in the case of a dilute gas [20][21][22][23][24][25][26] whose microscopic behavior can be described by a Boltzmann-like equation for the system's effective quasi-particles [27,28]. On the other hand, at strong coupling the fluid-gravity correspondence [29] provides an adequate framework to study the effects of spacetime gradients in a strongly coupled fluid.
In general, dissipation is included directly at the level of the equations of motion 3 , which in the absence of conserved charges, correspond solely to the conservation of energy and momentum ∇ µ T µν = 0 , where ∇ µ is the covariant spacetime derivative in a curved 4-dimensional spacetime described by a metric g µν , and T µν is the expectation value of the system's energy-momentum tensor operator. We shall consider here matter described by a relativistic quantum field theory giving the equation of state, P = P (ε), where ε and P are the local energy density and pressure of the fluid, respectively. The equation of state gives rise to the speed of sound in the fluid, c s = dP/dε. The basic idea of the gradient expansion is that the macroscopic degrees of freedom in the long wavelength, low frequency limit are only the local energy density, ε, 4-velocity, u µ , metric, g µν (in curved spacetime), and their gradients. In fact, the energy-momentum tensor can be generically decomposed as T µν = ε u µ u ν + P ∆ µν + π µν + ∆ µν Π (2.2) where the flow obeys u µ u µ = −1, and ∆ µν = g µν + u µ u ν is a local projection operator transverse to the flow. Note that such a decomposition inherently assumes that there is a well defined local rest frame (for examples of quantum field theories in far from equilibrium conditions without a local rest frame see Ref. [33]). Dissipation generally appears due to a nonzero shear stress tensor π µν = ∆ µναβ T αβ , (2.3) which is transverse to the flow 4 , symmetric, and traceless due to the definition of the tensor projector The last term in (2.2) denotes the dissipative contribution to the energy-momentum tensor with non-vanishing trace, Π, called the bulk viscous pressure. In terms of (2.2), one can 3 For recent discussions including attempts to formulate dissipative hydrodynamics in terms of an effective action see, for instance, [30][31][32]. 4 Note we use the Landau frame, i.e., uµT µν = −ε u ν [18].
show that the conservation of energy and momentum become Dε + (ε + P + Π)θ + 1 2 π µν σ µν = 0 , (ε + P + Π)Du µ + ∇ µ ⊥ (P + Π) + ∆ µ ν ∇ α π αν = 0 , (2.5) where D = u µ ∇ µ is the comoving derivative, ∇ α ⊥ = ∆ αβ ∇ β is the derivative transverse to the flow, θ = ∇ µ u µ is the scalar expansion rate, and σ µν = 2∆ αβ µν ∇ α u β is the shear tensor. The energy conservation equation can be written in terms of the equilibrium entropy density, s = (ε + P )/T , as follows In the gradient expansion approach, since only ε and u µ are the hydrodynamical variables, the dissipative components π µν and Π must be expressed solely in terms of derivatives of these quantities. To first order in gradients, this can be easily done and one finds where η is the shear viscosity and ζ is the bulk viscosity, respectively. Note that in this case the second law of thermodynamics in (2.6) imposes that η, ζ ≥ 0. If one uses the expressions above for the dissipative contributions in T µν , the conservation equations represent the relativistic extension of the well-known Navier-Stokes (NS) equations [18] 5 .
In kinetic theory, the transport coefficients η and ζ are proportional to their corresponding mean free paths, 6 . One can now see how the power counting scheme adopted in the gradient expansion works. Since is a microscopic scale and ε and u µ are taken to be slowly varying functions of time and space, one can associate with their gradients a characteristic (macroscopic) length scale ∼ 1/L macro such that /L macro 1. Therefore, terms such as those in (2.7) are taken to be of order 1 in the so-called Knudsen number K n ∼ /L macro 7 . Clearly, the continuous description of the system as a fluid hinges on the assumption that K n is sufficiently small. However, given that dissipation only appears at order 1 in this expansion, one may also entertain the case in which K n is still sufficiently small to ensure a well-defined continuous description but the flow is such that higher order terms may be taken into account. Nevertheless, one should keep in mind that the radius of convergence of the gradient series is limited by the first nonzero non-hydrodynamical quasinormal mode, as recently shown in [38] in the context of strongly coupled gauge theories and discussed earlier in [39] in the context of kinetic theory. 5 Another way to understand how dissipation appears is to notice that, for instance, in this NS fluid the inclusion of πµν breaks the time reversal invariance present in the ideal fluid equations of motion. However, it is possible to find nontrivial fluid patterns involving second order gradients where πµν is nonzero but time reversal is not broken -see [34,35]. 6 Note that the mean free path for bulk viscosity is different than that for shear viscosity [36]. However, for simplicity, we shall denote any mean free path here by . 7 We remark that the Knudsen "number" is actually a field since it depends on the spacetime coordinates.
Moreover, in general one may consider several types of Knudsen numbers associated with different properties of the flow, see for instance [37].
For instance, in the early stages of a ultrarelativistic nucleus-nucleus collision [40], the local energy density and flow are expected to have sizable gradients and corrections of second order in K n may be relevant at that stage of the QGP evolution. Also, as emphasized in [37], in collisions involving smaller systems such as proton-nucleus collisions at the LHC, the need for higher order Knudsen number corrections may be even more pressing. Therefore, it is reasonable to ask what are the expressions for π µν and Π including O(K 2 n ) terms. Generalizing the previous analysis involving 2nd order terms in a conformal fluid done in [41], Romatschke proposed in [9] the following expansion for the dissipative parts of a non-conformal relativistic fluid in curved spacetime valid at O(K 2 n ) and where R λ µσν is the Riemann tensor, R µν = R λ µνλ is the Ricci tensor, and R = g µν R µν is the Ricci scalar [42]. Moreover, we have also defined the vorticity tensor Ω µν = 1 2 ∇ ⊥ µ u ν − ∇ ⊥ ν u µ and the usual notation B µν = ∆ µν αβ B αβ for the traceless, symmetric, and transverse part of a second rank tensor B µν . Eqs. (2.8) and (2.9) can then be used in the conservation equations (2.5) to define the equations of motion for a non-conformal fluid in a curved spacetime valid at second order in gradients 8 .
One can see that besides the speed of sound squared c 2 s = dP/dε that is already present in ideal hydrodynamics and the coefficients η and ζ that appeared at 1st order, there are now altogether 15 new transport coefficients that appear at second order in the gradient expansion. Following [43], one may distinguish these coefficients by separating out those that are of thermodynamical origin and those that are not. The set of coefficients κ, κ * , λ 3 , λ 4 , ξ 3 , ξ 4 , ξ 5 , ξ 6 can be determined via Kubo formulas involving only equilibrium quantities and Euclidean two-and three-point functions of the energy-momentum tensor components and, thus, they are of thermodynamical origin being suitable for lattice calculations [43]. However, the other coefficients η, ζ, τ π , τ * π , τ Π , λ 1 , λ 2 , ξ 1 , and ξ 2 are associated with quantities that define the dissipative properties of the theory (for instance, η is proportional to the imaginary part of a retarded Green's function) such as σ µν and are, thus, of dynamical origin [43]. 8 Due to the fact that spacetime covariant derivatives generally do not commute in a curved spacetime, even when the metric is not dynamical (though still nontrivial) quantities such as R λ µνα , Rµν , and R are expected to appear in the equations of motion for ε and uµ. Also, we note that the expressions in (2.8) and (2.9) are in agreement with the corresponding terms (in flat spacetime) at O(K 2 n ) found in [25].
There are several interesting points regarding these second order terms. First, as discussed in [9,[43][44][45][46] not all of these coefficients are independent and we shall get back to this point in Section 5 where we compute several of these coefficients in the holographic model defined in Section 4. Second, for a conformal plasma only the terms that transform homogeneously under a Weyl transformation 9 are present [41] and this implies that κ * , τ π * , λ 4 , ζ, τ Π , ξ 1,2,3,4,5,6 vanish in the conformal limit (therefore, in this case, Π = 0). Third, even though κ, κ * , ξ 5 , and ξ 6 do not contribute to the equations of motion in flat spacetime, they do contribute to the Kubo formulas for the other coefficients relevant to flat spacetime hydrodynamics and should then be taken into account, as discussed in [41,43]. Moreover, in a non-conformal fluid all of these coefficients may be nontrivial functions of the temperature, i.e., η = η(T ), especially near a phase transition (even if of the crossover type). This shows how challenging numerical 2nd order hydrodynamics can be if all of these temperature dependent transport coefficients are taken into account. Clearly, as long as these 2nd order gradient terms can be taken as small corrections, their effect should be under control. However, it is not clear at the moment if this is indeed the case for the type of event-by-event hydrodynamic simulations fed by the complicated initial conditions that describe the early stages of a heavy ion collision [40].
Furthermore, unfortunately the equations of motion defined by (2.8) and (2.9) cannot be directly implemented in numerical hydrodynamic codes because, for instance, they are linearly unstable against small fluctuations around the static equilibrium (see Appendix C). In fact, in Section 6 we propose another second order theory that is more suitable for numerical investigations using the current relativistic hydrodynamic codes. This theory can be considered as a type of "UV completion" of the 2nd order theory in (2.8) and (2.9) in the sense that it possesses a well defined (and causal) UV behavior (at least in the linear regime) but its long wavelength, low frequency asymptotic hydrodynamical solution is identical to the one obtained by (2.8) and (2.9) with the same transport coefficients. We shall discuss these points in detail in Section 6.
3 Holographic calculation of the 2nd order coefficients κ and τ π In this section we discuss how we are going to evaluate the 2nd order hydrodynamic transport coefficients κ and τ π , originally defined in [41], using the gauge/gravity duality [4][5][6]. First, we shall consider the case where the action in the bulk corresponds to 5-dimensional pure gravity (with a negative cosmological constant) and then later we will generalize this discussion for the case where the bulk action contains a dynamical scalar field. 9 Under a Weyl transformation, gµν → e −2ω gµν , where ω is an arbitrary (positive-definite) scalar function.
In a conformal plasma, T µν → e 6ω T µν while the temperature transforms as T → e ω T , see [41].

Renormalized pure gravity action
In order to compute these transport coefficients, it is sufficient to consider only the small frequency ω, small momentum q limit of the (xy, xy)-component of the retarded propagator of the stress-energy tensor of the boundary quantum field theory, which can be written 10 in momentum space as follows [9] G xy,xy where κ * is a second order hydrodynamic coefficient which is non-vanishing only for nonconformal fluids (see Eq. (2.8)), being related to κ by the following constraint [43][44][45][46] According to the gauge/gravity duality dictionary, the stress-energy tensor of the plasma is sourced in the partition function of the boundary quantum field theory by the boundary value of a classical metric perturbation placed over an asymptotically AdS bulk. In (3.1), it was assumed that the xy-component of the metric perturbation has no dependence on the x-and y-directions, such that the 4-momentum of its Fourier mode is given by k µ = (−ω, q x , q y , q z ) = (−ω, 0, 0, q). From (3.1) and (3.2), one obtains the following Kubo formulas for κ and τ π (valid for both conformal and non-conformal fluids) 11 In order to compute κ and τ π from (3.3) and (3.4) via holography, we need to evaluate the renormalized on-shell bulk action and extract from it the retarded graviton propagator by following the prescription proposed in [48], which was later justified and generalized in [49][50][51]. In the case of pure gravity, the regularized action is defined by the sum of the Einstein-Hilbert (EH) action, the Gibbons-Hawking-York (GHY) action 12 [52,53] and the counterterm action 13 [41,54]

5)
10 This retarded Green's function is given by G xy,xy x),T xy (0, 0)] . 11 We thank G. Moore and K. Sohrabi for discussions about these coefficients. 12 The GHY action needs to be added to the EH action in order to properly define the variational problem with Dirichlet boundary conditions for the metric tensor when the spacetime manifold has a boundary. 13 The counterterm action is obtained through the holographic renormalization procedure [56][57][58][59]72].
where G 5 is the five dimensional Newton's constant and the cosmological constant enforcing the existence of asymptotically AdS 5 geometries with radius L as solutions of Einstein's equations is given by Λ = −6/L 2 . The metric induced at the boundary is given by γ M N = g M N −n MnN , wheren M is an outward directed unit vector orthogonal to the boundary. In a coordinate chart where the boundary of the asymptotically AdS space is at the value u = 0 of the radial coordinate, this is given byn the metric induced at the regularizing boundary surface 14 u = can be simply written as (3.6) The extrinsic curvature of the boundary of an asymptotically AdS space is given by (see [55] for a review) where the prime denotes a derivative in the radial direction. The boundary sectional curvature tensor and the boundary sectional curvature scalar are defined, respectively, by [56] where R µν (γ) and R(γ) are the Ricci tensor and the Ricci scalar evaluated using the induced metric on the boundary surface.
Let us now consider a small perturbation in the metric, h M N (u, t, z), placed over a diagonal and isotropic gravitational background 15 , g M N (u), which is assumed to be asymptotically AdS. Since we are only interested in the (xy, xy)-component of the retarded propagator of the boundary stress-energy tensor, as discussed in [54], if one fixes the gauge defined by the subsidiary condition h M u = 0, one only needs to consider h xy (u, t, z) = 0 and set to zero all the other components of the metric perturbation since the linearized equation of motion for the xy-perturbation decouples from the other components of the metric perturbation in this gauge. Therefore, we consider the following disturbed line element where we defined φ(u, t, z) = h x y (u, t, z) and, with the sign convention used in (3.9), we are considering the tt-component of the metric to be −g tt (u), with g tt (u) > 0.
14 Strictly speaking, this quantity diverges at u = 0. Therefore, in order to regularize quantities of interest, we introduce an ultraviolet cutoff 1 for the radial coordinate near the boundary, which must be taken to zero at the very end of the calculations after all the divergent terms in the on-shell gravity action have been canceled. 15 The index (0) refers to the undisturbed background.
Let us now explicitly show that the EH action for the disturbed metric, up to second order in the metric perturbation φ, can be written as the action for a massless scalar field in the undisturbed background g (0) M N plus total derivatives which shall contribute to the final expression for the regularized on-shell boundary gravity action 16 (3.10) Since we are assuming that the background metric is a solution of Einstein's equations, we can make use of these equations to write down two useful relations for our purposes. The first one comes from contracting the background metric with its equation of motion The second useful relation comes from substituting (3.11) into the xx-component of the Einstein's equations for the background metric Substituting (3.11) into (3.10), plugging the result into the EH action and integrating by parts the term proportional to φ∇ 2 (0) φ, we obtain, up to O(φ 2 ) where V 4 = ∂M 5 d 4 x is the 4-volume of the boundary, u = u H is the position of the background black hole horizon in the radial coordinate 17 and where we used relation (3.12). Substituting (3.14) into (3.13), we obtain, up to O(φ 2 ) (3.15) 16 We thank R. Critelli for discussions concerning this derivation. 17 If the background has no event horizon (or some kind of infrared wall), then one must take uH → ∞.
From (3.15) we see that, as stated before, the bulk term of the disturbed EH action up to second order in the metric perturbation φ corresponds to the action for a massless scalar field in the undisturbed background g (0) M N , as is well-known. Hence, the linearized equation of motion for the mixed metric perturbation φ is just the massless Klein-Gordon equation in a curved background Defining the Fourier transform 18 as one finds in momentum space Integrating by parts the bulk piece of (3.15) and making use of (3.16), we put the EH action on-shell up to O(φ 2 ) where, by following the prescription proposed in [48] for calculating the retarded propagator of the metric perturbation, we discarded the horizon contribution coming from the radial integration and took into account only the boundary contribution for the on-shell EH action.
Now we add (3.19) to the contributions coming from the disturbed GHY and counterterm actions evaluated up to O(φ 2 ) to find the following expression for the total regularized on-shell boundary action in momentum space 19 18 For the sake of notation simplicity, when necessary, we distinguish the position-space scalar field, φ(u, t, z), from its Fourier transform, φ(u, ω, q), only by their arguments. Notice also that, since φ(u, t, z) is real-valued, we obtain the following reality condition in momentum space: φ(u, −ω, −q) = φ * (u, ω, q). 19 This is a rather lengthy calculation but it can be straightforwardly done with the help of a symbolic mathematical software such as Wolfram's Mathematica [60]. We used the EDCRGTC code [61] written for Mathematica in order to deal with the lengthy tensor manipulations involved in the course of these calculations.
One can formally solve (3.18) in terms of the components of the undisturbed background metric, g M N (u), and the boundary value of the metric perturbation, ϕ(ω, q), by employing a perturbative expansion in ω T and q T . Such a solution up to O(ω 2 , q 2 ) is derived in details in Appendix A and the results near the boundary read . The renormalized on-shell boundary gravity action is defined by For completeness, let us also review some useful formulas for calculating the pressure, the entropy density, and the shear viscosity of the boundary plasma using the gauge/gravity duality. From (3.1), we see that once we have extracted the retarded propagator from the renormalized on-shell boundary action by following the prescription proposed in [48] S ren = lim one can obtain the pressure as follows Since the perturbation-independent part of the on-shell boundary action gives −P V 4 [41], we can also calculate the pressure by using the following alternative formula One can evaluate the entropy of the plasma by using the Bekenstein-Hawking's relation [62,63] where the "area" of the horizon is given by 28) and the entropy density reads From (3.1), one also obtains the following Kubo's formula for the shear viscosity Application: Thermal N = 4 Supersymmetric Yang-Mills (SYM) In this subsection we review the calculations for SYM at finite temperature [41] whose gravity dual is defined over an AdS 5 -Schwarzschild background with metric components The relation between G 5 and the number of colors of the strongly coupled SYM plasma reads [64] By substituting (3.31) and (3.32) into (A.5) and (3.29), one obtains, respectively Now we use (3.31), (3.32), and (3.33) to calculate (3.22) and, by substituting the result into (3.20) and then evaluating (3.23), we obtain the renormalized on-shell boundary action up to O(φ 2 , ω 2 , q 2 ) for SYM and by using (3.25) (or also (3.26)), we obtain the pressure of the conformal stronglycoupled SYM plasma The first order hydrodynamic transport coefficient η for the SYM plasma is found by From (3.33) and (3.37), one obtains the famous ratio which is actually valid for a broad class of gravity duals [65,66]. The second order hydrodynamic transport coefficients κ and τ π for the SYM plasma are found by substituting These results for the conformal SYM plasma were originally obtained in [41].

General holographic formulas for κ and τ π in Einstein+Scalar gravity duals
In this section we derive general formulas for κ and τ π which are valid also for systems where the metric couples to matter fields in the bulk (we shall focus here on the case where there is a scalar field in the bulk).
For geometries corresponding to solutions of the field equations which take into account the backreaction of matter fields in the bulk, the pure gravity regularized action (3.20) shall feature, in general, temperature-independent divergences as one takes the limit → 0. For instance, this happens in Einstein+Scalar models [67][68][69][70][71]. For this kind of system, the general procedure of holographic renormalization is discussed in [57]. Here, since we are only interested in evaluating hydrodynamic transport coefficients 20 , instead of dealing with a more complicated regularized action, we are going to apply a physical prescription which will allow us to obtain quite general formulas for κ and τ π through a simple procedure. This prescription is based on three main facts: 1. The equation of motion for the xy-component of the metric perturbation depends only on the background metric [54,65] and, therefore, the solution (3.22) remains valid also for Einstein+Scalar actions.
3. The coefficients κ and τ π η, which appear in the gradient expansion, must vanish at sufficiently low temperatures (as one approaches the vacuum).
Let us first discuss the coefficient κ in (3.3). We begin by defining the following regularized quantity (see eqs. (3.20) and (3.24)) and then we take the following expression, which is free of ultraviolet divergences where T high is a temperature that is sufficiently large so that we are near the ultraviolet fixed point and the temperature-dependent part 21 of κ (T high ) approaches κ SYM (T high ), u high = u H (T high ), and κ 0 is a constant to be subtracted (generally by numerical inspection) in order to ensure that κ(T min ) = 0 where T min is the lowest temperature considered in our numerical calculations (T min ∼ 10 MeV). From (3.32) and (3.39), one finds Analogously, for τ π in (3.4) we begin by defining and then we evaluate the UV finite expression

46)
21 Which is independent of the ultraviolet cutoff as mentioned before. where and Ω 0 is a constant that is subtracted to ensure that τ π (T min )η(T min ) = 0.
Eqs. (3.42), (3.44), and (3.46) can be used to compute the transport coefficients κ and τ π for Einstein+Scalar actions. In the next section we employ them to numerically evaluate κ and τ π in a holographic model for (2 + 1)-flavor QCD at finite temperature and zero baryon density proposed by Gubser and Nellore [70].

Thermodynamics of the Einstein+Scalar holographic model
In this section we shall briefly review a bottom-up holographic model [70] which is built to provide a quantitative description of the thermodynamical properties of (2 + 1)-flavor QCD as seen on the lattice 22 . In the next section we will apply the results of the preceding sections to compute all of the first order transport coefficients and most of the second order ones, presenting estimates for these coefficients that may be relevant for the hydrodynamic simulations of the QGP. The Einstein+Scalar holographic model has been successfully used in several other works to understand the temperature dependence of quantities evaluated near the crossover phase transition of the QGP such as the bulk viscosity [75], the expectation value of the Polyakov loop [76,77], the energy loss of heavy and light quarks [78][79][80], the electric conductivity [81], and the Debye screening mass [82].
The bulk action for the Einstein+Scalar model considered in [70] is given by where Φ is the scalar field and V (Φ) is the scalar potential, which we shall specify below.
We note that this model cannot be used to describe aspects of the QCD phase transition that are directly related to chiral symmetry 23 . For the purpose of solving the equations of motion, we use the following Ansatz for the metric (in the so-called Gubser gauge) [70] (see also [82] for more details) where the scalar field itself is considered as the radial coordinate. In order to have a black brane background, we require that h(Φ) has a simple zero at Φ = Φ H . We also impose that the metric is asymptotically AdS 5 and that Φ is associated with a relevant operator 22 For holographic studies of the deconfinement phase transition see, for instance, [73,74]. 23 For holographic bottom-up models that deal with effects of chiral symmetry see, for instance, [83][84][85].
in the gauge theory [75], which implies that Φ → 0 at the boundary. The entropy density is given by the Bekenstein-Hawking formula (3.29) and the corresponding temperature of the black brane is given by (A.5) The scalar potential V (Φ) is chosen in order to mimic the thermodynamics of (2 + 1)flavor QCD at zero baryon chemical potential as calculated on the lattice [86], represented by c 2 s (T ) = d ln T d ln s , the square of the speed of sound in the plasma as a function of the temperature. The potential we use is with γ = 0.606, b 2 = 0.703, b 4 = −0.1, b 6 = 0.0034 (and the asymptotic AdS 5 radius fixed as L = 1). The temperature scale is chosen in order to match the minimum of the speed of sound c 2 s computed holographically with that found on the lattice (we take this minimum to be at 143.8 MeV). The thermodynamics of (2 + 1)-flavor QCD is also used to fix G 5 by matching lattice data for P/T 4 ; this gives G 5 = 0.5013. In Figs. 1 and 2 we compare the lattice results with the results of this holographic model for c 2 s and P/T 4 as a function of the temperature T .
The model is able to describe the temperature region near the minimum of the speed of sound but the agreement does not persist at very high temperatures, which is expected since the model remains strongly interacting in this case while QCD is asymptotically free. Moreover, even though this model does not have the correct (hadronic) degrees of freedom at low temperatures, the temperature dependence of the thermodynamical quantities do follow lattice data even for T ∼ 130 MeV. The model may then be useful precisely in the temperature region T ∼ 130 − 450 MeV where a purely hadronic description is not adequate and the temperature may not be high enough to warrant a simple description using perturbative QCD 24 . 24 We note, however, that non-perturbative weak coupling approaches such as the one pursued in [87,88] do a fairly good job at describing the thermodynamic properties of QCD found on the lattice at high temperatures. The motivation for finding a holographic description of the strongly-coupled QGP relies on the fact that holography is not only able to describe the thermodynamics near the phase transition but it also allows for the direct calculation of non-equilibrium properties, such as transport coefficients, within the same setup.

Holographic calculation of the transport coefficients
Coefficient κ The transport coefficients that we will compute in this section appear in the second-order gradient expansion equations in (2.8) and (2.9). We first compute the coefficient κ using  Table 1.
(3.42). The first step is to fix a u high at which the geometry computed by solving the equations of motion has already reached its AdS 5 asymptotics. For numerical integration, it is necessary to keep the cutoff in (3.42). A value for the cutoff that is too small leads to truncation errors due to the subtraction of two small numbers in (3.42), which gives an artificial numerical divergence. Thus, to rule out errors introduced due to truncation, one should compute the integral in (3.42) with a range of choices for and search for a reasonable range that does not introduce numerical divergences in the integral (which gives a lower bound for ) and also satisfies < u high (which gives an upper bound for ). We found that the optimal region for numerical calculations is 10 −3 < < 10 −1 and we have chosen in this work = 10 −2 for the calculation of κ. We have chosen u high = 0.201, which corresponds to T high = 7.813 T c (where T c = 143.8 MeV). We have checked that at this high temperature both the thermodynamics as well as the transport coefficients have matched their conformal plasma limits.
Proceeding as discussed in the previous paragraph, we determine κ as a function of u H (T ) for the chosen u high . Using the conformal result at strongly coupling (3.43) and the value of G 5 determined in the previous section, we can determine the dimensionless ratio κ/T 2 . In Fig. 3 we show the numerical results for κ/T 2 as a function of T . Moreover, we see that κ/T 2 approaches the conformal limit from below rising monotonically with T . Our results are consistent with the lattice results in [89] where the authors obtained κ/T 2 ∼ 0.36 (15) for T = 2 − 10 T c for a pure glue SU (3) plasma (in this case T c is the critical temperature for the first-order deconfinement phase transition). However, in our model we are able to obtain an estimate for the behavior of κ near the crossover transition of (2+1)-flavor QCD.
For further use, we also present a fit to our calculated points. We use as a fitting function the seven-parameter fit where a to g are dimensionless fit parameters and T c = 143.8 MeV, as remarked before. A five parameter fit using parameters from a to e already yields good results -the two extra parameters f , g are used to provide a closer match to the points. This function corresponds to a modified Fermi-Dirac distribution. In Table 1 we present the parameters for the fit. These fit parameters provide a good description of our numerical data as shown in Fig. 3. With the same remarks as in the preceding section one can evaluate τ π using (3.44) and (3.46). The procedure is similar to the evaluation of κ. The same u high was chosen while in the present case = 2 × 10 −2 -the integrals in (3.46) are more complicated than the integrals in (3.42) and also more sensitive to the choice of the cutoff.
In Fig. 4 we show the numerical results for τ π η/T 2 (this is a convenient choice since, by combining the conformal results (3.37) and (3.40), we see that τ π η ∼ T 2 for a conformal plasma) as a function of T . We see the same general behavior as seen for κ/T 2 or P/T 4 , with a marked increase near the crossover region. The transport quantity τ π η/T 2 approaches its conformal limit (∼ 0.255) from below. We were able to fit these data points using the same parametrization as used for κ/T 2 in Eq. (5.1) -the fit parameters are shown in Table  2. For convenience, we show in Fig. 5 the quantity τ π T . One can see that this quantity approaches its conformal value, (2 − ln 2)/(2π), for T > 300 MeV while it displays a peak near the transition.

Coefficient κ *
We can now use our results for κ to evaluate several other transport coefficients of second order non-conformal hydrodynamics. Since these coefficients involve derivatives of κ with respect to the temperature, up to third order (in the case of the coefficient ξ 4 in (5.14), see below), to avoid discretization errors in the computation of the numerical derivatives we will use the parametrization given by Eq. (5.1) with the parameters displayed in Table 1. The numerical results for the coefficient κ * , defined by (3.2) and computed this way, are  Fig. 6. One can see that κ * /T 2 → 0 as T → ∞ -this can also be checked directly from the expression used for the fit, Eq. (5.1). We note that this is in agreement with the fact that this coefficient vanishes for a conformal plasma. Also, this coefficient has a very sharp dependence with the temperature near the phase transition (following the behavior displayed by c 2 s ) and κ * < 0 for all the temperatures used here.

Coefficient ξ 5
Another second order transport coefficient that we can directly evaluate is ξ 5 , which is determined by the following constraint equation [43][44][45][46] In conformal hydrodynamics one finds that ξ 5 = 0. We evaluated ξ 5 /T 2 as a function of the temperature using the equation above and the result can be found in Fig. 7. One can see that ξ 5 /T 2 has a broad peak in the phase transition around T ∼ 150 − 250 MeV and decreases at high temperatures.

Shear and bulk viscosities
For any isotropic holographic model with an effective gravitational action with at most two derivatives, the shear viscosity satisfies the ratio η/s = 1/(4π) [65] and this is the case of our model 25 .
We can now compute the bulk viscosity ζ using the results of [50,75]. This transport coefficient has attracted some attention recently due to its interplay with shear viscosity 25 We note, however, that the shear viscosity does not obey this ratio if higher order derivative corrections are included in the action [91][92][93] or if the plasma is not isotropic [94][95][96][97][98].
which is defined in terms of the retarded propagator of the spatial trace of the boundary stress-energy tensor T a a (t, x), Holographically, this coefficient is computed considering fluctuations of the xx-component of the metric 26 , h xx . The equation of motion for the perturbation ψ ≡ h x x = e −2A(φ) h xx is given by where the prime denotes a φ-derivative.
As usual, in order to apply the real time prescription for the holographic computation of retarded correlators, we consider the infalling wave condition at the horizon with the normalization condition ψ(φ → 0) = 1 at the boundary. The real time prescription implies that the imaginary part of the retarded correlator is given by where F(ω, φ) is a conserved flux in the radial direction which can then be conveniently evaluated at the horizon 27 where we used Eqs. (4.3), (3.38), and also the relation A (φ H ) = −V (φ H )/3V (φ H ), which can be derived using Einstein's equations for the background metric. Therefore, in order to compute ζ/s from (5.10), one only needs to evaluate (5.6) in the limit of zero frequency C = lim ω→0 ψ(φ → φ H ). We show the numerical results 28 for ζ/s in Fig. 8. The bulk viscosity displays a peak near the phase transition but its magnitude is still smaller than η/s.
For completeness, we also provide a fit to the numerical results of Fig. 8. The form of the function suggests a resonance-like fitting function. With this in mind we used a five-parameter trial function where a to e are fit parameters (and T c = 143.8 MeV). The first term of Eq. (5.11) describes the resonance-like peak of Fig. 8 while the second term describes a smooth background away from the peak. Using the parameters in Table 3, we obtain the fit shown in Fig. 8, which gives a good description of our numerical results.
At this point we have then directly computed 6 transport coefficients (besides matching lattice QCD thermodynamics): η/s, ζ/s, τ π , κ, κ * , and ξ 5 . Among these coefficients, only η/s was found to be a constant with T -all the other coefficients displayed some nontrivial behavior near the crossover phase transition.
Next, we use these results to give our best estimate for the temperature dependence of 6 other coefficients: λ 3 , λ 4 , ξ 3 , ξ 4 , ξ 6 , and τ Π . 27 We expand h(φ) around φ = φH to obtain h(φ → φH ) ≈ h (φH )(φ − φH ). We also make use of (4.4) to obtain h (φH ) = 4πT /e A(φ H )−B(φ H ) . 28 In order to solve Eq. (5.5) for the perturbation ψ we used a part of numerical code created by the authors of [50], which is available at https://www.princeton.edu/physics/research/high-energy-theory/ gubser-group/code-repository/.  Table 3. Estimates for the coefficients ξ 3 , ξ 4 , ξ 6 , λ 3 , and λ 4 Let us now examine three other second order transport coefficients of non-conformal hydrodynamics, ξ 3 , ξ 4 , and ξ 6 , which satisfy the following constraints [43] where the second order coefficients λ 3 and λ 4 are given by the following Kubo's formulas involving Euclidean 3-point functions [43,47] In order to compute ξ 3 , ξ 4 , and ξ 6 using the constraints (5.12) to (5.14) it is necessary to evaluate λ 3 and λ 4 . However, the holographic computation of 3-point functions is far more involved than the computation of 2-point functions. In fact, for the strongly coupled SYM plasma λ 3 has been evaluated explicitly by computing the Euclidean 3-point function G xt,yt,xy E (p, q), yielding λ 3 = 0 [90]. Since κ * = 0 in a conformal theory, from (5.15) we obtain that for a strongly coupled SYM lim pz,qz→0 In order to evaluate λ 4 one should in principle compute the Euclidean 3-point function G tt,tt,xy E (p, q). However, there is a shortcut which makes use of the constraint (5.12). In a 4-dimensional CFT, we know that c 2 s = 1/3 and we also know that κ * = 0 and that ξ 3,4,5,6 = 0, since these are coefficients of non-conformal hydrodynamics. Then, from Eq. (5.12) we deduce that in a CFT λ 4 = 0, which agrees with [9]. Then, from Eq. (5.16), we conclude that in a strongly coupled CFT lim px,qy→0 The evaluation of the full 3-point functions required to compute the coefficients λ 3 and λ 4 in the effective model of Einstein+Scalar gravity where the metric is only known numerically is far beyond the scope of this work. In order to fully determine them it is necessary to compute the full bulk-to-boundary and bulk-to-bulk propagators -and it is very difficult to compute these functions in terms of a numerical metric such as the one used in this work.
Thus, in this paper we will resort to a sort of "hybrid CFT/non-CFT" approximation. For the evaluation of the second order coefficients λ 3 and λ 4 from Eqs. (5.15) and (5.16) we shall mix the CFT 3-point functions (5.17) and (5.18) with the full non-conformal results for κ and κ * . In this case, the resulting approximations for λ 3 and λ 4 give With these approximations, and the full non-conformal results for κ, κ * , and c 2 s , we can approximately evaluate ξ 3 , ξ 4 , and ξ 6 from Eqs. (5.12) to (5.14). This sort of approximation provides our best estimate for these coefficients given the current lack of knowledge about 3-point functions in non-conformal holographic plasmas that display a crossover phase transition. In Figs. 9 to 11 we show the results for ξ 6 /T 2 , ξ 3 /T 2 , and ξ 4 /T 2 as functions of T -  all of these coefficients vary rapidly near the phase transition. We note that the approximations done here for these coefficients constitute the first deviations from the ultraviolet conformal regime and, therefore, they are much more reliable at high temperatures.
A lower bound estimate for τ Π Ref. [103] derived, using the asymptotic causality condition, a relation among the transport coefficients τ π , τ Π , η, and ζ ζ The computation of the transport coefficient τ Π directly from its retarded Green's function, as was done for the shear coefficient τ π , is beyond the scope of this paper. However, we note that one can use the relation (5.20) to obtain a lower bound for the coefficient τ Π that can still be useful for hydrodynamic modeling of the QGP 29 . The result for this lower bound in our holographic model is shown in Fig. 12 (this was computed using directly the fitting functions and Eq. (5.20)). This is the smallest value of τ Π T in our model that is still consistent with causality and linear stability [103]. One can see that this coefficient displays a peak near the phase transition, as was the case for τ π T , but it becomes very small at high temperatures, as expected due to conformal invariance. The results in Fig.  12 also admit a fit using the following parametrization, with the corresponding fit parameters a to d being given in Tab. 4. Therefore, in this section we presented results (computed within different levels of approximations) for 12 transport coefficients that appear at second order in the gradient 29 A similar idea was used in [104] to estimate the coefficient τπ in a hadronic gas with Hagedorn resonances given the result found in this case for η/s [105].  Table 4.
expansion of a non-conformal plasma that has thermodynamic properties similar to those found for the QGP on the lattice: η/s, ζ/s, τ π , κ, κ * , ξ 5 as well as ξ 3 , ξ 4 , ξ 6 , λ 3 , λ 4 , and τ Π . However, the equations of motion obtained from (2.8) and (2.9) are not in a form suitable for numerical implementation. In the next section we use the gradient expansion to find another 2nd order theory, similar to Israel-Stewart theory [106], that can be readily used in phenomenological studies of the QGP hydrodynamical evolution in heavy ion collisions. 6 Israel-Stewart-like 2nd order hydrodynamics for a non-conformal relativistic fluid It is known that relativistic NS theory leads to acausal propagation of sound and shear linear disturbances around an equilibrium state at rest and that in the case of a moving background fluid these disturbances are unstable [107,108]. It can be shown that the 2nd order theory in Eqs. (2.8) and (2.9) is linearly unstable above a certain critical wavenumber even for a fluid at rest, as demonstrated in Appendix C (see also [109]). Note also that the inclusion of these spatial gradients cannot modify the theory at very large momenta k where the asymptotic causality conditions are defined [103]. However, it is not the purpose of hydrodynamics to accurately describe small wavelength phenomena -this is beyond the scope of this effective theory. Nevertheless, it is desirable that for a system that may be coupled to gravity causality (and stability!) are preserved.
In this section we use the 2nd order theory in (2.8) and (2.9) to construct a relaxationtype theory (in curved spacetime) that is similar to that considered by Israel and Stewart [106] and also to those that appear naturally within kinetic theory using the moments method [25]. The main idea, already pursued in [41] in the case of a conformal fluid, is to write a simple relaxation-type theory that reduces to the gradient expansion theory in (2.8) and (2.9) in its asymptotic hydrodynamical limit. In this section we shall use the same procedure generalized to the case of a non-conformal fluid.
First, the terms involving comoving derivatives of σ µν and θ on the right-hand side of (2.8) and (2.9) are transferred to the left-hand side of those equations using the lowest order substitutions σ µν → −π µν /η and θ → −Π/ζ. Moreover, we use the leading order expression Ds = −sθ+. . . in (2.6) to simplify some of the terms (note that we are neglecting terms of third order in gradients that would appear in this general procedure). The same leading order substitution is done in the remaining 2nd order terms on the right-hand side of the equations with the exception of the term ∼ τ π η θσ µν in (2.8) where only the substitution σ µν → −π µν /η is done. This choice is made to make it explicit that the combination Dπ µν + 4θπ µν /3 is the correct combination that survives in the conformal limit [41] (being, thus, homogeneous under Weyl transformations). One can then show that this leads to and τ Π (DΠ + Πθ) + Π = −ζθ + ξ 1 η 2 π µν π µν + These are nonlinear, coupled partial differential equations of relaxation-type for the new dynamical variables π µν and Π in curved spacetime that require (independent) initial conditions in order to solve Eqs. (6.1) and (6.2) together with the conservation equations (2.5) 30 . These equations are similar to those found in Israel-Stewart theory [106] and they possess most of the terms found in kinetic theory in flat spacetime [25] 31 . Furthermore, note that the asymptotic, leading order solution of these equations necessarily reduce to the gradient expansion in Eqs. (2.8) and (2.9) up to O(K 2 n ).
However, for phenomenological applications in heavy ion collisions, the terms containing spacetime curvatures are negligible and can, thus, be dropped. Moreover, note that 30 This theory is qualitatively different than that in the gradient expansion -the dissipative parts of the energy-momentum tensor have their own differential equations and, thus, its initial conditions are not determined by the initial conditions for the hydrodynamic variables ε and u µ . 31 Ref. [25] used the Boltzmann equation and a completely distinct power-counting scheme to deal with the gradients in comparison to the one used to derived our equations for π µν and Π. For instance, according to the power-scheme of [25], the vast majority of their terms of order O(K 2 n ) do not appear in our equations. Thus, perfect agreement among these approaches should not really be expected.
while in our holographic model η/s is a constant, in general one should keep the term involving D ln(η/s) to make sure that this theory reduces to the correct gradient expansion theory if the asymptotic limit π µν → −ησ µν is taken. For the same reason, one should keep D ln(ζ/s) and, in fact, in our model since ζ/s is a not a constant its comoving derivative is not zero (note also that the conservation equations imply that, to lowest order, D ln(η/s) , D ln(ζ/s) ∼ −θ).
This leads us to the following (reduced) set of equations that can be used in hydrodynamic simulations of the QGP and One can show that the hydrodynamic theory described above is linearly stable and causal according to the criteria of [103] and, thus, it should be suitable for implementation in modern numerical viscous hydrodynamic codes such as [99,[110][111][112][113][114]. We stress that the terms involving D ln(η/s) and D ln(ζ/s) in the equations above are needed to recover the correct 2nd order gradient expansion and should not in principle be neglected in numerical simulations.
Among the 13 transport coefficients left in the equations above, results for 8 of them have already been presented in this paper while we have not yet discussed the coefficients λ 1 , λ 2 , ξ 1 , ξ 2 , and τ * π . The coefficients λ 1 and λ 2 have been studied at weak coupling in [115] and at strong coupling in [29,41,90]. The SYM values of these coefficients computed at strong coupling via holography are 32 We are not aware of any calculation of these coefficients in a non-conformal strongly coupled plasma with a crossover transition. However, within the phenomenological "spirit" of this section and since we currently lack a better way to compute them, one may take the SYM expressions above for the non-conformal case at hand. This would imply that λ 1 /T 2 ∼ s/T 3 and λ 2 /T 2 ∼ −s/T 3 . Thus, in this case these coefficients would display the same sharp rise observed by the entropy density near the phase transition.
Knowledge about the coefficients ξ 1 and ξ 2 is much more scarce. These coefficients only appear in non-conformal fluids and very little is known about them at strong coupling. An exception is the strongly coupled non-conformal plasma studied in Ref. [120] constructed via dimensional reduction of a higher dimensional pure gravity action. In this case, these coefficients can be extracted using the fluid/gravity correspondence and they read [9] 33 Also, in this theory one finds [9] τ * π = −3τ π In the absence of a better estimate for these three coefficients above in the non-conformal plasma proposed in this paper, it may be useful in hydrodynamic simulations of the QGP to use the expressions in (6.6) and (6.7) hoping that they get at least part of the nonconformal dynamics near the phase transition. Notice, however, that these expressions contain ( 1 3 −c 2 s ) and that this is an ubiquitous factor in non-conformal transport coefficients -for instance, the bulk viscosity of our model is proportional to this factor [75,121]. Thus, it is reasonable to assume that these expressions may describe the temperature behavior of these coefficients in our model as well. Moreover, we would like to remark that our calculation for the transport coefficients are qualitatively consistent with the results found in [122] for the large T expansion of the 2nd order transport coefficients for the nonconformal plasma dual to the Chamblin-Reall background [123]. Therefore, after this long discussion, the hydrodynamic theory described by Eqs. (6.3) and (6.4) together with the corresponding conservation equations (2.5) may be a good starting point for phenomenological applications of relativistic non-conformal hydrodynamics for the strongly-coupled QGP formed in heavy ion collisions 34 . To facilitate the use of our results in hydrodynamic simulations, we have provided a guide for all the relevant formulas for the 13 transport coefficients of this second order theory in Appendix B.

Conclusions
In this paper we used the gauge/gravity duality to determine the transport coefficients of a non-conformal, strongly interacting non-Abelian plasma that displays a crossover transition similar to that found for the QGP determined via lattice calculations. The 5-dimensional 33 The transport coefficients in the theory [120] are known analytically and, even though their numerical values are different than the ones found in this paper (their theory is different than ours), qualitatively they possess the same features found here -ξ5, ξ6 have the same signal as ours and would also display a peak where c 2 s has a minimum. However, their τπT is a constant while ours has a peak near the phase transition. 34 Note that the thermodynamics of the model is very similar to lattice data and, thus, in hydrodynamic simulations one may as well just use directly the lattice data for the thermodynamical quantities such as c 2 s or s. gravity dual model involves the metric coupled with a dynamical scalar field and its simplicity and capability of describing several nontrivial features of the QGP have motivated us to pursue the calculations of the several transport coefficients shown in this paper.
We first obtained holographic formulas for the transport coefficients κ and τ π present in the second-order gradient expansion of relativistic hydrodynamics in curved spacetime. Our method to compute these coefficients could also be applied in the case where the gravity dual possesses other fields besides the scalar field, such as the case of an Ein-stein+Scalar+Maxwell model with at most two derivatives in the action. Also, besides the well-known result for η/s, we were also able to directly compute five other coefficients that appear at second-order in the derivative expansion: τ π , ζ/s, κ, κ * , and ξ 5 . Apart from η/s, all of these coefficients displayed nontrivial behavior near the crossover transition. In particular, τ π T , ζ/s, and ξ 5 /T 2 display a peak near the transition while κ * /T 2 is similar to c 2 s (though it is negative) in that it displays a minimum at the crossover transition. On the other hand, κ/T 2 rises monotonically with T until it reaches its conformal limit (i.e., its value in SYM). Our values for τ π T only deviates from the SYM result (2 − ln 2)/(2π) at low temperatures. The coefficients κ, κ * , and ξ 5 only contribute directly to the equations of motion in a curved spacetime.
Our ζ/s is in general smaller than that found in other works [124,125] as it is clear in 35 Fig. 13, though it is similar in magnitude to the pQCD calculation in Ref. [36]. Thus, at least according to our calculations, cavitation [126][127][128] induced by a large ζ/s in the phase transition is not likely to occur in the QGP (a similar conclusion was reached in [129] using a kinetic theory-derived bulk relaxation time in a Bjorken expanding fluid). However, it could be that the other coefficients that appear in the bulk equation, together with the shear-bulk coupling terms such as the π µν Π term in (6.3), may in the end take the evolving plasma towards cavitation. This is an interesting possibility that can be checked in numerical simulations.
We used these calculations to provide estimates for the other coefficients ξ 3 , ξ 4 , ξ 6 , λ 3 , λ 4 , and τ Π . We found that ξ 3 /T 2 and ξ 4 /T 2 are negative and have a minimum near the transition while ξ 6 /T 2 is positive and displays a peak. Note that ξ 6 is only relevant in curved spacetime while ξ 3 , λ 3 , ξ 4 , λ 4 do affect the motion of the fluid in flat spacetime. In fact, ξ 3 and λ 3 are related to the vorticity tensor Ω µν whose role in hydrodynamic simulations has not yet been investigated in detail 36 . Moreover, ξ 4 and λ 4 also have not been investigated in hydrodynamic calculations and, thus, we hope the results of this paper may serve as motivation for a detailed investigation of their effects. In this paper we have used the asymptotic causality condition [103] to obtain the lowest possible value of the coefficient τ Π associated with bulk viscosity relaxation. In this case, this lower bound for τ Π T displays a peak near the transition (though its value at the peak is relatively small, 35 We thank J. Noronha-Hostler for making this plot available to us. 36 Note that for (0+1) purely Bjorken hydrodynamics this term disappears even in a non-conformal  [36] towards low temperatures, the blue squares shows the calculation using the PHSD model [124], while the black points shows the hadronic calculation from [125].
in agreement with the small value of ζ/s found here) 37 .
We have used the 2nd-order gradient expansion equations to construct an Israel-Stewart-like theory in flat spacetime, shown in Eqs. (6.3) and (6.4), which gives equations of motion that preserve causality and are linearly stable around thermal equilibrium 38 . These equations are similar to those used in current viscous hydrodynamic codes but they include additional 2nd order terms that are usually not taken into account.
This simplified theory still contains 13 transport coefficients and we have presented in 37 We remark that the coefficients τπ and τΠ, as defined via the gradient expansion, do not necessarily correspond to relaxation time coefficients. Clearly, relaxation time coefficients require at least relaxation equations for πµν and Π, such as those in Israel-Stewart theory. In fact, according to their definitions in the gradient expansion, these coefficients do not even need to be positive. For instance, [130] has found an example of a gravity dual in which τπ < 0 as defined via the gradient expansion. This, however, was shown to not generate instabilities as it would have been the case if that coefficient were indeed a measure of shear relaxation. As discussed in [39], only the coefficients extracted from the poles of retarded correlators do have the meaning of shear or bulk relaxation time coefficients, which is not generally the case for the coefficients defined via derivatives of retarded Green's functions (such as in the gradient expansion). 38 It is important to remark that the procedure used here to find relaxation equations using the gradient expansion has some ambiguities. Clearly, the set of relaxation equations obtained this way is not uniqueit only corresponds to one of the possible sets of equations that have the 2nd order gradient expansion as their asymptotic solution. At this level, this can be viewed as a type of UV completion procedure of the gradient expansion equations that is consistent with the asymptotic causality condition. Alternatively, this general procedure can be illustrated via a simple classical mechanics example. The differential equations x + γẋ + x = f (t) and γẋ + x = f (t) have the same asymptotic solution xasymp(t) ∼ f (t) for large times tγ 1 though their transient (short time tγ ∼ 1) behavior can be very different.
this paper either direct calculations or leading estimates for all of them. In the case of the transport coefficients τ * π , λ 1 , λ 2 , ξ 1 , and ξ 2 much more work needs to be done to obtain their exact temperature dependence in our holographic non-conformal model with a crossover phase transition. At the moment, our best estimate for their temperature behavior consisted in using the known expressions for λ 1 and λ 2 from SYM and the expressions for τ * π , ξ 1 , and ξ 2 from the (related) non-conformal model of [120]. This (admittedly uncontrolled) approximation must be taken with care, as emphasized in the main text. However, given our current ignorance regarding these coefficients, we believe that it still may be of phenomenological interest to heavy ion collisions to use these expressions and investigate their consequences. In particular, the direct shear-bulk coupling term τ * π may be very relevant in hydrodynamic simulations, as emphasized in [131] and [132,133].
Also, regarding the complete evaluation of the 2nd order transport coefficients that appear in the gradient expansion of non-conformal strongly-coupled hydrodynamics, the fluid/gravity correspondence [29] provides a way to compute all the coefficients. However, when the background metric is only known numerically, as it is in our case with a crossover phase transition, the actual implementation of the fluid/gravity approach becomes much more challenging. In this context, it may be useful to consider a simpler analytical model that still possesses a phase transition. For instance, one may consider the finite temperature holographic model of Ref. [134] where the background metric and scalar field are known analytically and the system displays a 1st-order deconfinement phase transition. In this case, it should be possible to carry out the fluid/gravity procedure and find expressions for all the 2nd order transport coefficients. This would allow for a complete study of entropy production in 2nd order non-conformal hydrodynamics [9] for a theory that displays a phase transition.
We remark that very similar equations of motion for a non-conformal plasma have been derived from the Boltzmann equation in [131] and, in that paper, the authors also gave explicit formulas for several 2nd order transport coefficients. It would be interesting to compare the results of hydrodynamic simulations computed using the strongly coupled transport coefficients of this paper with those obtained using the kinetic theory-derived coefficients of [131] 39 . We have gathered in Appendix B the fitting functions that describe the temperature dependence of all the transport coefficients in this 2nd order Israel-Stewart theory to facilitate their use in current hydrodynamic codes.
This brings us to an important point concerning the equations of motion of stronglycoupled relativistic hydrodynamics in the light of the gauge/gravity duality. As discussed in [39,135], the fact that the non-hydrodynamic modes of the xy − xy energy-momentum tensor retarded correlator possess comparable real and imaginary parts at zero momentum [136] implies that, strictly speaking, the effective theory that should be able to describe the hydrodynamical sound and shear modes as well as the lowest set of non-hydrodynamic modes is not of relaxation-type such as in Israel-Stewart theory (as obtained from the Boltzmann equation). Ref. [137] has recently proposed a way to describe the approach towards hydrodynamics in strongly-coupled SYM that involves equations of motion that are qualitatively different than those in (6.3) and (6.4) since they involve a second order, homogeneous differential equation for the part of the shear stress tensor associated with the two lowest quasinormal modes π µν QN M (also, in their effective approach nonlinear terms in π µν QN M are not taken into account). It would be interesting to investigate whether their effective theory is also applicable in the case of the non-conformal plasma studied in this paper.
A Solution for the metric perturbation up to O(ω 2 , q 2 ) In this Appendix we formally solve Eq. (3.18) in detail in terms of the coefficients of the undisturbed background metric, g M N (u), and the boundary value of the metric perturbation, ϕ(ω, q), up to O(ω 2 , q 2 ). In order to accomplish that one must specify the boundary conditions for the metric perturbation, φ(u, ω, q), at the boundary and at the horizon.
We first consider the boundary. The asymptotic form of Eq. (3.18) near the boundary u = for any asymptotically AdS background reads 40 where k 2 = −ω 2 + q 2 and I n (ξ) and K n (ξ) are the modified Bessel functions of the first and second kinds, respectively. Now we take which is a constant in the radial direction. Therefore, in the case of a massless scalar field, one can safely impose the Dirichlet boundary condition as follows where ϕ(k) denotes the boundary value of the metric perturbation prescribed by the Dirichlet boundary condition. Now that we have specified the boundary condition (A.3), let us discuss which kind of condition we must impose on the metric perturbation at the horizon. As discussed in [48][49][50][51], in order to obtain the retarded propagator we must single out the solution of the equation of motion (3.18) which is regular at the horizon and corresponds to a wave being absorbed by the horizon. We begin by assuming that the g tt component of the background metric has a simple zero at the horizon u = u H such that one can write with G(u H ) and F (u H ) finite. Notice that the Hawking temperature of the black brane gives and the infalling wave mode at the horizon is obtained by setting C 2 (ω) = 0 in Eq. (A.6). Therefore, one is motivated to separate the infalling behavior of the solution and take the following Ansatz for the metric perturbation with f (0, k) = 1 and f (u H , k) being regular. We considered only even powers of q in the series expansion (A.8) due to spatial isotropy. Substituting Eq. (A.7) into Eq. (3.18) we obtain the differential equation for f (u, k), Now we solve Eq. (A.9) order by order in the momentum expansion (A.8).
Zeroth orderf 0 At lowest order, one obtains the differential equation for f 0 (u), O(ω 0 , q 0 ) : ∂ u −g (0) g uu f 0 = 0. (A.10) 41 We only keep the dominant terms in each order in radial derivatives of the metric perturbation. where we just performed two indefinite integrations: the lower limit in the integral present in Eq. (A.11) has been conveniently chosen to lie at the boundary; however, we can choose any other value of the radial coordinate to be the lower limit of this integral and each possible different choice would just redefine the integration constants c 1 and c 2 . We immediately fix these constants by using the boundary condition f (0, k) = 1 and the regularity condition for f (u H , k). The former fixes c 1 = 1 and the latter fixes c 2 = 0; therefore, Integrating once, we obtain, = c 2 , (A.14) and, then, integrating again, we find Applying the boundary condition f (0, k) = 1 we obtain c 1 = 0. The horizon regularity condition implies 16) and, therefore, . In Eq. (A.19), the factors (4πT (u H − λ)) −1 and c 2 / −g (0) (λ)g uu (λ) diverge at the horizon while the factor −g (0) (ξ)g tt (ξ) diverges at the boundary; therefore, the boundary condition and horizon regularity fix c 1 = 0 and 42 Thus,

B Summary of the transport coefficients for hydrodynamic simulations
The equations of motion for the Israel-Stewart-like theory in flat spacetime are (6.3) and (6.4) together with the conservation equations (2.5). The 13 transport coefficients in this theory, namely η/s, ζ/s, τ π , τ Π , τ * π , λ 1 , λ 2 , λ 3 , λ 4 , ξ 1 , ξ 2 , ξ 3 , and ξ 4 were discussed already in the main text but here we gather the formulas that describe them in a single place with the intention to facilitate their use in hydrodynamic simulations.
Besides η/s = 1/(4π), the shear relaxation coefficient τ π η/T 2 is described by the function in (5.1) with parameters in Table 2. ζ/s is described by the fitting function in Eq. (5.11) with the parameters in Table 3. The lower bound for the bulk relaxation coefficient τ Π T is described by (5.21) with parameters in Table 4. The coefficients λ 3 = −λ 4 = 2κ * , where κ * in (3.2) is defined in terms of the coefficient κ. One can find a fit for κ/T 2 in (5.1) with parameters in Table 1. The coefficients λ 1 and λ 2 can be found in Eq. (6.5). Moreover, the coefficients ξ 1 , ξ 2 , and τ * π are found in Eqs. (6.6) and (6.7), respectively. The final coefficients ξ 3 and ξ 4 are given by Eqs. (5.13) and (5.14) and can be computed using the results for the previous coefficients discussed above. Furthermore, in all of these fitting functions and tables, the parameter T c is equal to 143.8 MeV.
Given that the thermodynamic properties of the model are very similar to those found on the lattice [86] when T ∼ 130 − 450 MeV, if one wants to use the transport coefficients computed in this paper in hydrodynamic simulations of the QGP formed in heavy ion collisions one may just use directly an interpolation for the lattice data when computing c 2 s and the entropy density s needed in the evaluation of the transport coefficients.

C Linear instability of the gradient expansion at 2nd order
In this Appendix we investigate the linear stability properties of a fluid described by the 2nd order gradient expansion theory in (2.8) and (2.9) (in flat spacetime) around static equilibrium. In a linear analysis, the relevant linear terms involving the dissipative part of the energy-momentum tensor are π µν = −ησ µν + ητ π Dσ <µν> , Π = −ζθ + ζτ Π Dθ . (C.1) We follow [103,107,108,138] and consider linear perturbations around a static background. In order to investigate the stability of the sound channel, it is sufficient Note that in the limit when τ π,0 , τ Π,0 → 0 one finds that k sound c diverges, which is in agreement with the fact that NS theory is stable against small perturbations in a fluid at rest [107]. Clearly, for a moving fluid the stability properties become more involved but it is possible to show that the same type of problems that appears in NS theory also appear in this case [109]. However, we remark that hydrodynamics is only expected to be valid in the low frequency, large wavelength limit. Nevertheless, the instability found here at finite wavenumber motivates the search for a UV completion of this theory (for instance, the one discussed in Section 6) that is linearly stable and can, therefore, be safely used in numerical simulations.
A larger critical wavenumber, k shear c , appears in the shear channel. Linear stability of shear modes can be studied by choosing a flow disturbance of the kind u µ shear = (1, 0, 0, 0)+(0, 0, δu y (t, x), 0) while the other relations in (C.3) remain valid (see [138]). This leads to the following dispersion relation For a recent study involving the first nonlinear corrections to the stability analysis and the propagation of waves in relativistic hydrodynamics see, for instance, [139].