Effective holographic models for QCD: Thermodynamics and viscosity coefficients

A finite temperature extension of the effective holographic models for QCD (EHQCD), proposed in Ref.[1], is investigated in the present work. EHQCD models are characterized by two parameters, the conformal dimension of the relevant operator that deforms the CFT and the associated coupling. We find that black hole solutions appear at temperatures higher than some temperature $T_{min}$ and can be categorized in two classes: large and small black holes. A large black hole is thermally stable and it is therefore interpreted as the gravity dual of a non-conformal plasma. A small black hole, on the other hand, is thermally unstable. We show that thermodynamic quantities such as the entropy density $s$, specific heat $C_V$, and speed of sound $c_s$ are sensitive to the model parameters. We investigate perturbations of the black hole solutions and calculate the viscosity coefficients of the corresponding dual non-conformal plasma. For the shear viscosity, we confirm that the ratio $\eta/s$ is given by the universal result $1/4\pi$. For the bulk viscosity, the ratio $\zeta/s$ varies with the temperature, displaying a rapid growth close to $T_{min}$, and it is sensitive to the model parameters. We compare our results for the thermodynamic quantities with the lattice $SU(N_C)$ results and find that they are compatible as long as the coupling is fixed appropriately as a function of the conformal dimension. We also compare our results for the viscosity coefficients against the JETSCAPE results that are obtained from the analysis of experimental data on heavy ion collisions.


Introduction
The emergence of the Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [2][3][4] about two decades ago opened a new window to investigate the very specific properties of the quark-gluon plasma (QGP) produced for the first time at the Relativistic Heavy Ion Collider (RHIC) [5][6][7][8]. Due to its collective (macroscopic) behavior, this plasma may be described by relativistic hydrodynamics. Moreover, the QGP phase lies in the strong coupling regime of QCD and the perturbative methods of quantum field theory are not useful. The AdS/CFT correspondence establishes a duality between strongly coupled conformal field theories in 4d flat space and weakly coupled gravitational or string theories in AdS 5 × M with M a compact manifold. The most notable example of this correspondence is the duality between N = 4 Super-Yang-Mills theory in 4d flat space and Type IIB Supergravity in AdS 5 × S 5 (see, for instance, Refs. [9][10][11] and the references therein). Motivated by the initial success of the AdS/CFT correspondence, phenomenological models were proposed to investigate Quantum Chromodynamics (QCD)-like theories. Since QCD is a non-conformal theory (except in the extreme ultraviolet (UV) region), a deformed CFT in which the conformal invariance has been broken by a scalar operator may be considered as a reasonable approximation for QCD. This kind of approach can be seen as a phenomenological bottom-up model for holographic QCD. To mention a few of previous works, Gubser and Nellore [12] built a holographic QCD model, where the CFT is deformed by a relevant operator, constrained by the equation of state of QCD. Viscosity coefficients were later obtained in Refs. [13,14]. In turn, Gursoy, Kiritsis, and Nitti [15] built another holographic QCD model, where the CFT is deformed by a marginal operator. Moreover, the authors of [15] found the appropriate IR behavior of the gravitational background that leads to confinement and to a linear spectrum. As a consequence, the model of [15], known as improved holographic QCD (IHQCD), leads to a glueball mass spectrum in agreement with the results of lattice QCD. The unconfined phase in the IHQCD model was investigated in Ref. [16] and the viscosity coefficients were obtained in [17].
The holographic QCD models described above have specific properties that make them unique. Considering the ultraviolet (UV) properties of the Gubser-Nellore model [12] and the confined infrared (IR) characteristics of the IHQCD model built by Gursoy, Kiritsis and Nitti [15], a simple holographic (hybrid) model was recently proposed and explored in Ref. [1]. The model parameters in [1] are the conformal dimension ∆ = 4 − of the relevant operator that deforms the CFT and the associated dimensionless couplingφ 0 . The model leads to a glueball mass spectrum in good agreement with lattice QCD and it was also realized that the mass spectrum is not sensitive to the variation of the conformal anomalous dimension as long as the couplingφ 0 is fixed as a function of (for details, see Refs. [1,18]).
In relation to the QGP, its transport coefficients (for example, the viscosity coefficients) must be finite, even if the plasma approaches a perfect fluid. Thus, the knowledge of their viscosity coefficients shall provide a better understanding of the properties and the evolution of the QGP. One way to calculate such coefficients is to consider the low-energy limit of the real-time spectral functions and then, using the Kubo formula, obtain the desired viscosity coefficients. This procedure is followed by some of the research groups working on lattice SU (N c ) theories (see, for instance, Ref. [19][20][21]). However, as it is well known, transport properties are still a very challenging problem in lattice SU (N c ) theories and there are large uncertainties [22] (for a discussion on the bulk viscosity, see for instance Ref. [23]). Using the holographic QCD approach, based on the AdS/CFT correspondence, one may be able to extract the viscosity coefficients of the dual plasma. One of the methods used is turning on a source for the desired response; for example, considering metric perturbations of the dual 5d asymptotically AdS black hole. Imposing incoming boundary conditions at the black hole horizon and Dirichlet boundary conditions at the AdS boundary one is able to calculate the corresponding retarded Green's functions. Then, using the Kubo formula it is possible to extract the shear and bulk viscosities of the dual plasma.
In this paper, we investigate the thermodynamic properties and viscous coefficients of a non-conformal plasma in an EHQCD approach, where the CFT is deformed by a relevant operator dual to a 5d scalar field (the dilaton). In this approach, the gravity dual of the non-conformal plasma is given by a 5d asymptotically AdS black hole arising from Einsteindilaton gravity with a specific dilaton profile that matches the conformal dimension of the 4d operator in the UV and satisfies the confinement criterion in the IR. An important consequence of the confinement scale, present in the EHQCD model, will be the emergence of a minimum temperature T min for the existence of black hole solutions. Moreover, at fixed temperature T > T min there will be two types of black hole solutions, namely the large and small black holes. The large (small) black hole will be characterized by the entropy that increases (decreases) as a function of the temperature. The large black hole will be thermally stable and therefore interpreted as the gravity dual of the non-conformal plasma.
We investigate the thermodynamics and the transport properties of the non-conformal plasma and find that they are sensitive to the two relevant parameters of the EHQCD model, namely the conformal dimension ∆ = 4 − of the scalar operator that deforms the CFT and the corresponding dimensionless couplingφ 0 . We obtain the entropy of the non-conformal plasma using the Bekenstein-Hawking formula for the black hole entropy whereas the plasma free energy will be obtained by the use of the first law of Thermodynamics. Expanding the action up to second-order we are able to extract the viscosity coefficients. For the shear viscosity, we rewrite the perturbations in terms of gauge-invariant variables describing the transverse and traceless sector. For the bulk viscosity, we change our strategy and rewrite the metric and equations of motion using the dilaton field as an independent variable. In such a procedure the dilaton field plays the role of the holographic coordinate, allowing us to simplify the calculations. For the thermodynamics, we compare our results with those found in lattice SU (N c ) theories [24]. Interestingly, we find that the thermodynamic quantities can be fit for any value of the conformal anomalous dimension as long as the dimensionless couplingφ 0 is fixed appropriately as a function of . This is consistent with our previous results at zero temperature [1,18], where fitting the glueball mass spectrum also led to fixingφ 0 as a function of . For the viscous coefficients, we compare our results with the phenomenological constraints found by the JETSCAPE collaboration, from a model-to-data analysis of heavy ion collision experimental data [25].
The present paper is organized as follows. In Sec. 2 we describe the extension of the EHQCD models to finite temperature and the physical interpretation of the model parameters. Sec. 3 deals with the thermodynamic properties of the EHQCD models, where we describe the two black hole solutions appearing at finite temperature and discuss their stability. In section 4 we investigate the viscosity coefficients considering the metric perturbations dual to the stress-energy tensor components associated with shear and expansion. Sec. 5 is devoted to comparing our results for the thermodynamics with the results found in lattice SU (N C ) theories and comparing our results for the viscous coefficients with the phenomenological constraints found by the JETSCAPE collaboration. We present our conclusions in Sec. 6, and the asymptotic analysis of the EHQCD models in Appendix A.
2 Effective holographic QCD models at finite temperature In this section, we summarize the general properties of the effective holographic QCD models (EHQCD) and present the ansatz for the black hole solutions that allow for a description of the non-conformal plasma at finite temperature. In sections 3 and 4 we will describe the thermodynamic and transport properties of the non-conformal plasma. For more details on EHQCD models at zero temperature, see the discussion presented in Refs. [1,18].
EHQCD models arise as asymptotically AdS solutions in 5d Einstein-dilaton gravitational theory. These 5d backgrounds are the gravity duals of 4d deformations of a CFT by a relevant scalar operator [12]. The Lagrangian of the deformed CFT is given by where φ 0 is the source, and O is the scalar operator dual to the dilaton field. The physical motivation behind this approach is the possibility of interpreting the Lagrangian (2.1) in terms of the Yang-Mills Lagrangian in the limit of a large number of colors, N c . The profile of the 5d dilaton field is constrained in the UV (close to the AdS boundary) and in the IR (far from the AdS boundary). The UV asymptotics is fixed by the relation m 2 = ∆(∆ − 4) between the 5d dilaton mass and the conformal dimension of the 4d scalar operator. For the IR asymptotics, we follow [15] where it was found a dilaton field that, far from the boundary, grows quadratically in the radial direction and leads to an asymptotically linear spectrum for the scalar and the tensor glueballs. Moreover, the authors of Ref. [15] showed that such a quadratic behavior also guarantees the key requirement for color confinement, namely that the warp factor in the string frame presents a global minimum [26]. The starting point in building holographic models for QCD is the five-dimensional action for Einstein-dilaton theory, where Φ is the scalar (dilaton) field, V (Φ) is the associated potential, and σ = M 3 p N 2 c is the effective 5D gravitational coupling. M p is the effective Planck mass and N c is a finite constant parameter associated with the number of colors in the four-dimensional dual field theory. The Einstein-dilaton equations of motion obtained from this action are where G mn stands for the Einstein tensor and T mn is the energy-momentum tensor defined by with S M = σ d 5 xL Φ being the matter action due to the scalar field. Let us notice that Eqs. (2.3) and (2.4) are expressed in general forms and then they apply to any spacetime geometry. However, in order to investigate holographic QCD models at finite temperature we consider the following black hole metric where f (z) is the horizon function (or blackening function), ζ 1 (z) is a global (inverse) scale factor for the 5D metric related to the warp factor through the relation ζ 1 = e −A 1 , with A 1 being the warp factor 1 . The black hole solution is characterized by the presence of an event horizon, z h , where the horizon function vanishes. With this, the holographic coordinate z is restricted to the interval 0 ≤ z ≤ z h . A confining thermal solution can also be investigated by setting f (z) to 1 for 0 < z < ∞. This is the trivial extension of the zero temperature solution dual to a confining vacuum. To describe the physics of a dual non-conformal plasma, in this work we focus on a black hole solution corresponding to a nontrivial f (z) that interpolates between f (0) = 1 at the AdS boundary and f (z h ) = 0 at the event horizon. After plugging (2.6) into (2.3), the equations of motion can be written in the form where the prime indicates the total derivative with respect to z. The coupled equations (2.7) may be solved by following different approaches. The traditional method is to build a dilaton potential that interpolates between the UV and the IR regime and then solve the equations for the fields Φ, ζ 1 , and f . This is the strategy considered in the pioneer works of [12] and [15]. As described in our previous work [1], an alternative strategy for solving the equations (2.7) is to presuppose that we know either the dilaton field Φ (models A) or the scale factor ζ 1 (models B). One can build simple interpolations in any of those scenarios that satisfy the UV and IR constraints. The potential is then reconstructed by the use of the 2nd equation in (2.7). Any of those two scenarios are usually known as the potential reconstruction method and have been successfully used in previous works, see for instance Refs. [27][28][29][30].
In the models of type A, where the dilaton profile is known, the (inverse) scale factor ζ 1 is calculated by solving numerically the first equation in (2.7). Then, solving the last equation in (2.7), we get an expression for the function f (z) in terms of the warp factor. The two integration constants arising in the blackening function are fixed by considering two boundary conditions. Namely, a condition imposed at the AdS boundary, where f (0) = 1, while the other one is given by the regularity condition at the horizon, where f (z h ) = 0. Thus, the blackening function results in In the present work we use subscripts to distinguish the scale factor ζ1 from the bulk viscosity ζ.
where the constant C h , which guarantees the regularity condition f (z h ) = 0, is given by . (2.9) The (inverse) scale factor ζ 1 is a non-negative monotonic function of z and hence the constant C h is non-negative. This in turn gives rise to a monotonically decreasing function f (z) that goes from f (0) = 1 at the AdS boundary to f (z h ) = 0 at the horizon.
In this work, we describe the finite temperature extension of the EHQCD model of [1]. We restrict ourselves to the model of type A1, where the dilaton profile is given by (2.10) As described in [1], in the zero temperature case the EHQCD background obtained from the dilaton profile (2.10) gives rise to a glueball spectrum in agreement with lattice QCD results [31] and the results from improved holographic QCD models [15]. As stressed in [1], the glueball spectrum does not change significantly when we consider the conformal anomalous dimension, , in the interval ∈ [10 −3 , 10 −1 ]. One of the aims of this paper is to investigate if the relevant thermodynamic variables are sensitive to the variation of the conformal anomalous dimension. The parametersφ 0 and Λ are related to the source φ 0 , and vacuum expectation value G by (2.11) The parameter C = Λ 2 appears as a coefficient in the large z behavior Φ(z) = Cz 2 and it is related to color confinement. The parameter Λ, with conformal dimension 1, plays a role similar to Λ QCD and fixes the units of the EHQCD model. The parameterφ 0 is the dimensionless version of the source (coupling) φ 0 . Thus, the relevant parameters of the model areφ 0 , , and Λ. Note that the parameter Λ controls the breaking of conformal symmetry, in the limit Λ → 0 the dilaton field vanishes, and conformal symmetry is restored. In this paper, we will work most of the time with dimensionless quantities so there will be no need of fixing Λ. In other words, we work in units where Λ = 1.
It is worth mentioning that, although the parameters andφ 0 are in general independent, as shown in [1], at zero temperature the glueball spectrum is a physical constraint to fitφ as a function of . A similar scenario occurs in the finite temperature extension considered in this work, where initially we take andφ 0 as independent parameters and will find that the thermodynamic properties become a physical constraint that fixesφ 0 as a function of . Finally, it is also worth pointing out that in the limitφ 0 → 0 a massless mode arises in the scalar sector [32].

Thermodynamics of the non-conformal plasma
In this section we investigate the thermodynamic properties of the black hole solutions arising in our EHQCD model, which will be interpreted in terms of the 4d non-conformal plasma.

Temperature and entropy
After finding the blackening function (2.8), it is possible to obtain the Hawking temperature of the black hole using the well-known relation According to the AdS/CFT correspondence, the black hole temperature obtained from the gravitational solution should be equal to the temperature of the dual quantum field theory. In general, the temperature is determined numerically after setting the model parameters.
In Appendix A we write the asymptotic behavior of the temperature for small z h (close to the boundary), see Eq. (A.5), and for large z h (far from the boundary), see Eq. (A. 19).
From the asymptotic behavior, we conclude that the temperature is non-monotonic and hence must has a minimum in the intermediate region T = T min at z h = z hc . In the following, we consider the case where the critical temperature for deconfinement in the non-conformal plasma is T c = T min . We shall see below that there are two black hole solutions for each temperature T > T c , although, as we will show later, one of them is unstable.
In order to investigate the effects of finite temperature in the dual field theory, we need to calculate some of the thermodynamic variables, such as entropy, free energy, and specific heat. In this paper, we follow a method commonly used in the literature (see, for instance, [13,16], for an earlier study see [33]) in which the entropy density of the 5d black hole, obtained from the Bekenstein-Hawking area formula, is identified with the entropy density of the dual field theory. From the knowledge of the entropy density and the use of the first law of Thermodynamics, one then obtains the free energy density which, in turn, allows us to calculate the other thermodynamic quantities.
To obtain the entropy density, we consider the transverse components of the metric (2.6) and calculate the event-horizon area through where V 3 represents the volume of the transverse Euclidean space. We may now calculate the entropy, using the Bekenstein-Hawking formula, where G 5 is the five-dimensional (5D) Newton constant, which is given in terms of the model parameters by G 5 = 1/(16πσ). Finally, by substituting (3.1) into (3.3) we find the thermodynamic relation for the entropy density s = S/V 3 , where the relation M 3 p = 1/(45π 2 ) has been used in the last equality [16].

The free energy and the trace anomaly
So far we have obtained the temperature and entropy density. Now we use the first law of Thermodynamics 2 , dF = −s dT, (3.5) in order to calculate the Helmholtz free energy density F . We may integrate this equation to calculate the free energy density of the black hole. Note that the entropy density and temperature are explicit functions of z h and therefore we may rewrite the integration of (3.5) by using z h as the independent variable, where z ha is an arbitrary initial value for the horizon radius with the corresponding to the free energy F a . One might consider the prescription where the free energy density F a is set to zero in the limit z ha → ∞, as done in Ref. [16]. Since in this limit the horizon function f (z) tends to unity, such a prescription is equivalent to setting the free energy density of the (confining) thermal solution to zero. In this work, however, we consider an alternative scheme where we integrate from a finite value of z h , namely, from z ha ≡ z hc , so that the integral representation of the free energy becomes Furthermore, set the free energy at the critical temperature T (z hc ) = T c to zero. This assumption guarantees that the free energy of the physical black hole solution, namely the large black hole with z h < z hc , satisfies the physical constraint F ≤ 0. This constraint results because in that regime the entropy density and the temperature are monotonically decreasing functions of z h . Interestingly, the free energy density of the small black hole solution, corresponding to the regime z h > z hc , will also satisfy the constraint F ≤ 0 because an inversion of the monotonic behavior of the temperature will compensate the inversion of integration limits. In particular, the free energy density of the small black hole solution will not vanish in the limit z h → ∞, where the small black hole solutions reduce to the confining thermal solutions 3 .
After integrating by parts, we can rewrite (3.7) in the form Note that the l.h.s is precisely the definition of the energy density, which is given by Therefore, the energy density may be calculated directly from the formula We will be interested in investigating finite temperature effects on the trace anomaly of the energy-momentum tensor which is given by T µ µ = ρ − 3p, where the pressure is p = −F . The trace anomaly can be written in the integral form In arriving to this equation we have considered the difference T µ µ (z h ) -T µ µ (z * h ) and used the property that the trace anomaly vanishes in the limit z * h → 0, assuring that conformal symmetry is restored.
The behavior of the trace anomaly at zero temperature was investigated in Ref. [1]. In the next subsection, we are going to investigate the effects of the temperature on such an important quantity. It is worth mentioning that, even though in some asymptotic regimes the analytic solutions can be obtained, as described in Appendix A, in the following analysis the thermodynamic variables shall be determined numerically,

Numerical results
Let us turn the attention to our numerical results for the thermodynamic variables. We will describe the thermodynamic variables and coordinates in units where Λ = 1. For different choices of Λ, for example, the values considered in [1], the thermodynamic variables and coordinates must be properly rescaled. For example, the temperature T in units where Λ = 1 is actually the dimensionless ratio T /Λ for a different choice of Λ.
In the following analysis, we consider the parameters andφ 0 as independent of each other. One way of describing the evolution of the thermodynamics quantities with the model parameters is fixing one of the parameters and varying the other. Consider, for example, setting the conformal anomalous dimension to = 0.07 (the value used in Ref. [13]) and varying the dimensionless couplingφ 0 in the region 0 ≤φ 0 ≤ 6.
The minimum temperature for the existence of black hole solutions, T c , is obtained by solving the equation ∂ z h T = 0. On the left panel of Fig. 1, we show that T c monotonically increases withφ 0 . A plot of the temperature as a function of z h , both in units where Λ = 1, for the few particular valuesφ 0 = {0, 1, 2} is displayed in the right panel of Fig. 1. The two branches, one of them corresponding to the large black hole regime (solid lines), and the other corresponding to the small black hole regime (dashed lines) are clearly seen. The inset figure displays a zoom of the region where the temperature reaches a minimum. We conclude that the results are sensitive to the value ofφ 0 , namely the critical temperature T c increases with the increasing ofφ 0 .  The numerical results for the dimensionless entropy s/(N 2 c T 3 ) as a function of the dimensionless temperature T /T c are displayed in Fig. 2. We observe that the dimensionless entropy is even more sensitive to the value of the parameterφ 0 , it decreases with the increasing ofφ 0 . As before, continuous lines represent the large black hole regime, while dashed lines correspond to the small black hole regime. Note that the caseφ 0 = 0 reaches the conformal limit, indicated by the dashed thin horizontal line, faster than the other cases.
So far, we have two black hole regimes. In order to discover which regime is physically relevant, we must investigate the thermodynamic stability of the corresponding solutions. A standard method for investigating the stability of a thermodynamic system is to calculate the corresponding specific heat, defined by the relation Stable systems are characterized by a positive C V . On the left panel of Fig. 3 we display the results of the dimensionless specific heat, , as a function of the dimensionless temperature T /T c for different values ofφ 0 . Solid lines represent the large black hole solutions, for which the specific heat is always positive (stable regime) whereas dashed lines represent the small black hole solutions, for which the specific heat is always negative (unstable regime). The figure also shows the dependence of the specific heat on the parameterφ 0 . For the large (small) black hole, increasing the value ofφ 0 leads to a decrease (increase) of the value of C V /(N 2 c T 3 ) for a fixed temperature.  Another relevant thermodynamic variable is the speed of sound, which may be calculated using the formula Our numerical results for c 2 s as a function of the temperature are displayed on the right panel of Fig. 3 for different values ofφ 0 , as indicated, and = 0.07. Interestingly, in the region of unstable black holes (dashed lines), c 2 s becomes negative, meaning that the speed of sound becomes imaginary. This is related to the fact that the specific heat is always negative in this region. To be more specific we may rewrite the speed of sound in a convenient form, (3.14) Since the entropy density is always positive, c.f. Fig. 2, a negative value of C V gives rise to a negative value for c 2 s , hence the results found for the small black holes on the left and right panel of Fig. 3 (dashed lines) are consistent. We then conclude that the instability of the small black hole is characterized by a negative specific heat and a purely imaginary speed of sound. Our numerical results also show that at T = T c , c 2 s crosses the horizontal axis (c 2 s = 0) whereas C V diverges. This is consistent with our formula (3.14). Note that close to T = T c , c 2 s increases rapidly with T and far from T c it varies slowly. Last but not least, we conclude that the speed of sound is also sensitive to the variation ofφ 0 . For the large (small) black hole solution c 2 s decreases (increases) for increasingφ 0 . Note also that c 2 s converges to 1/3 for the (physical) large black holes in the limit of very high temperatures, which is consistent with the restoration of conformal symmetry. Using our formula (3.7) for the free energy density, we can evaluate the pressure P = −F . Our numerical results for the dimensionless pressure 3p/(N 2 c T 4 ) as a function of the dimensionless temperature T /T c are displayed in Fig. 4 for = 0.07 and different values ofφ 0 . Note that the sensitivity to the parameterφ 0 depends on the region of interest. In the branch describing small black holes (dashed lines), the pressure is less sensitive toφ 0 far from the critical temperature T c . Meanwhile, in the branch of large black holes (solid lines), the pressure is more sensitive toφ 0 for temperatures much higher than T c . Our numerical results for the energy density ρ as a function of the temperature are displayed on the left panel of Fig. 5, where we plot the quantity ρ/(N 2 c T 4 ) for = 0.07 and different values ofφ 0 . As the figure shows, the energy density is also sensitive toφ 0 , it decreases with the increasing ofφ 0 for the large and small black hole regimes (solid and dashed lines respectively). We also observe that the energy density increases rapidly with T close to the critical temperature T c and varies slowly in the regime of high temperatures. The right panel of Fig. 5 displays our results for the dimensionless trace anomaly, , as a function of the dimensionless temperature T /T c for = 0.07 and different values ofφ 0 . We can see that the trace anomaly decreases with the increasing ofφ 0 . Note that for the large black hole regime (solid lines) the trace anomaly has a peak close to the critical temperature, then, it decreases and goes to zero in the limit of very high temperatures recovering conformal symmetry. For the small black hole regime (dashed lines) the trace anomaly becomes negative as T increases and goes to zero in the high temperature (conformal) limit. The main conclusion of the above analysis is as follows. When the conformal anomalous dimension is fixed, all the thermodynamic variables describing the non-conformal plasma in EHQCD are sensitive to the variation of the dimensionless couplingφ 0 . We did a similar analysis for the case where we fixφ 0 and found that the thermodynamic variables are sensitive to the variation of .
Other important lessons that can be learned from the analysis presented in this section are the following. Holographic QCD models that describe CFT deformations consistent with confinement at zero temperature lead to the formation of non-conformal plasmas only above a certain temperature T c . In other words, there is a minimum temperature for gluon deconfinement due to the CFT deformation. Moreover, above T c the black hole solution splits into two branches: a large stable black hole characterized by a positive specific heat and a real speed of sound and a small black hole characterized by a negative specific heat and an imaginary speed of sound. In our bottom-up framework, CFT deformations have been described by an effective 5d Einstein-dilaton theory. Interestingly, qualitatively similar results can be obtained considering a top-down approach, e.g. [34,35] 4 Viscosity coefficients of the non-conformal plasma There are at least three methods for calculating transport coefficients of a fluid arising from a strongly coupled theory within the framework of holography. One of them is through the hydrodynamic limit of the black hole quasinormal modes, where one calculates the dispersion relations of field perturbations on the gravitational side of the duality [10,36]. These relations are then compared to the ones obtained from the linearized hydrodynamic modes. In general, the dispersion relations arising from the relativistic hydrodynamics are constructed by using the gradient expansion method (see, for instance, Refs. [37][38][39][40]). Notice that in order to calculate the viscosity coefficients we need to consider first-order hydrodynamics.
A second method consists of the direct evaluation of the retarded Green's functions for the components of the stress-energy tensor associated with the transport coefficients.
Using the AdS/CFT correspondence, one identifies the dual 5d gravitational perturbations. Then, the action is expanded up to second-order in the perturbations to read off the retarded Green's function. Finally, Kubo's formula is used to calculate the desired transport coefficients. This is the approach used in the seminal work [41] for the shear viscosity and [14] for the bulk viscosity.
A third method is based on the fluid/gravity correspondence [42], where the full stressenergy tensor can be obtained from the 5d metric by solving the full non-linear Einstein equations via a gradient expansion. In this paper, we will follow the second method described above to calculate the shear and bulk viscosities of the non-conformal plasma arising in our EHQCD model.

Shear viscosity
In the case of the shear viscosity (η), considering the direction of propagation k µ = (ω, 0, 0, q), the relevant retarded Green's function is given by where T x 1 x 2 is one of the energy-momentum tensor components and ω (q) is the frequency (wavenumber) of the perturbation. From the holographic dictionary, we know that the source of the shear viscosity is related to the metric perturbation h x 1 x 2 , which couples to the component T x 1 x 2 of the energy-momentum tensor, relevant for calculating η.
In the following, we develop a general procedure to calculate the shear viscosity in the holographic model we are working with. Let us start with a general ansatz for the black hole metric in holographic QCD, which we write as where f , ζ 1 , and ζ 2 are functions of z only. We use this ansatz for the metric in order to compare with previous results in the literature and check the consistency of our procedure and results. The corresponding background equations are As expected, these relations reduce to (2.7) for ζ 1 = ζ 2 .
In order to calculate the shear viscosity, the next step is to consider a perturbation on the background black hole metric, g mn → g mn + h mn . To calculate the shear viscosity we need to consider the component h x 1 x 2 only, and this perturbation decouples naturally from the others. The relevant equation governing this sector is obtained from where G (1) x 1 x 2 and T (1) x 1 x 2 are first-order contributions to the Einstein and energy-momentum tensor expansions in h x 1 x 2 . Since x 3 is the direction of propagation of fluctuations in the transverse space, the differential equation for h x 1 x 2 (t, x 3 , z) may be written as However, this equation is not gauge-invariant. As discussed in Refs. [10,36] we may write it in terms of a gauge-invariant master field, the so-called Kovtun-Starinets (KS) master variable, defined by Z T = g x 1 x 1 h x 1 x 2 . In doing so, we obtain Considering the Fourier transform the fundamental equation becomes At this point, it is interesting to observe that considering ζ 1 and ζ 2 appropriately, equation (4.8) reduces to previous results presented in the literature. For example, in the conformal case, the dilaton field is zero and one may compare this equation with Eq. (6.6) of Ref. [10]. It is also interesting mentioning that equation (4.8) takes the same form as Eq. (2.11) of Ref. [14], setting q = 0. Finally, following the procedure implemented in Ref. [43], we may write this equation in terms of the Regge-Wheeler-Zerilli master variable 4 . To find the retarded Green's function (4.1) we now expand the action up to secondorder in the fluctuation, h x 1 x 2 . Then, we write the resulting on-shell action in terms of the gauge-invariant field Z T in the form where the Lagrangian is given by (4.11) As discussed in Ref. [14], it is possible to include the contribution of the Gibbons-Hawking-York surface term to the Lagrangian L (2) by adding a non-trivial scalar function G(z), such that we rule out contributions of the form Z 2 T and Z T Z T . In turn, as G(z) is an arbitrary (unknown) function, this will generate an ambiguity. However, the profile of G(z) may be fixed by phenomenology, eliminating the aforementioned ambiguity. Then, adding the contribution of this function to the Lagrangian we get Nevertheless, the additional term of the "improved" Lagrangian will not contribute to the imaginary part of Green's function. As we shall see, the shear viscosity coefficient depends on the imaginary part of Green's function; then, the explicit form of G(z) does not matter for our calculation. The next step forward is to promote Z T to a complex function, and this is because, in general, an arbitrary solution of Eq. (4.8) might be complex, depending on the boundary conditions of the problem. Thus, we rewrite the Lagrangian using the master variable and its complex conjugate, The form of this Lagrangian provides the same equation of motion (4.8). It is also possible to rewrite the last Lagrangian in a form where the presence of a surface term is evident, where At the boundary, J must be related to the retarded Green's function. To calculate the imaginary part of the Green's function (4.1) we need to define the number flux of gravitons, F, which is related to the imaginary part of J through As described in [14], the quantity F represents the number flux of gravitons in the radial direction and it is the conserved charge associated with the U (1) symmetry of (4.13). We shall use this expression in the last part of this section.
Let us now focus on the case we are dealing with, where ζ 1 = ζ 2 . To calculate the shear viscosity we follow the procedure implemented in Ref. [14]. First, we solve Eq. (4.8) in the limit of zero frequency and wavenumber, i.e. (ω, q) → 0, where the resulting differential equation has an exact solution, given by Here C 1 and C 2 are integration constants, and in the limit ω → 0 we must set C 1 = 0 in order to guarantee the regularity condition at the horizon, where f (z h ) = 0. Without loss of generality, we also impose the Dirichlet boundary condition Z T = 1 at the boundary, so that C 2 = 1. Now consider the case of small ω (with q = 0) where we solve Eq. (4.8) perturbatively, considering ω as the expansion parameter. In that case, the solution takes the same form as Eq. (4.17), where C 1 (ω) depends on the frequency. Then, expanding this solution close to the horizon we obtain the approximate solution The next stage is to solve the differential equation (4.8) close to the horizon, for arbitrary ω and q. To do so, we use the ansatz Z T = f β , where β takes the values: where f (z h ) = −4πT with T being the black hole temperature. We choose the solution associated with β 1 because it represents waves falling into the black hole. This condition is also related to the retarded Green's function. Then, considering the first subleading term close to the horizon, Z T is given by 20) or, expanding f (z) around z = z h , where F 0 and A − are constants. Considering q = 0 and expanding the leading term around ω = 0, we get This approximate solution must be equal to Eq. (4.18). Thus, we identify the corresponding coefficients: (4.23) Note that C 1 vanishes in the limit ω → 0, as expected. Having found the asymptotic solution for Z T near the horizon in (4.21) we can evaluate the (conserved) graviton number flux F in the limit z → z h . Plugging the result (4.21) for Z T into the number flux formula (4.16) we get (setting q = 0) .

(4.24)
The imaginary part of the retarded Green's function takes the form .

(4.26)
Just like the entropy density, the shear viscosity coefficient η depends on the (inverse) scale factor evaluated at the horizon. In Fig. 6 we display our numerical results for the shear viscosity, normalized as η/(N 2 c T 3 ), as a function of the dimensionless temperature T /T c for different values of the dimensionless couplingφ 0 , for a fixed value = 0.07 of the conformal anomalous dimension. We observe that close to the critical temperature the shear viscosity increases rapidly with the temperature whereas far from the critical temperature it varies slowly. We conclude that the shear viscosity is also sensitive to the value ofφ 0 , for fixed . This ratio is expected to change when higher curvature terms are added in the fivedimensional action, see for instance [44]. We would like to remark that, in contrast with Ref. [14], we did not need to redefine our radial coordinate z in terms of the dilaton field Φ. This is true whenever we express the perturbation in terms of a gauge-invariant master field. In conclusion, the master field Z T (z) is equivalent to the field H 12 (Φ) used in Ref. [14].

Bulk viscosity
Following the conventions of Ref. [14], the retarded Green's function related to the bulk viscosity is given by where T i i is the trace of the spatial part of the energy-momentum tensor. The holographic dictionary then indicates that the source of the bulk viscosity is related to the trace of the metric perturbations.
The bulk viscosity has been investigated previously in holography (see, for instance, Refs. [14,17,[45][46][47][48][49][50][51][52][53][54]). This transport coefficient is associated with the expansion scalar term ∂ µ u µ of the fluid stress-energy tensor and associated with conformal symmetry breaking. Here we follow the procedure implemented in Ref. [14] adapted to our case. The gauge where the dilaton field becomes the holographic coordinate simplifies dramatically the problem of calculating this transport coefficient. Such a choice is possible whenever the dilaton is a monotonic increasing function. Thus, within this gauge the metric (4.2) becomes It is worth pointing out that the metric (4.29) reduces to (2.6) considering ζ 2 = Φ (z)ζ 1 .
In this gauge, we do not need to care about the perturbation of the scalar field, which naturally shall be coupled to the trace of the metric perturbations. We restrict ourselves to the zero spatial momentum case (q = 0); thus, the perturbed metric depends on the time and dilaton field only, and the metric components are given in explicit form by where (4.31) Here ζ 2 , ζ 1 and f are functions of the dilaton field, whereas λ (λ 1) is a small parameter introduced to control the expansion. In turn, the background equations are: (4.32) Plugging (4.31) into the Einstein equations (2.3), and then performing a Fourier transform, we get the corresponding perturbation equations: (4.33c) As can be seen from the equations, once we solve the differential equation (4.33c) we may automatically calculate H 00 and H 55 . Using the background differential equations (4.32) we may rewrite Eq.(4.33c) in the reduced form (4.34) To obtain the bulk viscosity we need to calculate the retarded Green's function related to these metric perturbations. Thus, analogously to what has been done with the shear viscosity, we substitute the metric (4.31) in the action (2.2) and then expand up to second-order in λ. The resulting on-shell action may be written as The Lagrangian is given by where To get this Lagrangian we have used the equality The matrices are given by whereas L Φ and L t can be written as In this case, G is a symmetric matrix whose elements may depend on the dilaton. Such a function has the same role as the one introduced in Eq. (4.13). Thus, it will not contribute to the imaginary part of the retarded Green's function (4.28). In analogy to what was done in the previous section, we write down the improved Lagrangian where we have promoted the functions to complex ones and also defined The next step is to rewrite this Lagrangian in an analogous form to Eq. (4.14), where a surface term is evidenced. The result is Therefore, the flux number of gravitons associated with rotationally invariant perturbations, i.e., the imaginary part of J, is given by As observed in Ref. [14], the flux number F is proportional to the Wronskian of the complexified solution H 11 . Therefore, in order to calculate F, we only need the solution of H 11 . Moreover, asymptotic solutions may be obtained considering the asymptotic behavior of the background. Close to the horizon, H 11 behaves like H 11 ∝ f α , where f is a function of the dilaton field Φ. Plugging this ansatz into Eq. (4.35) we find two possible solutions: , (4.45) where Φ h is the location of the horizon. In terms of the temperature, defined by the first solution in (4.45) represents waves falling into the black hole horizon, while the second solution represents waves coming out from the horizon. Thus, we choose α 1 and we may write H 11 (Φ) close to the horizon as where B − is a constant. Expanding the last function around ω = 0, we get Note that for ω = 0, the function H 11 becomes regular (constant) at the horizon. On the other hand, close to the boundary we consider the ansatz H 11 ∝ Φ β , plugging in Eq. (4.33c) we get two solution: β 1 = 0 and β 2 = (4 − 2 )/ . Thus, the leading term in the asymptotic solution close to the boundary must be a constant, which can be fixed as Now we are able to find an expression for the flux number F. Plugging (4.47) in (4.43) we get the simplified expression (4.50) Note that F depends on the (inverse) scale factor, its derivative, and the constant |B − |.
The imaginary part of the retarded Green's function takes the form (4.51) To calculate the bulk viscosity we follow the procedure of Ref. [14] and find that (4.52) To write the last result we have used the entropy density relation (3.4).
A few comments are now in order. The bulk viscosity depends on the value of the inverse scale factor (and its derivative) evaluated at the horizon. Moreover, in contrast with the shear viscosity case, the bulk viscosity depends on the constant |B − |. It is worth mentioning that in the high-temperature regime the constant reduces to |B − | = 1. Thus, in this region, the result (4.52) is in agreement with the formula obtained in Ref. [49]. However, there is a subtlety when one compares both results, see the discussion in Ref. [50]. This is also true in the adiabatic approximation [17,49] (see also [50]). Considering |B − | = 1, we display our numerical results for the bulk viscosity to entropy density ratio ζ/s in Fig. 7. The solid lines represent the stable large black holes while dashed lines represent the unstable small black holes. We conclude that in the case |B − | = 1 the ratio ζ/s is less sensitive to the value ofφ 0 . Note that the bulk viscosity has a sharp rise close to the critical temperature. This result was previously reported in the QCD literature, see for instance Ref. [55] where a semi-analytic study is presented.
On the other hand, we may calculate the constant |B − | using the boundary conditions of the problem, the Dirichlet condition at the boundary H 11 (Φ = 0) = 1, and the incoming wave condition at the horizon given by (4.47). Following this procedure, we calculate numerically |B − | as a function of Φ h . The results are displayed on the left panel of Fig. 8. Meanwhile, the right panel of Fig. 8 shows the bulk viscosity using |B − | obtained numerically forφ 0 = 2 and = 0.3, solid lines represent the large black hole regime, while dashed lines small black hole regime. It is worth pointing out that depending on the combination of  Knowing the asymptotic form of the inverse scale factor we are able to find an expression for the bulk viscosity to entropy density ratio in the regime of large black holes. In this regime the derivative of the horizon function may be neglected, this means that the coefficient of H 11 (Φ) in Eq. (4.34) is zero in the limit of zero frequency. Thus, we may solve the resulting differential equation getting the asymptotic solution, (4.53) Hence, the boundary condition H 11 (0) → 1 reduces the solution to c 1 = |B − | = 1. Plugging this result and the leading term of the asymptotic warp factor (A.2) into (4.52), we get From this expression, it is easy to see the role of the parameter Λ. Notice that the bulk viscosity vanishes in the limit Λ → 0 where the dilaton field (2.10) vanishes, recovering conformal symmetry. We obtained numerical results for (4.52) considering |B − | = 1 and |B − | obtained numerically, an overlap of both results is displayed in Fig. 9, where the blue line represents the case |B − | = 1, while the red line represents the case where |B − | is obtained numerically. As it can be seen from the figure, both results are in agreement in the region of high-temperatures for both the large black holes (solid lines) and small black holes (dashed lines).

Comparison with lattice SU (N c ) theories and heavy ion collisions
In this section we compare our results for the thermodynamics quantities against the results from lattice SU (N c ) theories [24,56,57] and we also compare our results for the viscosity coefficients against the results from the JETSCAPE collaboration [25], obtained from a model-to-data analysis of heavy-ion collision experimental data. For this comparison, we will only be interested in the (physical) stable black hole solution, namely the large black hole.

Thermodynamic quantities
Let us start this subsection by describing the critical temperature for the formation of a non-conformal plasma in our EHQCD model. In our approach, we identify the critical temperature with the minimum temperature for the existence of black hole solutions.
Considering the parameters fixed by the glueball spectrum, obtained in Ref. [1], we calculate the critical temperature and observe its dependence on the conformal anomalous dimension . The critical temperature is obtained by solving the equation ∂ z h T = 0, and it is represented by a black dot in the plot displayed on the left panel of Fig. 10. The plot describes the variation of the temperature (in MeV) as a function of z h , whereφ 0 and Λ (in MeV) were fixed for a given from a fit to the glueball spectrum, as done previously in Ref. [1]. For the particular case shown in the figure, we have = 0.1, φ 0 = 5.6 and Λ = 743 MeV. As described in the previous sections, there are two branches: one corresponds to the large black hole (solid line) and the other corresponds to the small black hole (dashed line).
In turn, the right panel of Fig. 10 shows the critical temperature (in MeV) as a function of the conformal anomalous dimension . For each value of the parametersφ 0 and Λ (in MeV) were fixed appropriately from the glueball spectrum [1]. The interval considered here is ∈ [10 −3 , 0.5], extending the fit to the glueball spectrum performed in [1]. As the figure shows, the critical temperature is a slowly growing function of the conformal anomalous dimension and varies from T c ≈ 263 MeV to T c ≈ 274 MeV. It is worth pointing out that the present results for the critical temperature are very close to the recent results found in lattice SU (N c ) theories [56,57], namely T c ≈ 0.595 √ σ ≈ 262 MeV where σ ≈ (440 MeV) 2 is the phenomenological value for the string tension. It is remarkable that a holographic model with the metric coupled to a single scalar field provides such a result. This analysis also suggests an alternative approach for fixing the parameter Λ. We will see later in this section that for a given value of the anomalous conformal dimension we can use the lattice results for the trace anomaly to fix the dimensionless couplingφ 0 . Then we can also use the lattice result for T c to fix Λ in MeV units. In conclusion, the parametersφ 0 and Λ can be obtained as a function of considering either the glueball spectrum at zero temperature, as done in Ref. [1], or the thermodynamics of the non-conformal plasma at finite temperature. Next, we compare our results for the pressure of the non-conformal plasma against the results obtained in lattice SU (N c ) theories [24]. Our results for the pressure, properly normalized, as a function of T /T c are displayed in the left panel of Fig. 11 (solid lines) and compared against the lattice results (dotted lines with error bars). As the figure shows, the pressure in our EHQCD model is sensitive to the value of the parameterφ 0 , once the parameter is kept fixed. We present results for two values ofφ 0 : one of them being fixed with the glueball spectrum at zero temperature, as in [1], and the other being fixed by matching the maximum value for the dimensionless trace anomaly (ρ − 3p)/(N 2 c T 4 ) with the lattice result 0.381 found for SU (8) [24]. Although displaying a qualitative agreement with the lattice results, we observe that the values of the pressure in the case of parameters fixed by the glueball spectrum, i.e., for { ,φ 0 } = {0.3, 2} (orange dashed line), are far from a quantitative agreement. This apparent shortcoming was also observed in other holographic models for QCD, see for instance [58]. In turn, the results for the pressure in the case the parameters were fixed by the trace anomaly condition, i.e., for { ,φ 0 } = {0.3, 0.9} (orange solid line), are in quantitative agreement with those obtained from QCD on the lattice in the limit of large N c . We have found in this section that our EHQCD model with parametersφ 0 and Λ fixed as a function of from a fit to the glueball spectrum provides at finite temperature a good result for the critical temperature for deconfinement and pressure that displays a qualitative agreement with lattice SU (N c ) theories. However, for the pressure we could not obtain a quantitative agreement; the pressure in our EHQCD model is lower than the lattice result when the parameters are fixed by the (zero temperature) glueball spectrum. We also found a trace anomaly that has a qualitative behavior similar to that obtained in lattice SU (N c ) theories but presents a peak that is lower than the corresponding lattice result. Since the main goal of holographic QCD models is to provide an effective description of QCD in the limit of large N c , the fact that we found a discrepancy with the lattice SU (N c ) results for the thermodynamic quantities suggests that some more ingredients may be needed in our EHQCD approach.
A phenomenological solution for the discrepancy was considered in Ref. [59], where additional parameters were added to control the curve of the pressure in the region of high temperatures and the height of the peak of the trace anomaly close to the critical temperature. Our analysis suggests the following alternative approach. Once the parameter was fixed, we can fix the dimensionless couplingφ 0 imposing the condition (ρ − 3p)/(N 2 c T 4 ) = 0.381, which corresponds to the maximum value of the dimensionless trace anomaly obtained in the lattice SU (8) theory [24]. The running ofφ 0 as a function of is displayed on the left panel of Fig. 12 (solid line) and compared to the case whereφ 0 was fixed from a fit to the glueball spectrum (dashed lined). In turn, the parameter Λ can be fixed by matching the critical temperature for deconfinement T c with the value 262MeV obtained in lattice SU (N c ) theories in the limit of large N c . The dependence of Λ on is displayed on the right panel of Fig. 12 (solid line) and compared to the case where Λ was fixed from a fit to the glueball spectrum (dashed line).
We therefore reach the following conclusion. In our EHQCD model, when one fixes the conformal anomalous dimension there are two possible methods for fixing the parameterŝ φ 0 and Λ. The first method uses a fit to the lattice results for the glueball spectrum at zero temperature to fixφ 0 and Λ. This is the method we followed in our previous work [1]. The parameters fixed in that way are displayed as dashed lines in Fig. 12. The second method, found in this work, uses the results for some thermodynamics quantities describing the gluon plasma in lattice SU (N c ) theories to fixφ 0 and Λ, namely the maximum for the dimensionless trace anomaly to fixφ 0 and the value for the deconfinement temperature to fix Λ. The parameters fixed by following the second method are displayed as solid lines in Fig. 12.
Note that the behavior ofφ 0 ( ) is qualitatively similar for both methods (solid and dashed lines). As regards the behavior of Λ( ), in the first method we find that Λ is a slowly growing function of whereas the second method yields Λ as a slowly decreasing function of . On physical grounds, since the parameter Λ represents a dynamical mass gap similar to Λ QCD it is expected that it should be independent of the parameter . The fact that we find a slow variation for Λ( ) suggests that some further improvement would allow us to satisfy this physical condition.
We would like to remark that the numerical errors for the thermodynamic quantities in our model are very small. This is because the only differential equation that is solved numerically is the first equation in (2.7) that allows us to find the (inverse) scale factor ζ 1 as a function of z for a given dilaton profile Φ(z). Considering a moderate numerical precision, the uncertainty in ζ 1 is of order 10 −6 . The entropy density is inversely proportional to ζ 1 (z h ) 3 and therefore has an uncertainty of the same order. All the other thermodynamic quantities were obtained from the entropy density.

Shear and bulk viscosity
We end this section by comparing our results for the viscosity coefficients in the regime of large black holes against the results obtained by the JETSCAPE collaboration [25] (see also [60]), from a model-to-data analysis of by e heavy-ion collision data. The universal result for the shear to entropy density ratio η/s reproduced in this work, is displayed by the solid black line on the left panel of Fig. 13, while the dashed blue lines represent the upper and lower bounds found by the JETSCAPE collaboration. As it can be seen  Left: Shear viscosity to entropy density ratio as a function of the temperature. The black solid line represents the holographic result. The region bounded by dashed blue lines corresponds to the interval constrained by the JETSCAPE collaboration [25]. Right: Bulk viscosity to entropy density ratio ζ/s as a function of T forφ 0 = 2, and = 0.3. The black line was drawn by taking |B − | = 1. The red line represents the results for ζ/s by taking the value of |B − | determined numerically from the corresponding critical temperature T c = 0.269 GeV.
The bulk viscosity to entropy density ratio ζ/s is displayed on the right panel of Fig. 13. The sharp rise close to the critical temperature is in agreement with previous calculations on the lattice [23] (see also [61]). Using Eq. (4.52) with |B − | = 1, the maximum value of the ratio ζ/s ≈ 0.0695 (black line). Meanwhile, considering |B − | obtained numerically, the maximum value of the ratio is ζ/s ≈ 0.125 (red line), for the set of parameters { ,φ 0 , Λ} = {0.3, 2, 759 MeV} fixed by a fit to the glueball spectrum. Compared against the lattice results of [23], our results are within the error bars. In turn, our results are a bit larger than the results of [61]. However, our results fit very well into the expected region presented recently by the JETSCAPE Collaboration [25,60] enclosed by dashed blue lines. Compared against the holographic models of Refs. [14,17], the result for |B − | = 1 is of the same order, but for |B − | obtained numerically our result is larger than the one obtained in these papers, see also Refs. [17,[51][52][53] for additional discussions. Hence, we conclude that our results are in agreement with those results obtained previously in the literature.
We would like to remark that our result for the shear viscosity to entropy ratio η/s = 1/(4π) is exact whereas our result for the bulk viscosity to entropy ratio ζ/s is obtained numerically using the formula (4.52). According to this formula, the bulk to entropy ratio ζ/s only depends on powers of the inverse scale factor ζ 1 , its derivative, and the coefficient B − . For a moderate numerical precision, we find that the uncertainty in the scale factor ζ 1 is of order 10 −6 and the uncertainty in the coefficient B − is of order 10 −4 . Thus, the uncertainty in ζ/s is approximately of order 10 −4 . This uncertainty is very small compared with the uncertainties found in lattice QCD.

Conclusions
We have described the finite-temperature extension of the effective holographic models for QCD (EHQCD) proposed in [1] in order to describe the physics of a non-conformal plasma. At zero temperature the EHQCD model provides a spectrum of scalar and tensor glueballs in agreement with results obtained in lattice QCD [31]. The finite-temperature extension consists in embedding a black hole solution into the gravitational Einstein-dilaton theory, which is equivalent to creating a thermal state in the dual quantum field theory. We calculated some of the relevant thermodynamic variables, which are required to investigate the stability of the black holes. We showed that the large black holes are thermally stable, while the small black holes are thermally unstable. The unstable black hole is characterized by a negative specific heat and an imaginary speed of sound. We also showed that all the relevant thermodynamic quantities are sensitive to the variation of the model parameters, namely the dimensionless couplingφ 0 and the conformal anomalous dimension . These two parameters characterize the breaking of conformal symmetry in EHQCD. Interestingly, we found that the pressure and the trace anomaly display qualitative behaviors that are similar to the ones found in lattice SU (N c ) theories. In particular, the trace anomaly displays a peak near the critical temperature for deconfinement. In the limit of very high temperatures, all thermodynamic quantities reach the corresponding conformal values.
We also investigated the viscosity coefficients associated with the transport properties of the non-conformal fluid. We recovered the universal result for the ratio between the shear viscosity and the entropy density η/s = 1/(4π). This result was found by writing the relevant equations for the metric perturbations in terms of a gauge-invariant variable, calculating the retarded Green's function, and extracting the shear viscosity through Kubo's formula. Our numerical results show that the shear viscosity and the entropy density behave almost identically with the temperature, rising rapidly close to the critical temperature and then varying slowly in the region of high temperatures. We verified that the shear viscosity is sensitive to the variation ofφ 0 for fixed , and vice-versa, while the ratio η/s remains constant, as expected for this class of holographic models arising from Einstein-dilaton gravity. To study the bulk viscosity we followed a different approach, by adopting a gauge where the dilaton field plays the role of the holographic coordinate. After calculating the retarded Green's function we were able to find an expression for the bulk viscosity by using Kubo's formula. The numerical results indicate that the holographic bulk viscosity increases sharply close to the critical temperature. This result is in agreement with previous predictions from lattice QCD and other holographic models of QCD. We also showed that the bulk viscosity is sensitive to the parameters of the modelφ 0 and . It is worth mentioning that our results were obtained using Model A1 proposed in Ref. [1]. However, we have also considered Model A2 proposed in that work and obtained equivalent results. We decided not to present the results of Model A2 in this paper to avoid redundancy. We believe that Models B1 and B2 would provide equivalent results to those obtained from Models A1 and A2.
Finally, we compared our results on the thermodynamics against the data available from lattice SU (N c ) theories. We also compared our results for the viscosity coefficients against those found by the JETSCAPE collaboration. Regarding thermodynamic quantities, we found that our results for the pressure and trace anomaly are in qualitative agreement with the results found in the literature for lattice SU (N c ) theories. However, we found that the value of the parameters adjusted to fit the glueball spectrum does not provide a quantitative agreement with the thermodynamics of lattice SU (N c ) theories. We found, however, that a quantitative agreement with lattice SU (N c ) theories is possible if one fixes the model parametersφ 0 and Λ, for a given value of , by using the lattice results for the maximum value of the dimensionless trace anomaly and the critical temperature for deconfinement, respectively. We concluded that the results for the viscosity coefficients provided by our EHQCD model are consistent with the phenomenological constraints obtained by the JETSCAPE collaboration from the model-to-data analysis of the heavy-ion collision data. Although the shear viscosity did not always belong to the region bounded by the JETSCAPE collaboration, the bulk viscosity fits very well in the region of parameters considered by the collaboration.
A possible extension of this work would be by including flavor degrees of freedom in order to investigate chiral symmetry breaking. This task shall be reached by adding a non-Abelian SU (N f ) L × SU (N f ) R gauge symmetry (dual to the chiral currents) and a bifundamental scalar (dual to the chiral condensate). Another interesting direction would be investigating the role of a non-minimal coupling in the phase diagram of QCD in the same line of Refs. [28][29][30][62][63][64]. From the gravitational point of view, a natural next step would be the investigation of the quasinormal modes of the black hole solutions found in this work. This would allow us to describe the melting of scalar and tensor glueballs in a non-conformal plasma.
where the ellipsis stands for higher power on z h . Notice that the leading term in (A.4) corresponds to the AdS contribution, while the subleading term reveals the deformation introduced by the nontrivial dilaton field.
One may also write f (z) in terms of z h by plugging (A.4) into (A.3). It is clearly seen that the horizon function reduces to unity in the limit of zero z, as expected.
Let us now calculate the asymptotic expressions of the thermodynamic variables close to the boundary. The asymptotic expansion of the temperature is obtained by plugging (A.4) and (A.2) into (3.1), Again, the leading term is due to the AdS warp factor and the subleading term is the deformation generated by the dilaton field. It is worth pointing out that in the above analysis the natural independent parameter is the coordinate z h . However, in Thermodynamics one usually uses the temperature as being the fundamental degree of freedom. With this in mind, we invert the asymptotic expression (A.5) to get z h as a function of the temperature in the form (A.6) With this relation, we may express all thermodynamic variables as a function of the temperature. Let us then apply the procedure to the entropy density. Once again the leading term corresponds to the AdS warp factor contribution, which is equivalent to recover conformal symmetry, while the subleading terms correspond to the deformation from such symmetry.
To get an asymptotic expression for the free energy we write the integral representation (3.7) in the form where z hc is the value where the temperature reaches its minimal value, T (z hc ) = T c . Equation (A.9) indicates that we may split the free energy in two parts, the first one corresponding to large black holes,z h ∈ [z h , z hc ], and the second one related to small black holes,z h ∈ [z hc , ∞ . To guarantee the validity of the following analysis we rewrite the free energy of the large black holes in the form where z h < z h * < z hc .
In the large black holes regime, the main contribution is expected to come from the first integral in Eq. (A.10), and then we may use the asymptotic expressions for the temperature and entropy density to evaluate it. Thus, the result for F large is given by where O(z h * ) stands for the result of the integral evaluated at z h * , which is a subleading constant term in the limit of zero z h . Hence, as a function of the temperature the free energy becomes Now one may fix the factor (M p ) 3 by comparing F with the corresponding value obtained by using the Stefan-Boltzmann approximation of pure Yang-Mills, F YM = − π 2 45 N 2 c T 4 . Thus, from the leading term of (A.12) we get which is the same value obtained in Ref. [16]. It is worth pointing out that the subleading term in (A.12) may be of the same order as the leading term if the value of ∆ − is small enough, such that T 4 ∼ T 4−2∆ − . In turn, the asymptotic form of the trace anomaly is given by (πT ) 4−2∆ − + · · · (A.14) As it can be seen, the leading term of the trace anomaly goes like T 4−2∆ − . In the particular case where the leading term of the dilaton is linear in the UV, i.e., for ∆ − = 1, this expression reduces to the result obtained in Ref. [65], with the trace anomaly going like T µ µ ∼ T 2 .

A.2 Far from the boundary
So far we have dealt with the asymptotic analysis close to the boundary, and from here on we perform the asymptotic expansion of all relevant quantities far from the boundary, for large z.
Far from the boundary, the dilaton field behaves like the zero temperature asymptotic form Φ = C z 2 + · · · (A.15) Thus, plugging this expression in the Einstein equation (2.7) we get the asymptotic form for the warp factor, ζ 1 = e −A 1 (z) , Now, we may rewrite (2.9) in the form Additionally, we may split up the second integral in the intervals z hc ≤z ≤ z h * and z h * ≤ z ≤ z h . Hence, by plugging (A.16) in (2.9) and expanding the result in the region z h z hc , where z hc is the value where the temperature reaches its minimum, the integration constant may be approximated by This approximate expression confirms the linear behavior observed in our numerical results, as shown by the dashed lines in the left panel of Fig. 2.
In turn, analogously to what we have done in the large black holes regime, we may invert the series (A. 19) in the region of large temperatures. Hence, we get The entropy density in this region is given by By plugging (A.20) in the last result we get the entropy density in terms of the temperature, s σ = e − 2(πT ) 2 C 4π 5/2 T 3/2 C 3/4 + 15C 1/4 π 1/2 4T 1/2 + · · · . (A.22) The leading term is exponentially suppressed, this behavior can be seen in the right panel of Fig. 2 with dashed lines. It is also worth mentioning that the entropy density is always positive. Following the same procedure, it is easy to show that the free energy is given by F σ = e − 2(πT ) 2 C C 1/4 (πT ) 1/2 + 17C 5/4 16(πT ) 3/2 + · · · (A.23) The exponentially suppressed leading term is also observed in our numerical results, as it is shown by the dashed lines in the left panel of Fig. 4.