A holographic model for QCD in the Veneziano limit at finite temperature and density

A holographic model of QCD in the limit of large number of colors, $N_c$, and massless fermion flavors, $N_f$, but constant ratio $x_f=N_f/N_c$ is analyzed at finite temperature and chemical potential. The five dimensional gravity model contains three bulk fields: a scalar dilaton sourcing ${\rm Tr}F^2$, a scalar tachyon dual to $\bar qq$ and a 4-vector dual to the baryon current $\bar q \gamma^{\mu} q$. The main result is the $\mu,T$ phase diagram of the holographic theory. A first order deconfining transition along $T_h(\mu)$ and a chiral transition at $T_\chi(\mu)>T_h(\mu)$ are found. The chiral transition is of second order for all $\mu$. The dependence of thermodynamical quantities including the speed of sound and susceptibilities on the chemical potential and temperature is computed. A new quantum critical regime is found at zero temperature and finite chemical potential. It is controlled by an AdS$_2\times R^3$ geometry and displays semi-local criticality.

1 Introduction The phase diagram of QCD, as a function of temperature T and chemical potential µ, corresponding to baryon density or some other conserved charge like isospin, displays a rich structure [1]. Particularly interesting and important features of the phase diagram are the nature of the chiral phase transition, the location of the chiral critical point and its properties. All these have been extensively studied both with effective chiral models [2,3] and other approaches reviewed e.g. in [4], and holography [5][6][7]. Since first principle lattice methods [8][9][10] are currently still limited to small values of µ/T , the model studies provide important complement. However, the location of the critical point is very dependent on the details of the models [4]. In addition to temperature and density, typically also external perturbations of the chiral symmetry, i.e. finite quark masses, are present and provide further dimensions to the phase diagrams. For example, in the limit of vanishing quark masses, the finite temperature phase transition of two flavor QCD at zero chemical potential is expected to be of second order. At low temperatures, the chiral transition at finite µ is expected to be of first order [2]. At intermediate temperatures and densities the first and second order transition lines were conjectured to meet at a tricritical point. The finite quark mass softens the singularity at the second order phase transition which becomes a smooth crossover. The line of first order transitions is unaffected in the presence of small external perturbation, and the tricritical point becomes the critical endpoint for this line. The fate of the existence of the critical point of QCD in the (µ, T )-plane at the physical value of quark masses is ultimately determined by the form of the critical manifold in the multidimensional space of parameters µ, T, m q [11]. One can imagine several possibilities to occur. Indeed, it can be that the existence of a critical point near the chiral limit implies that the critical point does not exist at physical masses. Depending on the shape of the critical manifold, a variety of other possibilities can be imagined.
Effective field theories utilizing holographic methods, motivated by the AdS/CFT correspondence [12][13][14], have become a major tool in the analysis of strongly coupled theories both in elementary particle and condensed matter physics [15,16]. A class of bottom up models for QCD-like theories, which captures the entire renormalization group evolution of the corresponding quantum field theory from weak to strong coupling has been developed in [17][18][19]. A particular application of this framework is the determination of the vacuum and finite temperature phase diagrams of the associated quantum gauge theories [18,[20][21][22][23][24][25]. The framework has been extended to account for the dynamics of chiral symmetry breaking in the presence of flavors [26][27][28][29][30][31]. In order to consider effects coming from the backreaction of flavor to color, holographic models with dynamics close to that of QCD in the Veneziano limit were explored and developed [31][32][33][34]. In this work we consider adding finite chemical potential 1 in order to determine the phase diagram in the (T, µ)-plane by computing the pressure p(T, µ; m q = 0) in the phases where chiral symmetry is intact or spontaneously broken.
Concretely, we consider the holographic model for equilibrium QCD with N f massless quarks at the limit N f → ∞, N c → ∞ and fixed ratio x f = N f /N c . For a thorough discussion of the fundamentals of this type of bottom-up holographic model for QCD in the Veneziano limit (V-QCD) at zero or finite T but zero density, we refer to [32,33]. Here we only outline the features arising when we allow also finite density and chemical potential in V-QCD. According to the holographic dictums to add baryon density we must turn on a source for the five-dimensional gauge field A a . The dynamics of the baryon number gauge field A a is determined by its appearance in the tachyon DBI action, which can be schematically written as − det (g ab + κ ∂ a τ ∂ b τ + w F ab ). (1.1) Here κ and w are couplings, F ab = ∂ a A b − ∂ b A a , and τ is the tachyon, sourcingqq. To turn on a uniform constant density, the Ansatz A a = Φ(z)δ a0 should be made, where z is the coordinate of the 5th dimension and the only non-zero component of F ab is F z0 = ∂ z Φ(z). The action contains only the derivative of Φ and the finite density arises as the integration constantñ of the equation of motion of the cyclic configuration space coordinate Φ. The three bulk fields λ, Φ, τ correspond to the three arguments in p(T, µ; m q ), and we will consider only the case m q = 0 in this paper and denote the pressure simply by p(T, µ). As in [33] we find that there are two types of m q = 0 solutions: those with vanishing tachyon (chirally symmetric) and those with nonzero tachyon (breaking chiral symmetry spontaneously). To determine the pressure, the strategy is therefore to find black hole solutions with one or two scalar hair (corresponding to the dilaton and tachyon scalars) and a non-trivial charge density. Such solutions, when they exist, compete also with finite temperature but zero charge solutions without a black hole. The reason is that these zero charge solutions always have a constant Φ = µ and therefore correspond to saddle points with finite chemical potential but zero charge density. Such solutions are expected to dominate at small enough temperature and chemical potential, and we identify them with the "hadron gas" vacuum phase with zero pressure. Increasing the charge density, we have the possibility of a trivial or non-trivial tachyon field. The latter possibility describes a "deconfined" but chirality breaking plasma, while the former corresponds to chirally symmetric plasma. To determine which of these two dominates, one solves numerically for the coupled equations of motion of the fields, and finds pressures p s (T, µ) and p b (T, µ) corresponding, respectively, to the solutions with intact or spontaneously broken chiral symmetry. Equality of pressures, temperatures and chemical potentials then defines the phase boundary on the T, µ plane.
The main outcome of this work is the phase diagram shown in Fig. 1 which was obtained for the theory with x f = N f N f = 1, namely for the same number of massless flavors and colors. The chiral transition is of second order for small µ and of first order for larger µ down to T = 0 with a tricritical point in between. There is also a tentative deconfining transition at T h (µ) between the chirality breaking plasma and the "hadron gas" phase discussed above. This phase boundary is determined by the condition p b (T, µ) = p low = 0. 2 To motivate this in the field theory, note that the degrees of freedom of the low temperature phase are the Goldstone bosons of the spontaneously broken chiral symmetry, and their number is ∼ N 2 f . On the other hand, the number of degrees of freedom in the high-temperature phase is ∼ 2N 2 c + 7 2 N c N f . As we consider only the case x f = 1 we obtain p low /p high ∼ 2/11 ∼ 0. The relative weight of the low-temperature degrees of freedom grows with x f , and ultimately at some x c 4, in terms of the free energy, they become indistinguishable from the high temperature ones. This signifies the quantum phase transition from a confining gauge theory to the one whose long-distance behavior at zero temperature is governed by a nontrivial and stable infrared fixed point. We leave the study of the finite temperature and density phases in the limit x f → 4 for a further investigation. All these features of this phase diagram agree on the general expectations. There is, however, a surprise at low temperatures: There exists a new quantum critical regime at T = 0, with exotic properties which realize the symmetries of the associated geometry, that is AdS 2 ×R 3 . This exists on the T = 0 segment of the chirality breaking plasma as well as on the T = 0 line of the chirally symmetric plasma. The presence of the AdS 2 × R 3 geometry in the holographic solution indicates that there is a scaling symmetry of the time direction which does not act in the spatial directions. Such symmetries have been called semilocal. This is an unexpected symmetry in a theory at finite density, but it is natural and generic in the holographic context [36], and appears even in simple black holes as the Reissner-Nordström black hole [37]. The physics in this critical regime is similar to that of a theory with zero speed of light: all spatial points decouple in the IR.
It is well known that such AdS 2 solutions are highly unstable as AdS 2 has a rather restrictive Breitenlohner-Freedman bound. The instabilities associated to the fields we consider can be understood in terms of the physics of the phase diagram and we describe them in detail in section 3.4. However, there can be further instabilities associated with other operators which we have not included here. It is possible that such quantum critical points play an important role in the appearance of color superconductivity and color flavor locking at high density.
There are many technical obstacles one has to cross before obtaining the final numerical results for the phase diagram: First, to find the relevant charged black-hole solutions one has to guarantee that the metric function f (z) vanishes at the horizon z = z h . As the horizon is a singular point of the equations, the numerical evolution must start close to the horizon with the appropriate boundary conditions. Second, the UV quark mass will be fixed to zero in order to have exact chiral symmetry. To implement this, we must solve the entire coupled set of equations of motion and tune the boundary conditions so that the leading term of the tachyon field at small z (near the boundary) is ∼ z 3 , instead of a linear one corresponding to a finite quark mass. This requires high numerical precision in the solution of the non-linear equations of motion. The third difficulty is that the quantity to be computed is a function of two variables, p(T, µ). The numerics is correspondingly parametrised by two parameters, the value of the dilaton at the horizon λ h and the integration constantñ. These parameters cannot be continuous ones, but one can determine, say, T (λ h ,ñ) as a function of λ h for fixed values ofñ, and vice versa. Proceeding in this way one obtains p(T, µ) on two grids on the T, µ plane (see Fig. 21). Fourth and final issue is that one has to guarantee, using scaling properties of the equations of motion, that all the physically dimensionful quantities are expressed in the same units.
All of these considerations make the numerical problem at hand challenging. In this paper we focus on the details of introducing the chemical potential to the model, limit ourselves to one set of potentials chosen from [33] and to one value x f = 1. This allows us to show that the method works, produces interesting results and motivates further studies. We have released the numerical code which has been used to compute the results presented in this paper [38].
In Section 2 we specify the model and give the equations of motion and their scaling properties. In Section 3 we find all solutions with constant scalars as they are critical end points of flows. They correspond to AdS 5 and AdS 2 geometries and we analyze their RG stability. In Section 4 we discuss the horizon expansion required for initialising numerical solution and the physical values of the parameters λ h ,ñ of numerical integration. The main numerical results for the pressures, the transition temperatures T h (µ), T χ (µ) and sound velocity are shown and discussed in Section 5.1. The Appendices contain a detailed discussion of the numerical solutions, examples of the computed values of T and µ and a detailed presentation of the chiral phase transition line on the plane of numerical parameters λ h ,ñ.
2 Action and the equations of motion

Definition of the action
The action of the model for vanishing chemical potential has been discussed thoroughly in [32][33][34]. We focus here on the additional terms needed to describe the finite baryon density. The action of the model is, in standard notation [33], where the Lagrangian is The metric Ansatz is The functions b and f of the metric, the dilaton λ = e φ , the tachyon τ and the bulk density Φ depend only on the extra dimensional coordinate z. The Gibbons-Hawking counterterm is implied. The Lagrangian (2.2) is parametrized in terms of the potentials V g (λ), V f (λ, τ ), κ(λ) and w(λ) which are chosen to satisfy two basic requirements. First, in the ultraviolet, i.e. in the weak coupling limit, the model should reproduce the known perturbative behaviors of the corresponding field theory. Second, in the deep infrared the model should lead to the generation of a dynamical wall shielding the singular behavior as λ → ∞, which is responsible for confinement in the absence of the tachyon.
In the numerical study in this article we will take the gauge coupling constant function w in the DBI action to be proportional to the other model function κ as 3 Here the scale L A ∼ L UV appears in order to match the dimensions correctly. It can be formally eliminated from the Lagrangian (2.2) by rescaling Φ. In Appendix F we shall find that L A ≈ L UV (x f = 0) = 1. Note that if one expands the Lagrangian (2.2) in F ab and writes it in the form − 1 4e 2 F 2 , one can identify a dimensionless coupling (2.5) Explicitly, the potentials are [33] V g (λ) = 12 (2.6) where the function V f 0 (λ) is given by The scale L UV has a nontrivial dependence on x f , L UV = L 0 (1 + 7 4 x f ) 1/3 , which is determined by matching the pressure to the Stefan-Boltzmann limit at µ = 0 [33]. The function κ(λ) is given by The numerical factors appearing in (2.6) and (2.7) simply provide the equivalence with the known perturbative behavior in the weak coupling limit. This matching is obtained via the definition β(λ) = dλ db/b , (2.9) and recalling that the 2-loop beta function for the coupling λ = N c g 2 (µ)/(8π 2 ) of the boundary theory is in the Veneziano limit. Analogously, the numerical factors appearing in κ(λ) in Eq.
(2.8) are obtained by first defining 11) and then relating to the quark mass anomalous dimension with the scheme independent coefficient γ 0 defined by The actual numerical value of the quark mass (and the condensate qq ) is fixed by the UV expansion of the tachyon (remembering that the energy dimension of τ is −1): To have exact chiral symmetry one must find solutions for which m q = 0, and achieving this, is one of the technically most demanding tasks of this model (for details, see Appendix B). The behavior of the potentials at large values of the fields λ and τ is determined by requirements of a confining spectrum and breaking of the chiral symmetry in the deep infrared [18,20,21,32]. To fix the last remaining parameter we chooseμ = − 1 2 . This choice, according to [33], leads to regular thermodynamics at zero chemical potential.
With these definitions, the numerical results in this paper are given for the potentials (2.6)- (2.8) and for x f = 1 case only. Of course the above choice for the potentials and κ is not unique but other possibilities exist as discussed in [32,34]. The definitions presented above are taken in this paper to provide for a benchmark study of this model, and focused analyses of other potentials and other values of x f , in particular approaching the conformal region at x f ≈ 4, are left for future studies.
As a final remark here, we emphasize that the duality between classical gravity and field theory can be derived in the string theory framework only in the strong coupling limit. In our case, the matching to the scheme independent perturbative results in the weak coupling limit has to be regarded as a model assumption, to be judged on the basis of its consequences. Among these, an immediate and important one is that one can describe thermodynamics up to arbitrarily high T and µ and identify solid known behaviors. Actually it is quite nontrivial that this matching can be carried out and the correct running of the quark mass and the condensate implemented using the DBI action. The model is then an effective theory extending weakly coupled results at large T, µ to the strongly coupled domain.

Φ equation of motion
The fermionic part of the action, given by depends only onΦ so that Φ is a cyclic coordinate. Since both L f andΦ have energy dimension 2, we have a dimensionless constant of integrationn: From this one solveṡ where we have also introduced the dimensionless density factor The factor K defined above contains the density effects in this holographic model and will appear repeatedly in what follows. After the bulk fields λ and τ have been determined from their equations of motion and the Einstein's equations, Φ(z) can be computed by integrating Eq. (2.16): with the constraint that the field Φ vanishes at the horizon z = z h , from which µ is determined.

Equations of motion for other bulk fields
Using the previous results for Φ, differential equations for b, λ, f, τ can be derived. They are for b(z) (2.22) and for τ (z) For λ we also have the first order equation (2.24) It turns out useful to define the quantity Using this in τ = 0 case the equations can be written in more compact form as follows: First we have from the above definition Treating this as a function of λ and b, the three remaining equations of motion are then The energy unit of the solutions is determined by fixing the small-z behavior of the dilaton to the perturbatively known field theory behavior, i.e. b 0 λ(z) = −1/ ln(Λ 0 z) with Λ 0 = 1. To do this accurately enough, one has to go to extremely small values of z and it is better to use ln(z) or actually ln b as the coordinate. For numerics we thus write the equations in the A = ln b basis changing z to A via the relation Then we have for q(A), primes denoting derivatives with respect to A, The remaining equations of motion are for λ(A) and for τ (A) This has the formally notable consequence that the A-equations are not autonomous; there is explicit A dependence. The consequence of this will become explicit when we consider the scaling properties of the solutions in the following section; see Eq. (2.42). The equation (2.22) for f can be integrated once: using (2.16). Then f (z) is obtained by one more integration, with integration constants determined by f (0) = 1, f (z h ) = 0. Actually we are most interested in the charged black hole temperature, for which one obtains . (2.37)

Scaling properties of equations of motion
Numerical solutions have to be transformed to the required standard form by using scaling properties of the equations. A thorough discussion is given in Appendix B, and we summarize the main points in the following. The quantities which are not mentioned will remain unchanged and all bulk fields are taken to be either functions of z or A.
For the z equations (2.20)-(2.23) one performs the following scalings: • The boundary value of f (z) must be set to 1 so that the boundary metric is pure AdS, f (0) = 1. This is achieved by scaling In order to keep b 2 /f and K(z) in (2.17) invariant, this requires that further Note that also the integration constantn is scaled.
• The unit of energy can be changed by z → Λz, together with which leave the equations of motion invariant.
For the A equations (2.31)-(2.34) the corresponding scalings are: • Scaling of f to f 0 = f (∞) = 1 requires that q 2 /f be constant, so that Note that the density factor K(A) in (2.35) is not affected by this scaling.
• The scaling corresponding to z → Λz is The invariance of the density factor K(A) and Eq. (2.33) then demand that

Constant Scalar Solutions and IR Stability
To gain intuition on what to expect at zero temperature and finite chemical potential we now consider some special solutions of the equations of motion derived in Sec. 2. We need to determine the fixed point solutions with translational symmetry since flows between different such solutions categorize the various RG flows of the boundary theory.
In general the fixed point solutions with translational symmetry are AdS p solutions either with fixed scalars or hyperscaling violating solutions when the scalars run off to infinity, [36,39]. We have not found hyperscaling violating asymptotics in this theory. The other remaining scaling solutions must then have constant scalars. These solutions will be the non-linear generalization of AdS 5 Reissner-Nordström black hole (the so-called DBI black hole), and solutions with scaling AdS 2 regions in the IR, at extremality.
To search for these, we turn to the equations of motion and make the following replacements: where a zero sub-or superscript indicates the constant value of the appropriate quantity in the fixed point and The two classes of solutions are distinguished by whether the scale factor A is constant or not. If it is constant we obtain AdS 2 type solutions while if it is non-trivial we obtain AdS 5 type solutions.

AdS 5 and the DBI Black-Hole Solution
For constant scalars many of the equations become quite simple, and often can be decoupled. For example, the equation governing the warp factor, A = log b is just which has two independent solutions, A(z) = − log z or A constant. The first matches the AdS 5 result in these coordinates. This is the solution one anticipates as a UV fixed point in the dual theory. It will turn out to be the charged DBI black hole, which becomes the AdS 5 Reissner-Nordström solution in the limit of small gauge coupling. We can systematically insert this solution into the remaining equations of motion.
where 2 F 1 is the hypergeometric function. For consistency, the blackening function must be compatible with the constraint equation (2.24), given by (3.8) Note that this equation effectively governs the constant term in (3.6), or equivalently the near boundary value of the blackening function. To leading order inn the solution consistent with the above constraint is 4 which is exactly the form one would anticipate in AdS 5 with charged branes.
The remaining equations of motion, those for the dilaton and tachyon, contain algebraic constraints for various parameters of the theory. Specifically, the dilaton equation implies while the tachyon equation needs either in order to be satisfied. These constraints have a simple interpretation. The set of equations are simply the requirement that all the potentials are extremized at the appropriate value of (λ 0 , τ 0 ). Evidently, the same must be true for the gauge kinetic function w(λ, τ ). For the V-QCD potentials of interest, specifically those from Section 2, it turns out that the extremization condition in 3.10 can never be realized and the only possibility is the vanishing dilaton, λ 0 = 0. The gauge kinetic function w and the flavor potential V f are of the general form of Eq. (2.4) and Eq. (2.6), respectively: so ∂ τ w = 0 and the tachyon constraint reduces to (3.14) Therefore, the flavor potential is extremized in the τ direction for either τ 0 = 0 or τ 0 = ∞. Moreover, it can be explicitly checked that when the dilaton is zero there is no location in the parameter space (x f , τ 0 ) for which the combination V 0 f κ 0 diverges. Accordingly, one finds that in this V-QCD setup DBI black hole solutions exist at all x f so long as λ 0 = 0 and τ 0 = 0 or ∞.

AdS 2 Solution
There exists another simple solution to the constant scalar warp factor equation of motion (3.3). This is the constant solution A = A 0 . In this case, the Maxwell equation is satisfied by a constant electric field of the form giving a potential 5 which is the correct form for a gauge field in AdS 2 . The equation for the blackening function is and has the general solution The AdS 2 solution is simply the one in which C 1 = C 2 = 0, and we identify the AdS 2 radius, L 2 , as All the rest of the equations simply give constraints that determine when this solution can be realized. The "zero energy" constraint says that while the dilaton equation of motion requires and the tachyon equation forces In the following section we will investigate these constraints in more detail, to determine whether or not they can be realized in V-QCD models of interest.

A closer look at the AdS 2 solution
We can summarize the AdS 2 requirements succinctly by recalling the definition of the effective potential, Eq. (2.25), in the language of this section: in which case the AdS 2 constraints are simply The zero energy constraint V 0 eff = 0 shows that the volume form on the R 3 factor is just so one can think of this condition as an expression describing the size of the R 3 , as determined by the values of the potentials at the fixed point. For the class of potentials of immediate interest, this relationship fixes the volume of R 3 in terms of (λ 0 , τ 0 , x f ,n). Solving the zero energy constraint forn allows one to rewrite the extremization conditions like Note that these expressions depend only on x f , λ 0 , and τ 0 , and that the notation asks one to differentiate the potentials first, then evaluate the result at the constant scalar solution.
Finding simultaneous solutions to these equations provides the parameter space on a two-parameter plane in which the AdS 2 solution can be realized. For the class of potentials used in V-QCD (3.13), this constraint is again trivially satisfied for τ 0 = 0 or τ 0 = ∞. For vanishing τ 0 , it is easy to find solutions to the constraint numerically for the V-QCD potentials in Section 2. They appear in figure 2. Interestingly, there is a region at low x f where there are two solutions for constant (positive) dilaton. This behavior may be an artifact of the parametrization of the potential w(λ). The second fixed point is not expected, but we also find that it plays no role in the phase diagram.
In the case of the divergent tachyon, τ 0 = ∞ it is clear that V 0 f = 0. One can carry out the same analysis as in the τ 0 = 0 case to search for allowed AdS 2 solutions in V-QCD, carefully navigating the somewhat subtle limits implied by this solution. For  Allowed (n, λ 0 ) (left) and (x, λ 0 ) (right) values for the AdS 2 solution with vanishing tachyon. At left, the black dots mark the location of x f = 1 along each branch, and correspond precisely to the AdS 2 solutions found numerically and shown in figure 7. The black dashed line marks the Banks-Zaks limit at x f = 11/2. From the right plot we find that when x f 2.865 the constant dilaton solution becomes negative and is thus excluded as a fixed point candidate.
finiten but vanishing V 0 f one finds that a divergent tachyon implies an electric field (3.15) and AdS 2 radius of the form The extremization condition (3.26) becomes and the numerical results for the potentials in Section 2 are shown in figure 3. As before there are two branches of solutions-the smaller of which terminates at some finite value of x f within the Banks-Zaks limit at x f = 11/2.

Stability of the AdS 2 Region
The AdS 2 solutions can be a priori endpoints or starting points of RG flows. To determine exactly what happens we must do a scaling analysis of the perturbations around them.
We perturb the background AdS 2 metric like In this background, the IR is approached as r → 0 while the UV as r → ∞. Here L 2 is the AdS 2 radius as given by (3.19), C 0 controls the volume of the R 3 factor, and the other constants parametrize the fluctuations in the obvious way. Without loss of generality, we set C 0 = 1 in what follows. All fluctuation amplitudes are taken to be small.
The background fields are perturbed as well, The program is to insert these perturbation Ansätze into the equations of motion, linearize the equations about the fluctuations, and subsequently determine the scaling exponents and the fluctuation amplitudes that describe a given perturbation.
Operationally, one first sets all the fluctuations above the background proportional to the same power, which is to say The linearized fluctuation equations then reduce to a coupled set of homogeneous linear equations in the amplitudes of the fluctuations Importantly, the radial Anstze under investigation leaves a residual gauge freedom related to reparametrizations of r. Practically, this means that fixing B 1 constitutes a gauge choice, and the linear system consists of 5 independent equations. Requiring that the system have a non-trivial solution is equivalent to requiring that the determinant of the matrix of coefficients, M vanish for all r.
In this case, one finds that the determinant is of the form which vanishes for α * = {0, −2, −1, 1} and for the α = α * such that g(α * , λ 0 , E) = 0. The former correspond to "universal" modes, while the latter are "non-universal" in the sense that they depend on the details of the various potentials. Of the universal modes, we find that there are two types of IR relevant (α < 0) modes in the fluctuation spectrum, with exponents α * = {−2, −1}. That they correspond to relevant operators in the IR is clear from the fact that when α * < 0 these modes grow as r → 0.
To better understand the non-universal modes, it is useful to write them in terms of the effective potential (3.23). Note thatn can be easily related to the boundary value of the electric field (E in this section) via (3.15). The effective potential also turns out to govern the properties of two of the four non-universal exponents, while the other two are The superscripts signify the fact that these modes correspond to perturbations of the appropriate scalars as we will see below.
These exponents have a few noteworthy features. First, all of the exponentsuniversal or not -can be pairwise summed to give α + + α − = −1, which is the correct structure for modes in AdS 2 , in these coordinates. Moreover, we see from (3.39) that there is a BF-like bound signaling the onset of an instability when ∂ 2 . For the V-QCD potentials employed for numerical studies, these non-universal exponents are plotted in figures 4, 5 and 6 as functions of x f for both branches of the AdS 2 fixed point. Evidently, while the BF-like bound is never exceeded in the fluctuations corresponding to α λ , the fluctuation characterized by α τ realizes an analogous instability around x f ∼ 2.4 in the vanishing τ 0 case. When the tachyon is divergent, the equations of motion require E 2 w 2 0 = 1 and thus α τ ± saturates to {−1, 0}. It will turn out that the fluctuations described by α τ are appropriately named, as they correspond to fluctuations of the tachyon alone.
The full description of the perturbation is given by the exponent α * , which contains information about the dimension of the dual IR operator, and the amplitudes of the  various modes that are activated by this fluctuation. The following cases are pertinent for the two conjugate solutions: • If the operator is UV relevant then both perturbations vanish in the UV boundary.
• If the operator is IR relevant then both perturbations blow-up in the IR regime.
• If the operator is UV irrelevant then one perturbation vanishes and one blows up in the UV boundary.
• If the operator is IR irrelevant then one perturbation vanishes and one blows up in the IR regime.
The amplitudes are easily obtained by solving the linear system provided by a given α * , and in general depend on one undetermined (but non-vanishing) amplitude and a choice of radial gauge which can be fixed via B 1 . The results are listed in Appendix A.
We conclude this section by assessing the RG stability of AdS 2 solutions. The one that appears at small values of λ, denoted by a blue line in figure 2 has dilaton and tachyon perturbations that render it IR unstable. This explains the fact that it plays no role in the phase diagram we describe in this work.
The other AdS 2 solution that corresponds to the purple line in figure 2 has dilaton and other perturbations that are IR irrelevant but the tachyon perturbation is IR relevant in the non-tachyonic black-holes. This is as expected as we need to tune m q = 0 to reach this solution in the IR. Once we turn on m q = 0 we will avoid it and end up in the tachyonic black hole. On the other hand in the tachyonic case, the dilaton perturbation is IR irrelevant and the tachyon one is marginal. However it does not correspond to an extra parameter in the theory as τ = ∞ is a singular point in field space.

Numerical solution 4.1 Preliminaries
The equations of motion admit two types of solutions at finite temperature and chemical potential, which we call black hole and thermal gas solutions. The thermal gas solutions have no horizon in the IR. In this case the temperature is identified as the inverse of the length of the compactified time coordinate, while Φ = const. = µ. The blackening factor is trivial, f ≡ 1, and the z-dependence of the other fields is exactly the same as for the solutions at T = 0 = µ, which were constructed in [32]. When 0 < x f < x c , the dominant vacuum was found to have a nonzero tachyon field (and therefore broken chiral symmetry). The thermodynamics of the corresponding thermal gas solution is trivial: the pressure is independent of T and µ and will be normalized to zero here. Likewise, the condensate, which signals chiral symmetry breaking, will be nonzero but T independent.
The nontrivial task on which we concentrate in this article is the construction of the black hole solutions. The equations we have to solve numerically are the Einstein's equations (2.31) and (2.33), the equations of motion for λ, equation (2.32), and the equation of motion for τ , equation (2.34). Their solution forn = 0 has been discussed in detail in [33]. The numerical solving with given initial conditions as such is very simple using NDSolve of Mathematica. The main issue is the correct initialization and subsequent processing of the solutions via the scalings described in section 2.4.
An important general feature is that there are two types of black hole solutions: • The solutions with τ = 0 which describe the hot and dense matter in a chirally symmetric phase; these are expected to dominate the free energy at large T or µ.
• The solutions with τ (A) = 0. These will describe a chirally broken phase, expected to dominate at small T or µ. These solutions are parametrized by the value of the quark mass Since we are interested in solutions with exact chiral symmetry, we need to restrict to m q = 0. This is a technically very demanding task (see Appendix C) and necessitates going to very small values of z ≈ e −A , up to A ∼ hundreds. This is one of the reasons for using A as a coordinate. The details of the numerical solution and the associated scaling properties are discussed in detail in Appendix B.
In the numerical computations we choose the unit of number density so that L A = 1. In section F we shall actually fit that L A ≈ 0.97.

Initialization: expansion around horizon
For thermodynamics one needs solutions with a black hole. To generate them numerically, one has to start the integration at the horizon, which we place at A = A h such that f (A h ) = 0. Because of the singularities due to the 1/f terms in Eqs. (2.31)-(2.34) one cannot start the integration precisely at the horizon. Instead, one first writes the values of the fields at a small distance from the horizon by expanding in as which are then inserted to the equations of motion. Here and in the following the subscript h denotes quantities evaluated at the horizon. Then one expands in and demands that the divergences and the constant term vanish. Note that the input here is that in Out of the leading terms in (4.2)-(4.5) one can choose f h = +1 as the magnitude of f (A) will anyway be fixed by the scaling (2.41) to the boundary value f (A → ∞) = 1.
The dilaton value at the horizon λ h will remain as a parameter, closely associated with temperature. The second parameter, closely related to the chemical potential, isn. However, in the numerics it turns out to be more practical to use instead which is invariant in the scaling of (2.42) and (2.43). The tachyon value at the horizon will be fixed by the quark mass, τ h = τ h (λ h ,ñ; m q ). Including the terms up to O( ) for q and up to O( 2 ) for the other fields in (4.2)-(4.5) is sufficient to ensure that the values of these parameters in the resulting numerical solution match with their input values to a high precision. The remaining first-order derivative terms will be fixed by demanding that the 1/A (i.e. 1/ ) singularities cancel. Canceling the divergent 1/A term of (2.31) gives with the understanding that the potentials V g , V f , and w are evaluated at the horizon. Canceling the divergent 1/A term of the λ equation (2.32) gives and canceling the 1/A term of the τ equation gives In (4.7) and (4.8) we again have the important quantity, The remaining expressions are too complicated to be reproduced here but can be found in [38]. From the algebraic derivation of the initial conditions to the numerical integration of the system of differential equations (2.31)-(2.34), we treat the whole problem in Mathematica.

Observables
It is thus easy to produce some numerical solutions for the functions q(A), λ(A), f (A), and τ (A) with Mathematica, given λ h andñ, but an essential and nontrivial part of the numerical work is to transform the solutions to a standard form satisfying in z coordinates f (0) = 1 and that the scale of the UV expansions equals one (see Appendix B).
In A coordinates these conditions become The former is implemented by scaling f as in (2.41), the latter by scaling A as in (2.42).
To achieve this one determines the scaling factor Λ(λ h ,ñ) so that the asymptotic limit (4.12) holds. We start from a numerical solution having A h = 0, then according to (2.42) the value of b at the horizon in the scaled solutions is simply given in terms of the scaling factor by b h = exp(A h )/Λ = 1/Λ. From the standard configurations so obtained one then computes the temperature as the black hole temperature and the chemical potential using (2.16) and (2.19), otherwise the configurations as such are not of interest for this calculation. The procedure is described in detail in Appendix B. Summarising, from the numerical integration of equations of motion, for given (λ h ,ñ), one obtains the following quantities: (4.13) From these we obtain the entropy density using the basic formula (4.14) To obtain the 4d physical quark number density note first that, when deriving the Φ equation of motion from the fermionic part of the action, one has, for solutions of equations of motion, At z h one has to keep the value Φ(z h ) fixed to zero so that δΦ(z h ) = 0 and the upper limit does not contribute. Since S = −Ω/T and δΦ = dµ the fermionic contribution given by the above integral is the n dµ term in the free energy, and therefore the correct normalization of n is where we used the definition (4.6). This expression also gives a physical interpretation of the parameterñ of the integration of the equations of motion: Next we discuss what values of (λ h ,ñ) are possible and how the pressure is integrated from dp = s dT + n dµ. To compute the pressure we have to integrate over T and µ and these one-dimensional integrals are most simply carried out by converting them into integrals over λ h at fixedñ or vice versa, see Section 5.1.

Physical region in the λ h ,ñ plane
Forñ = 0 one found (see, e.g., [33], Fig. 7) that chirally symmetric solutions, i.e. the ones with zero tachyon, existed only for 0 < λ h < λ * , with λ * given by the extremum of the effective potential in (4.10), and chirally broken solutions with nonzero tachyon existed only for λ h > λ end with 0 < λ end < λ * . This is in harmony with the expectation that large T and chiral symmetry are associated with small coupling, small λ h , and strong coupling leads to chirality breaking. The introduction ofñ extends this pattern in an interesting and subtle way, exhibited in Fig. 7.
For smallñ the above pattern remains unchanged, only the curves λ end (ñ) and λ * (ñ) slowly decrease. Also thermodynamically the situation is only slightly modified: as we shall soon see along this curve the symmetric and broken phase pressures are equal and there is a continuous chiral phase transition. At aboutñ = 8 the λ end (ñ) curve has a discontinuity and breaks into two branches, λ end ≡ λ χb and λ χs so that λ χb is always the lower limit of the broken phase while λ χs exists within the symmetric phase. The branching signals the onset of a first order chiral transition; the broken phase at λ χb can be in thermal equilibrium at λ χs but so that the density jumps when crossing the phase boundary.
The computation of the curves in Fig. 7 is mostly numerical, but parts of the boundary of the symmetric phase can be found analytically. Note first that (4.7) Tachyonic chiral symmetry breaking solutions exist only above the blue curve λ end (ñ). This curve has a discontinuity atñ =ñ cr ≈ 8 at which it breaks into two branches, λ end ≡ λ χb and λ χs . Belowñ cr the symmetric and broken phases are in thermal equilibrium along λ end , aboveñ cr the states on λ χb and λ χs are in equilibrium. The deconfining transition line further in the IR is plotted in Fig. 22. The dashed lines are V eff = 0 and V eff (λ h ) = 0 at τ = 0 (see (4.10)).
implies that This with L A = 1 is the curve V eff = 0 in Fig. 7. Further, due to the interpretation β(λ) = λ (A) one usually expects that λ(A) monotonically decreases from its value λ h = λ(A h ) towards λ(A = ∞) = 0, and in particular that λ h < 0. From (4.8) this would imply However, deep in the IR the interpretation of λ (A) as a negative beta function need not be valid and solutions with signs opposite to those in (4.20) are also possible. This is confirmed by numerical computation and the real boundary is given by finding where T = 0 or where the scale factor Λ(λ h ,ñ) diverges. Requiring that both V eff and V eff vanish has two solutions marked AdS 2 , since actually the geometry at these points is asymptotically AdS 2 × R 3 in the IR. The lower AdS 2 point disappears at larger x f .
The boundary of the broken phase marked λ end in Fig. 7 is discussed in some detail in Appendix D. It is a lower limit for possible values of λ h . Forñ there is an upper limit, but there is no upper limit for λ h .
The significance of various parts of the physical region is also described by plotting curves of constant T and µ as in Fig. 8. Actually we show there the result only for the discrete values ofñ used in the pressure integration. Particularly interesting is the behavior of the µ = constant curves. Extrapolating them one sees that clearly asymptotically µ = 0 in the upper part of the T = 0 curve. States here have T = µ = 0 and thus represent vacuum. In the vertical part of the T = 0 curve correspondingly µ = ∞. This is also some special state. All the µ = constant curves end at the AdS 2 point, where thus all the exactly T = 0, µ finite symmetric phase thermodynamics resides.
A third important quantity is the dimensionless scale factor Λ(λ h ,ñ). It varies a lot as a function of λ h . The main part of the variation is contained in Λ(λ h , 0) which at the boundaries of the phase space, λ h → 0 and λ h → λ * (0), can be fitted by

Constant parameter curves
A straightforward way to generate the data necessary for solving the thermodynamics would be to compute black hole solutions in both the symmetric and non-symmetric branch of the solutions on a sufficiently dense lattice in the physical region of the (ñ, λ h ) -plane. In order to carry out the pressure integrals without accumulating large cumulative errors, a reasonably accurate continuum interpolation of the observables is needed. Since the observables as a function of (ñ, λ h ) are mostly smooth but have lines of zeroes and divergences, detailed in the following section, it would be necessary to use at least a somewhat sophisticated interpolation algorithm with an adaptive local grid size in two dimensions, or alternatively a very large amount of computing power for the brute force approach of simply a very dense uniform lattice. However, we have been able to avoid constructing a full 2D interpolation of the solutions by considering a grid of 1D interpolations, for which well-established adaptive algorithms are readily available. The two primary interpolations are curves withñ as a Figure 9. Constant values of the scale factor Λ, normalised to its value atñ = 0, on theñ, λ h plane for tachyonless solutions for the discrete values ofñ used in the pressure integration, other boundary curves as in Fig. 7. Curves of constant s and n can be inferred from this (see text).
constant and λ h as the variable, and those with λ h as a constant andñ as the variable. We compute the interpolations for a number of values ofñ and a number of values of λ h . Figures 18 and 21 in appendix D show images of these curves of the (µ, T ) -plane.We can then compute the pressure integrals along each of these, for both the symmetric and non-symmetric branches, fixing the constants as described in the next section.
In the next section, we also describe some additional one-parameter curves which have a thermodynamically interesting role that we have computed.
At least in the specific case handled in this paper, the constant parameter curve method has allowed us to extract all the thermodynamic features of interest without resorting to full 2D interpolation. However, in the case of a transition between two regions of the same branch of solutions, such as happens at small x f for some of the potentials explored in [33], a complete 2D interpolation may become necessary to extract the phase transition line.

Computation of pressure
According to the holographic dictionary, the pressure can in principle be computed by evaluating the on-shell action. However, this is numerically very challenging in this kind of model with corrections decaying only logarithmically near the boundary. Therefore, we use instead the usual thermodynamic formulas.
We first review how the pressure is computed by integrating s(T ) = p (T ) forñ = 0 since this is how the constant of integration is fixed in [33] and will be fixed here, too. One has, see Fig. 10, where the subscripts b and s denotes quantities in chirally broken and symmetric phases, respectively. What matters for the phase structure is the difference of the integration constants p b (∞) and p s (λ * ). This is simply fixed by requiring that pressure be the same for the two phases at λ h = λ end [33]. The outcome is plotted in Fig. 10. At this temperature there is a second order (both p and s ∼ p are continuous) chiral phase transition. The broken phase pressure vanishes at λ h = 3.19 at the temperature T h = T b (3.19, 0) = 0.14. This is the deconfining transition. At higher λ h or smaller T the dominant phase is the thermal gas phase with vanishing thermal pressure. For quantitative correctly normalised results one will need both the energy unit Λ 0 , which is implicit in formulas involving T and µ, and the constant 4G 5 . The former is fitted by the value of the critical temperature T χ (0) = 0.148Λ 0 . Taking T χ = 0.15 GeV, we fix Λ 0 = 1 GeV.
For 4G 5 normalisation to the T 4 Stefan-Boltzmann term at T → ∞ gives [33], see Eq. (F.14), In general, we wish to obtain the pressure by integrating dp = sdT + ndµ. All the quantities on the RHS are numerically known as functions of λ h ,ñ, see Appendix C. Note that we can write dp as where all quantities are functions of λ h ,ñ.
The differential dp can now be integrated either over curves of constantñ from λ b to λ t , 6) or over curves of constant λ h fromñ b toñ t , To test the path dependence of the integral, one can choose an arbitrary rectangle within the physical region in Fig. 7 for either of the phases. This is mapped to a foursided region on the grid in Fig. 21. One now integrates numerically around it using Eqs. (5.6) and (5.7) and checks whether the integral is zero. This is indeed what we find to a great accuracy.
This proof of path independence is a very impressive confirmation of the validity of our numerical computations. All the quantities included in (5.6) and (5.7) are the results of lengthy numerical solutions of Einstein's equations and it is striking to see that when they are put together as above, the outcome is path independent to a very good numerical precision.
The pressures of the two phases p s (λ h ,ñ) and p b (λ h ,ñ) can now be computed by fixing the relative integration constant by demanding that p s (λ end , 0) = p b (λ end , 0) and integrating to the point (λ h ,ñ) along any convenient path. These can trivially be converted to p s (T, µ) and p b (T, µ). Three-dimensional plots of pressure vs T, µ are numerically rather noisy and we shall focus on the main question: phase structure and phase transition lines.

Phase structure
We have discussed thoroughly the chirally symmetric and broken phases with pressures p s (T, µ) and p b (T, µ). Furthermore, as a model for the low T system we shall use the thermal gas phase, for which the metric Ansatz is like that in (2.3) but with f (z) = 1 and T is introduced by compactifying the imaginary time region, otherwise the equations of motion are as before. Note that from f = 1 and the equation of motion (2.22) it follows that one must haveñ = 0 so that also n = 0. Thus alsoΦ = 0 so that Φ = µ is constant. The property n = 0 for a low T chirally broken phase is consistent with the fact that this model contains no baryons in the spectrum of singlet states.
The pressure of the thermal gas phase is p low = 0. We identify this with the hadron gas phase of the field theory. This is well justified in the case of pure SU(N c ) gauge theory, for which the pressure of the plasma phase is ∼ N 2 c . In the case here (V-QCD) the degrees of freedom in the low temperature phase are the N 2 f Goldstone bosons of the spontaneously broken chiral symmetry, while the high energy degrees of freedom are the deconfined partons, 2N c + 7/2N f N c . The ratio of the number of degrees of freedom at low and high temperature is then x 2 f /(2 + 7/2x f ) which at x f = 1 is 0.18, and we expect that taking p low = 0 provides still a useful guide towards the location of the deconfinement phase boundary. However, as x f increases, the uncertainty associated with this approximation grows. At x f 4 the ratio becomes unity signalling the transition to a different vacuum phase as one enters the conformal window.
The stable phase has the smallest free energy or, equivalently, the largest pressure. For phase equilibrium one needs both kinetic, thermal and chemical equilibrium, i.e., the same pressure, temperature and chemical potential for the two phases. The outcome of the analysis has already been shown in Fig. 1.
Consider first the most reliable prediction of the model: the chiral transition at T χ (µ). At µ = 0 ofñ = 0 this took place at the point λ end in Fig. 10. This pattern continues whenñ is increased up to the valueñ =ñ cr ≈ 8, the chiral transition takes place along λ end (ñ) much as atñ = 0.
Atñ =ñ cr ≈ 8, called the critical point in Fig. 7, something novel happens. Note that this is also the point where λ end and the line V eff = 0 intersect. Suddenly the numerically determined value λ end , the lower limit of values of λ h for which tachyonic solutions with m q = 0 exist, has a derivative discontinuity. Whenñ is increased, it even leaves the region of the symmetric solutions. The question is what happens to the phase transition line.
At largeñ it is clear that the transition must be between disconnected points of the (ñ, λ h ) -plane, since the symmetric and non-symmetric phases do not overlap there. A brute force search of the transition would entail taking a point in both branches, say (ñ s , λ s ) in the symmetric phase and (ñ b , λ b ) in the broken phase, and requiring There are four unknowns and three equations, which could be expected to yield a oneparameter phase transition curve. This is somewhat non-trivial to implement in our numerical scheme, where we have up to this point been able to avoid generating full 2D interpolations in the space of solutions. After a visual inspection of the λ h -constant andñ -constant curves in the (T, µ, p) -space, it however seems that the λ end -curve comes very close to just touching the symmetric phase. This leads to the hypothesis that the λ end -curve continues to be the phase transition curve even at largeñ.
Testing this hypothesis requires computing the solutions and pressure along the λ end -curve, finding a curve of equal temperature and chemical potential from the symmetric phase and verifying that the pressures are also the same. Parametrizing λ end withñ b , the equations to be solved are Now we can take a fixed value ofñ s for which a curve as a function of λ h has been computed. This reduces the problem to a 2D problem of findingñ b and λ s such that (5.9) holds, which is numerically tractable as we are now indexing solutions along fixed curves where an interpolation of the observables has been pre-computed. For each fixed value ofñ s we obtain one point on the line of transitions, and using a large number of n s values gives an accurate representation of the whole transition line. Once this has been done, we can compare the pressures between the two phases. They are equal to within numerical accuracy, confirming our above hypothesis. We denote the curve of these equilibrium points in the symmetric phase by λ χs .
The transition is then between points on two curves, parametrized as λ χs (ñ s ) = (λ s (ñ s ),ñ s ) and λ χb (ñ s ) = (λ end (ñ s ),ñ b (ñ s )). (5.10) Symmetric phase states along λ χs are in equilibrium with those on λ χb for same values of the parameterñ s . The splitting of the line simply implies that the transition is of first order. This is again a striking manifestation of the ability of the holographic method to describe subtle effects. The equations (5.9) give directly the transition temperature T χ (µ). It is also of interest to plot the temperature as a function of n, which splits the phase transition line in two with a mixed phase in between. Theñ components in equation (5.10) express the change inñ; the change in physical density is The outcome of the computation is shown in Fig. 11 which plots both T χ (µ) and T χ (n). The discontinuity in λ end atñ = 8.0 corresponds to the critical values µ cr = 0.34, T cr = 0.0804, n cr = 0.0103. (5.12) Beyond these critical values in the 1st order region there is a density jump and the line splits in two leaving in between a mixed phase. Concretely, choosingñ s = 12.2, very close to the AdS 2 point in Fig. 7, one finds that λ s = 1.243 andñ b = 14.67, the physical densities are n s = 0.0180, n b = 0.0172 and the common pressure value is 0.0015. A symmetric phase point very close to the AdS 2 point in Fig. 7 is thus in equilibrium with a broken phase point within the bend of the curve λ χb . Below T χ (µ) the chirally broken phase is the stable one. Its pressure is positive but starts decreasing when T is further lowered, λ h increased, just as happens in Fig. 10 at µ = 0. Computing the pressure of the broken phase for arbitrary values ofñ, one finds that it vanishes for the values of λ h plotted in Fig. 22. The corresponding temperature T h (µ) is plotted in Fig. 11, see also Fig. 1. We interpret this as a first order deconfining transition between the chirally broken phase and the zero-pressure low T thermal gas phase. The transition temperature T h (µ) decreases monotonically with increasing µ; our numerical accuracy does not permit to make definite statements about the limit T → 0. Note that this corresponds to very large values of λ h , see Fig. 23.
One observes that at µ = 0 the temperatures T χ (0) and T h (0) = 0.95T χ (0) are very close to each other. For reference, one may note that a very similar situation with T h (0) = 0.94T χ (0) was observed in [40] in a completely different Schwinger-Dyson equation model for QCD thermodynamics. There the conclusion was that the chiral and deconfinement transitions probably coincide. These behaviors are most transparently understood on the basis of underlying exact and approximate symmetries and related order parameters [41,42].
From the computed pressure, we can determine the interaction measure − 3p = T s + µ n − 4p = (T ∂ T + µ∂ µ − 4)p, which is shown close to the chiral transition region in the symmetric phase in Fig. 12. Consider the curve for µ = 0, the structure of which is described in Fig. 9 of [33]. The analogous curve for QCD is plotted, e.g., in Fig. 3 of [43]. The V-QCD curve plotted in Fig. 12 starts by decreasing above T = T χ (0) but then changes direction and passes through a maximum at T ∼ 4T χ (0) with a QCD-like decay above that. This large T maximum can be interpreted [33] as a crossover transition. When x f is increased into the conformal region at x f > x c ≈ 4, this crossover is the only structure in p/T 4 which remains. It is now apparent that increasing µ does not change this overall pattern qualitatively. In particular, the large T decay is independent of the chemical potential.

Order of transition
The chiral transition was above numerically observed to be of second order at µ = 0, for the potentials used here. There are also potentials which lead to a 1st order transition, as concretely shown in [33]. It is commonly accepted that the chiral QCD transition at N f ≥ 3 is of first order [44], even though this has not been conclusively established with lattice Monte Carlo computations, say, for N f = N c = 3, x f = 1 [45]. It is useful to see how our gauge/gravity duality model, valid, in principle, for N c 1 relates to the general effective theory arguments.
The order parameter for the effective theory of the QCD chiral transition is a complex . . , N f , x is the d = 3 dimensional spatial coordinate. The potential term in the action is To study the phase transition one should compute the effective potential of the theory.
In the 1-loop approximation this was carried out, for m = 0, in [46]. Much information can already be obtained from the beta functions of the couplings in d = 4− dimensions: if there is an infrared stable fixed point, zero of the beta function away from g 1 = g 2 = 0, the transition probably is of second order. If the couplings run to infinity, the transition is of first order. In the computation of [46], the color and hence the value of N c is hidden in the color contraction in qq . Opening up these color interactions in the 1-loop computation in full is an impossible task, but in the large N c limit a single tr is always one quark loop and thus suppressed by a factor 1/N c . Thus we expect that in the above effective potential g 1 ∼ 1/N 2 c and g 2 ∼ 1/N c . According to [46] the β-functions of the two couplings in (5.13) (scaled by a factor π 2 /3) have a fixed point at with the eigenvalues , − so that the fixed point is unstable, the flows are plotted in [46]. This is true also at large N c , N f , indicating a first order transition in this limit, too. However, one may argue that when N c = ∞, the term with g 1 in (5.13) should be entirely neglected. Then only the β function for g 2 remains and it has an infrared stable fixed point at This indicates a 2nd order transition. The two arguments are compatible if the latent heat of the 1st order transition is ∼ 1/N c . Another way to say this is that as N c → ∞ and g 1 /g 2 → 0, the Hermitian model becomes equivalent to the O(2N 2 f ) model that is known to have a second order phase transition.
One should also remember that the expansion cannot give any definite answer. A good example of this is another standard model transition, the electroweak phase transition. There the expansion method also leads to a first order transition [47] while a numerical computation leads to a first order transition for small Higgs masses, m H 75 GeV, while at larger Higgs masses there is only a cross over [48].

Quark number density
The quark number density as a function of µ is plotted in Fig. 13 for some values of T . As can be seen from Fig. 1, T = 0.1 crosses the phase diagram at a 2nd order transition and T = 0.05 at a 1st order one. The figure has the expected structure, the jump in n is relatively small at T χ and hardly visible at T h ; the transitions are in this sense weakly first order. When plotted for T > T χ (0), the fixed-T curves contain just the symmetric phase with monotonically increasing n(µ). At large µ the stable phase is always the symmetric phase (solid red) which can exist as a metastable phase (dotted red) even below the transition. Below the chiral transition the broken phase is stable (blue) and at the lowest µs the stable phase is the thermal gas phase with n = 0 (black, shown for T = 0.1). At the first order chiral transition there is a jump in n, as in Fig.11. There is also a numerically hardly visible jump between the broken and the thermal gas phases. Above T c = T χ (0) only the symmetric phase exists. Right: The quark number susceptibilities χ 2 (continuous) and χ 4 (dashed) per N 2 c including the normalisation factor 4/(45π). The limit of χ 2 /T 2 at large µ is 1 3 , that of χ 4 is 2/π 2 .
It is common to characterize the µ dependence by the susceptibilities at µ = 0: In the ideal gas limit and for L A = 1, These (per N 2 c ) are also plotted in Fig. 13, now also including the normalisation term (see Appendix F). One is approaching the ideal gas limit but rather slowly.
It may also be useful to express critical quantities in physical units. We have, e.g., taking N c = 3, from (5.12) a not totally unreasonable value.

Polyakov line
The basic difficulty in the study of thermodynamic deconfinement is that there is no symmetry and thus no order parameter associated with deconfinement. For N f = 0 the Polyakov line, trace of path ordered exponential of A 0 over the periodicity range 0, 1/T in imaginary time, signals breaking of Z(N c ) symmetry and separates low and high T phases. Even though it is not an order parameter for finite N f , it is a gauge invariant measurable observable, which in lattice Monte Carlo studies varies together with the chiral condensate qq .
Constructing the gravity dual of the Polyakov line is very complicated, but we can model it in a simple way following [49]. The idea is to start from a duality determination of a string tension in a thermal ensemble, interpret this as dF/dz = dF/dT × dT /dz, compute from here dF/dT , integrate F (T ) by choosing the integration constant by physical arguments and finally plotting L = exp(−F (T )/T ). Here one can start from the determination of the spatial string tension σ s [50], determined from Wilson loops with sides in spatial directions, Apart from the string tension all quantities here are known and the evaluation, using (4.7), gives . In the symmetric phase we shall, somewhat arbitrarily, fix the integration constant in the integration of F (T ) so that F s (T χ (µ)) = 0. Normalising L to 1 at large T (actually, due to limitations of numerics, at T = 10T c ) one obtains the red large T curves in Fig. 14 for µ = 0, 0.2, 0.4. For the broken phase we shall enforce continuity at T χ (µ) by demanding that also F b (T χ (µ)) = 0. This is reasonable in the 2nd order range, but in the 1st order range L should be discontinuous. However, in this simple model for the Polyakov line we do not have the means of computing the discontinuity. A discontinuity for µ = 0 at T h is also visible, a very small similar discontinuity exists also for µ = 0.2 which crosses T h in Fig. 1.
It would be very valuable to derive a theoretically better founded gravity dual for the Polyakov line. There is no order parameter for deconfinement but this operator anyway is and will be used in lattice Monte Carlo studies.

Sound speed
Now that we have two thermodynamic variables we have three second derivatives of p(T, µ). Out of the standard quantities C V and C p are complicated to compute in the present framework, but it so happens that the formula for the sound speed squared where all quantities are to be taken at fixedñ, can be directly evaluated. Rather fortunately, our method of computation makes it trivial to take into account the extra condition among the fluctuations of pressure and energy density in (5.21), they are to be taken at fixed n/s and, due to (4.17) this is just fixedñ. In particular, forñ = 0 the derivative should be taken in the direction of T , as is usually done, though in this direction the volume density of entropy varies. The most interesting region is that near the phase transitions. Fig. 15 shows c 2 s plotted vs T /T χ (0) (numerically T χ (0) = 0.1484) for fixed µ = 0.2, 0.4, 0.6, i.e. as one moves vertically in the T direction in Fig. 11. At very large T , outside the figure, c 2 s approaches the conformal value 1/3. For µ = 0.6 one does not cross any phase transition and c 2 s approaches some fixed value at T = 0. For µ = 0.4 one hits the chiral transition at T = 0.30T χ (0) and c 2 s jumps downwards and continues almost constant to T = 0. At µ = 0.2 the jump corresponding to the chiral transition takes place at a higher temperature. Note that it is consistent to have a jump in c 2 s even at a 2nd order transition since c 2 s is a second derivative. The jump would be largest at µ = 0 (not shown), where the broken phase only exists over a small range from T χ (0) = 0.1484 to T h (0) = 0.1417.
For the broken phase Fig. 15 shows c 2 s plotted vs T for fixedñ. The range of variation at eachñ corresponds to the region between the curves λ end and λ p b =0 in Fig. 22. Beyond that there is a brief supercooled region and the system becomes unstable when c 2 s goes to zero. According to the phase diagram in Fig. 1 the temperatures are less than 0.15, which is obtained forñ = 0. One sees that the numerical value of c 2 s is very small for the larger temperatures. The largest value for eachñ is always obtained when λ h is at λ end . When λ h grows further into the broken symmetry region, c 2 s decreases and ultimately reaches zero. Even here the conformal value 1/3 will be reached in the limit T = 0 only whenñ → 23.99. In Fig. 15 the curve forñ = 16 does not reach T = 0, only an extremely small value.

The T → 0 limit
The T → 0 limit is particularly interesting. The end point at T = 0 can be reached by putting f = 1 in the metric Ansatz [31] or by imposing T ∼ f (z h ) = 0 in the thermal Ansatz. The flow properties of the latter case were studied in Section 3. For full thermodynamics one needs a numerical computation of p(T → 0, µ) in both symmetric and broken phases. The computations in the broken tachyonic phase are not yet complete, but the essentials are as shown in Fig. 11 along the T = 0 axis.
At T = 0 a 1st order phase transition occurs at µ = 0.507. At this µ the symmetric and broken phases are in equilibrium. The symmetric phase equilibrium point will approach the AdS 2 point at λ h = 1.108,ñ = 12.295 in Fig. 7. According to Fig. 11 one has a T = 0 quantum phase transition at µ = 0.507 between a chirally symmetric phase at n s = 0.0184 to a chirally broken phase at n b = 0.0175. According to (4.17) this would correspond to T = 0 values s s = 0.0188, s b = 0.0092: entropy density jumps The numerical value µ = 0.507 for the T = 0 phase transition follows concretely from the broken side by plotting µ along λ end as in Fig. 23. It thus corresponds in the physical region to asymptotically large values of λ h , τ h and, depending on the potentials, also those ofñ. In the symmetric phase, Fig. 8 shows that µ = 0.507 plays no special role among all the constant µ curves approaching the AdS 2 T = 0 point. Note also that the curves marked T = 0 in Figs. 7 and 8 correspond to T arbitrarily close to 0, T = 0 is only at the AdS 2 points. How all values of µ in the full 5d computation are mapped to the single AdS 2 point is so far incompletely understood. We plan to return to this question as well as to the full construction of the chirally broken tachyonic T = 0 states in a later paper.

Conclusions
We have analysed in this paper a non-fine-tuned gauge/gravity duality model for hot and dense QCD in the limit of large number of colors and flavors. The model contains 5-dimensional gravity with AdS 5 symmetry on the boundary, a scalar dilaton for confinement and asymptotic freedom, a scalar tachyon for quark mass and condensate and the zeroth component of a bulk 4-vector for chemical potential and quark number density. The potentials of the model are constructed so that one obtains the correct QCD beta function and mass running in the weak coupling region and color confinement in the strong coupling region.
Tuning the quark mass to zero, the main result of this paper is the phase diagram and a description of dynamical chiral symmetry breaking when temperature or density is decreased. Chiral symmetry corresponds to solutions with vanishing tachyon, which automatically leads to the vanishing of both m q and the condensate. Spontaneous chiral symmetry breaking corresponds to solutions with non-zero tachyon, which are constructed such that m q = 0 but nevertheless the condensate is nonvanishing. By explicit calculation of the pressures of chirally symmetric and broken phases we find that the broken one dominates at small temperature and chemical potential, T < T χ (µ). The transition in between is of second order for µ µ cr = 0.34 and of first order for µ µ cr .
When T is further decreased below T χ (µ), the system ultimately goes to another phase at some T h (µ) with a non-zero tachyon but without black holes: the thermal gas phase. We use this as a model for the low T hadron phase. In lattice Monte Carlo simulations one normally finds that the chiral transition (as identified by variation of the quark condensate) and deconfinement transition (as identified by the Polyakov line or energy density discontinuities) coincide; there is effectively just one transition line between a quark-gluon plasma phase and a hadronic phase. Here we find that these lines are separated.
The numerical effort needed to obtain the results presented here is extensive and we have thus limited ourselves to a quantitative study of one set of potentials and the case N f = N c . With improved numerical techniques many further questions can be addressed. A set of potentials which describes all T = 0 QCD physics in quantitative detail has been identified [31,34]. Computing also its thermodynamics with good accuracy would make it possible to correlate zero and finite T properties reliably. For example, how does the requirement of linear Regge trajectories in particle spectra affect the thermodynamics? How does the quark condensate behave as a function of density? Further, the approach to the conformal limit at x f = x c ≈ 4 has been computed for µ = 0 in [33] and it will be interesting to study also the T, µ phase diagram in this limit.
There is a surprise in the phase diagram at low temperatures. Our analysis suggests that there is a new quantum critical regime with exotic properties at T = 0 which realizes the symmetries of the associated geometry, AdS 2 × R 3 . This exists both on the T = 0 segment of the chirality breaking plasma as well as the T = 0 line of the chirally symmetric plasma.
The presence of the AdS 2 × R 3 geometry in the holographic solution indicates that there is a scaling symmetry of the time direction which does not act in the spatial di-rections. Such symmetries have been called semilocal. This is an unexpected symmetry in a theory at finite density, but it is natural and generic in the holographic context [36] and appears even for simple black holes like the Reissner-Nordström black hole [37]. The physics in this critical regime is similar to that of a theory with zero speed of light: all spatial points decouple in the IR.
The local RG pattern of such AdS 2 solutions is fully compatible with the phase diagram we derived. It is an interesting question to determine physical implications of this scaling regime as well as its potential experimental signatures.

Acknowledgments
We thank Blaise Goutéraux, Misha Stephanov and Aleksi Vuorinen for discussions. This

A Fluctuation modes around the AdS 2 point
In this appendix we compute explicitly the amplitudes for the fluctuations of the AdS 2 region discussed in Sec. 3.4. The amplitudes are easily obtained by solving the linear system provided by a given α * , and in general depend on one undetermined (but nonvanishing) amplitude and a choice of radial gauge which can be fixed via B 1 . For the V-QCD AdS 2 fixed point, the fluctuations are as follows: This is a relevant perturbation that couples the metric and the gauge field like This mode corresponds to one sort of finite temperature perturbation to AdS 2 . It is the AdS 2 black hole studied in [37].
There exists another finite temperature perturbation, which is again characterized by a relevant mode. In this case the perturbation has Φ 1 , D 1 = 0 and B 1 = −D 1 and which corresponds to a shift in the chemical potential and a metric of the familiar form Again, this is a black hole in AdS 2 , related to the one obtained from the α * = −2 fluctuation by a radial coordinate transformation. More specifically, under the transformation the metric (A.4) becomes  When C 1 = 0, then the volume form on the R 3 changes by a factor of and when D 1 = 0 then one obtains a shift in the time coordinate These are the conjugate modes to the α * = −1 finite temperature fluctuation above.
The last universal mode is the irrelevant perturbation conjugate to the α * = −2 mode. This perturbation couples all of the fields except for the tachyon, and is partially responsible for driving the system away from the AdS 2 fixed point. The amplitudes are somewhat complicated, but take the form and That this mode interpolates between the IR and UV solutions is suggested by the fact that for λ 1 = 0 the spatial part of the metric acquires a non-trivial radial dependence as per (A.11). For the non-universal exponents, one finds a simple perturbation This is gauge equivalent to a mode consisting of only a tachyon fluctuation. Finally, there exists a somewhat more complicated perturbation where where α λ is given by (3.39). This is a perturbation which couples the dilaton fluctuations to the metric fluctuations, leaving the spatial part of the metric unchanged. Note that the above expressions for the perturbation amplitudes hold for generic values of the constant scalars λ 0 and τ 0 . In the special case of the divergent tachyon, many of the amplitudes simplify as in this case E 2 w 2 0 = 1. To wit, (A.11) becomes

B Numerical solution of the equations of motion B.1 Definitions
The purpose of this appendix is to clarify the technical details of the numerical solutions to the equations of motion, the scaling properties of the results, and their dimensional analysis in the A = ln b -coordinates. This is essential for extracting the physics out of the numerics. All the details are built in the numerical code SolveFiniteTTachyons deposited in [38]. In this treatment, the solutions in the A-coordinate system are considered primary, and the z-system is just an auxiliary coordinate system used to relate the results to known holographic formulae. The treatment extends [33], but we shall rewrite it explicitly in the form that the actual numerical code [38] uses, and keep each stage of the equations dimensionally consistent. We first define the notation: the fields q 1 (A 1 ), f 1 (A 1 ), λ 1 (A 1 ) and τ 1 (A 1 ) are the fields produced by numerical equation solving, expressed as a function of the coordinate A 1 . These will be referred to as level 1 solutions. The level 1 coordinate A 1 is the one in which the numerics is defined, and the horizon sits at A 1,h = 0. Level 2 solutions are obtained after f -scaling (see next section) and level 3 solutions [33] are the final ones with the fields, observables and the coordinate in the units corresponding to the desired UV boundary conditions.
We shall consider at first only V-QCD at µ = 0, and then devote a separate section to the µ = 0 case.
We define the coordinate z by with the boundary condition z(A = ∞) = 0. Notice that this is defined with the final scaled level 3 fields, and so we have precisely one system of z-coordinates, which we never scale.

B.2 The f -scaling
We generally want the function f to asymptote to 1 in the UV (z → 0 or equivalently A 1 → ∞) in order to have the standard Minkowski coordinate system with c = 1 on the boundary. When the boundary conditions are set at the horizon, this is not in general guaranteed. Fortunately the equations of motion are invariant under a combined scaling of f and q, such that if f 1 , q 1 are solutions, then also the pair with no change to the other fields or the coordinate A 1 = A 2 , is a solution for any value of f scale , although with different boundary conditions. Choosing f scale = 1/ f (A 1 = ∞) gives us the desired solution. From now on, fields and coordinates with the subscript 2 denote the numerical solutions scaled in such a way. We will call these the level 2 solutions. In [38], this scale factor appears as fscale and is explicitly used in generating the scaled solutions. The solutions produced by SolveFiniteTTachyons are level 1 in this notation, whereas SolveAndScaleFiniteTTachyons produces functions that are level 2 in this notation. The scaling itself is carried out in ScaleSolution, which can also be used to convert the level 1 solutions produced by SolveFiniteTTachyons to level 2 solutions.
The solution produced by this scaling no longer corresponds to the initial conditions set in the numerics. Specifically, if the original equation solver was started with the condition then the new solution corresponds to with the initial conditions for the other fields unchanged. The code [38] sets the initial conditions 6 where the first is chosen arbitrarily, since the magnitude of f is anyway set by f -scaling, and the second was derived in (4.7). The post scaling boundary condition then simply is with the rest of the fields unchanged.

B.3 The Λ-scaling
The UV -expansion (see (4.12) and Appendix A in [33]) is , (B.14) whereÂ 0 is a constant of integration. Here L UV is the asymptotic value of −q 2 (A) at large A (or equivalently A 2 ). Using these, we find Since we want to find a solution where Λ = Λ 0 , we write this in the form 16) and observe that the solution with the shifted coordinate A = A 2 −Â has the required asymptotics. Since the equations of motion are invariant with respect to shifts in A, this is also a solution of the equations, although with different boundary conditions. These are the level 3 solutions. We further denote eÂ = Λ/Λ 0 ≡ Λ scale . This is the factor that appears in [38] as Λscale, although in some places it is (inaccurately) denoted as simply as Λ. This factor is dimensionless. Using (B.16) converges somewhat slowly for practical purposes due to the O(λ) = O(A −1 ) corrections. We speed up that convergence by consideringÂ as a function of A max as given by (B.16), where A max is the limit up to which the numerical solution has been computed. From the numerical process, we know the derivatives of the fields, so we can computeÂ (A max ) and derive the formulâ which cancels the O(λ) corrections. The value ofÂ computed by this method is returned by SolveAndScaleFiniteTTachyons, which uses ScaleSolution to scale the level 1 solutions to level 2 and to derive Λ scale . We now write the transformation equations explicitly: Especially note that the horizon value of, for example, λ(A h ) = λ 2 (A 2 (A h )) = λ 2 (A h + A) = λ 2 (−Â +Â) = λ 2 (0). In other words, the horizon value of any field h in the solution with the correct asymptotics, is the same as the value of the original function h 1 coming from the numerics, evaluated at A 1,h = 0. The conformal factor of the metric appears in several physical observables. In level 3 coordinates it is simply In addition, the derivatives of fields in the z-coordinate system often play a role, and we observe that for example 24) and especially at the horizon An identical result holds for any field.
Since it is possible in this way to eliminate the need to explicitly shift the fields, and thus the need to keep track of one extra variable, the actual numerical code [38] does precisely this. In the code A always refers to A 2 = A 1 , the horizon is always at A 2 = 0, and the fields used to compute the physical observables are level 2.
Since the equations of motion are invariant under shifts of A without any corresponding change in the fields, the initial conditions for the fields themselves at horizon when expressed in terms of the A -coordinates are not changed by this scaling. The relation between A and z is what changes. Note however that once we introduce the gauge field Φ corresponding to a chemical potential, this changes since Φ explicitly breaks this shift invariance. We will return to that later.

B.4 Physical observables at µ = 0
Using the previous results, we can work out the formulas for physical quantities used in the code. The temperature is where we used (B.25), (B.7) and (B.6) . In the code, T is returned by TemperatureFromSols.
Note that both f scale and Λ scale are dimensionless, with the function q carrying one dimension of length, giving the correct unit of 1/length = energy. Also note that from (B.6) one sees that the unit of length in q ultimately comes from the potential, which is proportional to 1/L UV 2 . This shows that Λ scale is the dimensionless factor which tells the relation between the 4D boundary units and the units of the potential. The entropy density comes from This is returned in the code by s4G5FromSols. Note that b(A h ) is dimensionless, so the entropy density picks up its units from the gravitational constant G 5 .
The quark mass is expressed as where A max is the maximum A to which the equations of motion have been solved. Except for the appearance of Λ scale , the shift between A and Similarly as with the determination ofÂ, the A −1 corrections to m q are rather large at easily reachable values of A max . As before, we can consider m q as a function of A max and take its value at another point A b < A max . Using from the above that m q (A) = m q (1 + kA −1 ) for some unknown coefficient k, we can cancel the O(A −1 ) corrections: Since we know the derivatives of the fields from the numerical process, we can go further and take the limit In practice the finite difference method of (B.29) is slightly more stable and converges only very slightly slower, and that method is therefore used in the code by default. The function QuarkMass computes the mass from the solutions with this method.

B.5 Chemical potential
In (2.2), we introduce a zero component Φ of a gauge vector field in the bulk to model a chemical potential in the boundary theory. It turns out that the UV asymptotics are not affected by this addition, and so we will want to do similar scalings as in the zero chemical potential case. However, the full structure of the solution with respect to scaling the UV-variables does change, since there are new terms in the equations of motion, of the form where we have written (2.33) without substituting the solution of the Φ equation of motion.
If we have a solution with subscripts 1, including Φ 1 which has as yet undefined transformation properties, we have From this we conclude that if f 1 , q 1 , λ 1 , τ 1 , Φ 1 solve the equations of motion with UV asymptotics corresponding to the level 1 fields, then the corresponding level 3 functions solve the equations of motion with the correct UV asymptotics, if the function Φ 1 is replaced with Φ, such that It is apparent by inspection that the rest of the equations of motion are also invariant under this substitution. We naturally call Φ(A) a level 3 gauge field. The addition of a new field of course also adds a new pair of initial conditions. The fundamental physical constraint (2.19) in A-coordinates requires that we set Φ h = 0, so the remaining initial condition is determined by the derivative of Φ at the horizon, Φ h . Now given a level 1 solution, corresponding to the initial condition Φ 1,h , the scaled solution clearly corresponds to the initial condition The field Φ is a cyclic coordinate: its equation of motion is d dA which we can immediately integrate to the form (2.15): Different values ofn correspond to different initial conditions for the Φ field. Evaluating this at the horizon for a given solution or a set of initial conditions gives us the value ofn corresponding to that solution. Specifically, using the standard boundary conditions for starting the numerics we have in terms of the level 1 solution On the other hand, applying known scaling properties of the fields to the lhs of the same expression for the level 3 solution leads tô Since the level 3 solutions were the final ones, this gives us the scaling property ofn. We can solve (B.38) to yield an explicit expression for Φ in terms ofn and the other fields (forΦ, see (2.16)): Since Φ appears in the equations of motion always in the combination L 2 A Φ , we could entirely eliminate the choice of L A at this stage by rescalingn →nL 2 A . Therefore we can set L A = 1 without loss of generality in the numerics. Plugging the resulting formula into the equations of motion gives us the equations (2.31)-(2.34) on which [38] is based on. Solving the highest derivatives from those leads to the form in TachyonEquationsOfMotion. Since A 1,h = 0,n 1 matches with the scale invariant quantityñ of (4.6):n With this substitution and using the scaling properties it is apparent that, once the equations have been solved and subjected to the f -scaling to yield λ 2 , f 2 , τ 2 and q 2 , we can write Φ as This function is returned in the code by APrimeFromSols. It is the fully scaled form, but expressed as a function of the coordinate A 2 , which has not been shifted, i.e. it is in the same coordinate system as all the other functions returned by the code. It is level 2 in the same sense as the rest of the level 2 functions: the coordinate system is such that the horizon is at zero, but the units are such that it needs no further factors of Λ scale . It can be used fully consistently with all the other output functions, but note that if for some reason the coordinate system would be shifted again, Φ would then be scaled again according to (B.35). An interesting point is that the dependence on the actual numerical solution is precisely the same as for the entropy density s in (B.27). Thus one has a physical interpretation for the input parameterñ, it is simplyñ = 4πn(λ h ;ñ)/s(λ h ;ñ), where one also inserted the constants given in (4.17).
The other observable is of course the chemical potential itself. The holographic formula for it is for which we need to integrate (B.38). The correct boundary condition is that Φ(A h ) = 0, yielding This, and also the function Φ(A 2 ), is returned in the code by AAndMuFromSols, with ∞ replaced by the upper limit A max of the range for which the equations have been solved. In addition, the code uses L A = 1, but any other choice can be implemented by simply scaling n and µ.
The quark mass m q can be computed from the formulas presented in Appendix B given the initial conditions at horizon. However, for computing physical results, we For µ the smallest values ofñ are 0.001, 1. Note that T develops a minimum around λ h = 0.3 forñ > 9.5. The T, µ derived from here is in Fig. 18. are interested rather in finding a class of solutions corresponding to a predetermined value of m q , in this paper specifically m q = 0. This is in principle a simple problem of numerical function inversion, but it is complicated by the fact that the inverse is multivalued and that computing values of m q (λ h ,ñ, τ h ) takes a considerable amount of time (of the order of 1 second per point on a single core of an Intel i7 level processor).
The main task involves determining τ h (λ h ,ñ; m q ), given a fixed pair (λ h ,ñ). We will denote m q (τ h ) ≡ m q (λ h ,ñ, τ h ). As discussed in previous papers [32,33], m q (τ h ) may have several zeroes, corresponding to different Efimov vacuums, of which the most stable is the one with the largest τ h . This means that it is not enough to find a zero, but rather we have to be able to bracket an interval containing the last zero before the asymptotic rise of m q (τ h ) at large τ h , or alternatively deduce that no zero at finite τ h exists. In addition, since this needs to be done in at least tens of thousands of points on the (λ h ,ñ) -plane, the search must be fully automated and reliable enough to not need manual checking of the solutions.
The function τ hFromQuarkMass in [38] does this with a heuristic method that will be briefly described here. We omit some details, for which we invite the interested reader to look into the the code itself.
1. First we find a point where the solution exists and m q < 0, starting from an initial guess, by default τ h = 1. If the solution does not exist at all at the initial guess, τ h is multiplied by 2 to form a new guess. This is repeated until a point τ h,exist where the solution exists is found. After that, a point where m q (τ h ) < 0 is found by a binary search between ]0, τ h,exist ]. We denote that point by τ h,min . If this point is not found in a predetermined number of bisections (default is 80), we conclude that a chiral symmetry breaking solution does not exist for this pair (λ h ,ñ).

2.
We look for τ h,max such that m q (τ h,max ) > 0 by progressively doubling τ h until τ h,max is found.
3. We now have two points τ h,min < τ h,max such that at least one root of m q lies between them. Starting a numerical root finder in this bracket with Brent's method would be guaranteed to find a root, but unfortunately there is no control over which root. We need to start looking for zeros of m q in this interval. This is complicated by the fact that the distances between the zeros in τ h become exponentially larger toward increasing τ h . The heuristic we use determines an initial step length by the formula ∆τ h = τ h,min (( . Then m q (τ h,i ) is computed at points τ h,min +n∆τ h,min , n = 0, 1, 2, and we form the unique parabola passing through all of these points. The distance between its two roots is used to provide a local estimate of the distances between zeroes, which is used to determine a new step length ∆τ h . We then compute m q at intervals of ∆τ h until we find a zero, that is, m q changes sign during a step. 4. Once a zero is found, we update ∆τ h to be the distance between τ h,min and the zero divided by small safety factor (default = 5).
5. We continue to iterate with step length ∆τ h , and whenever a new zero (a change of sign in m q ) is found, we update the step length ∆τ h to the distance between the two latest zeroes divided by the safety factor. This iteration is terminated, and we take the last zero found, τ h,last , as the correct root, when the following conditions hold simultaneously: (a) m q has not decreased from the last step, since we know that asymptotically m q grows.
, where k is a heuristically determined number, typically a few hundred. This is the main condition used to ensure that the search goes on for long enough to reach the region of asymptotic growth.
Once the iteration described above completes, we are reasonably confident that τ h,last and τ h,last − ∆τ h bracket the largest zero, and simply use a standard root finder implementing Brent's method to find the precise location of that root.

D Numerical results for T and µ
As everywhere in this paper, numerical results in this appendix are all computed for the potentials (2.6)-(2.8) withμ = − 1 2 and for x f = 1. First, Fig. 16 plots T and µ as functions of λ h for fixed values ofñ for the tachyonfree chirally symmetric solutions. In Fig. 17

E Thermodynamics along λ end
The lower limit λ end (ñ) of the physical region of the tachyonic solutions on theñ, λ h plane plays an important role in the thermodynamics. We shall here analyse its properties.
The lower limit λ end (ñ) arises when one tries to determine what values of τ h are possible so that after integrating towards the boundary m q = 0 is obtained. One finds   that τ h = τ h (λ h , m q = 0,ñ) is a monotonically growing function of λ h (see, e.g., Fig.  5 of [33]) which starts at some λ h = λ end (ñ). As long asñ 8 one has the normal situation in which τ h (λ end ) = 0 but ifñ 8 the value τ h (λ end ) > 0, (see Fig.22). An accurate plot of λ end (ñ) has been presented in Fig. 7, in Fig. 22 it is shown together with the location of the deconfining transition. As discussed in Section 5.1, the chiral transition takes place along λ end . The extreme situation is that when τ h is so large that V f decouples due to the e −aτ 2 factor. The largeñ limit of (4.19) is then simplỹ n max = V g (λ h )κ(λ h ).
(E.1) Figure 22. Left: Behavior of τ h (λ h ,ñ, m q = 0) for largeñ. Note the jump in the beginning of the curve. Right: The physical region of tachyonic solutions λ h > λ end compared with the values of λ h at which the broken phase pressure vanishes. This is where the deconfining transition takes place, the corresponding temperatures T h (µ) are plotted in Fig. 1. The chiral transition takes place along λ end (see Section 5.1).
The upper limit of λ h moves to infinity and the large λ h limit of (E.1) is n max = 23.99 + 12.34 In the T = 0 limit it seems that s ∼ b 3 h goes to a finite limit even though T = 0. In fact, s b varies only very little along the chiral equilibrium curve.

F Large scale behavior
Since our model has asymptotic freedom built in it, we can at large T fit the magnitude of the pressure to the ideal gas limit  Large T means small λ h and we thus want to compute the pressure p(λ h ,ñ) at very small λ h and finiteñ in the approximation In this section, the limit λ h → 0 is always implied and the argument λ h (equivalent to z h ) is often omitted. Thus here V g = 12, κ = 1 and, for x f = 1, L ≡ L UV = (1 + 7 4 x f ) 1/3 = 1.401 and The pressure is obtained by first doing the pressure integral (5.6) along λ h atñ = 0 and then at fixed λ h the pressure integral (5.7) fromñ = 0 to someñ. The former is simple and gives (F.4) The latter becomes in the approximation (F.2) To evaluate this we must work out A 0 (z) and µ = µ(z h ,ñ) from (2.18) and T = T (z h ,ñ) from (2.37) in the approximation (F.2). Here z h is equivalent to λ h due to (F.2) and the order of arguments is irrelevant.
Noting that 4 3 , −x 6 y 2 ) (F.6) one finds and This form of T shows explicitly that T vanishes atñ = L 2 A V 2 g − V 2 f = 10.457 L 2 A , i.e., at the physical region boundary in Fig. 7 (where L A = 1). By taking the ratio one sees that µ/(πT ) is essentially determined byñ so that it grows monotonically from 0 to ∞ at the physical region boundary.
Inserting these exact forms to (F.5) and integrating one finds that so that the final total pressure in the limit of λ h → 0,ñ finite becomes where we also used (4.17). This further implies that = T s − p + µn = 3p in this UV corner of parameter space. Note the mixed notation, p = p(z h ,ñ) is given directly by the above equations, but if we want p(T, µ) we must solve z h = z h (T, µ) andñ =ñ(T, µ) from the exact expressions (F.7) and (F.8).
To compare with (F.1), consider first the limit T → ∞, µ = constant. Taking the ratio of (F.7) and (F.8), expanding inñ and inverting the series one obtains (F.11) and 1 z 3 Inserting this to (F.10) gives The two parameters G 5 and L A can be fixed by the magnitudes of the T 4 and µ 2 T 2 terms. Comparing the T 4 terms of (F.1) and (F.13) gives first [33] (F.14) Using this the µ 2 T 2 terms agree if, inserting (F. 3) and x f = 1, For completeness, the µ 6 /T 2 term in (F.13) is Thus p starts falling below p idea , the non-expanded result is in Fig. 25. The comparison of (F.1) and (F.13) can also be carried out in the limit T → ∞, µ/T = constant. Fixing the normalisation at µ = 0 we simply have  whereñ( µ T ) is to be determined by inverting the ratio of (F.7) and (F.8) numerically; the small-µ/T terms were given in (F.11). The result is plotted in Fig. 25. The value of L A was fixed so that the ideal gas µ 2 term was correctly reproduced. Now one sees that the good agreement extends to large values of µ, at µ = 4T the deviation is 3%. One can work out analytically the asymptotic limit at µ T which corresponds tõ n →ñ max = L 2 f . It depends on V g = 12 and V f and its numerical value is 0.572. There is no obvious constraint leading to the value 1, but it is nevertheless rather close to this value.
It may be useful to compare holographic and perturbative QCD predictions for other quantities, too. At T = 0 [31,34] finds that the holographic and perturbative QCD results for the correlator of vector flavor currents agree in the UV if see Eq. (C.10) in [34] with w 2 = L 4 A κ 2 = L 4 A and W 0 = V f /x f . This matches exactly with the combination of (F.14) and (F.15). Note that both the pressure at large T, µ and the vector correlator at large momentum depend only on the combination L 2 A κ(λ h = 0), here we have assumed κ(0) = 1. These quantities can be fixed separately using the scalar correlator. Combining the result in Eq. (C.21) of [34] and (F.15) one finds that This modified value of κ would affect the normalisation of τ and consequently that of the chiral condensate, but not the results in this article. As another example of comparisons of weak coupling and holographic computations one may also compare this result with an analogous analysis of the finite temperature correlators of the energy momentum tensor in the UV [52]. Using the thermal normalisation (F.14), the holographic result of the shear correlator is too small by a factor 4/9 with respect to the perturbation theory one, for the bulk correlator the results agree. This result illustrates the fact that this action cannot describe all the phenomena in the weak coupling region at the same time, an action with higher derivatives is needed.