Towards a Holographic Quark Matter Crystal

We construct the gravity dual of $d=4$, $\mathcal{N}=4$, SU($N_\rm{c}$) super Yang-Mills theory, coupled to $N_\rm{f}$ flavors of dynamical quarks, at non-zero temperature $T$ and non-zero quark density $N_\rm{q}$. The supergravity solutions possess a regular horizon if $T>0$ and include the backreaction of $N_\rm{c}$ color D3-branes and $N_\rm{f}$ flavor D7-branes with $N_\rm{q}$ units of electric flux on their worldvolume. At zero temperature the solutions interpolate between a Landau pole singularity in the ultraviolet and a Lifshitz geometry in the infrared. At high temperature the thermodynamics is directly sensitive to the Landau pole, whereas at low temperature it is not, as expected from effective field theory. At low temperature and sufficiently high charge density we find thermodynamic and dynamic instabilities towards the spontaneous breaking of translation invariance.


Introduction and summary
Quantum Chromodynamics (QCD) at non-zero quark density is notoriously difficult to analyze. The only first-principle, non-perturbative tool, namely lattice QCD, is of limited applicability due to the so-called sign problem [1]. It is therefore useful to construct toy models of QCD in which interesting questions can be posed and answered. The gauge/string duality, or holography for short [2,3,4], provides a framework in which the construction of such models is possible. The goal is not to do precision physics but to be able to perform first-principle calculations that may lead to interesting insights applicable to QCD (see e.g. [5] for a discussion of the potential and the limitations of this approach). In the case of QCD at non-zero temperature, the insights obtained through this program range from static properties to far-from-equilibrium dynamics of strongly coupled plasmas (see e.g. [6] and references therein).
The simplest four-dimensional example of holographic duality can be obtained from the supergravity solution for a collection of N c D3-branes. In this setup one finds an equivalence between four-dimensional, N = 4 super Yang-Mills (SYM) theory with gauge group SU(N c ) living on the stack of D3-branes and string theory on the near-horizon limit of the geometry sourced by the D3-branes [2]. Since the matter in these theories is in the adjoint representation of the gauge group, in order to consider a non-zero quark density new degrees of freedom in the fundamental representation must be included. On the gravity side this can be done by adding N f so-called flavor D7-branes to the D3-brane geometry [7]. For conciseness we will refer to these new degrees of freedom as 'quarks' despite the fact that they include both bosons and fermions. Placing the theory at a non-zero quark density N q then corresponds to turning on N q units of electric flux on the flavor branes [8]. This flux sources the same supergravity fields as a density N q of fundamental strings dissolved inside the flavor branes. We will thus refer to N q as the quark density, as the electric flux on the branes, or as the string density interchangeably. Note that N c and N f are dimensionless integer numbers, whereas N q is a continuous variable with dimensions of (energy) 3 .
If N f N c and N q is sufficiently small then there exists an energy range in which one can study this system in the so-called 'probe approximation' (see [9,10] for early references and [11,6] for reviews), meaning that the gravitational backreaction of the flavor branes and the strings on the original D3-brane geometry can be neglected. However, this approximation inevitably breaks down both at sufficiently high and at sufficiently low energies.
At high energies the probe approximation breaks down because it ignores the positive β-function of the gauge theory. One way to see this is to note that the β-function is proportional to N f /N c and hence the logarithmic running of the coupling is a small correction to the physics over energy ranges that are not exponentially large in N c /N f . However, at sufficiently high energies the coupling constant eventually diverges because the theory develops a Landau pole. In this region the backreaction cannot be ignored and the probe approximation ceases to be valid. Nevertheless, this high-energy regime can be correctly described holographically [12] by including the backreaction of the D7-branes [13].
At low energies the probe approximation breaks down because the backreaction of the charge density always dominates the geometry sufficiently deep in the infrared (IR) no matter how small N q is. Changing the value of N q simply shifts the energy scale at which this happens. Therefore the inclusion of backreaction is not an option but a necessity in order to identify the correct ground state of the theory.
In this paper we will find the fully backreacted supergravity solutions for the D3-D7 system at non-zero temperature T and non-zero quark-density N q . As we will explain in Sec. 2, we will distribute, or smear, the flavor branes [13] (see [14] for a review) and the strings [15] over the compact part of the geometry. We will focus on the case in which this geometry is an S 5 , but our results are also valid (with the sole modification of some numerical coefficients) if this is replaced by any other five-dimensional Sasaki-Einstein (SE) manifold. On the gauge theory side this corresponds to replacing the maximally supersymmetric SU(N c ) gauge theory by an N = 1 quiver theory.
The solutions that will be our main focus correspond to charge densities between zero and a certain critical value. Solutions with charge densities above this critical value display different, somewhat pathological IR physics, and will only be discussed in Appendix F. In the rest of the paper, including the rest of this section, we will only discuss sub-critical solutions. These are pictorially summarised in Fig. 1, in which each solution on the gravity side, or equivalently each Renormalization Group (RG) flow on the gauge theory side, is represented by a curve running from top to bottom. Each curve corresponds to a different value of the  charge density. A point on a curve corresponds to a given energy scale. At zero temperature a solution is described by an entire curve. Cutting off a curve at different points would correspond to different solutions with the same charge density and different non-zero values of the temperature. In the limit T N 1/3 q the solutions were constructed perturbatively in N 1/3 q /T in [16,17]. In this limit, however, all the relevant IR physics is hidden behind the horizon.
The structure of the solutions can be simply understood as follows. In the ultraviolet (UV) all solutions approach the neutral solution of [13], since the geometry is dominated by the backreaction of the D3-and the D7-branes. In other words, in this limit the temperature and the quark density only produce subleading corrections. The asymptotic geometry is singular, as it corresponds to the presence of a Landau pole in the gauge theory. This introduces a physical scale in the theory, Λ LP , that can be compared, for example, with T or N 1/3 q . Despite the presence of the LP singularity, holographic renormalisation can be straightforwardly implemented and physical quantities can be computed unambiguously [12]. This means that, just as in Quantum Electrodynamics, the presence of the Landau pole is no impediment for the extraction of sensible IR physics.
At zero temperature and non-zero charge the non-compact part of all solutions approaches a five-dimensional Lifshitz geometry with dynamical exponent z = 7 in the IR. In this regime the geometry is dominated by the backreaction of the strings and the D3-branes, with the D7-branes producing only a subleading effect. We emphasize that the precise IR theory to which the theory flows depends on the value of the charge density and other parameters. In other words, although in all cases the dynamical exponent is the same, other features are different. A simple example is the number of active degrees of freedom at low temperature, as measured by the entropy density. The Lifshitz symmetry implies that this must scale with temperature as s = f T 3/z = f T 3/7 , (1.1) but the T -independent function f is not fixed by the scaling properties of the IR solution.
We will see in Sec. 5 that this function depends on N q and other parameters. When the charge density vanishes the zero-temperature solution reduces to the supersymmetric solution of [13]. This configuration flows from the Landau pole geometry in the UV to an IR that differs "only logarithmically" from AdS 5 , meaning that observables exhibit conformal invariance up to a logarithmic dependence on the energy scale. This flow is represented by the upper diagonal straight line in Fig. 1. The solutions with non-zero temperature that can be obtained by adding a horizon at different points along this line were constructed in [12].
The logarithmic dependence of the log-AdS solution is a consequence of the running of the coupling constant caused by the presence of the flavor branes. It is possible to perform a scaling of all the charges in such a way that this effect disappears while the backreaction of the strings is retained. In this limit there exists a flow from an AdS 5 geometry in the UV to a Lifshitz geometry in the IR. This flow was found in [15] and is represented by the lower diagonal straight line in Fig. 1.
We will show in Sec. 2.3 that, up to simple rescalings, the full set of zero-temperature solutions can be reduced to a set parameterized only by the dimensionless combination where s is the string length and we have omitted purely numerical factors (see (2.60) for the exact expression). Solutions with non-zero temperature depend on N q and where we have omitted purely numerical factors (see (3.8) for the exact expression). The factors of N c , N f and s in these equations simply reflect a convenient choice of units on the gravity side. In particular, s can be traded for a function of Λ LP , N q , N c , N f and the value of the coupling constant at some reference scale. This shows that N q and T are truly gauge-theoretic parameters that can be thought of as a dimensionless charge density and a dimensionless temperature, respectively. Moreover, the s factors will always cancel in dimensionless gauge theory quantities and, as we will see, the scaling with N c of physical quantities suffers from an ambiguity. We will come back to this point in Sec. 6. We can now understand the general structure of the results. For large values of N q , which can be thought of as a large-charge density limit, the supergravity solutions transition directly from the LP geometry to the Lifshitz geometry. Instead, for small values of N q the solutions exhibit an intermediate region in which they display the physics of the log-AdS region. Adding temperature simply cuts off the flows at different scales, possibly hiding the Lifshitz and/or the log-AdS regions behind the horizon.
The main result of our analysis is the phase diagram of the system. This is schematically summarised in Fig. 2. We will refer to this figure as a "phase diagram" despite the fact that it is a slight abuse of terminology, since some of the solutions are unstable and we have not identified the putative stable solution that would be thermodynamically preferred. Regions  I, II and III are all locally thermodynamically, as well as dynamically, unstable. We will now summarise the properties of these regions. Region I describes the high-temperature behavior of the system. In this region the solutions are locally thermodynamically unstable because the specific heat at constant charge, C Q , and the charge susceptibility, χ, are both negative. The system is afflicted (at least) by two dynamical instabilities associated to a negative speed of sound squared, c 2 s , and to a negative charge diffusion constant, D. The behavior of C Q and c 2 s in this region is the same as in the high-temperature limit of the neutral solution of [12] (in which χ and D are not defined). In the neutral case C Q and c 2 s become negative at exactly the same point because they are related through c 2 s = s/C Q , with s the entropy density. It is remarkable that this feature extends to charged solutions in an almost charge-independent manner. In other words, the top curve in Fig. 2 is almost a straight horizontal line. As explained in [12], the hightemperature behavior follows directly from the properties of the LP. Since our main interest is in IR features that are independent of the UV completion of the theory, we will not discuss Region I further.
Region IV is the only stable region in the phase diagram. Local thermodynamic stability is guaranteed by the fact that both C Q and χ are positive. Our analysis shows no sign of a dynamical instability either, since c 2 s and D are both positive. We also do not see any indication of a global instability since the pressure (free energy) at any point in Region IV is higher (lower) than that of the neutral phase at the same temperature. We emphasize that this does not prove the global thermodynamic stability of Region IV, since there could be other phases that we have not identified that have a lower free energy.
Region III is particularly interesting since it is the low-temperature, high-charge density region. This region is locally thermodynamically unstable since C Q > 0 but χ < 0. It is also dynamically unstable since c 2 s < 0. Interestingly, Region III includes subregions where the charge diffusion constant can have either sign. Roughly speaking, D is positive in all of Region III except in a small subregion at sufficiently low temperature and at a charge density higher than those shown in Fig. 2. As we will see, this is correlated with the negative sign of the enthalpy density, E + P , which itself indicates another type of dynamical instability (see Fig. 13(bottom-left) and Sec. 6). Region III also includes two different subregions from the viewpoint of how the pressure (equal to minus the free energy) compares to the pressure of the neutral phase at the same temperature. These are separated by the dashed curve in Fig. 2.
In the subregion to the left (right) of the curve the pressure of the charged phase is higher (lower) than the pressure of the neutral phase at the same temperature, indicating that the charged phase is thermodynamically (not) preferred with respect to the neutral phase. Thus bubbles of the neutral phase can form inside the charged phase in the region to the right of the dashed curve. This instability against finite charge fluctuations is not particularly relevant in this context since the region where it is present is also unstable against infinitesimal charge fluctuations. Finally, at higher charge densities than those shown in Fig. 2 the quantity E +3P becomes negative. We will come back to this in Sec. 6.
Region II is locally thermodynamically unstable since C Q > 0 but χ < 0. The speed of sound c 2 s is positive. However, in part or in all of Region II the charge diffusion constant is D < 0, indicating a dynamical instability. As we mentioned above, the sign of D is correlated with that of E + P . Our current numerical results suggest that c 2 s is always positive in Region II, thus D would be negative throughout the entire region. However, with our current resolution it is difficult to definitively exclude the possibility that c 2 s becomes negative in some small subregion of Region II. If this happens then D will turn positive. Just like Region III, Region II also includes subregions that are or are not thermodynamically preferred with respect to the neutral phase, as indicated by the dashed curve.
The negative values of c 2 s and/or D that we have encountered imply dynamical instabilities in the hydrodynamic sound and charge diffusion channels, respectively, towards the spontaneous breaking of translation invariance. This suggests that the putative, stable phase in the corresponding regions may be a crystalline phase. We will come back to this point in Secs. 5 and 6.

Gravitational description
We study holographically a family of theories consisting of SYM with N c colors and N f massless flavors at non-zero temperature T and non-zero charge density N q . The string dual to this configuration is given by the geometry sourced by a stack of type IIB N c D3-branes coupled to a set of N f D7-branes with N q units of electric flux on their worldvolume. The temperature is represented by the presence of a black brane horizon. We work in the approximation in which the system can be described by means of type IIB supergravity with sources. The regime of validity of this approximation will be discussed in Sec. 3.5.
Each D7-brane extends along the directions parallel to the D3-branes and along the holographic radial direction, and wraps a three-cycle inside the five-dimensional compact part of the geometry. We smear the D7-branes in these compact directions in the most symmetric way compatible with supersymmetry. In this way we recover some of the isometries that would be broken by a single D7-brane (or by a collection of overlapping D7-branes) and reduce the problem to solving ordinary (as opposed to partial) differential equations.

Ten-dimensional ansatz
Four-dimensional, N = 4 SYM theory with N c colors is holographically dual to supergravity solutions sourced by N c D3-branes. In the supergravity description this is encoded in a flux of the self-dual Ramond-Ramond (RR) five-form through an appropriate five-dimensional compact manifold M 5 : Here ω 5 is the volume form of M 5 , whose total dimensionless volume we denote as V 5 . Quantization requires that the D3-brane charge is related to the number of colors through Note that, unlike in Refs. [18,19,20,21], we follow [12] and work with a RR charge quantized in units of N c instead of g s N c (for a comparison between both normalizations, see e.g. Sec. 4.1 of [22]). The latter choice is convenient in situations in which there is a natural factorization of the dilaton φ of the form e φ = g s eφ. This is the case, for example, if the gauge theory is conformal, since this means that the dilaton is constant and one can simply normalize it so thatφ = 0 everywhere. If the gauge theory is not conformal but approaches a fixed point in the IR or in the UV then it is natural to normalize the dilaton so thatφ = 0 at the corresponding fixed point. In contrast, in the solutions that we will consider the dilaton will run from zero to infinity and there will be no natural factorization into a constant piece and a running piece. We will therefore work with the full dilaton, which is related to the running YM and 't Hooft couplings through In the simplest setups the metric supported by the F 5 flux is AdS 5 × M 5 , with M 5 a Sasaki-Einstein (SE) manifold. The radius L of these two spaces is related to the D3-brane charge through Q c = 4L 4 .

(2.4)
If M 5 = S 5 then the gauge theory is N = 4 SYM; otherwise it is a non-maximally supersymmetric theory. For example, if M 5 = T 1,1 then the gauge theory is the Klebanov-Witten quiver. In general the rank of the gauge group, the radius L and the five-dimensional effective gravitational coupling (see (2.39) below) are related through In order to add flavor to any of the theories above it is convenient to view the SE manifold as a U(1) fibration over a four-dimensional Kähler-Einstein (KE) base. This geometric construction is naturally equipped with an SU(2) structure characterized by a real one-form, η, and a real two-form, J, which is the Kähler form of the KE manifold. These satisfy the relations dη = 2J , dJ = 0 , They also close on each other under Hodge dualization on the SE manifold: We normalize the curvature of the SE manifold so that R ab = 4 g ab . Making use of the fibration structure the most general ten-dimensional metric compatible with the symmetries that we wish to preserve takes the form (in string frame) where t and x are the gauge theory Minkowski directions and r is the holographic radial coordinate. Note that we have allowed for (i) a squashing between the fiber and the base of the SE manifold, i.e. G b = G f , since this will be produced by the backreaction of the flavor branes, and (ii) for the possibility that G tt = G xx , since we will consider solutions with non-zero temperature and/or non-zero charge density that will therefore break Lorentz invariance.
As explained in Sec. 1, the addition of flavor and a quark density on the gauge theory corresponds on the gravity side to the addition of D7-branes with an electric Born-Infeld (BI) field on their worldvolume of the form with B the Neveu-Schwarz (NS) potential and the BI potential. These objects act as sources for both the RR fields and the NS field strength H = dB, and hence they modify their equations of motion. For the RR fields, through the Hodge-duality relations that these fields obey, this also leads to a modification of their Bianchi identities, and therefore to a modification of their very definition in terms of gauge potentials. The full action in the presence of the most general set of D7-brane sources is discussed in Appendix A, to which we refer the reader for additional details. In this work the general equations simplify because F ∧ F = 0. Under these circumstances the Bianchi identities take the form In these equations Q f is the D7-brane charge, related to the number of D7-branes through with V 3 = J ∧ η the dimensionless volume of the three-dimensional submanifold M 3 ⊂ M 5 wrapped by any of the D7-branes. Eq. (2.11c) is satisfied by (2.1). Eq. (2.11a) is the "violation" of the usual Bianchi identity for F 1 that expresses the fact that D7-branes are magnetic sources for F 1 . This violation is already present in the case of flavor without strings discussed in [13] (see [14] for a review) and it is solved by taking For later reference, we note that in our conventions the tension of a D7-brane is (2.14) with κ 2 the ten-dimensional gravitational coupling. The second term on the right-hand side of (2.11b) is a violation of the usual Bianchi identity for F 3 which was not present for the solutions with strings but without D7-branes studied in [18]. Indeed, this term is only present when the system contains both strings and D7-branes [16], since it comes from the magnetic components of C 6 sourced by the term contained in the Wess-Zumino (WZ) part of the D7-branes' action. Through the duality relation these magnetic components give rise to the second term on the right-hand side of dF 3 . We thus see that, although it is the string density represented by F that sources C 6 , the presence of the D7-branes is necessary since the coupling between F and C 6 is supported on their worldvolume. This coupling implies that C 6 must contain a term of the form where B(r) is a function that will be determined below. We thus take the following ansatz for the dual three-form: The second term in this expression is the one implied by (2.17), whereas the first one, as we will see, is related to the string density [23,15]. This first term can also be interpreted as describing a density of baryon-vertex-like D5-branes [24] wrapping the compact manifold M 5 . In writing this ansatz we have already made use of the fact that we will seek solutions with H = 0, i.e. we have discarded the first term on the right-hand side of (2.11b). To see the relation between the first term in (2.18) and the string density we turn to the equation of motion for the B-field: with S the total supergravity-plus-sources action. This may be written as the contribution of the Dirac-Born-Infeld (DBI) part of the D7-brane action. Two observations are important. First, there is no explicit contribution from a variation of the WZ term of the D7-branes because this has combined on the right-hand side of (2.20) with other terms into the gauge-invariant, modified field strengths F n , whose definition is given in (A.11). Second, the left-hand side of (2.19) is not the same as the left-hand side of (2.20), since we have moved terms between both sides of the equation when going from (2.19) to (2.20).
As anticipated above, we will solve Eq. (2.20) with B = H = 0, which means that its right-hand side must vanish. At first sight it may seem surprising that the presence of strings does not automatically lead to a non-zero H, but the non-linearities of supergravity imply that the H sourced by the strings can be exactly cancelled by the H sourced by the products of RR forms in (2.20). Requiring this cancellation fixes the BI field on the D7-branes to To close the circle we note that this is automatically a solution of the equation of motion for the BI field obtained by varying the full action with respect to A. The reason is that dA always appears together with B in the gauge-invariant combination (2.9). This means that the so-called electric displacement, i.e. the momentum conjugate to A, is given by and therefore that the equation of motion for A, is automatically implied by the exterior derivative of (2.19). A straightforward calculation shows that the electric displacement is given by The string density in the 123-directions is obtained by integrating this expression over M 5 , which finally yields the relation between the string density and Q st :

Five-dimensional effective theory
In the previous section we have written a ten-dimensional ansatz that includes all the solutions of interest. The RR forms are given by (2.1), (2.13) and (2.18). The NS B-field vanishes. The BI field is given by Eq. (2.22). The ten-dimensional metric takes the form (2.8). The functions that we must solve for are the metric components, B, and the dilaton φ, all of which depend only on the radial coordinate r. For the analysis of the solution it is convenient to reduce the ten-dimensional system to an effective five-dimensional action. This exercise was done in a general setup in Ref. [25], from where we extract the final result, truncated to the fields of interest for us. We begin by parameterizing the ten-dimensional string-frame metric as ds 2 10 = e φ/2 e − 10 3 σ ds 2 5 + L 2 e 2σ−2w ds 2 KE + L 2 e 2σ+8w η 2 , (2.28) and the effective five-dimensional metric as ds 2 5 = g tt (r) dt 2 + g rr (r) dr 2 + g xx (r) dx 2 3 . (2.29) The RR fluxes F 1 and F 5 thread only internal directions so, upon dimensional reduction, they simply produce terms in the scalar potential. In contrast, in order to dimensionally reduce the RR three-form we must expand it in terms of five-dimensional forms that we denote G n (not to be confused with the ten-dimensional RR field strengths G n defined in Appendix A).
Computing the Hodge dual in (2.18) we see that F 3 has the following components: Imposing now the Bianchi identity (2.11b) we see that G 2 and G 1 can be written in terms of gauge potentials C 0 , C 1 (not to be confused with the ten-dimensional RR potentials C n defined in Appendix A) as We will see below that, in the gauge C 0 = 0, the vector C 1 is directly related to the function B in Eq. (2.18). We find it convenient to dualize the three-form G 3 to a vector field in fivedimensions.
Since G 3 appears in the five-dimensional action coupled topologically to the NS potential B [25], Therefore we define the dual vector A 1 (not to be confused with the BI potential A) and its field strength F 2 (not to be confused with the ten-dimensional RR field strengths F n defined in Appendix A) through the relation In summary, the dimensional reduction of F 3 gives three gauge potentials C 0 , C 1 , A 1 with field strengths G 1 , G 2 , F 2 . In terms of these the final result for the five-dimensional effective action is [25] and We emphasize that this five-dimensional action is a useful way to encode the equations of motion of the theory within the ansatz that is of interest to us, but it is not a consistent truncation. The reason is that only solutions of the five-dimensional action for which can be lifted to a solution of the original ten-dimensional equations. We have kept H in the five-dimensional action despite the condition that it must vanish for the uplift to exist because its equation of motion (2.43c) is non-trivial even after setting H = 0. The five-dimensional Newton's constant in (2.35) is related to the ten-dimensional one through 1 2κ 2 We have written the DBI part of the action in a form that vanishes identically if F = 0, so that the contribution of this part of the D7-branes' action in the neutral case is contained entirely in the term linear in Q f of the scalar potential This potential can be derived from the superpotential The equations of motion for the differential forms are Note that the equation for G 1 follows from the exterior derivative of (2.43b). Eq. (2.43a) is simply the Bianchi identity for G 3 -see (2.34) -and it is solved by In terms of the vector potential this implies where we have used the fact that in our solution B = 0.
To treat the remaining equations it is convenient to work with the momentum conjugate to the massive vector C 1 in the gauge in which C 0 = 0. Taking the ansatz one can define the conjugate momentum as and similarly We have chosen the normalization in such a way that this momentum coincides exactly with the function B(r) appearing in (2.18 This automatically solves the equation of motion for A t that follows from the action (2.35): where the electric displacement is which reproduces the ten-dimensional result (2.25) upon using the relations between κ, κ 5 , L and Q c .

Scalings
The action (2.35) enjoys several scaling symmetries that will allow us to work with appropriate dimensionless variables. To see this let us recall the length dimensions of the different dimensionful variables in our ansatz (in this section we set the temperature to zero): In particular, note that B rt , C t , A t are dimensionless. In order to work with dimensionless variables we make the following replacements where we recall that L is a length scale in the solution related to the D3-brane charge through (2.4). Upon this replacement the action scales homogeneously as with S independent of L. The effective Newton's constant (the normalization in front of S) has dimensions of (length) −4 , in accordance with (2.56), since we rescaled all coordinates except the four Minkowski ones. We will use these redefined variables in our numerical code and, from this point onward, we will drop the over-bars in most places for notational simplicity. The action (2.35) and the length scale L are also invariant under the following rescaling of the five-dimensional fields: (2.59) Upon implementation of (2.56) and (2.59) the rescaled five-dimensional action becomes independent of L and Q f whereas, in terms of the rescaled fields, the right-hand side of Eq. (2.44) becomes proportional to the dimensionless combination This is a crucial observation because it implies that, up to trivial rescalings, the zerotemperature solutions only depend on N q . In practice this means that we can set Q f = 1 and Q c = 4L 4 = 1 in our numerical code, use Q st as the parameter labelling our solutions, and then replace Q st by N q and undo the rescalings above in the final answer to obtain the full dependence of the physics on N f , N c and N q . In Appendix B we provide the equations of motion for the rescaled functions of our ansatz. See the comments below Eqs. (1.2) and (1.3) concerning the s and the N c dependence of N q . There is a final scaling symmetry that leaves the action, L and Q f unchanged, and is given by Since this involves only dynamical fields it leads, via Noether's theorem, to a radially conserved quantity Notice that this quantity, evaluated at a horizon with the boundary conditions that the vector fields vanish there, equals the product of the Bekenstein-Hawking entropy and the Hawking temperature, therefore corresponding to the heat energy.

Numerical solutions
To solve the equations of motion given in Appendix B we resort to numerical methods. The approach we follow is to begin with the Q st = 0 black brane solutions of [12]. These solutions depend on one parameter, the radius of the black hole horizon r h , which is related to the temperature of the dual field theory. As we will explain below, adding a non-zero string density has an important effect in the IR part of the geometry if the temperature is low, but produces only a subleading correction in the UV. Therefore we can use the Q st = 0 solutions of [12] with large r h as seeds for solutions with non-zero but small string density. By taking small increments we can then span the (Q st , r h ) space and construct in this way solutions for different values of T and N q . We begin by fixing the gauge choice for the radial coordinate. We do so by specifying the functional form of the spatial component of the ten-dimensional Einstein-frame metric. In terms of the string-frame metric (2.8) we require that In Eqs. (2.28)-(2.29) we introduced one parameterization of the ten-dimensional string-frame metric that is particularly convenient for the purpose of reducing the system to five dimensions. In contrast, in order to perform the numerical integration of the equations of motion it is more convenient to use the parameterization in terms of which the condition (3.1) becomes This radial gauge fixing and parameterization are exactly those used in [12] for the neutral case. As we will review below, as r → ∞ the geometry approaches a hyperscaling-violating (HV) geometry with θ = 7/2 that encodes the physics of the Landau pole. Consequently, the dilaton diverges at r → ∞.
Using the equation of motion for h and the radial Einstein equation one can solve algebraically for c and c . Plugging this back into the remaining equations of motion we are left with five second-order differential equations for the functions b, f, g, φ and B which we do not show explicitly since they are not particularly illuminating. The bottom line is that we have reduced the problem to solving for five different functions subjected to ten boundary conditions. This is consistent with the fact that the boundary conditions in the UV and in the IR each depend on five parameters, as we will show now.

Boundary conditions in the UV
In the neutral case it was shown in Ref. [12] that the asymptotic UV solution is given by an HV metric with parameter θ = 7/2, supported by running scalar fields. Setting Q st = 0 in our five-dimensional action (2.35) this asymptotic solution takes the form where subleading corrections are of order 1/R 1/2 . In this expression we have changed the radial coordinate according to r L in order to make the HV-form of the metric explicit.
The inclusion of a non-zero charge density does not modify this asymptotic solution at leading order but produces only subleading corrections (except for the function B, see below). Returning to the parameterization of Eq. (3.2) the corrected solution takes the form In our numerical code we have constructed the expansion in r −4 up to order r −54 . All the nonleading coefficients in this expansion are determined in terms of the five unspecified constants κ b , κ f , κ f2 , κ φ , and κ B that will depend on the temperature and the charge density. We note that the only difference at leading order between the neutral and the charged configurations is in the asymptotic value of the function B: In the neutral case this tends to zero, whereas in the charged one it approaches a non-zero constant. Looking at (2.22) we see that the specific value of this constant ensures the physical requirement that A t vanishes at, and therefore that there is no electric flux through, the Landau pole. This should be contrasted with the situation in the lower-dimensional setup of Ref. [20]. In that configuration there is no Landau pole (since the field theory is asymptotically free) and both B and A t vanish asymptotically.

Boundary conditions at the horizon
The presence of a finite temperature in the field theory is captured by the existence of a black brane in the gravitational solution. The boundary conditions at the horizon are the usual requirement that the blackening factor in (2.8) has a simple zero at r = r h and the rest of the fields attain a finite value: In our numerical code we have gone to seventh non-trivial order in the distance to the horizon in these expansions. All the non-leading coefficients are determined in terms of the five In the presence of a non-zero charge density the temperature of the solution depends not only on r h but also on Q st . Nevertheless, r h can still be thought of as a proxy for T since the latter is still a monotonically increasing function of r h for fixed charge, as illustrated by Fig. 3. From this point onwards we will typically work with the dimensionless temperature as well as with the dimensionless radial coordinate introduced in (2.56):

Zero-temperature limit
Our numerical analysis below will show that the IR geometry of the zero-temperature solutions can be understood on general grounds. Indeed, we will see that it is given by a Lifshitz spacetime with dynamical exponent z = 7 and running dilaton and B fields. In other words, despite the fact that in order to construct the solutions with non-zero temperature we only impose regularity at the horizon, this boundary condition automatically approaches a Lifshitz boundary condition as T → 0. The physical reason for this is that, in the deep IR, the tension of the D7-branes is a subleading effect with respect to the backreaction of the electric field on the D7-branes. The IR-limit of the zero-temperature solution is therefore the same as in the D3-brane + strings system without D7-branes of [15,18], and it takes the form together with a rapidly vanishing B ∝ r 10 . With respect to the radial variable u in [18] we have performed the rescaling u = 2 5/8 17 5/24 11 5/12 r .
At zero temperature the blackening function is simply b(r) = 1. At low temperatures it is approximately given by (3.12) The dynamical exponent z = 7 can be seen by inspection of the Einstein-frame metric given within the square brackets of (3.10). 1 We have left undetermined the normalization of the time coordinate, encoded in the constant c t . This factor contains UV information via the RG flow connecting the UV-normalized time coordinate to the IR one of (3.10). The ten-dimensional solution (3.10) can be seen to arise as the uplift of an exact solution of our setup (2.35) in the following limit. 2 First one should perform a Legendre transform in the action to trade the radial derivatives of A t in favor of the charge density parameter Q st , with the help of the relation (2.22). Then one must take Q f → 0 keeping Q st finite. The result of this procedure is that the DBI action of the D7-branes becomes a Nambu-Goto action for a set of fundamental strings [8]. In other words, in this limit the tension of the D7-branes is subleading with respect to the effect of the charge density -see Eqn. (C.7). Once this limit is taken it is a simple exercise to check that the equations of motion for B and for the metric are compatible with setting B = 0 and G f = G b , meaning that the full symmetry of the five-dimensional compact manifold is recovered in the deep IR, i.e. the squashing vanishes in this limit. An exact scaling solution of the resulting equations was first found in [27] and studied in detail in [15,18].
The argument in the previous paragraph just shows that (3.10) is an asymptotic IR solution in the flavorless limit with an external charge density. However, in Appendix D we provide an analysis of the deformations away from this limit once the flavor effects are included. It is shown that the irrelevant deformations of the solution (3.10) (those that vanish in the IR) are still determined by five independent constants of integration, in agreement with the analysis in Sec. 3.2. In particular, the function B vanishes for r → 0 as with β 1 one of the integration constants.

Numerical integration
We have seen above that the UV and IR limits of the solutions that we are seeking depend on five integration constants in both instances. Our goal in this section is to show that there is a choice of these ten parameters that allows us to join smoothly the UV and the IR limits. As explained below (2.60) we can set Q c = Q f = 1 without loss of generality. In particular this makes the radial coordinate r dimensionless. In most places in this section we will not indicate this explicitly, but in the plots we will do so by writing r instead of r.
Since we set Q c = Q f = 1 our equations depend only on one external parameter Q st , which at the end can be replaced by N q through (2.60). The horizon radius r h enters the equations of motion via the IR boundary conditions and the range of integration. Therefore the solutions depend on r h and Q st , or equivalently on T and Q st .
Our integration strategy is simple. We use the UV and the near-horizon expansions (3.6) and (3.7) to shoot from r = 60 and r − r h = 10 −8 , respectively, at fixed values of r h and Q st . Imposing continuity of the functions and their derivatives at an intermediate radius gives us ten conditions that allow us to completely fix the five UV and the five IR free parameters for every value of r h and Q st , by an iterative Newton-Ralphson method. We declare continuity when the difference between the UV-and the IR-integrated functions is no larger than 10 −50 .
As an initial seed we use the neutral solution constructed in Ref. [12]. In solutions with a relatively large value of the horizon r h ∼ 1 we are able to include charge up to a critical value Q st = Q * st ≈ 10.16. Using these new charged solutions we then fix the charge density and start to decrease r h (equivalently, the temperature) until the Lifshitz scaling of Eq. (3.10) is observed. The smallest non-zero value of the charge for which we construct a solution is Q st = 1/100.
One outcome of the numerical integration is the value of the IR and UV parameters as a function of Q st and r h . These parameters encode much of the physics of the solutions. For example, we will see in Sec. 4 that the thermodynamic properties can be extracted from the UV parameters. Here we will see that they help elucidate the geometric properties of the solutions.
To see the emergence of the Lifshitz solution (3.10) for small values of r h we examine the behavior of the functions with respect to the horizon radius at fixed charge. First of all, comparing the metrics (3.2) and (3.10) and using the form of the blackening factor b in (3.12) we expect that, for sufficiently small r h , the constant b h in (3.7) will scale as (recall that (small horizon and finite charge) , (3.14) which is a large deviation with respect to the b h = 4/r h behavior of the neutral configuration. The Q st -dependent proportionality constant b 0 h in (3.14) cannot be obtained analytically since it depends on the relative normalization between the IR and the UV time coordinates, and hence on the entire numerical solution. For the dilaton we see from (3.10) that its value at the horizon should scale with the horizon radius as  Fig. 4. We see that for large values of r h all the curves converge to a single one, reflecting the fact that the leading UV behavior is controlled by the neutral solution. In contrast, for sufficiently small values of r h the charged curves deviate from the neutral one and approach straight lines with slopes that agree precisely with those predicted by Eqs. (3.14) and (3.15). For the dilaton the normalization also agrees with that in (3.15), whereas for the b h parameter we have performed a fit to (3.14). As expected, the value of r h at which the IR behavior sets At a qualitative level, one would expect the value of the dilaton at the horizon for a solution with non-zero temperature to be approximately the same as the value of the dilaton at the position r = r h for a zero-temperature solution. In other words, the scalar fields in a solution with non-zero temperature should be similar to those in a solution at zero temperature cut off at the corresponding value of the radial coordinate. This is confirmed by the similarity between Fig. 4(left) and The fit of the quantity b 0 h as a function of Q st is presented in Fig. 6. From this figure we observe that b 0 h diverges as Q st → 0 approximately as Q −9/2 st . At Q st ∼ 7 the behavior of the function changes and eventually starts increasing again with increasing Q st . This change of behavior at Q st ∼ 7 is also observed in other numeric parameters presented in this section. As we will see in the next section, at these values of the charge density a thermodynamic instability has appeared, but we have not found a direct relation between these two properties.  The two IR parameters of the metric, namely f h and g h , are shown in Fig. 7, where we plot the ratio with respect to the horizon radius, such that the curves asymptote a constant at small values of r h .  Both in the log-AdS region and in the Lifshitz region the squashing of the S 5 disappears and f h and g h become equal to one another. For the neutral solution the two parameters approach 1 logarithmically, whereas in the charged case they approach 136 1/4 /11 1/2 ≈ 1.029. Moreover, at large r h both parameters must approach their value in the neutral solution, where squashing is present. The value of r h at which the curves transition from the squashed behavior at large-r h to the non-squashed behavior at small r h depends on Q st .  Finally, Fig. 8(left) shows the horizon value B h . We see that at small r h the slope is that determined by the leading term in (3.13), i.e. B h ∼ r 10 h , whereas at large r h the horizon value of B approaches the negative constant corresponding to the leading term in (3.6). Since this constant is negative we cannot observe the behavior in the logarithmic plot we present. In the same plot, the curve corresponding to Q st = 10 is not present since it is negative for all values of r h . This is again a consequence of the generic change of behavior in the IR and UV parameters that occurs when Q st ∼ 7, as commented above.
To summarise so far, we have seen that the large-and small-r h values of the IR parameters of our numerical solutions match precisely those predicted by the LP and Lifshitz geometries, respectively. This strongly supports the fact that our configuration interpolates between the UV LP geometry of the neutral solution and the IR Lifshitz one. Moreover, the behavior of the dilaton at intermediate values of r h for small values of Q st supports that in this situation the transition between the LP and the Lifshitz geometries proceeds through an intermediate log-AdS region.
We now turn to the UV parameters in our numerical solutions. We expect their large-r h values to match those of the neutral one. This is confirmed by Fig. 8(right), Fig. 9 and Fig. 10.
In particular, κ B and κ f2 decrease as r −12 h and r −8 h , respectively. In contrast, the small-r h values of the UV parameters are not predicted by any general argument. The reason is that, even if the Lifshitz geometry emerges in the IR, the UV is still governed by the LP geometry of (3.6), with the free parameters entering in the subleading terms. Therefore the Q st -dependent values of these parameters for small r h are a prediction of our numerical solutions. One particular value that will play an important role is that of κ b . This approaches a negative constant at small r h if Q st is small, but for Q st 7, as in the Q st = 10 curve of Fig. 9, the parameter κ b tends to a positive value at small r h . In Fig. 10(right) the curve for Q st = 10 goes to a constant at r h = 0 approximately equal to κ f 2 = 17Q 2 c . The presence of an IR Lifshitz geometry in our solutions can be further verified from the scaling of the different functions. For any ψ(r) this scaling behavior is extracted as r ψ /ψ.   changes we evaluate this expression at the horizon. Thus we define and similarly From the Lifshitz solution (3.10) we expect to find that in the small-r h limit  The value of R φ extracted from our numerical solutions is shown in Fig. 11. As expected, we see that R φ approaches 4 and 6 at large and small r h , respectively, for all solutions. Moreover,

Regime of validity
We are now ready to determine the regime of validity of the solution. Since we are only interested in parametric dependences, in this analysis we will ignore all purely numerical factors. Moreover, we will focus on the zero-temperature limit, since turning on a non-zero temperature would simply cut-off the zero-temperature solutions. Similar discussions can be found in [18,12]. We must require that both supergravity and the smeared DBI action for the D7-branes be valid. We begin with supergravity. The first condition that we must impose is If this is not satisfied then string loop corrections become important and degrees of freedom not included in supergravity, such as D-branes, become light. Since Fig. 5 shows that the dilaton is a monotonically increasing function of r, we expect that the most stringent implication of (3.20) is obtained from the UV behavior of the dilaton. Using the asymptotic form of the solutions given in Secs. 3.3 and 3.1 the dilaton condition becomes UV: As expected, the UV result coincides with that obtained for the neutral solution in [12].
The other two conditions that we must require in order for supergravity to be valid are that the curvature of the string-frame metric be small in string units, and that the curvature of the Einstein-frame metric be small in Planck units. 3 The result can be understood simply by considering the two asymptotic forms of the solution in the IR and in the UV. The reason is that, as shown in Fig. 12, the curvature as measured by the square of the Ricci tensor, Ric 2 = R mn R mn , exhibits two simple behaviors separated by a rapid crossover around r ∼ 1. We will see below that the only exception to this statement occurs when the charge density N q becomes parametrically small. In this situation the flow follows the neutral flow and penetrates arbitrarily deeply in the log-AdS region, where eventually the supergravity + DBI description ceases to be valid [12]. This situation will not be of interest to us since, as we anticipated in Sec. 1, all the interesting physics takes place at values of N q of order 1.
As in the case of the dilaton, the two characteristic behaviors of the curvature seen in Fig. 12 at small and large r are controlled by the IR and the UV asymptotic solutions of Secs. 3.3 and 3.1, respectively. We have checked that the Ricci scalar, R = G mn R mn , the square of the Ricci tensor, Ric 2 = R mn R mn and the square of the Riemann tensor, R mnpq R mnpq , behave in the same way parametrically. Since presumably the same is true for other curvature invariants it suffices to examine Ric 2 . Through explicit calculation we find that the conditions that the string-frame curvature be small in string units take the form UV: IR: These r-dependences match the r 1 and the r 1 behaviors shown in Fig. 12, respectively. The curvature of the Einstein-frame metric in Planck units turns out to be proportional to 3 Note that in our conventions 4 p ∼ 4 s because we are not factoring the dilaton into a constant times a position-dependent part.
the string-frame curvature in string units, both in the IR and in the UV: It follows that the condition (3.20) and the requirement that the string-frame curvature be small in string units imply that the Einstein-frame curvature are small in Planck units. In fact, the Einstein-frame curvature is non-divergent in the IR, as expected for the Lifshitz times a constant-size sphere solution (3.10). 4 Therefore we will ignore the Einstein-frame curvature in the following. Since Eqs. (3.21) and (3.23) must be valid at the transition point r ∼ 1, it immediately follows that we must have the hierarchy Under these circumstances the IR conditions (3.22) and (3.24) are automatically satisfied at the transition point for values of N q = O(1), and the supergravity description is valid over the range where the lower bound comes from (3.24) and the upper bounds come from (3.21) and (3.23).
We now turn to the constraints imposed by the requirement that the Abelian DBI action for the D7-branes be valid [29,30]. The first requirement concerns the characteristic distance between nearby D7-branes, and it consists of two complementary conditions. On the one hand, this distance must be small in macroscopic terms in order for the distribution to be treated as continuous. This simply implies that N f 1. On the other hand, this distance must be large in string units, since otherwise strings stretching between nearby D7-branes would become light and the non-Abelian nature of the DBI action would become important. Since all the D7-branes wrap the η-fiber in the internal geometry of the metric (3.2) we must consider their separation in the KE base. The characteristic size of this manifold is with G st KE the coefficient in front of ds 2 KE in the metric (3.2). Since inside the four-dimensional KE base the branes are co-dimension two objects, one may effectively think of them as points in a two-dimensional space of volume ∼ 2 . The volume available to each of the branes is therefore 2 /N f . As a consequence, the typical inter-brane distance is 2 /N f . The requirement that this distance is large in string units is thus (3.29) Using the UV and the IR asymptotics we obtain UV: IR: The UV condition implies the hierarchy which is more stringent than (3.26). Similarly, the IR condition implies that which is more stringent than the IR part of (3.27). Note that the hierarchy (3.32) does not determine which of the two upper bounds on the right-hand side of (3.27) is more stringent. The second requirement for the DBI action to be valid is that the effective coupling between open strings be small. In the absence of smearing this coupling would be e φ N f . However, in the presence of smearing not all the N f branes but only the fraction contained in a volume of string size can participate in a characteristic process involving open strings. As argued above, this fraction is N f 2 s / 2 . The requirement that the effective string coupling (squared) be small is therefore This follows automatically from (3.20) and (3.29), so it does not yield any additional independent conditions. Putting together the various constraints coming from supergravity and from the DBI action we conclude that, in order for the gravity-plus-branes description to be valid, the number of colors and the number of flavors must obey (3.32). The region of validity in terms of the r coordinate is parametrically large and is given by

Thermodynamics
We will now compute the thermodynamic quantities necessary to study the equilibrium properties of the theory. This goal is greatly simplified by the way in which we wrote the fivedimensional action (2.35). The key point is that, since we traded G 3 = Q st dx 123 in favor of F 2 , none of the fields in our solution has purely spatial components turned on. Combined with the radial dependence of all fields this implies that the on-shell Lagrangian density is given by Since this is a total derivative we can write the on-shell Euclidean action I 5 as where τ is the Euclidean time and the only contribution comes from the evaluation at the LP. The relation between the functions in the five-dimensional effective metric (2.29) and those in the ten-dimensional metric (3.2) is g tt = −b g xx , g xx = 2 · 2 2/3 L r  Using the expansion (3.6) in these formulas one easily obtains the asymptotics for the fivedimensional functions. Substituting the result in (4.2) one finds that the on-shell Euclidean action diverges as We must therefore regularize the integral by introducing a cutoff for the radial coordinate, add the necessary boundary terms and take the r → ∞ limit. We first add the standard Gibbons-Hawking piece, written in terms of the trace of the extrinsic curvature at a constant-radius slice. Using the UV expansion (3.6) we get As shown in Ref. [12], counterterms associated to the asymptotically HV metric characterising the LP can be obtained via analytic continuation from the standard counterterms for asymptotically AdS spacetimes, in an analysis similar to the one in [31]. In the present situation the action (2.35) is more complicated than the one in [12] due to the inclusion of massive and massless vector fields, which also enter via the DBI action, not just via a Maxwell-like action. However, as stated above, the near-LP behavior of the solutions is dominated by the neutral, supersymmetric sector of Eqn. (3.4). In particular, the divergences are totally determined by the metric and the scalars, so we can simply take the same counterterm as in [12], which is proportional to the superpotential (2.41). Explicit evaluation of this term gives Adding up the three pieces we obtain a finite result in the r → ∞ limit: As usual, we will identify this result with the grand canonical free energy density 5 (4.8) By varying I GH and I W (I 5 gives a vanishing contribution) with respect to the induced metric on the constant-radius slice we obtain two divergent expressions whose sum yields a finite result for the boundary stress tensor: The fact that the pressure and the free energy (4.8) differ only by a sign confirms that G is the free energy in the grand canonical ensemble. Recalling now that the quantity h defined in (2.62) is equal to T s and that it is radially conserved we obtain an expression for T s in terms of the UV data: It is now easy to check that (4.12) We will now see that the right-hand side of this expression is precisely of the form −µ Q st , with µ the quark chemical potential dual to Q st . Indeed, from the action (2.35) we can compute the momenta conjugate to the three different vectors occurring in (4.12). For A t and A t this leads to two constants whereas for the massive vector C t it leads to a function of r: (4.14) In this way, we can rewrite (4.12) as As usual in holography, the asymptotic values of the conjugate momenta are identified with conserved charges, whereas the asymptotic values of the vector fields themselves are identified with the dual chemical potentials. Therefore the right-hand side of (4.15) has the form However, due to our boundary conditions there is only one independent charge Q st . This is already apparent in the fact that both conjugate momenta in (4.13) are proportional to Q st . Moreover, using the asymptotic expansion (3.6) we see that it is also true for Π C since As a consequence, all three charges vary simultaneously when Q st changes, and the contribution to the change in free energy is of the form µ dQ st with an effective chemical potential given by In terms of this Eq. (4.15) can be written as which is nothing but the usual definition of the grand canonical free energy. Similarly, we can now define the free energy in the canonical ensemble as 6 20) which in terms of UV data takes the form Note that the third term inside the brackets comes from the term proportional to C t (r LP ) in (4.11) or (4.18). This term has been evaluated explicitly to using (2.49) to relate C t (r LP ) to B (r LP ) and (3.6) to evaluate B (r LP ). Note that (2.49), together with the near-horizon expansions, implies that at the horizon This immediately implies C t (r h ) = 0, which is nothing but the regularity condition for a vector field at the horizon. Analogous conditions must also be imposed on A t and A t in order to compute A t (r LP ) and A t (r LP ), which are obtained by integrating (2.46) and (2.52) from r h to r LP , respectively, with the boundary conditions that A t (r h ) = A t (r h ) = 0. At this point it is useful to comment on the dimensions of several quantities. We have already defined a dimensionless temperature T in (3.8). Similarly, we will work with a dimensionless entropy density given by where we recall that V 5 is the volume of the SE manifold and we have used the expression (2.5) for the five-dimensional Newton's constant. Since the dimensions of C Q and s are the same we will use the same factor to make C Q dimensionless: Since E and P (and F, G, etc) have the same dimensions as T s we use the product of the two factors above: We refer to Q st as the charge density because of the proportionality relation (2.27), but note that it has dimensions of (length) −1 . Therefore we use precisely as its dimensionless counterpart, as defined in (2.60). Since µ Q st must have the same dimensions as T s it follows that the chemical potential has dimensions of (length) −3 . We emphasise that this is simply an unusual convention with no physical implications. Our dimensionless chemical potential is therefore The necessary coefficient to make the charge susceptibility χ dimensionless then is just the ratio of that in Q st over that in µ, i.e.
Finally, we anticipate here that the dimensionless versions of the functions f, g to be introduced in Eqs. (5.10) and (5.12) are We close this section with plots in Fig. 13 of the energy density, E, the pressure, P , the enthalpy, E + P , and the chemical potential, µ. The sign of the enthalpy will play an important role when we discuss the charge diffusion constant below. From (4.10) we see that E + P ∝ −κ b . As we emphasised around Fig. 9(right), κ b is always negative except in a small region at low temperature and high charge contained within Region III of Fig. 2. In this small region, shown explicitly in the bottom, right-hand corner of the E + P plot of Fig. 13, κ b turns positive and E + P becomes negative. Although our numerical results suggest that the chemical potential becomes slightly negative at high temperature, our resolution is not good enough to establish this definitively.

Instabilities
We will now examine the local thermodynamic stability of the system. This is equivalent to the positive-definiteness of minus the Hessian in the grand-canonical ensemble, namely of the susceptibility matrix where all derivatives with respect to T and µ are taken with the other one constant. The system is stable if and only if The first condition is the positivity of the charge susceptibility. This can be conveniently rewritten in terms of the canonical free energy F as The second condition can be related to the specific heat at constant charge: Indeed, using the chain rule and the equality of the crossed derivatives of F we find that In Fig. 14(left) we show a contour plot of the specific heat of our solutions as a function of T and √ 2 N q , while in Fig. 14(right) we show slices at several different values of the charge. We see that C Q is positive for T 0.8 but becomes negative above this value. This property is also illustrated by Fig. 15, where the negative slope of the entropy density curves as a function of T is evident at high T . In Figs. 14(right) and 15 we see that all curves converge to the same one for T 0.8. The common curve is the same as in the neutral case studied in [12]. As explained in that reference, the negativity of C Q in this region is an UV effect associated to the presence of the LP. Since in this paper we are interested in IR physics that is safe from LP effects, we will not elaborate further on this instability. In Fig. 16 we show a contour plot of the inverse charge susceptibility of our solutions as a function of T and √ 2 N q . We see that the region where χ −1 is negative includes roughly the same UV region as in the case of the specific heat. Presumably this is also a UV effect related to the presence of the LP. However, unlike the specific heat, the inverse charge susceptibility is also negative at arbitrarily low r h and T for charge densities √ 2 N q 0.8. The negativity of χ signals an instability towards charge clustering, which suggests that the putative, stable phase in this region may break translational invariance spontaneously. We will come back to this point in Sec. 6.
A further instability comes from the speed of sound of the system. For a system with both energy and charge density this is given by (see e.g. [36]) After some manipulations, c 2 s can be expressed in terms of our UV data as where˙= d/dQ st and = d/dr h . In Fig. 17(left) we show a contour plot of the result, and in Fig. 17(right) several slices at constant charge. As for C Q and χ, the region where c 2 s is negative includes the UV region r h 1, but also the IR region roughly defined by √ 2 N q 0.8 and 0 ≤ T 0.25.
At asymptotically low temperature we can use the fact that the far IR is controlled by a Lifshitz geometry to verify the negativity of c 2 s . The Lifshitz asymptotics imply that, at sufficiently low T , the entropy density must scale as T 3/z with a possibly charge-dependent coefficient, namely The function f (Q st ) can be easily extracted from the numerics and is shown in Fig. 18(left). Since where g(Q st ) is a function of the charge only, related to the chemical potential through Using this formula we extract the function g (Q st ) from the numerically computed chemical potential. The result is shown Fig. 18(right). We can now compute all the quantities of interest. In particular (5.14) Given that g (Q st ) exhibits a maximum at around Q st = 0.8, we infer that χ −1 changes sign, becoming negative, at that value of the charge. Using (5.12) and the standard thermodynamic relations for the energy and the pressure E = F + T s , we can compute (5.8) at T = 0 with the result Since µ is positive at low temperature (see Fig. 13) this result shows that at T = 0 the speed of sound squared and the charge susceptibility change sign at exactly the same value of the charge density.

Discussion
We begin by addressing the fact that we used factors of s and N c in our definitions of dimensionless quantities such as the temperature (3.8), the entropy density (4.24), the charge density (2.60), etc. The first observation is that the s factors cancel out in dimensionless ratios of physical quantities. For example, ignoring purely numerical factors we have that Assuming that s/T 3 is of order 1, we see in the first ratio that not only s cancels out but also that we seemingly get the scaling of the entropy density with N c expected from the number of color degrees of freedom in the theory. However, this scaling is actually ambiguous. The reason is that in the region where s/T 3 is of order 1, namely around the transition point between the LP and the Lifshitz geometries, so is the combination where λ is the 't Hooft coupling (2.3). This means that, parametrically, at that scale one can freely replace factors of N f and/or N c with powers of λ. The same ambiguity affects other physical quantities such as the second ratio in (6.1). Ultimately, this ambiguity arises from the lack of a privileged scale or a fixed point at which to anchor the value of the coupling. The second observation concerns the relation between dimensionful gauge theory quantities and the Landau pole scale. While the conventions we have chosen are convenient on the gravity side, from the gauge theory viewpoint it may be more natural to measure dimensionful quantities in units of the intrinsic scale of the theory, namely in units of Λ LP . In particular, it may be natural to construct the phase diagram by changing N q and T while keeping Λ LP fixed instead of s . In order to do this one must define the energy scale associated to the Landau pole. A convenient and well-motivated choice is to define Λ LP as the mass M of a string stretching from the IR bottom of a zero-temperature geometry all the way up to the LP [12]. Fig. 19 shows this mass, normalized to the mass M 0 in the neutral case, as a function of N q . The important message of this plot is that not only M/M 0 does not change parametrically as N q is varied but also that it changes rather smoothly. This means that, in order to refer all our dimensionful quantities to the LP scale instead of to s , we would need to do a rather mild N q -dependent rescaling of our results. This would slightly distort our phase diagram but it would leave unchanged all the qualitative conclusions, which we now proceed to discuss. One of the most interesting outcomes of our analysis is the presence of instabilities in our system, some of which are summarised in Fig. 2. The most relevant ones are those that are present at low temperature, since they are the least sensitive to the UV completion of the theory.
One of these instabilities is associated to the negative speed of sound squared, c 2 s < 0, that is present in Region III of Fig. 2. This indicates a dynamical instability for the hydrodynamic sound mode, whose dispersion relation takes the form Negativity of c 2 s implies that c s is purely imaginary. This indicates that, if perturbed by a small-amplitude fluctuation with sufficiently long wavelength, the system will develop an inhomogeneous profile that will initially grow exponentially in time as e |cs|kt . Presumably, when the non-linearities become important the system will settle down to an equilibrium, inhomogeneous configuration in which translation invariance is broken spontaneously. It would be interesting to identify the endpoint of this instability under dynamical evolution along the lines of [32,33]. Note that this is related, but not identical, to identifying all the possible equilibrium, inhomogeneous states of the system. For example, the latter may include states with domains of finite characteristic size, namely crystalline phases, as well as phase-separated configurations in which two semi-infinite, homogenous phases separated by an interface coexist. From a thermodynamic viewpoint one of these inhomogeneous states will be absolutely preferred, rendering the rest metastable. However, the lifetime of some of these metastable states may be very long in general, and in particular this will be the case in our large-N c limit. Moreover, to which final state the initial unstable configuration will evolve is a dynamical question that depends on the "landscape" of configurations above and beyond the purely equilibrium ones.
Similar considerations apply to the instability associated to the negative charge susceptibility, χ < 0. As explained in Appendix E, this is related to the charge diffusion constant through In this expression σ is the electrical conductivity, which must be positive in order for the divergence of the entropy current of first-order hydrodynamics to be non-negative, i.e. in order for the second law of thermodynamics to hold. A negative value of D indicates a dynamical instability towards charge clustering (anti-diffusion), thus also generating inhomogeneities.
Our numerical results allow us to establish that this is the situation in a subregion of Region II in Fig. 2, where all the factors in (6.4) are positive except for χ, but they do not allow us to establish definitively whether c 2 s becomes negative in some small subregion of Region II. If this were the case then D would be positive in this subregion.
In contrast to Region II, D is positive in most of Region III, since there all factors are positive except for both χ and c 2 s . Our numerical resolution allows us to establish the existence of a subregion of Region III where E + P becomes negative. Whenever E + P < 0 in Region III one finds consequently, from Eq. (6.4), that D < 0, thus indicating a dynamical instability in the charge diffusion channel. This change of sign of the enthalpy is visible in Fig. 13 but not in Fig. 2 because it lies at charge densities larger than those shown in the latter figure.
It is interesting to note that the subregions where E + P < 0 suffer from a further hydrodynamic instability in addition to those associated to the negativity of c 2 s and/or D. Indeed, the transverse momentum density in a plasma (the density of momentum pointing orthogonally to the direction of propagation of the wave) obeys a diffusion equation with diffusion constant given by [36] D P = η E + P . (6.5) Since η must be positive for the same reason as the electrical conductivity σ, we see that D P becomes negative if E + P < 0, leading to an instability towards transverse momentum clustering.
A contour plot of the combination E + 3P is shown in Fig. 20. This quantity is interesting because it would control the backreaction of the fluid on the spacetime metric if we were to couple the fluid to dynamical gravity. As is well known in the context of Cosmology, the condition for accelerated expansion is E + 3P < 0. In our context this is not particularly relevant because the phase where this condition is satisfied is unstable. However, it is amusing to consider the question of whether a stable quark-matter phase with E + 3P < 0 may exist, and if so what would be the possible consequences if this phase were to be realised in e.g. the interior of a neutron star.
To conclude, we reiterate that some of the instabilities that we have identified suggest the possible existence of quark matter crystalline phases in our model, but establishing this definitively requires further analysis. We hope to report on these issues in the near future.
In the last line we have used the equations of motion (A.3). The duality relations to be imposed after the equations of motion have been obtained are with * * = 1 for odd forms. Then the equations of motion for G 7 and G 9 become the Bianchi identities for G 1 and G 3 (as usual G 5 is self-dual and its equation of motion is its Bianchi identity). Add now the sources piece coming from the D7-branes. The WZ term gives a linear coupling of the worldvolume fields to the RR forms Here Γ is a two-form that describes the distribution of the D7-branes and with A the BI field living in the worldvolume of the D7-branes. Note that The linear source in the WZ term modifies the equations of motion for the RR forms which now read (we denote by F n the RR field strengths in the presence of sources, as opposed to the unsourced G n ones) In order to impose the equivalent duality relations to (A.5), i.e.
we must modify the definitions of the RR field strengths to where γ is a one-form satisfying dγ = −Γ . (A.12) In section 2 we have written γ = Q f η , With this definition of the field strengths the equation of motion for the NS forms gets two extra contributions with respect to the one derived in Eqn. (A.4), first from the explicit F terms in the WZ action (A.6) and second from the implicit F terms in the definitions of the F n in (A.11): where in the last equality we have used the equations of motion (A.9). The equation of motion for the dilaton and the metric can be obtained from the same democratic-formulation action, adding the DBI piece. Coming back to the non-democratic formulation of supergravity, the equations of motion for the different F n and H fields, as well as the ones for the dilaton and the metric, can be obtained from the ten-dimensional action found in Einstein frame in [25], where a five-dimensional reduction over the compact manifold preserving the SU(2) structure of the Sasaki-Einstein manifold is given.

B Equations of motion for rescaled variables
In this section we write the equations of motion for our ansatz explicitly, for string frame metric. We give the scaled quantities described in Sec. 2.3, and for simplicity we omit all bars in the different variables. The scaled electric field satisfies and the B equation of motion is In order to write the remaining equations of motion let us define the combination In terms of this, the metric and dilaton equations are

C A probe in the asymptotic geometries
In this Appendix we provide more evidence supporting our claim that the zero-temperature IR geometry is dominated by the asymptotic solution (3.10) and the UV by the LP geometry (3.4). To this end we note that the addition of a (charged) probe D7-brane in these solutions does not have an important backreaction on the background configurations in the radial regimens of interest, i.e. the origin for the IR geometry and the boundary for the UV one. An analogous discussion for a lower-dimensional case was given in Ref. [20], and we provide here a slightly different (but equivalent) argument. Consider a set of n f probe D7-branes with a worldvolume gauge field, a, turned on. The action describing this set of branes is (in string frame) where T D7 n f is the total tension, and we have assumed that there is a non-trivial potential C 6 in the background of the form given by Eqn. (2.17): We also allow for a non trivial NS potential, which we will take in the B rt directions to be aligned with da = a t (r)dr ∧ dt. For a massless probe brane the embedding profile is a constant in the transverse directions to the brane. Therefore, working in the static gauge we can write the probe action as where we have defined the following functions with V 3 the volume of the 3-cycle wrapped by the probe branes and the metric functions those of the generic metric (2.8). The function H 3 comes from C 8 via Hodge-dualization of F 1 = Q f η and H 4 has its origin in C 6 . The electric field on the worldvolume of the D7-branes has a first constant of motion given by δS D7 δa t = n q ⇒ a t + B rt 2π 2 2 , (C.5) whose backreacted version for the ansatz considered in this paper can be found in Eq. (2.22). Working at fixed charge density n q we can use the former expression to eliminate a t in favor of n q in the action by performing a Legendre transform. This results iñ

C.1 Charged flavorless setup and the IR geometry
Consider now the behavior of (C.6) in the limit n f → 0, but keeping n q finite. This is the limit described in Sec. 3.3. The Legendre-transformed action becomes and we recognize this as a smeared Nambu-Goto action where Ξ = n q dx 1 ∧ dx 2 ∧ dx 3 is a density of fundamental strings, extended in the radial direction and distributed on the spatial ones. The same limit can clearly be achieved asymptotically if one works in a radial regime of a geometry such that The solution (3.10) satisfies this criterium when r → 0, suggesting that near the origin this is a valid asymptotic solution. Actually, one should consider the first correction to the function B = 0 to ensure the requirement H 4 → 0 near the origin, i.e. that the leading correction to this function is vanishing when r → 0 (otherwise the probe approximation would fail). In Eq. (3.13) we show that indeed B → 0 in the backreacted setup.

C.2 Supersymmetric chargeless solution and UV geometry
From the probe argument (C.6) we would expect that the supersymmetric solution is valid asymptotically if near the Landau pole (n q +H 4 )/(H 1 H 2 ) → 0, provided the NS form vanishes. Plugging all the values from the asymptotic expansion (3.6) in the H i functions we obtain which near the Landau pole does not go to zero unless there is a precise cancellation between H 4 and the charge density. In other words, it is not possible to obtain a parametrically small quotient of the charge versus D7-tension effects (such that the charge is subdominant in the UV) by just going to sufficiently large values of r. Gravitationally, this is an effect of having the end of the geometry at a finite proper distance from any point in the bulk.
We argue now that the setup we have considered automatically tunes itself to avoid this issue: the function B goes to a constant value near the Landau pole that cancels exactly the contribution from the charge density, such that the backreacted version of (C.9) indeed goes to zero as one approaches the Landau pole. To see this take the full system of backreacted equations of motion, given in appendix B, and expand around the supersymmetric solution. In particular one can consider a perturbative expansion in which Q st = Q st and B = B, where is a book-keeping parameter to denote the inclusion of a small charge density in the supersymmetric setup. At first order in the equation of motion for B decouples and reads with all functions evaluated in the asymptotic solution (3.6). Close to the Landau pole we can approximate the equation by 11) and the particular solution is given at leading order by a constant By plugging this value in (2.22), which is the backreacted version of (C.6), we observe that at the Landau pole A t → 0, which allows us to construct the asymptotic expansion (3.6) around the supersymmetric solution.

D Scaling solution and IR geometry
In this Appendix we show that the configuration in Eq. (3.10) is an asymptotic solution in the presence of a finite number of flavors. One can perform an expansion of the full equations for small Q f , and solve order by order. At order O(Q 0 f ) we have Eq. (3.10) with B = 0 as a solution. To find the first order correction we expand the functions as where the metric components are Z = {tt, xx, rr, b, f } and the superindex (0) refers to the solution (3.10). The equation of motion for β(u) decouples and has the solution with β 1,2 integration constants. For the Q f -expansion to be well defined we need to set β 2 = 0 and we observe that as r → 0 the first-order correction β vanishes, and in particular the explicit solution to the non-homogeneous part of the equation goes like r 12 . The remaining first-order corrections satisfy a coupled system of equations, and the solutions are a combination of eight different modes which we group in pairs as with i = 1, 2, 3, 4. For fixed i and sign ±, the factors γ i,± Z and ϕ i,± are not all independent, but given in terms of just two free parameters, with one of these parameters corresponding to a gauge fixing of the radial coordinate.
We have grouped the expansion modes in pairs of fixed i. The two powers characterising each of the pairs add up to ∆ i,+ + ∆ i,− = −10. The first of these pairs is given by the values the ∆ 1,+ = 0 mode corresponding to a rescaling of time and the ∆ 1,− = −10 mode corresponding to turning on a temperature, with the corresponding non-vanishing coefficients given by where, once again, one of the two undetermined coefficients corresponds to a fixing of the radial coordinate. As seen from the negative value of the ∆ 1,− coefficient this is a relevant mode that modifies the IR geometry; since in this appendix we are interested in the zero temperature solutions, with the Lifshitz geometry in the IR modified just by irrelevant deformations, we take the coefficients multiplying the ∆ 1,− modes to be zero. The mode ∆ 1,+ = 0 is associated to a free γ 1,+ tt coefficient, corresponding to the choice of c t in Eq. (3.10). The remaining three pairs appearing in the solution are given by The corresponding coefficients γ i,± Z , ϕ i,± can be determined analytically. We give here only the relations between the two coefficients in the compact part of the manifold The modes ∆ 2,± and ∆ 3,± appeared already in the analysis made in [15,18]. In particular and γ 3,± f = γ 3,± b , which guarantees that these deformations do not contribute to the squashing of the five-dimensional compact manifold, preserving the full isometry group. However, γ 4,± f = −4γ 4,± b + 5 2 ϕ 4,± and these are the modes responsible for the breaking of the symmetry by means of the squashing of the compact manifold allowed in our ansatz (2.8), which is nothing but a consequence of the presence of flavor in the setup. Considering just the irrelevant deformations (those with ∆ i,± > 0) to have the Lifshitz solution (3.10) in the IR implies choosing the free parameters so that they cancel the ∆ i,− modes, allowing only the ∆ i,+ ones, with i = 2, 3, 4.
Thus, we have taken to zero all relevant deformations of the geometry, meaning that our perturbation holds always close to the origin, even if we relax the initial condition Q f 1 that lead us to the perturbation equations we just solved. In other words, we have just proven from the supergravity equations of motion with sources that there is a solution near the origin that asymptotes to the (3.10) solution, with corrections given by the 5 irrelevant deformations ∆ 2,+ , ∆ 3,+ , ∆ 4,+ , c t and β 1 .

E Charge diffusion constant
In this section we provide a short derivation of Eq. (6.4). We start with the expression derived on page 27 of [36]: where σ is the electric conductivity and α i , β i are thermodynamic derivatives defined in [36].
Using the expressions at the bottom of page 26 of that reference one can rewrite this as where χ ab is the susceptibility matrix Note that χ 33 = χ as defined in the first equation in (5.2). Using thermodynamic identities one can show that

F Super-critical solutions
In the main part of the paper we have limited our discussion to solutions with √ 2 N q below a critical value, √ 2 N q < √ 2 N * q 10.16. The reason is that all these sub-critical solutions share a common IR given by the Lifshitz asymptotics. For completeness we will present here the results for integrations with √ 2 N q > √ 2 N * q . To construct them, we follow the same integration strategy as before (see Sec. 3.4), i.e. we start with a solution with relatively large value of r h at fixed N q , and then decrease the radius of the horizon step by step using the previous solution as a seed. For low values of r h we observe a dramatically different behavior in the IR compared to the sub-critical case. In particular, the super-critical solutions do not flow to Lifshitz in the IR, but to a different scaling solution with negative specific heat.
The IR parameters of the expansion (3.7) extracted from the numerical super-critical solutions are shown in Fig. 21. We observe that for large values of r h all the curves converge to a single one, reflecting the fact that the leading UV behavior is controlled by the neutral solution. On the other hand, for small values of the horizon radius these solutions exhibit a scaling behavior different from the Lifshitz scaling observed for sub-critical solutions presented in Sec. 3.4. This IR scaling behavior takes the form e φ ∼ r 4 , e g ∼ r 2 , e f ∼ r −2 , B ∼ constant < 0 , (for r h → 0) .

(F.1)
This can be easily seen in Fig. 22, where we have plotted the functions B and R φ (r) = r ∂ r e φ e φ , R g (r) = r ∂ r e g e g , R f (r) = r ∂ r e f e f , (F.2) for √ 2 N q = 11 and r h / √ 2 = 0.5. The behavior of the blackening function b near the origin when r h → 0 is b ∼ constant > 0 , (for r h → 0) , (F. 3) indicating that Lorentz symmetry is restored in the deep IR for super-critical solutions in the limit in which the horizon disappears. The determination of the horizonless metric in the Irrespective of the actual value of this power c, the metric takes the asymptotic form of a hyperscaling-violating metric with parameter θ > 3, which is compatible with the null energy condition. This in turn implies a thermodynamic instability with negative specific heat, since the entropy scales with the temperature as s ∼ T 3−θ . Note that T diverges as r h → 0, but the entropy density tends to zero. For this reason r h , but not T , is a good proxy for the energy scale at which the solution is being probed. In particular, the limit r h → 0 can be thought as the limit in which the solution approaches the ground state of the theory. In this limit the super-critical solutions possess a curvature singularity at r = 0. The fact that Lorentz symmetry is restored in the IR geometry of the super-critical solutions in the limit r h → 0 suggests that the IR description is no longer dominated by the backreaction of the charge density, or equivalently by the fundamental strings dissolved inside the flavor D7-branes. This is further confirmed by the behavior of the two IR parameters of the metric, e f h and e g h . Since they scale differently the compact manifold remains squashed even at very small values of r h , as opposed to the situation for sub-critical solutions. Given that the squashing is caused by the D7-branes, this suggests that the tension of the branes is not subdominant in the IR for super-critical solutions. This is further supported by revisiting the conditions (C.8) for the particular scaling of Eq. (F.1). We have checked that the last of these conditions, namely H 4 ∝ B → 0, is not satisfied in the IR of the super-critical solutions. One finds instead that B approaches a negative constant. Since the source for this function is supported on the flavor branes, as explained in the paragraph leading to (2.17), we conclude that for these values of the charge density the flavor has an important effect in the IR.
The UV parameters appearing in the expansion (3.6), obtained from our integration in the super-critical cases, are presented in Fig. 23. Once again we observe that all curves approach the same one at large values of r h , since the neutral solution controls the UV, whereas at low values of the horizon radius there is a dependence on N q . We now turn to the thermodynamic properties of the super-critical solutions. In Fig. 24(left) we present the temperature as a function of the horizon radius, and in Fig. 24(right) we show the entropy as a function of the temperature. As observed in this figure, the temperature is a double-valued function of the radius of the horizon: at fixed charge there is a small and a large black hole corresponding to the same temperature, and temperatures below a certain (N q -dependent) critical one, T * (N q ), are not realized in our solutions. 10 10 The small black holes are described near the IR by the scaling solution discussed around Eqs. (F.1) and (F.3), and they are thermodynamically unstable since the entropy decreases with the temperature, making the specific heat negative. For the large black holes this also occurs if the temperature is sufficiently high, reflecting the presence of the Landau pole in the UV, as explained in Sec. 3.1. For a small range of temperatures just above T * (Q st ) the large black hole solutions possess a positive specific heat and therefore are not locally thermodynamically unstable against energy density fluctuations at fixed charge density. We have not examined whether other possible instabilities are present. In Fig. 25 we plot the free energy associated to several super-critical solutions. For comparison we also show the free energy of one sub-critical solution with √ 2 N q = 10. For super-critical solutions we see the cusp in the free energy typical of a first-order phase transi-tion. This suggests the existence of new additional solutions that may be locally stable and thermodynamically preferred at low temperatures, but we have not been able to obtain them with our numerical strategy.