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 $v_s=1/\sqrt{3}$. 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 thermodynamic derivative 1 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 ranges 1 The symbol s denotes the entropy density here.
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 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 (topdown) 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 Sec. 2 we introduce both the top-down and bottomup models we work with, and in Sec. 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 Sec. 4, which is followed by a thorough analysis of the stability of our solutions in Sec. 5. Conclusions are finally drawn in Sec. 6, while a number of computational details will be discussed in the Appendices of the paper.

Holographic models
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. 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

Top-down model
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 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 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.7). 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 (2.10) 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.12) 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 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 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.
out in the r coordinate, one must perform the shift 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.81). 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πG5 , such that (2.19) After defining and taking the renormalized expectation values (B.81), detailed in Appendix B.2, we then get

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.11), 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.12). 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 m 0 = (φ (0,0) ) 1/(4−∆) .
For convenience when obtaining the numerical solutions, we will perform the change to the u coordinate (2.13). The near-boundary expansions of the fields are given by Eq. (2.15), except for the scalar field, which now reads The map between the coefficients in both coordinates is given by (2.16), except for the scalar, which now takes the formφ 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.17).
It is convenient to define our thermodynamic variables (µ, T ) in units of the mass m 0 , as in (2.18) and (2.19). The normalization of the expectation value of the scalar operator reads in the (2.25) Taking the renormalized expectation values detailed in the Appendix of [28], we get for ∆ = 3: (2.26) 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 energymomentum 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 Fig. 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.21), κ 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.
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, 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].
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 : Λ κ m 0 = a 0 e −κ1/8 .

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 (non-interacting)  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 . 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 Fig. 5. 4 We observe that, for κ 1 < 0, when |κ 1 | is increased the speed of sound becomes larger at low values 4 These plots correspond to values of the chemical potential that are much below the regime of validity of the  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.37) that controls the contribution of the logarithmic terms in (2.21) increases with increasing |κ 1 |. When this happens, the logarithmic terms become large in magnitude. If the logarithmic terms in (2.21) 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 probe approximation used in [28]. For tr = 1 one should go at least to values µr > 150 before we reach the probe limit. 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 Fig. 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 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 Sec. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 8 we plot the speed of sound as a function of the chemical potential for a fixed temperature t r = 0.1 and   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 Fig. 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 Fig. 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) . (5.46) 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.
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 corresponding to making z 1 or z 2 zero at the horizon. When these solutions are taken to the boundary, a linear combination of them will be normalizable for the values of the frequency corresponding to the quasinormal modes. In the u coordinate, the expansion of the solutions at the boundary u → 0 is, to leading order, 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 Fig. 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 non-vanishing in the range we are interested 0 ≤ µ r ≤ 50 (right plot in Fig. 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 Fig. 12). We then fix the temperature to t r = 0.1 and increase the chemical potential µ r (right plot in Fig. 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 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/

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

A.1 Near boundary series expansions
The near boundary behavior for the scalar field is plus sub-leading terms that we do not put here.
Note that f (1) H can be fixed in terms of t r and A 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,t0) , ||M || < 10 −9 , then we can take this vector as initial guess on the next computation, i.e., X (µ0,t0) → X guess (µ0+δµ,t0) .

B Calculation of thermodynamic quantities 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.7) evaluates to Let us now use the fact that for our solutions is the extrinsic curvature and K α µ = g αβ K βµ , K = g µν K µν . Using also the simple result we can write Here, we defined γ µν = g µν as the boundary metric and used 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.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.11) 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 which will cancel out the volume divergence, we need to include also the Gibbons-Hawking term, The details of the holographic renormalization of bottom-up models can be found in [28]. In the following we focus on the top-down model, that present some small differences due to the more complicated form of the kinetic term and the potential for the scalar field.
From the near boundary behavior of the metric field, we note that it is necessary to add the following counterterm that will cancel out divergences due to the backreaction of the scalar field, Another counterterm may also be added, which will introduce non-trivial finite contributions to our QFT.

C.1 Equations for gauge invariant combinations
We will use radial gauge δg µr = δa r = 0. At zero spatial momentum fluctuations split in decoupled sectors according to their representation under the group of spatial rotations. There are three sectors: • Tensor: h ij − 1 3 δ ij δ kl h kl .
In principle we expect instabilities to be related to changes in the scalar, thus we will restrict the analysis to the scalar sector. We see that there are five components of the fields in the scalar sector. The equations of motion (Einstein, Maxwell, and the equation of motion for the scalar) include a second order (dynamical) equation for each mode plus three first order (constraints) equations. This adds up to eight coupled equations for the five modes. However, the actual number of independent dynamical modes is just two and the system can be reduced to two coupled differential equations (of second order). We will do this in the following.
In the radial gauge there are residual diffeomorphisms ξ M (x) and gauge transformations λ(x). The linear variations of the fields are where c 0 , c i are arbitrary functions of the frequency. We can construct a basis of two independent combinations of the scalar components that are invariant under these gauge transformations z 1 , z 2 ; these are the expressions given in (5.48). The equations of motion can be found in a straightforward way by taking radial derivatives of z i and using the equations of motion of the scalar modes. They take the generic form (5.49). The result with q = 0 is quite cumbersome, so we will give here expressions for the bottom-up models with q = 0 and canonical kinetic term, but generic potential. The off-diagonal coefficients vanish A 12 = A 21 = B 12 = B 21 = 0 and the diagonal ones take the values: (C.94)