Stiff phases in strongly coupled gauge theories with holographic duals

According to common lore, Equations of State of field theories with gravity duals tend to be soft, with speeds of sound either below or around the conformal value of υs=1/3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\upupsilon}_s=1/\sqrt{3} $$\end{document}. This has important consequences in particular for the physics of compact stars, where the detection of two solar mass neutron stars has been shown to require very stiff equations of state. In this paper, we show that no speed limit exists for holographic models at finite density, explicitly constructing examples where the speed of sound becomes arbitrarily close to that of light. This opens up the possibility of building hybrid stars that contain quark matter obeying a holographic equation of state in their cores.


Introduction
Astrophysical observations of neutron stars with masses up to two solar masses [1,2] imply that the Equation of State (EoS) relating the energy density ε and pressure p of the matter inside the stars should be very stiff [3]. The stiffness can be measured by the where v s can be identified as the speed of propagation of sound waves, naturally obeying the causal bound v s ≤ 1. According to our current understanding, the nature of this matter 1 The symbol s denotes the entropy density here.

JHEP11(2017)031
ranges from a relatively dilute gas of nuclei immersed in a sea of electrons in the crust of the star to dense nuclear and superdense neutron matter deep inside the star, expected to reach at least a few times the nuclear saturation density, n s ≈ 0.16/fm 3 , in the cores of the most massive stars. With the deconfinement transition of Quantum Chromodynamics (QCD) expected to take place around these densities, it is at the moment still unclear, whether quark matter should be present inside the stars or not.
There are a variety of nuclear matter EoSs that predict very high speeds of sound, some of them even exceeding the speed of light [4]. In all of these cases, the region of validity of the approach is, however, restricted to densities below (roughly) the nuclear saturation density, so that a straightforward extrapolation of the results to the large densities met in the cores of neutron stars is likely to suffer from uncontrollable systematic uncertainties (see [5] for a discussion of this topic). In particular, there is no hope of extending the description of these nuclear matter models to the quark matter phase, possibly relevant for the description of the stellar cores. At the same time, it is equally clear that approaches based on weak coupling expansions in the quark matter phase, such as perturbative QCD [6][7][8][9], cannot be used to describe the transition region, and therefore the standard approaches for the description of this regime typically include model calculations (see e.g. [10] and references therein) and interpolations between the low-and high-density regimes [11].
Considering the above difficulties, there is clearly room for alternative approaches to describing dense strongly interacting nuclear and quark matter. Such a novel approach could be provided by the gauge/gravity, or holographic, duality [12][13][14], which offers a way to relate problems in strongly coupled field theories in their large-N c limit to calculations performed in classical supergravity in a curved spacetime. An interesting observation pointing towards neutron star matter indeed behaving like a strongly coupled system can be seen from the so-called Taub inequality [15] (see also [16]), 2 which states that in a relativistic kinetic theory causality imposes the condition where ρ stands for the mass density. For instance, it is easy to check that degenerate fermionic matter satisfies Taub's inequality for any value of the chemical potential. The inequality clearly implies that ε ≥ 3p, which is saturated by conformal theories. As shown in [3], such an EoS is, however, too soft to support the heaviest observed stars, which clearly implies that one of the assumptions behind Taub's inequality must fail. The most likely culprit is the assumption of the validity of a quasiparticle description, which is far from being guaranteed for the matter found inside neutron stars. In fact, it may well be that the correct expansion point would be that of infinite (or very strong) coupling instead of a system of weakly coupled quasiparticles.
The holographic approach has already been used to describe both the confined [17][18][19][20][21][22][23][24] and deconfined [25][26][27][28] phases of QCD matter through the study of strongly coupled non-Abelian gauge field theories containing fundamental matter with a global U(1) baryon symmetry. In [27], we adopted the strategy of describing the low-density phase of QCD 2 We thank Luciano Rezzolla for drawing our attention to this inequality. matter using the Chiral Effective Theory (CET) results of [29], supplemented by the extrapolations provided in [5], and matching them with the EoS of N = 2 Super Yang-Mills theory at finite baryon density, corresponding to a D3-D7 brane intersection on the gravity side. While successful in providing a consistent description of dense QCD matter, this setup led to the prediction that the deconfinement transition would always be of such a strong first order type that the resulting hybrid stars become unstable as soon as even a microscopic amount of quark matter is generated in their cores. The reason for this behavior was found to be the soft nature of the holographic EoS, with v 2 s < 1/3, in comparison with the stiff low-density EoSs of [5].
The softness of the holographic EoS constructed in [27] came as no surprise; in fact, already in [30,31] it was conjectured that any field theory with a gauge/gravity dual can have a speed of sound at most as large as that of a conformal theory, i.e. v s ≤ 1/ √ 3. In [28], we, however, showed that this conjecture is generically not valid at finite density (even though it might hold in certain theories [31]), and more recently a violation of the bound has been proposed even at zero density through the introduction of multitrace deformations in the dual gauge theory [32]. However, in both cases the violation is not nearly large enough to allow for the existence of quark matter inside neutron stars, and the question remains, whether at least a moderate softness of the EoS of strongly coupled deconfined matter is a universal prediction of holography. We should also note that a bound on the speed of sound at fixed chemical potential has been proposed in [33], and it seems to hold in holographic models that reproduce thermodynamic properties of QCD computed using lattice techniques at small densities [34].
In the present work, we shall demonstrate that the speeds of sound obtained in gauge/gravity models can be arbitrarily close to the speed of light by considering several examples where this turns out to be the case. On the gravity side, the models consist of Einstein-Maxwell theory minimally coupled to a scalar field, which can be either charged or neutral. These models are dual to a strongly coupled gauge theory in its large-N c limit. The bulk gauge field is then dual to a global U(1) current on the field theory side, while the scalar field is dual to a relevant scalar operator. A relevant deformation breaking conformal invariance is introduced by turning on a coupling for the scalar operator. The first example we will study has a string theory (top-down) realization with a known field theory dual, while the rest of the cases considered form a family of bottom-up models. Interestingly, we observe that the simplest scenario including a quadratic potential for a canonically normalized scalar field does not lead to large enough values for the speed of sound. To reach higher values, it is necessary for the scalar field to possess self-interactions, which will be reflected in the properties of higher order correlators of the dual operator. This point should be a very interesting one to investigate further in the future.
Our paper is organized as follows. In section 2 we introduce both the top-down and bottom-up models we work with, and in section 3 we discuss a subtle issue related to the spontaneous generation of a scale in the top-down model. After this, we move on to presenting our main result, the EoS in both types of models, in section 4, which is followed by a thorough analysis of the stability of our solutions in section 5. Conclusions are finally drawn in section 6, while a number of computational details will be discussed in the appendices of the paper. We will use holographic models as a tool to study the EoS of strongly coupled gauge theories at finite density and temperature, although we will be more interested in low temperatures. The models will be chosen in such a way that the theory is well defined in the UV, in the sense that there is a fixed point at asymptotically large energies. If the theory was conformal, the EoS would be fixed by symmetry; here, this will be avoided by introducing a relevant deformation of the UV fixed point that breaks conformal invariance explicitly. We will consider two cases in parallel: a top-down model with a well defined string theory construction, and a family of phenomenological bottom-up models that allow a wider analysis while keeping the main ingredients of the top-down model.

Top-down model
The first case we are going to consider is a deformation of N = 4 SU(N c ) super Yang-Mills (SYM). The theory has a global SU(4) R SO(6) R R-symmetry group associated to rotations of the supercharges. N = 4 SYM contains vector bosons, fermions, and scalars, all in the adjoint representation of the SU(N c ) gauge group. They can be listed as There are three mutually commuting U(1) i=1,2,3 ⊂ SU(4) R in the R-symmetry group. We will study states with charge for the diagonal U(1) (equal charges for all of the U(1) i ). Since N = 4 SYM is a conformal field theory, we will also need to turn on additional couplings that break explicitly conformal invariance. We will do this by introducing a mass for the gauginos, i.e. we will add a term to the Lagrangian of the form As we are not adding similar mass terms for the scalars, this also breaks supersymmetry explicitly.
In the N c → ∞ limit and for very strong 't Hooft coupling λ Y M 1, the N = 4 SYM theory has a holographic dual description as type IIB string theory in a AdS 5 × S 5 geometry, at weak string coupling g s ∼ 1/N c and large curvature radius compared to the string scale L 4 /(α ) 2 ∼ λ Y M . The leading order behavior of the theory is thus captured by classical supergravity (SUGRA) in AdS 5 × S 5 [12]. Turning on a charge density and/or additional couplings in N = 4 SYM is realized in the holographic dual by turning on dual fields that modify the background geometry.
Rather than dealing with the full ten-dimensional SUGRA description of the theory, we will restrict to a subsector that admits a consistent truncation to a simpler five-dimensional theory. The truncation is explained in more detail in [28]. The action reduces to the one of Einstein-Maxwell theory coupled to two real scalars

JHEP11(2017)031
where e is the volume density and with the coupling constant g related to the AdS radius as g = 2/L. The bulk gauge field A µ is dual to the diagonal U(1) R-current J µ and sources for the current.
We introduce the complex field, in such a way that the action takes the form with a charge q = 2 and kinetic and potential terms For small Φ, Therefore, the canonically normalized scalar has a mass m 2 L 2 = −3 which corresponds to a field dual to an operator of conformal dimension ∆ = 3, the gaugino mass operator O = tr λλ, and the associated coupling m 0 . Therefore, the five-dimensional action of the truncated SUGRA subsector contains all the necessary ingredients for our analysis.

Bottom-up models
Taking the top-down model as a guide, we are going to consider a family of models with a gravity dual consisting of Einstein-Maxwell theory minimally coupled to a scalar. Thereby, we will be describing a subsector of the dual field theory including a global U(1) current J µ and a relevant scalar operator O. The usual large-N c and strong coupling limits are assumed to hold for the classical gravity approximation we take to be valid. In order to obtain different EoSs, we will allow for some freedom in the choice of the action for the scalar field. This means that in most cases the field theory dual, if it exists, is not known. We will use these models as an exploratory mean to determine whether holographic models can produce a stiff EoS, with the perspective of looking for proper holographic duals with similar properties in the future. One can in principle allow the kinetic term and the potential for the scalar to be generic functionals, although we will fix their form to be able to do explicit calculations. The five-dimensional action for these models will be as given in (2.5). For the bottom-up models we will take the charge to be zero q = 0, as eventually we would like to identify the U(1) symmetry with baryon symmetry, which is unbroken. For simplicity, we will fix the kinetic term to be canonically normalized K(Φ) = 1 and the potential to be of the form

JHEP11(2017)031
Here, we will allow the masses to lie in the interval 0 > m 2 L 2 ≥ −3, in such a way the scalar field will be dual to a scalar operator of dimension in the interval 4 > ∆ ≥ 3. We will study first the case with a purely quadratic potential V 4 = 0 and then the behavior when V 4 is changed.

Charged black hole solutions in the top-down model
We take an Ansatz for the metric of the form in such a way that it is asymptotically AdS at r → ∞. There is a black hole horizon at r = r H , where f (r H ) = 0. The scalar field and the time component of the gauge field are also turned on and depend only on the radial coordinate, i.e. Φ 0 = Φ 0 (r), A 0 = A 0 (r). The equations of motion and the near boundary behavior of the bulk fields are detailed in appendix A. If the fields were decoupled, their expansion at the boundary would take the form 3 (2.10) We can identify µ with the chemical potential in the dual field theory and φ (0,1) with the coupling of the dual operator. If it is nonzero, this amounts to introducing a relevant deformation that breaks explicitly conformal invariance in the dual field theory. In this case, φ (0,1) gives a mass to the gauginos. For ∆ = 3, φ (0,1) has dimension one, so we can in fact identify it with a mass scale m 0 ≡ φ (0,1) . The coefficients A 0(0,2) , f (0,4) , and φ (0,3) determine the R-charge density, energy, and the expectation value of the scalar operator (gaugino bilinear), respectively. The coefficient of the logarithmic term φ (1,3) is finally proportional to m 0 .
For convenience when obtaining the numerical solutions, we will perform the variable and gauge field redefinitions so that the AdS boundary is now at u → 0 while the horizon is at u → 1. This change implies that in order to correctly match with the holographic renormalization scheme adopted, carried out in the r coordinate, one must perform the shift 3 We use a notation where X (n,m) is the coefficient of (L 2 /r) m (log(r/L)) n in the expansion of the field X.
in such a way that the dependence on r H in the near-boundary series solution is absorbed. Such a shift has to be applied also to the boundary operators (B.18). In the u coordinate, the near-boundary fields read The map between the coefficients in both coordinates is Near the horizon we will impose regularity of the solution plus vanishing boundary conditions for the warp factor and the gauge field. Then, at leading order, we have where the subleading terms can be found in appendix A.2.
It is convenient to define our thermodynamic variables (µ, T ) in units of the mass m 0 so that a 0 = µ r α. We will also normalize the thermodynamic potentials and expectation values of the charge and scalar operators by the mass and a common factor N = L 3 16πG 5 , such that (2.17) After defining We have computed the solutions by means of the shooting technique, thoroughly explained in appendix A.3. We plot the results as a function of µ r for a fixed temperature t r = 1 in figure 1.

Charged black hole solutions in bottom-up models
In the bottom-up models, we will proceed in a similar manner to the top-down one. We take an Ansatz for the metric of the form given in (2.9), and fix q = 0 for simplicity and because the potential application to the physics of dense nuclear matter requires the U(1) symmetry to be unbroken. The equations of motion and the near boundary behavior of the bulk fields are detailed in the appendices of [28]. For ∆ = 3 the equations and expansions take a similar form as in the top-down model (2.10). For ∆ = 3 only the expansion of the scalar field at the boundary changes to We can identifyφ (0,0) with the coupling of the dual operator. If it is nonzero, this amounts to introducing a relevant deformation that breaks explicitly conformal invariance in the dual field theory. Similarly to the top-down model, we will introduce the mass scale For convenience when obtaining the numerical solutions, we will perform the change to the u coordinate (2.11). The near-boundary expansions of the fields are given by eq. (2.13), except for the scalar field, which now reads The map between the coefficients in both coordinates is given by (2.14), except for the scalar, which now takes the form

JHEP11(2017)031
Near the horizon we will impose regularity of the solution plus vanishing boundary conditions for the warp factor and the gauge field as in (2.15). It is convenient to define our thermodynamic variables (µ, T ) in units of the mass m 0 , as in (2.16) and (2.17). The normalization of the expectation value of the scalar operator reads in the general case (2.23) Taking the renormalized expectation values detailed in the appendix of [28], we get for ∆ = 3: (2.24) We will fix κ 2 = 0 in the following, since this parameter is irrelevant for the speed of sound. This selects V 4 = −2/3 as a special value, for which the logarithmic terms drop and the conformal anomaly vanishes, although there are still terms contributing to the trace of the energy-momentum tensor proportional to the expectation value of the scalar operator. For ∆ = 3, we on the other hand get We have computed the solutions using the same numerical methods as for the top-down model. The results are plotted as functions of µ r for a fixed temperature t r = 0.1 in figure 2.

Generation of a new scale in the top-down model
There are some subtleties entering the EoS of the top-down model that we shall presently discuss. In (2.19), κ 1 and κ 2 are the coefficients of finite counterterms. These terms are scheme dependent but once the renormalization scheme has been fixed, their values are related to physical quantities such as the expectation value of the scalar operator and the charge density. This implies that the theory is not completely determined by the bulk action of the gravity dual, but it is necessary to specify the value of the finite counterterms as well. From the point of view of the field theory, consider that in addition to the N = 4 SYM fields there is a decoupled scalar field ϕ and a Yukawa coupling Y ϕ between the scalar and the N = 4 SYM gauginos, In the large-N c limit, we can treat the scalar field as quenched, neglecting loop effects from the N = 4 SYM theory. Nevertheless, this coupling breaks conformal invariance (even though it is classically marginal) and will introduce a logarithmic dependence log(E/Λ) on the energy scale E in physical observables, such as scattering cross sections. In particular, a wave function renormalization of ϕ will show up in the kinetic term of the scalar field, having the same form as the finite counterterm associated to κ 1 . The scale Λ that appears inside the log depends on the scheme, but can be fixed by measurement. After this, the value of Λ will be different in different schemes, but physical quantities will naturally have the same values in each of them. On top of the scale appearing due to logarithmic terms, if the scalar field acquires an expectation value ϕ = m 0 , this will affect the N = 4 SYM theory as an explicit breaking of conformal invariance. Note that in principle the scale of explicit breaking m 0 and the scale that determines the running of the coupling Λ would be completely independent, if no further condition is imposed.

JHEP11(2017)031
To illustrate the above with an example, consider the computation of a one-loop contribution to the self-energy of a scalar field due to a loop of a fermion field of mass m 0 . There is a logarithmic UV divergence that in dimensional regularization in d = 4 − 2 dimensions becomes a pole as → 0. Depending on the scheme, removing this divergence leaves behind different finite terms, taking the forms Here, β ∼ Y 2 ϕ is a scheme-independent factor, M S and M S denote the usual (modified) minimal subtraction schemes with scale parameters Λ andΛ, and F S stands for a fixed scale scheme with an arbitrary finite term κ F S . The physical mass of the scalar M corresponds to the position of the pole in the propagator where m 0 is the bare mass. This can be viewed as fixing the arbitrary renormalization scales of the M S and M S schemes and the constant in the F S scheme, For a given scheme, changing the renormalization scale or the finite counterterm amounts to a change of the physical scale and thus a modification of the theory.
In the holographic calculation we fix the scheme of holographic renormalization by using L as the reference scale in the asymptotic expansion of the fields and m 0 in the definition of the finite counterterms. We could have chosen a different scale, say L , in such a way that Physical results would be unchanged as long as we appropriately identify the values of the finite counterterms in each scheme, We could also have changed the scheme by using a scale different from m 0 in the logs leading to a somewhat different relation between the finite counterterms in different schemes,

JHEP11(2017)031
This shows that an arbitrary scale can indeed be introduced through holographic renormalization.
Once we have fixed our renormalization scheme (for instance one could choose schemes where κ 1 = 0 or κ 2 = 0), different values of finite counterterms correspond to different values of physical quantities (i.e. renormalization group invariants). However, one can see that the effect of κ 2 is to add a term independent of the temperature or the chemical potential that shifts the value of the vacuum energy. It is therefore unimportant for thermodynamics, and a valid physical choice could be that the effective cosmological constant term in the dual field theory vanishes. A similar term appears in the D3/D7 model [35], where the counterterm is fixed by supersymmetry and gives a vanishing expectation value for the scalar operator [36]. In principle, the dual field theory is supersymmetric with soft-breaking terms (the gaugino mass in eq. (2.1)), and supersymmetry could then fix the value of the finite counterterm κ 2 at zero temperature and chemical potential. The value of κ 1 might also be fixed by similar considerations, but in the present context with a nonzero chemical potential, it is not known how or whether the finite counterterm is to be fixed. To our knowledge, this is still an open problem in holographic renormalization in general. This implies that giving arbitrary values to the finite counterterms might spoil the identification of the dual field theory corresponding to the consistent truncation. Nevertheless one could interpret models with different values as extensions of the original supersymmetric model with additional terms (similar to the one in eq. (3.1)) that introduce an explicit breaking of supersymmetry.
Compared to κ 2 , κ 1 has a more interesting and physical effect: it changes the argument of logarithms of α according to This means that a new scale Λ κ has been spontaneously generated in the dual field theory, and that its relative size in comparison with the scale of the explicit breaking of conformal invariance is controlled by κ 1 : In particular, for |κ 1 | sufficiently large, Λ κ can be pushed towards the UV. In figure 3 we plot Λ κ as a function of the reduced chemical potential for various negative values of κ 1 . When the red line crosses the other curves, Λ κ = µ r and the argument of the logarithm in eq. (3.9) becomes unity. Note that the hierarchy is not parametrically large with N or the 't Hooft coupling, so we will remain in the realm of classical supergravity. Furthermore, as the hierarchy is introduced through finite terms in the boundary action, there is no change in the bulk supergravity equations of motion and solutions; in particular the consistent truncation is unaffected. It is still the same subset of operators in the dual field theory that close under the OPEs.

Equation of State
If a weakly coupled quasiparticle description is possible for the system under study, it is appropriate to use kinetic theory to derive its Equation of State. In a relativistic theory causality then imposes a constraint, Taub's inequality [15] where ρ is the mass density. One can check for instance that for a degenerate (noninteracting) Fermi liquid τ = τ F ≥ 1, where τ F = 9 16 (µ 2 r −1) 3 2µ 6 r −3µ 4 r +µ 2 r +log 2 µ r + µ 2 r −1 −2 µ 2 r −1µ 3 r log µ r + µ 2 r −1 (4.2) and µ r = µ/m F where m F is the mass of the fermions.
In a strongly coupled theory the above condition may easily be violated. A simple example is the D3/D7 model [35] that is used to model flavor physics at strong coupling, and that contains quarks and squarks with a mass m q . The EoS is known analytically [37][38][39][40][41], and the pressure, energy density, and mass density at zero temperature read as functions of the chemical potential where λ is an unimportant constant factor. Defining the reduced chemical potential as µ r = µ/m q , one finds Taub's inequality is obviously violated, which indicates that the theory is indeed strongly coupled and that it possesses no good quasiparticle description. Note that there is, however, a (weaker) bound that constrains the Equation of State. Indeed, as long as τ ≥ 0 we will have a condition ε ≥ 3p .

Top-down model
We can compare the values of τ in the models we study here with those of the degenerate Fermi liquid and the D3/D7 model, see figure 4. We observe that τ < 0 for a range of values of the chemical potential (µ r 32 for κ 1 = −10; for even more negative values of κ 1 the curve of the top-down model goes further down). It thus appears that in these regions the EoS is stiffer than in a conformal theory, but how stiff can it be? In order to answer this question we would need to compute the adiabatic speed of sound (1.1). However, it is technically easier to work at fixed temperature and compute the isothermal speed of sound If the pressure has an analytic expansion in T /µ for T /µ 1 (as it will be the case in our models) and the entropy goes to zero at zero temperature, one can neglect the terms proportional to ∂sr ∂µr tr and the two speeds become the same. At non-zero temperature, the difference is suppressed by a factor of at least O(T /µ). Moreover, in many practical applications the temperature is taken to be zero as a good approximation. Therefore, we will study the isothermal speed of sound in the following and drop the label.
The behavior of the speed of sound v s in the top-down model is depicted in figure 5. 4 We observe that, for κ 1 < 0, when |κ 1 | is increased the speed of sound becomes larger at low values of the chemical potential, eventually becoming quite close to the speed of light, and the region where the speed of sound is large also grows. A possible way to understand this is to recall that the scale Λ κ defined in (3.10) that controls the contribution of the logarithmic terms in (2.19) increases with increasing |κ 1 |. When this happens, the logarithmic terms become large in magnitude. If the logarithmic terms in (2.19) dominate, the EoS becomes stiff but remains compatible with causality, as ε r ∼ p r . Therefore, there is no fundamental obstacle towards obtaining a stiff EoS for a large interval of chemical potentials, as long as a significant separation of scales is present.
An important issue to consider is the possibility that the theory might become unstable in the stiff regime. A necessary but not sufficient condition for thermodynamic stability is that the charge susceptibility be positive, In figure 6 we plot χ| tr=1 for different values of κ 1 . For large enough values of |κ 1 |, the susceptibility is positive and the theory is thermodynamically stable with respect to density fluctuations. There is a critical value κ 1 = κ c ≈ −6.44, for which the theory becomes 4 These plots correspond to values of the chemical potential that are much below the regime of validity of the probe approximation used in [28]. For tr = 1 one should go at least to values µr > 150 before we reach the probe limit. unstable at low values of the chemical potential. Therefore, the models with a large speed of sound are thermodynamically stable in the stiff regime. We will study their dynamical stability in section 5.

Bottom-up models
Moving again to the bottom-up models, we first consider the case without a quartic term in the potential V 4 = 0. The results are summarized in figure 7. We find that the speed of sound can be larger than the one in a conformal theory, and that larger deviations occur for operators of lower dimensions, close to ∆ = 3 for our allowed range. The left plot of figure 7 reflects this: there, we have fixed the temperature, computed the speed of sound as function of the chemical potential, and plotted the largest value we have found for each dimension of the scalar operator. This behavior holds for a range of low temperatures. The right plot of figure 7 shows the largest value of the speed of sound for a fixed dimension ∆ ∼ 3 as we vary the temperature. We see that the magnitude increases as we lower the temperature, but it seems to saturate at an absolute maximum. The maximum value is just slightly larger than the conformal value by some 3 %, while for phenomenological purposes it should be at least ca. 30 % larger. Next, we turn on the quartic term in the potential, i.e. let V 4 = 0. In figure 8 we plot the speed of sound as a function of the chemical potential for a fixed temperature t r = 0.1 and different values of V 4 . We observe that making V 4 more negative increases the value of the speed of sound, while making V 4 more positive has the opposite effect. It is possible to reach values of the speed of sound 20-40 % larger than the conformal value for V 4 ∼ −1.5 and µ r ∼ 0.6-0.75. The speed of sound seems to be growing further at lower values of the chemical potential. This shows that stiff phases are possible in generic holographic models.
However, in contrast to the top-down model, we find that in most cases there are violations of causality (v s > 1) or thermodynamic instabilities (v 2 s < 0) at small values of  the chemical potential, so there is likely a phase transition between the high temperature, zero density phase and the low temperature, non-zero density one. Nevertheless, as we show in figure 9, for any given temperature there is a range of values of V 4 where the speed of sound remains in the physical range 1 ≥ v 2 s ≥ 0. This happens around the special value of V 4 = −2/3, for which the conformal anomaly vanishes; we even observe that near the special value the speed of sound becomes very close to its conformal limit and almost independent of the chemical potential.  To inspect thermodynamic stability, we have plotted the charge susceptibility in figure 10. In the cases where v 2 s < 0 we also find that the susceptibility becomes negative, although this appears to happen at lower values of the chemical potential, so it may correspond to a different kind of instability. When the speed of sound becomes superluminal there can also be a small interval with χ > 0 for V 4 −1. In the window where 1 ≥ v 2 s > 0 for all values of the chemical potential we find that χ > 0, so these correspond to thermodynamically stable phases. Numerically, it seems that the values of V 4 for which v 2 s = 0 and χ = 0 at zero chemical potential coincide.

Stability
The aim of this section is to determine whether the stiff phases we have found are indeed local minima of the free energy in the space of homogeneous configurations. To this end, we will introduce a small time-dependent perturbation, expecting that if the equilibrium configuration is unstable we will witness the exponential growth of some of the modes. Otherwise, we expect the perturbation to oscillate and/or decay back to equilibrium. As we are considering only homogeneous configurations, we can suppress the spatial dependence. On the gravity side, this translates into studying a linear perturbation around the previously obtained background solution 5 Φ → Φ(r) + δΦ(r, t), A 0 → A 0 (r) + δA 0 (r, t), g µν → g µν (r) + δg µν (r, t) .
Since our background is stationary, we can expand in plane waves of a given frequency ω, Dynamical modes are normalizable and satisfy an ingoing boundary condition at the horizon. This is possible typically only for a discrete set of complex frequencies, the quasinormal frequencies ω n . If the imaginary part of the quasinormal frequency is negative or zero, Im ω n ≤ 0, the associated quasinormal mode decays in time or is oscillatory, and the background is stable. On the other hand, if Im ω n > 0 is positive, the quasinormal mode grows exponentially in time and the background is unstable. 5 We work in the δgrµ = Ar = 0 gauge. An important piece of information is that at zero chemical potential and high temperatures -µ r = 0, t r 1 -the model is known to be stable, as the quasinormal modes should approximate those of a probe scalar in an AdS black hole background, all of which are on the lower half of the complex frequency plane [42]. As the background changes continuously, a quasinormal mode has to cross the real axis to the upper half plane in order to develop an instability. Physically we expect that if the background becomes unstable there will be another stationary solution corresponding to the true vacuum of the theory.
In that case the crossing should happen at the origin of the complex plane. Therefore, the onset of the instability can be determined from the appearance of a quasinormal mode at zero frequency.
We study the appearance of a zero frequency quasinormal mode using the determinant method of e.g. [43,44] and standard techniques, most details of which can be found in appendix C. First, we introduce gauge invariant combinations of the fields under diffeomorphisms that preserve the condition g µr = 0. There are two independent scalar modes where h = δ ij h ij /3 is the trace of the spatial components of the metric fluctuation. If q = 0, the two modes are coupled with coefficients A ij , B ij that depend on the background fields. If q = 0, the off-diagonal components of A and B are zero and the two modes decouple. We impose that the solutions are ingoing at the horizon. There are two independent solutions z the expansion of the solutions at the boundary u → 0 is, to leading order, where we identify the coefficients of the non-normalizable (nn) and normalizable (n) solutions. We arrange the solutions in a matrix with constant entries at the boundary M depends on the frequency, and a normalizable solution exists when M (ω) has a zero eigenvalue, i.e. det(M (ω)) = 0. For the top-down model we have computed the determinant at zero frequency ω = 0 first for a zero chemical potential µ r = 0 starting at high temperatures and decreasing the temperature to values t r < 1 (left plot in figure 11). As the determinant never vanishes, the background is stable for µ r = 0, t r = 1. We then repeat the same calculation but keeping t r = 1 fixed and increasing the chemical potential µ r . We find that the determinant is nonvanishing in the range we are interested 0 ≤ µ r ≤ 50 (right plot in figure 11). Therefore, the theory remains dynamically stable in the regime where the EoS is stiff.
For the bottom-up model we do a similar stability analysis, we first compute the determinant at zero frequency ω = 0 at zero chemical potential µ r = 0 starting at high temperatures and decreasing the temperature to values t r < 0.1 (left plot in figure 12). We then fix the temperature to t r = 0.1 and increase the chemical potential µ r (right plot in figure 12). We find that the determinant is non-vanishing for values of V 4 where the speed of sound remains in the physical window, even when the speed of sound is close to the speed of light. Therefore, these models have sensible physical behavior and no obvious instabilities even in the regime where the EoS is stiff.

Conclusions
In the paper at hand, we studied the thermodynamics of cold and dense strongly coupled matter via simple holographic models. The models include the minimal ingredients of finite -20 -

JHEP11(2017)031
charge density and breaking of conformal invariance through a coupling for a relevant scalar operator of conformal dimension 4 > ∆ ≥ 3. We find that for some of these cases it is possible to find very stiff Equations of State, with the speed of sound almost reaching the speed of light. A simple stability analysis of the models furthermore showed no obvious thermodynamic or dynamic instabilities.
We observe that the simplest models possessing a quadratic action for the scalar field do not reach speeds of sound significantly larger than the conformal limit of v s = 1/ √ 3. In bottom-up models with a quartic potential, the speed of sound can on the other hand reach the speed of light if the quartic term has a negative coefficient V 4 < 0 with large enough magnitude. However, except for a small range of values around V 4 = −2/3, the isothermal speed of sound becomes superluminal or imaginary (indicating the presence of an instability) at low values of the chemical potential. Concerning the superluminal behavior, it should, however, be noted that when the chemical potential is of the same order or smaller than the temperature, one should rather consider the adiabatic speed of sound, which may affect to the range of values of V 4 , for which causality is respected.
In addition to the bottom-up models, we also studied a top-down model with a more complicated action for the scalar, determined by a consistent truncation of supergravity. The issues of superluminal or imaginary speeds of sound do not appear in this case, which suggests that adding higher powers of the scalar field to the scalar potential might ameliorate the behavior of these quantities also in the bottom-up models. On the other hand, a stiff EoS is achieved in the top-down model only when there is a large separation between the scale of explicit breaking of conformal invariance and another scale that is spontaneously generated due to logarithmic divergences. The conclusion seems to be that although there is no fundamental obstruction to achieving a stiff EoS, this may not be possible in the simplest models and/or for the most "natural" values of the parameters of the system. On the positive side, the conditions required to achieve a stiff and physically consistent EoS may prove to be quite restrictive and thus turn out to be useful in constraining possible holographic models of QCD.
An interesting question is if the increase in stiffness is due to some underlying physical mechanism involving microscopic degrees of freedom in the dual field theory, at least for the top-down model where the dual is known. However, microscopic fields are not gauge invariant and therefore not directly accessible using the duality. So far we can just make a broad qualitative statement, it appears that one needs a combination of explicit and anomalous breaking of conformal invariance, with a hierarchy between these scales such that the scale of anomalous breaking is the larger.
An obvious phenomenological application of our results lies in the physics of neutron stars, where a holographic quark matter EoS has previously been matched to nuclear matter EoSs in [27]. The fact that very stiff EoSs can be obtained from holography opens up the possibility to construct matched EoSs exhibiting a weakly first order or even a cross-over deconfinement transition, thus allowing for the existence of a macroscopic amount of quark matter in the cores of the stars. Recalling the ease, with which quantities such as neutrino emissitivities and transport coefficients can be computed in holography, this paves the way for very interesting astrophysical studies.

A Background solutions
Varying the action (2.5) of the top-down model (2.6) with respect to the bulk metric, gauge, and scalar fields yields the equations of motion (A.1)

A.1 Near boundary series expansions
The near boundary behavior for the scalar field is We shall assume the following series expansions for the other fields which upon implementing the equations of motion become plus sub-leading terms that we do not put here.

JHEP11(2017)031
tracking tangent lines. For a guess X, compute the Jacobian J of partial derivatives of the mismatch vector. The new vector X shall be The Jacobian is computed through finite differences, once the solutions in a neighborhood of the guess point (on each direction on the constants space) are known. In particular, as step in the Jacobian we will take 10 −10 . On each numerical integration, u hor = 1 − 0 , u boun = 0 , 0 being some sensitive cut-off; we use 10 −8 and u * = 1/2. As for the initial data X guess , a sensitive choice for mild reduced chemical potential and temperature is the solution inherited from the scalar field in probe approximation, wherein the geometry reduced to an AdS-RN [28], where α P , β P are obtained from integration of the scalar equation in this approximation, once φ P H is set. If the norm of the mismatch ||M || lies above some threshold fixed a priori, the iteration starts once again, but taking X as the new starting point and stops if otherwise. In our computations, we will fix the threshold to be 10 −9 . Our attempt to connect the model to neutron star physics implies that we will focus in regimes at which X AdS-RN works not very well, but luckily, thanks to the smoothness of the solutions, if for some choice X (µ 0 ,t 0 ) , ||M || < 10 −9 , then we can take this vector as initial guess on the next computation, i.e., X (µ 0 ,t 0 ) → X guess (µ 0 +δµ,t 0 ) .

B.1 On-shell action
For the holographic models we consider, one can write Einstein's equations in the form From the trace of these equations, we find that the Ricci scalar reads implying that the on-shell action (2.5) evaluates to Let us now use the fact that for our solutions Γ α µν = Γ r rν = Γ α rr = 0 , Γ r µν = − 1 √ g rr K µν , Γ α µr = √ g rr K α µ , Γ r rr = 1 2 g rr ∂ r g rr , (B.4) where K µν = 1 2 √ g rr ∂ r g µν (B.5) -25 -

JHEP11(2017)031
is the extrinsic curvature and K α µ = g αβ K βµ , K = g µν K µν . Using also the simple result ∂ r √ −g √ −g = Γ r rr + √ g rr K , (B.6) we can write Here, we defined γ µν = g µν as the boundary metric and used √ −g = √ g rr √ −γ. On the other hand, from Einstein's equations we obtain where we only focused on the nonzero components of the solutions. Solving now for V Φ and introducing the result in the on-shell action, one gets Finally, we use the equation of motion for the gauge field, We can then replace the q 2 term in the action by a derivative term and write the action as a total derivative: (B.11)

B.2 Holographic renormalization
In order to be able to read off the speed of sound, we need the energy density ε and pressure p, which can be read from the diagonal components of the expectation value of the stress energy tensor, T µν . We can decompose the line element (2.9) into its transverse and longitudinal components, We will now determine, which counterterms we need to consider in order to obtain finite one point correlation functions. Together with the cosmological constant term

JHEP11(2017)031
which will cancel out the volume divergence, we need to include also the Gibbons-Hawking term, 14) The details of the holographic renormalization of bottom-up models can be found in [28].