Massive holographic QCD in the Veneziano limit

QCD at finite, flavor independent quark mass is analyzed by using bottom-up holography in the Veneziano limit, where the backreaction of quarks to the gluon dynamics is fully included. The dependence on the quark mass of observables such as the bound state masses, the chiral condensate, the S-parameter, and the critical temperatures is studied. Many of the results are argued to be universal, i.e., independent of the details of the holographic model, and compared to explicit computations in the V-QCD models. The effect of adding four-fermion operators in QCD is also discussed.


Introduction and summary
QCD displays an interesting phase structure as the number of flavors N f and number of colors N c are varied. It is natural to discuss the phase diagram in the Veneziano limit [1,2]: as a function of the variable x which has become continuous in this limit. The "standard" expectation for the diagram is shown in Fig. 1 (at zero temperature and quark mass). We restrict to the interval 0 < x < 11/2 ≡ x BZ where the theory is asymptotically free. Various regimes can be identified: • The QCD regime 0 < x < x c where the infrared (IR) dynamics is similar to ordinary QCD (having N c = 3 and a few light quarks), with confinement and chiral symmetry breaking.
• The walking regime with 0 < x < x c and x c − x 1 where the coupling constant of the theory varies very slowly, i.e., "walks", over a large range of energies.
• The conformal window x c < x < x BZ where the theory runs to an IR fixed point (IRFP), which is located at such a small value of the coupling that the chiral symmetry stays intact.
The existence of the conformal window is solid in the Banks-Zaks (BZ) limit x → x BZ , because the value of the coupling is parametrically small and perturbation theory is trustable [3]. It is also credible that dynamics at small x 1 is similar to ordinary QCD -the ∼ 1/N 2 c corrections arising in the Veneziano limit are not expected to change the picture qualitatively. But the nature of the "conformal transition" at x = x c , and the behavior of the theories near the transition, is an open question.
The phase diagram of Fig. 1 with a walking regime is obtained in the Dyson-Schwinger approach [4]. The transition is then of the BKT type [5], associated with the so-called Miransky scaling law [6]. The existence of the walking regime is important, since theories in this region may have properties, which are desirable for technicolor candidates [7,8] (see also the reviews [9]). The location of the transition has been estimated by using different approaches [10]. However, it has also been suggested that the transition is discontinuous, such that the dynamics "jumps" and walking is absent [11]. There is an ongoing effort to clarify these issues by using firstprinciples lattice simulations [12,13], but as it turns out, obtaining reliable results in the transition region is difficult.
The holographic V-QCD models [14] also have the phase diagram of Fig. 1, including a BKT transition an walking. The "V" in V-QCD refers to the Veneziano limit of (1.1), where the models are defined. V-QCD is based on two building blocks, the improved holographic QCD (IHQCD) [15] for the gluon sector, and a method for adding flavor by inserting a pair of space filling D4 − D4 branes [16,17]. This means that the gluon sector is described by a five dimensional Einstein-dilaton gravity, and the flavor sector is described by a tachyonic Dirac-Born-Infeld action. The sectors are fully backreacted in the Veneziano limit. The structure of V-QCD at finite temperature and chemical potential has been studied in [18][19][20][21]. Two-point correlators and bound state masses were analyzed in [22,23].
There have also been numerous other holographic models addressing the phenomena expected in QCD at finite values of x. Walking gauge theories have been modeled by using as a starting point the traditional bottom-up models for example in [24][25][26]. Walking within the top-down framework has been found and studied in [27]. The conformal transition [28] was studied in a top-down setup in [29], and by using a tachyon-Dirac-Born-Infeld action in [30]. IHQCD has been used to study walking and IR conformal theories by modifying the holographic RG flow "by hand", without inclusion of dynamical fermionic degrees of freedom [31,32]. Walking dynamics and the conformal transition have also been studied in Dynamic AdS/QCD [33,34], which has partially similar ingredients as V-QCD.
Including effects due to finite quark masses in holographic models for QCD is important for various reasons. First of all ordinary QCD, which describes the strong interactions observed in nature, has finite quark masses. Therefore, proper understanding of the dependence of the model on the quark masses should lead to better models for ordinary QCD.
Knowledge on the dependence of quark masses is also useful when analyzing the data from lattice simulations, which are often carried out at unphysically large quark masses due to technical reasons, and therefore extrapolation in the quark mass is needed. For lattice studies which aim at uncovering the phase diagram of QCD as a function of N f , extrapolations in (typically flavor independent) quark mass are particularly important in the probably most interesting region near the conformal transition where lattice simulations are demanding. For example, the socalled "hyperscaling relations" [35] in the conformal window have turned out to be useful close to the conformal transition [36]. Holography may further help to analyze the physics of the transition, because it can provide a unified description of the phases in the vicinity of x = x c where most important observables can be computed rather easily for all values of the quark mass.
Finally, it is also important to ensure that the behavior of the holographic model is realistic for all values of the quark mass. Studying limiting behavior (such as the limit of large quark mass) may lead to nontrivial constraints on the model. Fixing the behavior of the model to agree with QCD in limiting cases where field theory computations are tractable, is also expected to improve the model for all values of the quark mass and other parameters.
In this article we study in detail the dependence of holographic QCD on a flavorindependent quark mass m q in the Veneziano limit. We argue that many of the basic results, such as the dependence of the energy scales on x and m q , are essentially universal, i.e., independent of the details of the holographic model (given some natural assumptions which will be specified below). Additional scaling results for meson and glueball masses, the decay constants, the chiral condensate, the S-parameter, and the critical temperatures are derived for V-QCD analytically at large and small m q and in the different regimes for x. Explicit results are computed numerically in V-QCD and are shown to agree with the expectations from the analytic studies.
The numerical analysis of [14,18,23] was mostly done at m q = 0, but some observables were already computed for a few values of m q and for a limited range of x. Here we extend these results to cover all relevant regimes (with 0 < x < x BZ ) on the (x, m q )-plane. Importantly, this extension does not require tuning the models or adding any new terms in the V-QCD action, as the value of the quark mass is determined through the boundary conditions of the model. Actually, we will use the choice for the V-QCD action with various potential terms in the action defined exactly as in [23].
We carry out a particularly detailed analysis of the dependence of the chiral condensate on m q , again comparing analytic formulas with explicit numerical results in V-QCD. It is then shown how this analysis is used to study the deformation of QCD by a four-fermion operator ∼ (qq) 2 . Understanding the effect of four-fermion operators is interesting for purely theoretical reasons, but also because such terms naturally arise in technicolor models due to the so-called extended technicolor interactions [9].
This article is organized as follows. In the remaining part of introduction we summarize our main results, and discuss the status and future of the exploration of V-QCD. In Sec. 2, we give a brief review of the V-QCD models. In Sec. 3, we discuss the universal scaling results for the energy scales of (holographic) QCD in the Veneziano limit, without explicitly referring to V-QCD. These results are confirmed numerically for V-QCD in Sec. 4, and extended to include the m q -dependence of the bound state masses and decay constants. In Sec. 5, we analyze the mass dependence of the chiral condensate, and prove the Gell-Mann-Oakes-Renner (GOR) relation for the pion masses for the fully backreacted case. The results of Sec. 5 are then used to study four-fermion deformations in Sec. 6. The mass dependence of the S-parameter and the pion decay constant is analyzed in Sec. 7. Finally, a finite temperature is switched on in Sec. 8, where the scaling laws for critical temperatures are derived and compared to numerical results. The Appendix provides technical details for the derivation of the results.
The article is long but many of its sections are largely independent, so a reader only interested in a specific topic may want to jump directly to the corresponding section. There are, however, exceptions: Sec. 4 requires studying Sec. 3 first, and Sec. 6 requires Sec. 5.

Summary of results
Let us then summarize the main results of this article. We identify three different regimes on the (x,m q )-plane, shown schematically in Fig. 2, where the dependence of the model on the quark mass is qualitatively different. In regime A, the quark mass is a small perturbation. Regime B is the "hyperscaling" region where the couplings constant walks, and the amount of walking is controlled by the quark mass 1 . Regime C refers to the limit of large quark mass (in units of the scale of the ultraviolet (UV) RG flow Λ UV ). The white regions in the plot indicate crossovers between the regimes A, B, and C -there are no phase transitions at finite quark mass and zero temperature.
In this article, we discuss in detail how the diagram arises quite in general in holographic 2 models. We show how scaling results for the energy and mass scales in the various regimes can be obtained by using relatively simple assumptions and straightforward analysis (explicit results are given in Eqs. (3.9), (3.10), and (3.14)). The most important findings are: • The physics at small m q , including regimes A and B of Fig. (2), is universal, i.e., qualitative features are independent of the details of the holographic model.
• At large m q (in regime C) the results are model dependent, and can directly be compared to field theory results which are tractable in the limit of large m q . We demonstrate that V-QCD reproduces most important features such as the decoupling of the massive quarks. Remarkably, we find that a specific subclass of the V-QCD models, which was found to be closest to QCD by analyzing the asymptotic meson spectra in [23], also works best at large quark mass. These are the models with "potentials I" below.
We compute numerically the meson and glueball masses in V-QCD, and show that they agree with the generic scaling results in all regimes.
Notice that Fig. 2 only shows changes in the dependence of the quark mass. There are also other features: In regime A, the theories close to x = x c are walking and therefore much different from the theories at low values of x, reflecting the behavior of the m q = 0 backgrounds, which are perturbed by the quark mass. Within regime B, the scaling laws change in a nonanalytic manner exactly at x = x c (see Eqs. (3.9) and (3.10)).
We also analyze in detail how the chiral condensate depends on the quark mass. The main results are the following: • In the QCD and walking regimes, there is an interesting spiral structure, which we call the Efimov spiral. We give an asymptotic formula [38] of the spiral at small m q in Eq. (5.14), and compare this to data in Fig. 9. We demonstrate how the spiral is related to the subleading Efimov vacua 3 of the theory and to the Miransky scaling law in the walking regime.
• The standard holographic proof of the GOR relation is extended to the case of V-QCD, which requires handling of the full backreaction and the logarithmic corrections which appear in near the UV boundary due to the RG flow. The relation is given in Eq. (5.25) and checked numerically in Fig. 11.
We use the results for the chiral condensate to check how Witten's method [40] for adding multi-trace deformations in holographic models by modifying the UV boundary conditions works for V-QCD. The resulting phase diagram for the case of a double-trace (four-fermion) deformation ∝ g 2 (qq) 2 at zero quark mass is shown in Fig. 13 in the (x, g 2 )-plane. For positive coupling g 2 , we see no change with respect to g 2 = 0, whereas negative g 2 induces an instability.
The mass dependence of the S-parameter, the pion decay constant, and related quantities is studied in detail.
• As any finite m q is turned on in the conformal window, the S-parameter discontinuously jumps from zero to a O (N f N c ) number. Except for this discontinuity, the dependence on x and m q is weak.
• We demonstrate a novel power-law scaling of the subleading terms in the Sparameter at small quark mass in the conformal window and in the walking regimes (see Figs. 19 and 20). It is argued that the power can be expressed in terms of the dimension of Tr F 2 at the IRFP as in Eqs. (7.28) and (7.30).
• We analyze the x-dependence of the S-parameter by writing it as a series over the contributions from the low-lying vector and axial vector meson poles (see Eq. (7.31)). It is shown that the increase of the S-parameter with x in the QCD and walking regimes can be attributed to slower convergence of the series.
Finally we analyze the dependence of the critical temperatures of the deconfinement transition and various crossovers on m q and x. The scaling laws of (8.5) - (8.8) and (8.10) - (8.11) are demonstrated numerically in Fig. 22. The second order chiral transition, which is found at zero quark mass in the walking regime [18], transforms into a crossover as finite m q is turned on. We also discuss how the thermodynamics approaches that of the YM limit (x → 0) as the quark mass is taken to infinity so that the quarks are decoupled.

Outlook
This article is part of an ongoing program [14,[18][19][20][21][22][23] for studying the properties of the V-QCD models, and in more general the structure of QCD in the Veneziano limit. There are several possible future directions to explore. As a continuation of this study, one could consider the case of flavor dependent quark masses, which would be interesting in order to construct more realistic models for ordinary QCD, where all quark masses are unequal. In holography, this means that the background solutions are nontrivial dependence on the flavor indices, and consequently are described in terms of a non-Abelian Dirac-Born-Infeld action. The precise definition of such a non-Abelian action is not known, so requiring the physics to be correct might lead to interesting constraints for it.
The study of the CP-odd terms of the V-QCD action is in progress at the moment [41]. These terms govern the physics of the axial anomaly and the theta angle of QCD. Similarly to the analysis at zero values of theta, requiring the regularity of the solutions in the IR, and among other things the correct behavior of the asymptotic meson spectra at finite values of theta result in constraints for the potentials appearing in the CP-odd terms, in addition to the constraints obtained in the YM limit [15,42,43]. After analyzing these constraints, we may make physically reasonable choices for these potentials, and in particular for their asymptotic behavior near the UV boundary and deep in the IR.
So far the potentials of the V-QCD have not been tuned to fit any nonperturbative QCD data, and therefore all predictions of the model are qualitative. But once the CP-odd sector of the action has been fixed, one can start fitting the parameters of the potentials both to experimental results for QCD (such as meson masses) and to lattice data (both for QCD at finite N f /N c [12] and for YM [44,45]). The hope is that the overall fit, together with the other constraints for the potentials, fixes the predictions of the model to a good accuracy for all relevant values of the parameters (such as x, the quark mass, temperature, and chemical potential). Consequently the model would give a effective description of QCD with real predictive power, rather than being just a toy model.

V-QCD
In this section, we will briefly introduce a class of bottom-up models for QCD, which we call V-QCD [14]. The V in the name refers to the fact that the models are defined in the Veneziano limit: V-QCD is based on two "building blocks". The first block (the glue sector of V-QCD) is IHQCD [15] which is a bottom-up model for Yang-Mills (YM) theory inspired by five-dimensional noncritical string theory. The second block (the flavor sector of V-QCD) is a framework for adding flavor via tachyonic Dirac-Born-Infeld actions [16,17]. This framework has been tested previously [46] in the probe (or 't Hooft) limit, i.e., without including the backreaction of the flavor branes to the background. However in V-QCD, and more generally in the Veneziano limit, the flavor and glue sectors are fully backreacted.

The V-QCD action
Let us then discuss briefly the dictionary of V-QCD. The most relevant fields are • The dilaton φ. The exponential e −φ is dual to the operator TrF 2 . We will denote λ = e φ below. As this notation indicates, its background value is identified as the 't Hooft coupling on the field theory side.
• The tachyon T ij which is a N f × N f matrix in flavor space. The combination T + T † is dual to the operatorq i q j whereas i(T − T † ) is dual toq i γ 5 q j . For the background solutions considered here we will take T ij = τ (r)δ ij so the flavor structure does not appear explicitly.
• The left-and right-handed gauge fields A L/R µ which are dual toqγ µ (1 ± γ 5 )q. They are also matrices in flavor space but we have hidden the flavor indices here. These fields evaluate to zero for the backgrounds considered in this article.
In addition, the scale factor A of the metric of the vacuum solution is identified as the logarithm of the energy scale on the field theory side. The blackening factor f in the metric may be either identically equal to one or a nontrivial function of r. Our convention will be that the UV boundary lies at r = 0, and the bulk coordinate therefore runs from zero to infinity (when f (r) ≡ 1) or up to a horizon at a finite value of r (when f (r) has a nontrivial profile). The metric will be close to the AdS metric in the UV: A ∼ − log(r/ ), where is the (UV) AdS radius. In the UV, r is therefore identified roughly as the inverse of the energy scale of the dual field theory. The action for the V-QCD model consists of three terms: where S g , S f , and S a are the actions for the glue, flavor and CP-odd sectors, respectively. Explicit expressions for the first two terms will be given below. As discussed in [14], only these terms contribute in the vacuum structure of the theory if the phases of the quark mass matrix and the theta angle vanish. The last term S a is important for the realization of the theta angle and the axial anomaly of QCD [17]. This term has been written down explicitly in [23], and it will be zero for all configurations discussed in this article. The glue action is that of IHQCD [15]. It includes five dimensional Einstein gravity and the dilaton λ = e φ : The flavor action is the generalized Sen's action [17,47] (see also [48]), (2.5) where the quantities inside the square roots are defined as with the covariant derivative The trace Tr is over the flavor indices -recall that the fields A L , A R as well as T are N f × N f matrices in the flavor space. It is not known, in general, how the determinants over the Lorentz indices in (2.5) should be defined when the arguments (2.6) contain non-Abelian matrices in flavor space. However, for our purposes such definition is not required: our background solution will be proportional to the unit matrix I N f , as the quarks will be all massless or all have the same mass m q . In such a case, the fluctuations of the Lagrangian are unambiguous up to quadratic order.
The form of the tachyon potential that we will use for the derivation of the spectra is This is the string theory tachyon potential where the constants have been allowed to depend on the dilaton λ. For the vacuum solutions (with flavor independent quark mass) we will take T = τ (r)I N f where τ (r) is real, so that The coupling functions κ(λ, T ) and w(λ, T ) are allowed in general to depend on T , through such combinations that the expressions (2.6) transform covariantly under flavor symmetry. In this article, we will take them to be independent of T , emulating the known string theory results.

Potentials and the holographic RG flow
In order to fully fix the action, the potentials V g (λ), V f 0 (λ), a(λ), κ(λ), and w(λ) need to be specified. It turns out [14] that V g (λ) must satisfy the same constraints as in IHQCD [15]. The other potentials will be subject to analogous constraints. We review the main idea here, and the details can be found in [14,23]. First, identification of the field λ as the 't Hooft coupling and the scale factor A as the logarithm of the field theory energy scale defines the holographic renormalization group (RG) flow and the holographic beta function for the coupling as in IHQCD, with the understanding that the fields are evaluated on the r-dependent background solution. In IHQCD, the dilaton potential V g can be directly mapped to the holographic beta function [15] at any value of r, that is, at any energy scale. In V-QCD, there is an additional field, the tachyon, whose background value is linked to the running quark mass. Therefore one may define a holographic gamma function, which controls the holographic RG flow of the quark mass [14]. The mapping between the beta and gamma functions and the potentials is, in general, more complicated than in IHQCD, but simplifies in the UV. The behavior of the potentials in the UV (where λ → 0) is then restricted by requiring that the holographic beta and gamma functions match with their QCD counterparts in the UV. For the UV structure to be consistent, all potentials are chosen to be analytic at λ = 0, and the series coefficients can be related to those of the perturbative beta and gamma functions in QCD. It turns out that the dilaton potential V g is consequently mapped to the perturbative beta function of YM theory as in IHQCD. Due to the backreaction, the beta function of QCD (in the Veneziano limit) is mapped to the combination V g − xV f 0 . This mapping leaves one undetermined parameter, the UV normalization of V f 0 , which we call W 0 . The gamma function of QCD fixes the UV behavior of the ratio a/κ. Here we will match the λ = 0 expansions of the potentials up to two loops for the beta function and up to one loop for the gamma function of QCD.
Using perturbative QCD to determine the UV behavior of the potentials may be surprising, since holography is in general not expected to work at small values of the coupling λ. The idea is, however, that by using this procedure correct boundary conditions for the more interesting IR dynamics are obtained. Notice also that the procedure can be seen as a rather mild generalization of what is usually done in bottom-up holography. For example, the bulk mass of the tachyon is typically required to satisfy the relation −m 2 2 = ∆(4 − ∆) (at least in the UV), where is the IR AdS radius and ∆ is either the dimension of the quark mass or the chiral condensate. Here this relation is effectively generalized to include loop effects, i.e., the perturbative anomalous dimension of the quark mass, which is roughly mapped to the first few coefficients in the expansion of the bulk mass of the tachyon at λ = 0.
In this article, we will carry out one more check which demonstrates that our UV boundary conditions make sense. Namely, we prove that the GOR relation holds even in the backreacted case and that the UV RG flow of the quark mass and the condensate, imposed by the matching to perturbative QCD, cancels in this relation, as it should (see Sec. 5.3).
The choice of the potentials in the IR is more relevant since it affects how the nonperturbative physics of QCD is modeled. Since we are working with bottomup models, there is a lot of freedom in choosing the potentials, and it is important to choose them such that the IR physics resembles that of QCD. The mapping to the beta functions is not useful at large values of the coupling because of scheme dependence. The asymptotic behavior of the potentials at large λ can however be constrained heavily by comparing to several different observables, most importantly the asymptotics of the spectra at large mass.
The IR behavior of the dilaton potential V g , i.e., its asymptotics as λ → ∞, can be fixed by requiring (among other things) confinement and correct asymptotic behavior of the glueball spectrum at large excitation numbers. The remaining parameters have been fitted to YM data [43]. Similarly, asymptotics of the meson spectra sets strict constraints on the large λ asymptotics of the other potentials in the V-QCD action [23]. The remaining degrees of freedom will be fitted to experimental and lattice data in future studies.
In this article, we will be using the choices "potentials I" and "potentials II" which are exactly the ones given in [22,23], and are defined explicitly in Appendix G. The potentials I reproduce more accurately qualitative features of QCD. This choice has the UV parameter W 0 set to a constant value 3/11, and assumes w = κ. The potentials II are included in order to study the model dependence of the results. This choice has the UV parameter W 0 set to a value which guarantees "automatically" the Stefan-Boltzmann normalization of pressure at high temperatures [18], and assumes w = 1. Most of the results derived in this article are qualitative (e.g., scaling laws of various observables) and therefore insensitive to the details of the potentials.

Background solutions
Let us discuss some general features of the background solutions of the V-QCD models, first restricting to the standard case, which has a phase diagram similar to what is usually expected to arise in QCD. Such a phase diagram is obtained if the potentials are chosen as discussed above.

Zero temperature
In this article we will mostly discuss solutions at zero temperature. In this case, the blackening factor f in (2.2) is trivial, f ≡ 1. To find the background, we consider r-dependent Ansätze for λ, and A. As pointed out above, we assume that the quark mass is flavor independent, and therefore take T = τ (r)I N f . We also set all other fields to zero, and look for solutions to the equations of motion (EoMs). The models are expected to have two classes of (zero temperature) vacuum solutions [14]: 1. Backgrounds with nontrivial λ(r), A(r) and with zero tachyon τ (r) = 0. These solutions have zero quark mass and intact chiral symmetry.
2. Backgrounds with nontrivial λ(r), A(r) and τ (r). These solution have broken chiral symmetry. As usual, the quark mass m q and the chiral condensate are identified as the coefficients of the normalizable and non-normalizable tachyon modes in the UV.
In the first case, the EoMs can be integrated analytically into a single first order equation, which can easily be solved numerically. The regular solution ends on an IRFP, where the dilaton approaches a constant value, and the geometry is asymptotically AdS 5 . In the second case, one needs to solve a set of coupled differential equations numerically. The regular solution ends in a "good" IR singularity, where both the dilaton and the tachyon diverge. Let us first recall what happens at zero quark mass. The ratio x = N f /N c is constrained to the range 0 ≤ x < 11/2 ≡ x BZ where the upper bound was normalized to the Banks-Zaks (BZ) value in QCD, where the leading coefficient of the β-function turns positive. The standard 4 phase diagram at zero quark mass has two phases separated by a phase transition at some x = x c within this range.
• When x c ≤ x < x BZ , chiral symmetry is intact. The dominant vacuum solution is in the first class with the tachyon vanishing identically and an IRFP. Therefore the geometry is asymptotically AdS 5 in the IR.
• When 0 < x < x c , chiral symmetry is broken. The dominant vacuum therefore is in the second class with nonzero tachyon even though the quark mass is zero.
The IR geometry is singular. The singularity is of the "good" kind [49], and therefore fixes all IR boundary conditions uniquely.
For potentials satisfying some reasonable requirements [14,23], the phase transition at x = x c is due to an instability of the tachyon at the IRFP [28,50]. That is, the solution with vanishing tachyon is unstable if the bulk mass of the tachyon at the IRFP −m 2 * 2 * , where * is the IR AdS radius, violates the Breitenlohner-Freedman (BF) bound [51]: where ∆ * is the dimension of the quark mass at the fixed point. Notice that when the bound is violated, ∆ * becomes complex. As a consequence of violating the BF bound, the phase transition at x = x c (which is only present at zero quark mass) involves BKT [5] or Miransky [6] scaling, for values of x right below the critical one. The order parameter for the transition, the chiral condensate σ ∼ qq vanishes exponentially, as x → x c from below. Here the constant K is positive. At finite quark mass, the BKT transition disappears: the background is always in the second class with finite tachyon, and the dominant vacua at all values of x are smoothly connected. Therefore the geometry ends in a "good" singularity in the IR. In particular, the IR geometry changes in a discontinuous manner (from an IRFP to the good singularity) when a small quark mass is turned on in the conformal window (x c < x < x BZ ). This discontinuity causes interesting behavior of observables which will be discussed in the following sections.
The models may also have unstable subdominant vacua when 0 < x < x c . We will discuss such vacua in more detail in Sec. 5.

Finite temperature
At finite temperature, one can first identify two types of background geometries [18]: 1. The thermal gas solutions which have the same (and therefore temperature independent) r-dependence as the zero temperature solutions, and in particular f ≡ 1.
2. The back hole solutions which have a nontrivial f (r) and end on a horizon in the IR.
As usual the temperature is given by (the inverse of) the length of the compactified time direction, and is equal to the Hawking temperature of the black hole for the second type of solutions. Both of these geometries further split into two classes, one having zero and the other having nonzero tachyon. Therefore there can be up to four qualitatively different competing saddle points. The finite temperature phase diagram has been studied in [18]. The following structure was found at zero quark mass: • For x c ≤ x < x BZ , the finite temperature phase is the tachyonless black hole.
A tachyonless thermal gas solution which has the same r-dependence as the zero temperature solution also exists at all temperatures but is subdominant. Therefore if the system is first prepared at zero temperature, any amount of heating makes the system jump to a different phase immediately. At small temperatures, the finite temperature backgrounds are obtained by deforming the zero temperature backgrounds, which are asymptotically AdS 5 , only very close to the IR end. Consequently, small temperature thermodynamics is that of AdS, p ∝ T 4 , and the zero temperature transition is of 4th order.
• For 0 < x < x c , the low temperature phase is the tachyonic thermal gas phase which is smoothly connected to the zero temperature solution and breaks chiral symmetry. The high temperature phase is the tachyonless black hole phase. There is always a first order ("deconfinement") transition separating the phases, but it is also possible that an intermediate, chirally broken tachyonic black hole phase appears. If this is the case, chiral symmetry is restored at a separate second order transition, which has higher critical temperature than the deconfinement transition. As x → x c from below, the critical temperatures go to zero following Miransky scaling.
At finite quark mass, all phases are tachyonic due to the UV boundary conditions, and chiral symmetry is broken. At low temperatures the dominant vacuum is the thermal gas phase. When the system is heated it undergoes a first order transition to the black hole phase. The phase structure is therefore similar for all 0 < x < x BZ . The dependence of the critical temperature on m q and x will be discussed in Sec. 8. In addition to the phase transition, we find several crossovers which are linked to changes in the zero temperature geometry.

Energy scales of (holographic) QCD in the Veneziano limit
In this section we shall discuss the dependence of observables (at zero temperature) on the quark mass and x = N f /N c for holographic models of QCD in the Veneziano limit and in general, i.e., not necessarily only for V-QCD. The results will not be proven rigorously, but we will sketch how they arise from rather natural assumptions. Comparison to quantitative results from V-QCD will be carried out in Sec. 4.
It is useful to discuss separately the regions with small (possibly zero) and large quark mass.
• At small quark mass, we shall demonstrate that under certain natural assumptions, a universal picture arises from holography, which is not dependent on the details of the model.
• At large quark mass, the m q -dependence of energy scales and most observables is most naturally derived from QCD directly. Large quark mass probes the UV structure in holography, where predictions may not be reliable. We shall see that the results from holography are model dependent in this region. We discuss what is needed to match with the known results in QCD.
In this section we will concentrate on the behavior of the various energy scales for a generic holographic model. In the remaining of this article we will analyze concrete observables such as bound state masses, the chiral condensate, and the S-parameter for V-QCD. Also the behavior of these observables at small quark mass is to large extent independent of the model details, but it is just much easier to work with the explicitly fixed V-QCD action. We will exclude the BZ limit (i.e., the limit x → x BZ ) from our analysis for the moment. In this region the theory is fully under perturbative control and holographic approach may not be useful.
First we need to give rough definitions for various energy scales at zero temperature, assuming a generic holographic model. There must be a field dual to theqq operator, which we call the tachyon, and the geometry must be asymptotically AdS in the UV. We take the UV boundary to lie at r = 0 as above. Because we want to keep the discussion generic, the definitions for the energy scales will be rather sketchy. Precise definitions in the case of V-QCD will be given in Sec. 4. One could take Eqs. • Λ UV is the scale of the UV RG flow in QCD. In holography, it can be identified as the scale of the UV expansions of the vacuum solution.
• Λ IR is the soft IR scale which governs the IR expansions (or is set by an IR cutoff).
• The quark mass m q is defined as the source for the tachyon.
• Λ τ is the scale of chiral symmetry breaking (whenever it is broken). In holography it can be identified as the energy scale (inverse of r) where the tachyon grows large (or becomes O(1), more precisely). As we shall see, its dependence on m q is similar to that of the constituent quark mass in nonrelativistic quark models.
Notice that we only discuss the dynamical scales appearing through the background solutions, and for example the UV AdS radius is assumed to be fixed.

Small quark mass
By small quark mass we mean here that m q /Λ UV 1, and m q may also be zero. In this case our results will be essentially universal, i.e., independent of the details of the holographic model. Naturally, several assumptions must be made on the model, in order to ensure that the physics resembles that of QCD in the Veneziano limit. Most importantly, there must be an IRFP and a conformal window, and the conformal transition at x = x c is assumed to be a BKT transition, which is naturally implemented in holography by the tachyon hitting its BF bound at the IRFP [28,50]. For simplicity we also assume that the upper edge of the conformal window lies at the QCD value, x = 11/2 ≡ x BZ .
The model assumptions can be made explicit in terms of the tachyon background in various ranges of the bulk coordinate r. We will only need the solutions for small tachyon, so that the tachyon EoM is approximately linear. In the UV (where r → 0), we assume the standard behavior determined by the dimension of the quark mass: where σ is proportional to the chiral condensate. Here logarithmic running of the quark mass (and the condensate) could also be included, but it is not important because it will not affect the leading scaling behavior. In the vicinity of the fixed point when the tachyon is small and the BF bound is satisfied, corresponding to x c ≤ x < x BZ in QCD, we find that where the precise x-dependence of the anomalous dimension of the quark mass ∆ * depends on the model. When the BF bound is violated, and x c − x is small where ν = Im∆ * π (x c − x)/K as x → x c from below. The coefficient K depends on the model, and it satisfies Now it is straightforward to fix the integration constants (σ, C m , C σ , C w and φ) and compute the mass dependence of the various energy scales by using the following recipes: • Both the normalizable and nonnormalizable terms of the tachyon solution are separately continuous 5 . Therefore an approximate tachyon solution for all r 1/Λ τ can be found by gluing together the results from the different regimes (3.1), (3.2), and (3.3). Because the bulk mass of the tachyon varies smoothly, the terms of the tachyon can only jump by a factor O (1) as we move from one regime to another.
• The IR behavior of the tachyon in the regime r 1/Λ τ is nontrivial and qualitatively different from the formulas given above. Therefore, the IR boundary conditions obtained by requiring continuity at r 1/Λ τ are taken to be generic, i.e., it is assumed that the conditions are not fine tuned to pick any specific solution for r 1/Λ τ .

Zero quark mass
Let us first recall how the energy scales depend on x at zero quark mass. The results for x < x c are obtained by using the above assumptions and recipes. The computation can be found in Sec. 10 of [14] and the scaling results will also be reproduced by the analysis of Sec. 5 below, where the chiral condensate is studied in detail (see also [29]). We will only list the results here. When x < x c , Λ IR ∼ Λ τ (i.e., the tachyon becomes sizeable where the asymptotic IR geometry starts), and Λ τ is not defined in the conformal window (x c ≤ x < x BZ ). Depending on the value of x = N f /N c , there are three regions with qualitatively different behavior: • In the QCD regime, meaning that 0 < x < x c and x c − x 1, there is only one scale, Λ UV ∼ Λ IR , which corresponds to Λ QCD .
• In the walking regime, which is found when x < x c and x c − x 1, the IR and UV scales are related through Miransky scaling: The chiral condensate satisfies • In the conformal window (x ≥ x c ) the tachyon is zero as chiral symmetry is intact. In this case one expects that the scales of the UV and the IR expansions are the same, Λ UV ∼ Λ IR (as there is no obvious mechanism which would separate the scales).

Small but finite quark mass
Let us then generalize these results to nonzero quark mass, first assuming a small mass (m q /Λ UV 1). The detailed analysis uses the strategy formulated above, and can be found in Appendix A. We will only discuss the results here.
One can identify, for m q /Λ UV 1, two regimes on the (x, m q )-plane, shown schematically in Fig. 2: A The regime where the quark mass is a small perturbation. This is possible for 0 < x < x c as the m q = 0 solution is continuously connected to the solution with finite m q . In order to determine the extent of regime A, we will perturb the solution at vanishing quark mass by adding a small m q . It is a small perturbation so long as the nonnormalizable term remains small for r 1/Λ UV . This is equivalent to m q |σ(m q = 0)/Λ 2 UV |, or m q Λ UV 1 (3.7) in the QCD regime, and in the walking regime. Apart from the perturbation caused by the quark mass, the energy scales behave as in the m q = 0 case discussed above. In particular, the Miransky scaling law (3.5) applies for x c − x 1.
B The "scaling" regime which involves walking of the coupling, and the amount of walking is determined by the quark mass. Continuity of the tachyon implies for x ≤ x c that (see Appendix A for details) (3.10) Notice that (3.9) is obtained from (3.10) by setting ∆ * = 2, which is indeed the value of the anomalous dimension as x → x c from above, and the scaling behavior is therefore continuous and both expressions are valid at x = x c . Moreover, (3.10) gives what is termed the "hyperscaling" relation for the chiral condensate (see, e.g., [35]) and often written in terms of the anomalous dimensions γ * = ∆ * − 1, Such scalings have been studied recently in a specific holographic model [34]. The behavior of Λ τ is the same as for m q = 0 and in the regime A, i.e., Notice that we have not discussed the m q -dependence of bound state masses and decay constants, because they are difficult to analyze without specifying the details of the model. In the regime A and at small x, however, it is clear that the lowest bound state masses (except for the light pions) and decay constants must be O (Λ IR ) because this is the only available scale (apart from the small perturbation due to the quark mass). The pions are expected to obey the GOR relation since this arises in holography quite in general. In the regime B, as well as in the regime A when x c − x 1, the most natural expectation is that the masses and decay constants are still given by the soft IR scale O (Λ IR ). This will be demonstrated for V-QCD below.

Large quark mass
Let us then discuss the limit m q /Λ UV 1 (regime C in Fig. 2). Since the quark mass is large, the tachyon grows large very close to the UV boundary where the dilaton is still small. Therefore regime C probes the limit of small dilaton and large tachyon in the holographic model. Due to the smallness of the dilaton it is not obvious that the holographic description is reliable for all observables in this regime. But in analogy to how the UV structure of V-QCD is fixed in order to guarantee correct boundary conditions for the IR dynamics (as explained in Sec. 2), it is important that the holographic model is as close to QCD as possible also at large quark mass, in order to have the best possible boundary conditions for the physics at small and O (1) quark masses. We will therefore analyze what are the possibilities in this limit.
Deep in the UV, for r 1/m q , the tachyon will have the standard form of (3.1) which already implies scaling results for σ and Λ τ , with similar assumptions as above. At r ∼ 1/m q it will grow large, after which nonlinear effects are important and a generic solution cannot be written down. In particular, the RG flow is not expected to approach the IRFP at any point, and the IR structure of the tachyon solution shown above is not relevant.
Instead, additional scaling results can be derived assuming that the holographic model implements some key features of QCD (as will be the case in V-QCD). These are the RG flow of the 't Hooft coupling in QCD, and the decoupling of the massive quarks at energy scales much smaller than m q . Due to the decoupling, the RG flow is that of full QCD for r 1/m q and that of YM theory for r 1/m q . Continuity of the 't Hooft coupling (rather than the tachyon) is required.
The collected results for the regime C in Fig Here σ will be connected to the properly renormalized chiral condensate (see Appendix C for details). The renormalization also involves scheme dependence which is important at large quark mass (see [46] for a discussion in the context of holography). The condensate is proportional to m 3 q as in (3.13) for generic schemes. Further, analysis of the RG flow leads to where b 0 (b YM 0 ) is the leading coefficient of the beta function for QCD in the Veneziano limit (YM theory at large N c ). That is, Let us then briefly comment on the size of the bound state masses in regime C. First, as the quarks are decoupled for energy scales smaller than the quark mass, the glueballs are expected to decouple from the mesons, and have a mass gap of O (Λ IR ). Recall that the meson states in QCD become nonrelativistic at large m q , and therefore their mass gap is ∼ 2m q , and the mass splitting of the low-lying states is much smaller than the gap. The example of V-QCD, which we shall discuss below, shows that obtaining such a mass gap in holography is nontrivial (see the analysis for V-QCD in Appendix B), and the gap can be either O (m q ) or O (Λ IR ) for actions which produce reasonable results at small m q . If an IR cutoff is placed at the point where the tachyon grows large (as in the dynamic AdS/QCD models [33]), mass gap ∝ m q is obtained. In the presence of such a cutoff, Λ τ ∼ m q is the only scale in the system. But in this case the splitting between the bound state masses is also expected to be O (m q ).
Finally let us comment on the structure of Fig. 2 near the BZ region. Namely, the boundary of the regime C bends toward higher m q /Λ UV as x grows. Similarly the upper boundary of the regime B bends down in the BZ region. This is a generic feature due to RG flow which becomes slower and slower in the BZ limit x → x BZ . Due to slowness of the flow the separation of the energy scales m q and Λ UV needs to be larger for the scaling results to apply. This phenomenon is discussed in more detail in Appendix A.
Actually the BZ regime could be identified as an additional scaling regime. Here this is not done, however, since from the point of view of holography this regime is not the most interesting one. The coupling constant is restricted from above by its small value at the IRFP, the theory is perturbatively soluble, and therefore holographic description is not expected to be useful.

Scaling of bound state masses in V-QCD
While we discussed above the generic behavior of (holographic) QCD in the Veneziano limit, we shall now derive explicit predictions for the V-QCD models defined in Sec. 2, and demonstrate that they agree with the results of the previous section.

Energy scales in V-QCD
The various energy scales can be defined explicitly in terms of the background solutions. For definiteness we will write down the definitions here. We will use the UV and IR expansions which can be found in Appendix D of [14] and in Appendix D of [23]. The UV expansions can also be found in Appendix C.3.
• First, Λ UV is the scale of the UV expansions. The precise definition is most conveniently given in terms of the UV expansion of the dilaton λ: (4.1) Recall that the coefficients of the potentials in the V-QCD action were matched to those of the QCD beta function in the UV. We used this mapping to write the coefficients in the expansion (4.1) in terms of the coefficients b i of the QCD beta function in the Veneziano limit, • The IR scale Λ IR can be defined analogously in term of the IR expansions.
For the standard geometry in V-QCD with an IR singularity the definition of Λ IR = 1/R is given in limit r → ∞ as [14] log λ = 3 2 In the presence of an IRFP the definition is modified to where δ = ∆ F F − 4 is the anomalous dimension of the TrF 2 operator at the fixed point.
• The quark mass is determined 6 in terms of the UV asymptotics of the tachyon (r → 0), where is the UV AdS radius and ρ = γ 0 /β 0 is the ratio of the leading coefficients β 0 and γ 0 of the beta function and the anomalous dimension of the quark mass in QCD, respectively. Notice that in V-QCD, the logarithmic running of the quark mass is therefore included, and is matched to agree with that of QCD.
• The scale Λ τ where the tachyon grows large is defined simply by 7 where τ and A are to be evaluated on the background solution. 6 It is well known that the quark mass can be defined only up to a constant. This constant is set to one for notational simplicity. 7 For potentials II slightly modified definition is used: the number on the right hand side is set to 1/10. This is necessary because the transition to the IR region where nonlinear terms in the tachyon are important turns out to happen when the tachyon is still smaller than one, and in the nonlinear region the tachyon grows relatively slow. Therefore using exactly (4.6) would lead to an energy scale which does not precisely reflect the change in the dynamics.

The background EoMs for V-QCD are invariant under the transformation
This transformation changes all energy scales defined above by the same number, and reflects the choice of the units of energy on the field theory side. Consequently, only the ratios of the above energy scales are a priori well-defined.
There is a small issue with the above definitions which will be visible in the explicit results below. Namely, the definition of Λ UV which is natural at generic values of x is not optimal at high x and in particular in the BZ limit. In this limit one would expect Λ UV to match the scale of the UV RG flow, but this turns out not be be the case. Also at zero quark mass there should be only a single scale and therefore Λ UV should equal Λ IR , but as it turns out, actually Λ UV becomes exponentially suppressed with respect to Λ IR . Because the RG flow of the coupling is controlled by the two-loop beta function in the BZ limit, a single scale Λ can be defined which has better behavior in this limit (see Appendix B of [52] This scale is related to Λ UV in the BZ limit by grows large in the BZ limit (see Appendix A for some more details). Notice that Λ UV will anyhow be used as a reference scale in the numerical analysis below -since the BZ region is not very interesting from holographic viewpoint, we have chosen to use the same definition of Λ UV as earlier literature, even though it has unnatural behavior in this region.

Flavor nonsinglet masses and decay constants: scaling results
Before going to the numerical results let us argue how the scaling laws for the masses and decay constants can be derived in V-QCD. The meson masses in each sector (for vectors, axial vectors, scalars, and pseudoscalars) may be defined into two classes: the flavor singlet and nonsinglet states. The former are singlets under the vectorial SU (N f ) transformation, whereas the latter are the other fluctuation modes, which correspond to quark bilinear operators involving the Hermitean traceless generators t a of SU (N f ). We restrict analytic considerations to the nonsinglet states, because the singlet mesons mix with the glueball states, which would lead to complications [23,53]. Only the main points are summarized here and the details can be found in Appendix B.
The fluctuation equations for the scalar nonsinglet mesons can transformed in the Schrödinger form as detailed in Appendices A and B of [23]. The masses for each sector are then given as eigenvalues of a Schrödinger equation with a certain potential term and the Schrödinger coordinate running from u = 0 (UV) to u = ∞ (IR). In order to find the behavior of the mass gaps and splittings one then needs to study the Schrödinger potentials V S (u).
Let us first take finite but small quark mass (m q /Λ UV 1), so that the regimes A and B are covered. We use here the V-QCD action, but results in this region are not sensitive to the details of the action. The Schrödinger potential for u 1/Λ IR is given by the background with small or walking dilaton and small tachyon, so that the geometry is close to AdS, and consequently V S (u) ∼ const/u 2 . For u 1/Λ IR the diverging tachyon creates the confining potential (V s (u) ∼ Λ 4 IR u 2 if the excitation spectrum is linear, m 2 n ∼ n with n being the excitation number). Therefore the only relevant scale is Λ IR and the (lowest lying) meson masses as well as the mass splittings between the modes are given by this scale: which is also consistent with the results from Dyson-Schwinger and Bethe-Salpeter approaches in regime A [54]. In units of Λ UV , by using the scaling laws from Sec. 3 we therefore obtain the same results as at vanishing quark mass [22,23] within regime A: m n ∼ Λ UV in the QCD regime and the masses go to zero obeying Miransky scaling, , in the walking regime. In regime B, we find the "hyperscaling relations" [34,35] for the low-lying meson masses: The sole exception to these scaling results is the pion mode (for x < x c ), which is massless at m q = 0 and obeys the GOR relation in regime A: so that it is lighter than the other meson states. The GOR relation will be discussed in more detail in Sec. 5. At the level of the Schrödinger formalism the absence of the mass gap in the pseudoscalar sector is reflected in the negativity of the Schrödinger potential in the UV region. In regime B, the pseudoscalar masses will also obey the scaling (4.10).
The various decay constants are slightly more difficult to analyze. In Appendix B it is argued that they are similarly of the order of Λ IR when the quark mass is small. At large quark mass (regime C) the meson masses and decay constants depends more on the details of the action. Actually only the form of the flavor action S f in the limit of large tachyon and small dilaton is relevant, as is shown in Appendix B. In this limit the functions a, κ, and V f 0 become almost constants and the action of (2.5) takes the form of the standard DBI action with the exponential Sen potential.
Interestingly, as shown in Appendix B, the choice V f (λ, τ ) = V f 0 (λ) exp(−a 0 τ 2 ) of potentials I, where the function a(λ) is set to a constant value a 0 , is a special case. This choice was motivated by the asymptotics of the meson trajectories at zero quark mass [14,23], and noticeably a(λ) is also constant in tachyon potentials obtained from string theory [16,17,55]. In Appendix B we show that this choice is essentially the only one which produces physically reasonable mass gap and splitting, if the value of a 0 is slightly modified from the value of potentials I which reproduces the correct dimension of the quark mass and condensate in the UV.
We will anyhow discuss the results for potentials I without this modification of a 0 , because such a modification was not introduced for our numerical studies. It is found that the mass gap for the mesons in all sectors (vectors, axials, scalars and pseudoscalars) is given by where is the UV AdS radius, and ξ = 1 would be required to match with the field theory result for nonrelativistic bound states 8 . We stress that this formula was obtained by fixing the value of a 0 to produce the UV dimension of the quark mass and the condensate, as was done for the potentials I used in the numerics. The last expression in (4.13) suggests that ξ = 1 could be obtained by tuning the value of W 0 . This is, however, problematic because the contribution involving W 0 should vanish in the probe limit x → 0, and also because W 0 would negative, which causes V f to have a node [14] at which the tachyon background equation becomes singular. The mass splitting is suppressed with respect to the gap as is the case in real quarkonia. For potentials I it scales as the inverse of the mass gap. The decay constants of the lowest meson states are exponentially suppressed in the quark mass and they are therefore decoupled. This is found for all mesons with masses below O (m q ). The decay constants of the states with masses ∼ m q are of the order of Λ IR .

Numerical results
Computing the energy scales and masses numerically is straightforward (but tedious) after the background has been constructed numerically [14,23]. The potentials I  Figure 3: Dependence of the energy scales on x and m q . Also the masses of the rho meson m ρ , the lowest singlet scalar m ss , and the pions m π are shown for reference. The blue solid curves is Λ UV , the red dashed curve is Λ τ , the thin dotted magenta curve is the rho mass, the thin dotdashed green curve is the pions mass, and the thin long-dashed brown curve is the mass of the lowest singlet scalar in each of the plots. See text for more details.
(given explicitly in Appendix G) are used here unless stated otherwise. We choose three reference values x = 1, 4, and 4.5, which lie in the QCD regime, walking regime, and conformal window, respectively, and plot the observables as functions of the quark mass. We also show plots where m q /Λ UV is fixed to 10 −6 and x is varied over the whole parameter space. These choices cover the most interesting structures of Fig. 2.

Energy scales
In Fig. 3 the dependence of the energy scales on x and m q is demonstrated. We have chosen to show the scales Λ UV and Λ τ in the units of Λ IR since this makes the details visible. Some of the bound state masses are also shown as thin lines.
The top-left plot shows the dependence of the scales on x at a tiny quark mass (m q /Λ UV = 10 −6 ). The data extend only up to x = 5.1 because in the BZ region it is difficult to do reliable numerics. The top-right plot shows the mass dependence in the running regime (x = 1). The bottom-left plot is in the walking regime (x = 4, which is close to x c 4.0830). The bottom-right plot is in the conformal window (x = 4.5).
The thick solid blue curves show the ratio Λ UV /Λ IR in each plot. In the top-left plot the crossover from the QCD-like regime to conformal window is clearly visible: As x → x c from below, the ratio first grows according to the Miransky scaling law (3.5) in regime A until the condition (3.8) no longer holds. Thereafter the ratio saturates to roughly Λ UV /m q in regime B as predicted by (3.9), and then one moves out of regime B at even higher x.
The dependence of Λ UV /Λ IR on m q also follows the predicted scaling laws. In the regime A (low m q /Λ UV in top-right and bottom-left plots) the ratio is constant as the quark mass is a small perturbation. In the regime C, i.e., at large m q /Λ UV , the ratio follows the power law of (3.14) as best seen in the top-right plot. In the regime B the ratio obeys the other power law of (3.9) and (3.10) as best seen at intermediate m q in the bottom-left plot. Notice that much of the solid blue curve was left out in the bottom row plots in order to make the details of the other curves better visible.
The thick dashed red curves show the ratio Λ τ /Λ IR . At small quark mass, including regimes A and B, the ratio is close to one as expected (except in the BZ limit). For the explanation of the divergence of Λ τ /Λ IR in the BZ limit see Appendix A, Eq. (A.15). At large quark mass (regime C), we find Λ τ ∼ m q as predicted in (3.13).
At high x (plots in the bottom row) the convergence toward this scaling law is quite slow due to the slow running of the quark mass. It could be demonstrated by continuing the plots up to much larger m q /Λ UV , but we have chosen not to do so in order to show the other details in the plots more clearly. Actually all variables vary slower and slower as functions of m q /Λ UV when x is grows, and therefore we have substantially increased the range of m q /Λ UV in the plots with higher x, but this is still not enough to demonstrate the large m q scaling convincingly. The change in the m q -dependence is due to the RG flow (see the end of Appendix A for some more details).
The dependence of σ on m q will be discussed in Sec. 5.

Flavor nonsinglet masses and decay constants
The flavor nonsinglet spectra can be computed as explained in [23]. The results for different values of x and m q are shown in Figs. 3, 4, and 5. First, the ρ and π masses are shown in Fig. 3 with thin dotted magenta and thin dotdashed green curves, respectively. As expected, the ρ mass is O (Λ IR ) at small quark mass (regimes A and B) and obeys the power law (4.13) in regime C. The pion mass is close to the ρ mass in regimes B and C, but in regime A it obeys the GOR relation instead, as best visible from the top-right and bottom-left plots at small m q .
In the top-left plot, m π /Λ IR is O m q /Λ UV at small x but then increases with x  in the walking regime (actually obeying the Miransky scaling law) until the ratio is O (1) at the point of the crossover near x = x c . Fig. 4 shows the masses of the lowest meson states (i.e., the mass gaps) in each sector normalized to the ρ mass. The choices of x and m q are the same as in Fig. 3. The masses of the lowest vector, axial, scalar, and pseudoscalar states are given by the solid blue, dashed red, dotted magenta, and dotdashed green curves, respectively (the brown curves give the scalar singlet mass gap which will be discussed below). These ratios are mostly constant and close to one as predicted by the above scaling arguments. The exception is the pion mass (dotdashed green curves) which obeys the GOR relation in regime A. Notice that in regime C all meson mass gaps should approach the same number (roughly 2m q ) as expected for nonrelativistic bound states, but we find instead that the axial and pseudoscalar gaps are larger than those of the vectors and scalars. The reasons for this are analyzed in Appendix B. Notice also that the lowest scalar states are lighter than vectors even at small values of x, which seems to be in conflict with QCD. Such details are, however, sensitive to the choice of the potentials in the V-QCD action, and can be changed by tuning the potentials.
Finally we plot the masses of the four lowest vector states in Fig. 5 in order to demonstrate the dependence of the mass splittings in the spectra on m q . The masses are given in units of Λ IR (left hand plot) and normalized to the lowest mass, i.e., the ρ mass (right hand plot). The splittings decrease with increasing m q in regime C which is in qualitative agreement with the bound states becoming nonrelativistic, but the power laws are not exactly correct (see Appendix B for details).

Scalar singlet masses
The singlet sector is qualitatively different from the nonsinglet sector because it also contains glueball states which mix nontrivially with the singlet meson states. Such mixing takes place in the scalar and pseudoscalar sectors. The scalar sector will be discussed in detail here while the pseudoscalar sector will be analyzed in a future publication [41].
Before going to the numerical results, let us discuss the generic features of the spectrum. The singlet mesons masses are expected to show similar 9 m q and x dependence as the nonsinglet mesons above, given by (4.10) and (4.13). In the singlet case there is, however, nontrivial mixing of the mesons with the glueball states, the masses of which should be independent of m q and therefore always characterized by Λ IR . In particular in regime C the masses of the mesons become much larger than the glueballs which suggest that the meson and glueball states decouple.
Let us then demonstrate these features numerically for the scalar singlet sector in V-QCD. The mass gap of the scalar singlets is shown by the long-dashed brown for all values of m q and x as seen from Fig. 3. The fact that the lowest glueball mass is suppressed with respect to the rho mass m ρ in regime C can also be seen in Fig. 4: the ratios m ss /mρ given by the brown curves, where m ss is the mass gap for the scalar singlet states, decrease with m q at large m q . The decoupling of the meson and glueball states is most clearly demonstrated by Fig. 6, where the masses of the five lowest scalar singlet states are plotted as a function of m q in the QCD regime (x = 1, top-left), in the walking regime (x = 4, top-right) and in the conformal window (x = 4.5, bottom). When m q /Λ UV 1 the spectrum has both glueball and meson states which are nontrivially mixed. As m q increases, the meson masses increase while the glueball masses stay constant, which leads to the crossing structure seen in the plots. This would be expected to happen at m q /Λ UV ∼ 1, but as was discussed above, Λ UV is not exactly the scale of the UV RG flow when x is large, and therefore the crossing structure shifts to higher values of m q as x increases. At large m q /Λ UV only the glueballs are left. Their decoupling from the mesons is demonstrated by the fact that the limiting values of their masses as m q → ∞ are independent of x, and in fact it can be checked that they match with the glueball masses obtained in the YM limit (x → 0) of V-QCD.
It was shown in [22,23] that V-QCD does not have a light "technidilaton" mode [8]  (which would be the lightest scalar singlet state) as x → x c . Both the singlet and nonsinglet scalars fluctuations do have critical [28] behavior in the near conformal region (see Appendix I in [23]) and the Schrödinger potential for the nonsinglet scalars is negative (as was also found in [32]), but it was shown numerically that this is not enough for a technidilaton to appear. The negative result is also seen in Figs. 4 and 6: the lowest singlet scalar does not become light with respect to the other states in the walking regime (x = 4). We do see, however, that it is lighter in regime A than in regime B.

Masses near the conformal transition
To conclude this section, let us add a few comments on the scaling of the bound state masses near the conformal transition. We plot the masses of the rho meson, the pion, and the lowest scalar singlet state in Fig. 7. The left hand plot shows the masses in units of Λ IR and the right hand plot shows the masses in units of Λ UV . Notice that in units of Λ IR the pion mass deviates from the masses of the other states at extremely small m q (regime A) where it obeys the GOR relation, but in units of Λ UV the pion mass obeys the same power law ∼ m q /Λ UV in regimes A and B.
Knowledge of the dependence of the meson masses on m q near the walking regime is important for the lattice studies which take place at finite quark mass, and aim to locate the conformal transition at m q = 0 [13,36]. As we have pointed out above, the regime B extends even to x < x c (see Fig. 2), and the crossover between regimes A and B moves to lower x as m q is increased. This suggests that studies at finite m q lead to an underestimate for x c . One should recall, however, that the scaling in regime B involves ∆ * = γ * + 1 which depends strongly on x. We show Reγ * , which controls the scaling exponents of the masses, as a function of x for both potentials I and II in Fig. 8. The kinks in the plots are located exactly at x = x c , and γ * drops rapidly right above the kinks. This supports the idea that x = x c can be located by extracting γ * from the meson masses on the lattice. In fact, recent lattice results for γ * in the conformal window report very low values [36,57] that are in apparent contradiction of the curves in Fig. 8. Recall however that the model has not been tuned to fit any QCD data yet 10 . It is actually not difficult to construct potentials for which γ * drops much more rapidly when x > x c .

Quark mass and the chiral condensate
Let us then discuss the mass dependence of the chiral condensate. Recall that the tachyon solution in the UV reads where σ can be identified as the chiral condensate 11 (the exact identification is studied in Appendix C.1), and ρ can be expressed in terms of the leading coefficients of the beta and gamma functions as ρ = γ 0 /b 0 = 9/(22 − 4x).
Notice that the analysis of previous sections was restricted to the standard, dominant vacuum. It is known [14,23], however, that there are subdominant "Efimov" vacua in the QCD and walking regimes (x < x c ) which quite in general appear in connection to the BKT transition (see, e.g., [29,30,38]). These vacua are mapped to different values of the chiral condensate on the (m q , σ)-plane: all possible regular vacua form a spiral structure, which will be called the "Efimov spiral" below. The results for the spiral structure will be used to analyze four-fermion deformations of QCD in Sec. 6.

Efimov spirals
Let us first review the structure of the subdominant vacua. Including the solutions with finite quark mass [14]: • When x c ≤ x < x BZ , only one vacuum exists, even at finite quark mass.
• When 0 < x < x c and the quark mass is zero, there is an infinite tower of (unstable) Efimov vacua in addition to the standard, dominant solution 12 .
• When 0 < x < x c and the quark mass is nonzero, there is an even number (possibly zero) of Efimov vacua. The number of vacua increases with decreasing quark mass for fixed x.
The infinite tower of Efimov vacua, which appears at zero quark mass, admits a natural enumeration n = 1, 2, 3, . . . where n is the number of tachyon nodes of the background solution as we shall demonstrate below (see also Sec. 10 and Appendix H in [14]). A generic feature of these backgrounds is, that they "walk" more than the dominant, standard vacuum, so that the scales Λ UV and Λ IR become well separated for all 0 < x < x c when n is large enough. It is possible to show that for any 0 < x < x c . The coefficient ν will be given below in Eq. (5.8) (see also Appendix F in [14]). In the walking regime, one finds that for any value of n. In particular, n = 0 corresponds to the standard solution discussed in the previous sections, and the relation (5.3) gives the standard Miransky scaling, whereas for n > 0 the scaling is even faster. We also found a similar scaling result for the free energies of the solutions as x → x c in [14], therefore proving that the Efimov vacua are indeed subdominant in this limit, and verified this numerically for all 0 < x < x c . In [23] it was shown that the Efimov vacua are perturbatively unstable (again analytically as x → x c − , and numerically for all 0 < x < x c ). When the BF bound is violated at the IRFP, the quark mass and the condensate are known to show an oscillating behavior for the (chirally broken) backgrounds where the coupling flows very close to the fixed point [14]. Let us first discuss how these oscillations arise from the tachyon EoM.
First, take a background at zero tachyon which reaches the fixed point as r → ∞. Then consider turning on an"infinitesimal"tachyon. It satisfies the linearized tachyon EoM where λ * is the value of the coupling at the fixed point and * is the IR AdS radius.
in terms of the dilaton and tachyon potentials. The BF bound is thus given by When the BF bound is violated we denote In this case the asymptotic infinitesimal tachyon solution is oscillatory, Let us consider next specific tachyon solutions τ m and τ σ , where either an infinitesimal quark mass m q or a condensate σ is turned on in the UV, respectively. These solutions are conveniently expressed in units of Λ UV , and will have the asymptotics of (5.9). We denote where the coefficients K i and φ i cannot be computed analytically but it is easy to extract them from numerical solutions. One can require that K i > 0 and −π/2 ≤ φ i < 3π/2. We are, however, interested in the solutions where the tachyon is small and finite. In this case the tachyon will eventually grow large when r ∼ 1/Λ τ ∼ 1/Λ IR , and drive the system away from the fixed point. The above formulas (5.10) and (5.11) then hold as approximations for 1/Λ UV r 1/Λ IR . The solution for r 1/Λ IR depend on the details of the model in the IR. However, in order to satisfy the boundary conditions imposed by the good IR singularity, the tachyon must have certain fixed normalization and phase when expressed in IR units: In general the IR asymptotics of the background depends on one parameter (e.g., T 0 for potentials I [14]) but as the fixed point is approached, the dependence on this parameter appears only through the ratio Λ UV /Λ IR whereas K IR and φ IR take fixed values (see Appendix I in [23]).
As the fixed point is approached, the result (5.12) must match with the sum of (5.10) and (5.11). Therefore one finds From here one can solve As u varies the equations (5.14) define a spiral on the (m q , σ)-plane, which has been studied recently at finite chemical potential in a different context (see [38]). Notice that ν does not need to be small. It can be verified numerically that the handedness of the spiral is such that its phase increases with increasing u (counter clockwise direction) on the (m q , σ)-plane. This means that As we shall show in Appendix C.1, this is also required in order for the standard solution (with nonzero and nodeless tachyon) to be dominant (for x < x c so that the BF bound is violated and the spiral exists).
Finally, let us point out some properties of the spiral as x → x c from below. Notice from (5.8 . The approximations (5.10) and (5.11) are valid for 1/Λ UV r 1/Λ IR , but they should join smoothly with the UV asymptotics of the tachyon in (5.1). At r ∼ 1/Λ UV we find that (5.17) These estimates, and similar estimates for the solution (5.11), can be matched with the UV asymptotic formulas if An analogous argument shows that for (5.12) to satisfy generic IR boundary conditions, We have found numerically that (with the convention (and φ m − φ σ > 0 such that (5.16) holds) for all potentials which we have studied.
Then the first node of the mass in Eqs. (5.14) in their regime of validity u 1 occurs at This node is identified as the "standard" solution, where the tachyon has no nodes (in particular no nodes appear in the regime of validity of (5.10), (5.11), and (5.12)). Notice that solving Λ UV /Λ IR from (5.21) results in Miransky scaling 13 of (3.5). Nodes at larger values of u, i.e., νu (n+1)π, with n = 1, 2, . . . are identified as the Efimov solutions, where the tachyon has n nodes.

Numerical results
As an example, we have computed m q and σ numerically for potentials II with SB normalization 14 for W 0 at x = 2. The results are shown in Fig. 9. The blue dots are our data and the red curves are given by Eqs. (5.14). The blue line on the right hand side is a power-law fit σ −Cm 3 q . The spiral structure is not well visible on the left hand plot because the distance of the curve from the origin decreases exponentially with increasing u in (5.14). In order to make the spiral structure visible, we have 13 In principle it could be possible to satisfy the boundary conditions with . We speculate that this happens in the model of [58] where no Miransky scaling was found.
14 It is much easier to extract the Efimov spiral numerically for potentials II than for potentials I. plotted the logarithms of (the absolute values of) σ and m q on the right hand side of Fig. 9. Since the action is symmetric under τ → −τ there is actually also another spiral which not shown in the left hand plot but can be obtained simply by a rotation of 180 degrees around the origin. Notice that when m q = 0 it is difficult to define σ unambiguously in practical calculations, because the vacuum expectation value (vev) solution of the tachyon cannot be separated from the subleading terms of the source solution. This means that we also have to specify more carefully how K m and φ m in (5.10) are defined in our numerical analysis: we pick a reference solution which has σ = 0 by definition, and consequently defines K m and φ m (see also Appendix C.1). It is natural to expect that the σ = 0 solution has a quark mass of O (Λ UV ) in the QCD regime. Therefore the reference solution was fixed to be the standard, dominant vacuum solution with m q /Λ UV = 1. This kind choice is important in order to avoid unnatural fine tuning effects.
The various coefficients in the solutions in (5.14) were determined as follows. Three solutions were chosen for the background for which σ and m q are very small such that (5.9) holds as a good approximation for at least two periods of oscillation. The first (second) solution was tuned to have approximately zero σ (m q ) and was used to fit the constants K i , φ i in (5.10) [in (5.11)]. The third solution was used to fit the constants K IR and φ IR in (5.12).
In the walking regime, extra care is needed in the choice of the reference solution having σ = 0. The chiral condensate is expected to have a node around such values of the quark mass where the normalizable and nonnormalizable terms in the tachyon solutions are nontrivially coupled. This happens for the standard vacuum at the crossover point between the regimes A and B, i.e., for for The spiral in the walking regime (x = 3.6, close to x c 3.7001) is show on the loglog scale in Fig. 10 (left). By studying numerically various quantities (for example the S-parameter, plotted below in Fig. 20) it is found that the choice m q /Λ UV = 3×10 −11 lies at the crossover and is therefore chosen as the reference point with σ = 0.
As the critical value x c is approached, the Efimov spiral becomes "squeezed", as seen by comparing plots in Fig. 9 (right) and 10 (left). In the walking regime, the consecutive solutions with vanishing m q and σ are close so that the spiral approaches a straight line in the log-log scale. This reflects our analytic results for the spiral (5.14) derived above in the limit x → x c − : in particular φ m − φ σ ∼ √ x c − x → 0. We have also repeated the analysis in the conformal window. In this case the scaling of (3.10) is expected to hold for small quark mass and that σ ∼ m 3 q for large quark mass. The results for potentials II with SB normalization for W 0 at x = 4 are shown in Fig. 10 (right). The red dots are the data while the blue curves are given by the power laws mentioned above.
Again we need to specify how σ is extracted in the numerical analysis because the tachyon vev solution cannot be separated from the subleading terms of the source solution. In the conformal window, however, the solution is simple because the values of σ grow fast with increasing quark mass. Therefore one can choose any reference solution "with σ = 0" at very small m q , and the results are essentially independent of the choice.
Notice that the plots in the walking regime (left) and conformal window (right) do not seem too different. Indeed the curves undergo a smooth transition at x = x c .
The coefficients of the power laws are continuous at the transition, and the nodes where σ = 0 or m q = 0 approach the origin of the spiral very fast as x → x c from below. Actually the rate of the approach is given by the Miransky scaling law.

Gell-Mann-Oakes-Renner relation
The GOR relation can be obtained as usual by combining the results from two computations. The first result is the expression for the pion mass at small m q , which is obtained by analyzing the fluctuation equations at small m q as done in Appendices D and E. One finds that 15 where κ 0 = κ(λ = 0) and W 0 = V f 0 (λ = 0). The second result is the relation between σ and the chiral condensate as m q → 0, which is obtained by deriving the renormalized on-shell action (i.e., the vacuum energy E) with respect to m q (see Appendix C.1): The combination is the GOR relation: It can be checked that the proportionality coefficient (here minus one) is correct for our normalization of f 2 π , which differs by the factor N f /2 from the standard normalization in chiral perturbation theory (see, e.g., [59]).
Notice that 1. One might expect that the logarithmic terms which appear in the UV expansion of the fields would result in correction that are only suppressed by 1/ log m q in (5.23). Remarkably, the logarithmic terms completely cancel, as they should, since the QCD the relation holds up to linear corrections in the quark mass. Why this happens is shown in Appendix D. The cancellation of these corrections is also consistent with the combination m q qq being RG invariant [60].
2. The derivative in (5.24) is nontrivial, since changing the quark mass also affects the geometry even at m q = 0 due to the full backreaction between the flavor and glue sectors, which may add contributions to the derivative (see Appendix C.1). 3. The relation is valid only in regime A: the correction terms in (5.25) become large at the crossover between regimes A and B.
We have tested the GOR relation numerically for potentials I at x = 1. Both sides 16 of the relation are plotted as functions of the quark mass in Fig. 11 (left) and good agreement is found. The subleading correction in Fig. 11 (right) is clearly linear in the quark mass as it should.

Four-fermion operators in V-QCD
After we have constructed the solutions on the (m q , σ)-plane it is straightforward to analyze the effect of multitrace deformations following the recipe of [40]. In the presence of such deformations, the coefficient of the source term, which was denoted by m q above, is no longer trivially related to the quark mass. Let us therefore denote this coefficient by α m in this section. Similarly, the coefficient of the normalizable term is denoted by β m .
One can now study the following extra terms on the field theory Lagrangian: Then the boundary conditions are obtained by setting the source term α m to the value [40] (see also [61,62]) and the vev is given by The possible vacua can then be identified by overlapping this condition with the curves of regular solutions on the (α m , β m )-plane (shown above in Figs. 9 and 10).
We are most interested in the case of double trace deformation, g 2 = 0 with other couplings equal to zero, since this operator becomes marginal at the critical point x = x c . Let us also set m q = 0. In this case The overlap plot is shown for the phase with the Efimov spiral (0 < x < x c ) in Fig. 12. Notice that the second branch of the spiral, obtained by reflection about the origin, which was omitted in earlier plots (for example Fig. 9) was now also included 18 . For this special case where g 2 is the only nonzero coupling, as the mapping (α m , β m ) → (m q , σ) in (6.3) and (6.4) implies, the change of the UV boundary conditions with respect to the standard case g 2 = 0 corresponds to changing the vertical axis to the line defined by (6.5) (examples are the red and dashed magenta lines in Fig. 12) while the horizontal axis is kept fixed. Notice that the solutions with g 2 = 0 lie on the vertical axis, and are denoted by the green line in Fig. 12. The green dot shows the (stable) standard solution, but there is also an infinite tower of additional intersection points near the origin, which are not visible as the spiral converges very fast towards zero. These intersection points give the Efimov solutions.
In order to draw the phase diagram at nonzero g 2 , we need to solve the free energy of each solution and find the dominant vacuum. We can start from the identity In Appendix C.1 we show that the conditions (6.3) and (6.4) are indeed consistent with (6.6) and that the higher order expectation values satisfy in agreement with the large N factorization of expectation values. Recall from Sec. 5 that the Efimov spiral can be parametrized as Notice that if only g 2 is nonzero, m q and σ still satisfy (6.8) (with α m (β m ) replaced by m q (σ), respectively), if the coefficients K i and φ i are redefined. This is seen by inserting the conditions (6.3) and (6.4), and is consistent with the change of boundary conditions simply corresponding to a new choice of axes in Fig. 12. By inserting this parametrization in (6.6) and integrating along the Efimov spiral we find that for the solutions with m q = 0 (see Appendix C.1) where E 0 is the free energy of the solution with α m = 0 = β m . Since sin(φ m −φ σ ) > 0, the dominant solution from those in the range of validity of (6.8) is that with the largest value of u. In the walking regime this can be seen explicitly. Namely, we argued in Sec. 5 that in the walking regime ν π (x c − x)/K → 0, and that the solutions are found at νu (n + 1)π, with n = 0, 1, 2, . . .. Therefore 11) which agrees with the scaling of the free energy found in [14] (see Sec. 10 and Appendix H there).
Also quite in general the solution with largest |σ| is the dominant vacuum. From (6.6) we see that the energy density is given in terms of the (oriented) area between the spiral and the horizontal axis. For clockwise oriented spirals the minimum energy is reached at the solutions furthest away from the origin.
It is then straightforward to construct the phase diagram, which is shown in Fig. 13 for the case of m q = 0 and nontrivial g 2 . For x < x c the diagram can be found simply by analyzing the solutions in Fig. 12. When g 2 > 0 the configuration is qualitatively similar to that at g 2 = 0: the dominant vacuum for x c < x is the vacuum with the standard tachyon solution having no nodes. When g 2 < 0, the dominant solution is that of Fig. 12 (right) having sizeable |σ|. As this solution is absent for g 2 ≥ 0, there is a discontinuity (a "zeroth order" transition) at g 2 = 0−. It is smoothly connected, though, to the dominant solution of g 2 > 0 through the limit g 2 → ±∞ (where the constraint (6.5) gives a horizontal line in the plots of Fig. 12). There is also a subdominant solution shown in Fig. 12 (left), namely the continuation of the standard, dominant solution at g 2 = 0 to negative g 2 . We will argue below that this solution is metastable for small |g 2 | and unstable for large |g 2 |.
Notice that the chiral condensate is a scheme dependent quantity, and the scheme dependence is important, in particular, at large m q (see [46] for a discussion of the scheme dependence in the context of holography). Therefore the dominant solution g 2 < 0 appears scheme dependent at least for small |g 2 |. We have, however, found the behavior of generic schemes σ ∝ m 3 q with a negative proportionality constant for all choices of potentials we have tried, and independently of the precise definition of σ. This is enough to guarantee that the phase diagram is that of Fig. 13. Actually, as stressed in Appendix B, the results at very large m q are essentially independent of the choices of the various functions in the V-QCD action.
In the conformal window, the situation is even simpler. For g 2 ≥ 0 there is only the solution with zero quark mass and the condensate. When g 2 < 0, there is again a solution with large |σ|, analogous to that shown in Fig. 12 (right). When g 2 > 0, the transition at x = x c is similar as at g 2 = 0, i.e., a BKT transition, since changing the value of g 2 only amounts to changing the axes in Fig. 12 without affecting the structure of the spiral.
One can also show that when a finite m q is turned on, the BKT transition and the chirally symmetric phase disappear, but the discontinuity at g 2 = 0 − remains.

Perturbative stability
Finally perturbative stability of the solutions with modified boundary conditions could be analyzed following [23]. Here we will only discuss which of the solutions are expected to be unstable, and will not prove the stability of any solution. Naturally, the modification of the boundary conditions for the background also implies that the boundary conditions of the fluctuations are similarly changed as the fluctuations must preserve the physical value of m q .
Let us first analyze any solutions with walking, i.e., u 1. Recall that the standard solution has been shown to be stable, while the Efimov solutions are unstable when g 2 = 0 [23]. The instability appeared in the scalar flavor singlet and nonsinglet sectors. In order to study stability of the other solutions, one should look at the scalar fluctuation equations. In the walking regime, they admit simple solutions in the UV and in the vicinity of the (approximate IR) fixed point. In fact, as argued in Appendix I of [23], the fluctuation equations (for sufficiently small mass of the fluctuation) take the same form as the EoMs for the background. This is true both in the flavor singlet and in the nonsinglet sectors. Therefore, the fluctuations in the vicinity of the (approximate) IRFP are analogous to (3.3): where ψ S is the radial wave function of any scalar fluctuation mode. The solution necessarily has a node if this approximation is valid for more than half a period of the sine function. In terms of variable u = log(Λ UV /Λ IR ) the node therefore appears for νu π (where we used the fact that Λ IR ∼ Λ τ whenever walking is present). Such a node of the wave function (say at zero momentum) implies an instability, for generic UV boundary conditions for the fluctuations. It is straightforward to prove this in the case of nonsinglet scalars for which the fluctuation equations can be cast into the Schrödinger form, and with the standard UV boundary conditions such that the fluctuation wave function is normalizable in the UV (see also the analysis in [63]). We shall sketch the proof here.
Denote the UV normalizable (but not necessarily IR normalizable) Schrödinger wave function by φ and the location of any of its nodes by r 0 . By studying the variation of with respect to the mass, we find that We see that all nodes move towards the IR as m 2 is lowered. But as m 2 → −∞ the solution for generic UV boundary conditions is φ ∝ exp(|m|r) and has no nodes for r 1/m. Therefore all nodes must disappear either by moving to r = ∞ or by merging with other nodes. However, the Schrödinger equation does not admit solutions with a double node for regular potentials V S , so merging of the nodes is not possible, and consequently all nodes must disappear by moving to the far IR. In particular, if there the wave function has a node when m 2 = 0, the node must move to r = ∞ at some negative value of m 2 . At this value, φ is IR normalizable: otherwise it could not have a node in the far IR if m 2 is slightly perturbed. Tachyonic normalizable mode marks the presence of instability. Putting the above observations together, we conclude that a node of the m 2 = 0 wave function marks the presence of an instability.
The above proof does not apply directly to our case because the UV boundary conditions are modified, and consequently the integral in (6.14) is divergent. The divergence can be regulated by introducing a UV cutoff at r = , but this results in extra counterterms on the right hand side of (6.14) and its negativity is no longer obvious.
The proof can be fixed, however, in our case (i.e., scalar fluctuations and node in the region of validity of (6.12)) when x → x c from below. This is because the node becomes well separated from the UV region: the location of the node satisfies log(r 0 /Λ UV ) ∼ 1/ √ x c − x as seen from (6.12), where ν ∼ √ x c − x. For generic boundary conditions, the right hand side of (6.14) is then dominated by the contributions to the integral near the node. The counterterms which cancel the divergence of the integral are essentially independent of x for any reasonable boundary conditions, and therefore negligible. The rest of the proof remains unchanged. Let us then study the solutions with m q = 0 as x → x c . Recall that the phase differences φ IR − φ σ − π and φ σ − φ m will be ∝ √ x c − x in this limit as we argued in Sec. 5. The standard solution (the green dot on the vertical axis in Fig. 12) is in the regime of validity of the approximations leading to (6.8). It is found at Based on the above analysis, the onset of the perturbative instability is also expected As the standard solution is stable, the critical value of u must be larger than that given in (6.15). Inserting this in (6.8), the critical value of g 2 is given by the ratio α m /(β m c σ ). The factors of √ x c − x cancel leaving a O (1/Λ 2 UV ) number, which must be negative given the handedness of the spiral. The critical value of g 2 is therefore of the same order as the value corresponding to the dashed magenta line in Fig. 12. The parts of the spiral which are closer to the origin from the critical points (near the intersection points with magenta dots) are unstable.
The above analysis cannot be applied in the QCD regime where x c −x is not small. This regime could be analyzed numerically. The natural expectation is that the critical value is still negative and O (1/Λ 2 UV ). This is supported by the fact that the Efimov vacua are unstable, which has been shown by a numerical computation [23].
To conclude this section, the phase diagram is that shown in Fig. 13: For g 2 > 0, the diagram is identical to that of g 2 = 0. The result is somewhat different from the result obtained in the gauged Nambu-Jona-Lasinio model, where the four-fermion operator is also slightly different, and the lower bound of the conformal window appears to increase with g 2 [64]. As any negative g 2 is turned on, the "standard" vacuum immediately becomes unstable and the dominant vacuum has much larger |σ|. For x < x c and g 2 < 0, the perturbation analysis suggests the the standard vacuum is metastable for small |g 2 | and becomes perturbatively unstable for larger |g 2 |, with the critical value being O (1/Λ 2 UV ).

S-parameter and current-current correlators
The vector-vector and axial-axial correlators have the structure and t a are the generators of SU (N f ) with the normalization Tr t a t b = δ ab /2. The factors 2/N f on the right hand side were added to ensure that Π V and Π A are proportional to N f . The numerical factor was chosen such that these factor are equal to one when N f = 2, which is the smallest number for which the flavor nonsinglet currents are defined. This will result in the standard normalization of the S-parameter. The sign convention is η µν = diag(−, +, +, +). Notice that Π L vanishes for zero quark mass because then ∂ µ J a (A) µ = 0. When the spectrum is discrete (i.e., 0 < x < x c or m q is finite), we may formally write the correlators as sums over the contributions from the meson states: (7.6) Depending on the choice of the holographic action, these series may not converge. This issue will be discussed below in the context of the S-parameter. The residues S 0 of the "spurious" q 2 = 0 pole in (7.4) and (7.6) must identical for the pole to cancel in the full correlator. At zero quark mass also S 0 must be related to the pion decay constant, S 0 = f 2 π ≡ (f (P ) 1 ) 2 , since Π L vanishes. The difference of the vector-vector and axial-axial correlators involves the quantity which is nontrivial only when chiral symmetry is broken. The expansion of D(q 2 ) at q 2 = 0 defines the S-parameter: Here S 0 is the same coefficient which appears in (7.4) and (7.6), and we will also study the higher order coefficient S 2 below.

Correlators and the S-parameter in V-QCD
Let us then recall how the correlators can be computed in V-QCD. As pointed out in Sec. 2, the vector currents are dual to the gauge fields in the DBI action (2.5).
Following the standard approach (see [23] for additional details), we carry out the fluctuation analysis writing down an Ansatz which separates the spatial and radial dependence of the fluctuation modes in momentum space. The spatial (radial) wave functions of the vector, transverse axial, and longitudinal axial modes are denoted by V, A and P (ψ V , ψ A , and ψ L ). The radial wave functions are IR normalizable and satisfy the UV boundary conditions Here ψ P is the radial pion wave function. The terms of the on-shell V-QCD action which are quadratic in the vector fields can then be written as 10) where projects to the transverse parts of the wave functions. To compute the vector-vector correlators, we need the precise dictionary given in terms of the couplings to field theory currents: 14) Applying the gauge/gravity correspondence with these couplings leads to the following expressions for the form factors: where I = V, A, L. Notice that for m q > 0 (so that Π L is nonzero) the wave functions ψ A and ψ L satisfy the same fluctuation equations as q 2 → 0. Consequently lim q→0 q 2 Π A (q 2 ) = lim q→0 q 2 Π L (q 2 ) (7.17) and this number equals S 0 of (7.4) and (7.6) which is therefore well defined. This equality ensures the cancellation of the "pole" at q 2 = 0 as pointed out above.  Figure 14:

S N c N f x c x BZ
The normalized S-parameter as a function of x for m q /Λ UV = 10 −6 (blue solid curves) and m q = 0 (red dashed curves). Left: potentials I with W 0 = 3/11. Right: potentials II with SB normalized W 0 .
By inserting (7.16) in the definitions (7.7) and (7.8) we obtain for the S-parameter This formula is, however, not convenient for high precision numerical computations, because the subleading terms at r = are only suppressed by logarithms of in V-QCD, so that extremely small values of would be needed to obtain accurate results. By an analysis of the fluctuation equations, it is possible to derive more convenient integral representations for S (see Appendix E). We find that where u is the Schrödinger coordinate defined in Appendix E, and that The formula for the S-parameter is well-known in the context of simple bottom-up models (see, e.g., [65]). In Appendix E we write it in a form which holds for very generic holographic models. We used the fact that ψ V | q 2 =0 = 1 (7.22) in order to obtain the last expression (7.21).

Numerical results for the S-parameter and f π
Using equations (7.16), (7.19), and (7.21) it is straightforward to compute the Sparameter, the pion decay constant as well as the coefficients S 0 and S 2 of (7.8) numerically (see [23] for additional details, and [26,66,67] for analysis of the Sparameter in other holographic models). Figure 14 shows the results for S at m q = 0 (dashed red curves) and at m q /Λ UV = 10 −6 (blue curves) for potentials I (left) and II (right). The numerical value of M 3 was fixed such that the asymptotics of the vector-vector correlator matches with perturbative QCD (see Eq. (C.10) in [23]). The most striking feature in these plots is the discontinuity of the S-parameter in the conformal window. When x ≥ x c , the S-parameter immediately jumps from zero to a O (N f N c ) number when any finite m q is turned on. The mechanism which leads to this discontinuity will be discussed in detail below (from the holographic viewpoint), but it appears rather natural: the S-parameter is O (N f N c ) whenever the geometry has the IR singularity, and vanishes only for zero quark mass in the conformal window where there is an IRFP instead. The result is also consistent with the analysis based on field theory at qualitative level [68]: the S-parameter is finite except for exactly zero mass in the conformal window.
There is, however, one striking difference [22,23] with respect to previous results: the S-parameter increases with x in regime A, whereas many earlier analyses [26,[69][70][71] suggest that the S-parameter is suppressed in the walking regime and may even vanish as x → x c . Recall, however, that the IR behavior of the potentials in the V-QCD action has not yet been fitted to QCD or lattice data, and such fits may affect the x-dependence of the S-parameter. By analyzing the form of the fluctuation wave functions (see Appendix G in [23]) in the integral formula (7.21) one can indeed check that it is dominated in the IR (for small m q and for all values of x). This is consistent with the analysis of Appendix I of [23], where the same is argued to hold for the meson masses and decay constants.
We have also computed the pion decay constant which is defined in terms of the residue of Π L at the pion mass whenever m q = 0 (see Eq. (E. 21) in Appendix E). The results for both potentials are given in Fig. 15, and they are also compared 19 to the constant S 0 . As expected the pion decay constant (blue curves) match with √ S 0 (dashed red curves) in the QCD regime. The ratio f 2 π /S 0 decreases fast with increasing x for both potentials when x x c , suggesting that the pion decouples.
The higher order coefficient S 2 (in units of Λ IR ) is also shown for potentials I in Fig. 16 (left) and the product S 0 S 2 in Fig. 16 (right). The dependence on x is similar as for the S-parameter when x x c . 19 Notice that S 0 is UV divergent whenever the quark mass is finite, as can be seen by inserting the UV expansions of the wave functions [23] in (7.16), and needs to be renormalized. The divergence is ∝ m 2 q in agreement with the one-loop field theory computation of Appendix F. At small quark masses it is irrelevant how the renormalization is done because the difference between all reasonable renormalization schemes is negligible due to the smallness of the coefficient in the divergent term. We have not tried to analyze the various observables in the BZ limit x → x BZ because this region is not the most interesting one from a holographic viewpoint. In general, however, the slow RG flow in the BZ limit causes that all ratios of energy scales to be typically ∝ exp[#(x BZ − x) −2 ]. Therefore even ratios which are expected to be close to one for generic values of x easily blow up in the BZ limit -in this case # in the above relation is small, but it is difficult to define the scales such that it would be exactly zero. In view of this, it is not surprising that only the dimensionless quantities S and S 0 S 2 approach finite values in the BZ region.
Let us then analyze the mass dependence of f π and S in more detail. for all values of x = 0 (apart from the discontinuity at m q = 0 which is only present in the conformal window and the fact that S varies slower as x increases which is due to the RG flow as discussed at the end of Appendix A). In particular the limiting value as m q → ∞ is independent of x. In this limit the S-parameter is expected to approach the value N c N f /12π from perturbative QCD (see Appendix F). Even though V-QCD is not expected to reproduce perturbative results in general, the limiting value in V-QCD is numerically close to the QCD number 1/12π 0.0265.
The dependence of f π on m q is demonstrated in Fig. 18. Again the different plots are in the QCD regime (x = 1, top left plot), in the walking regime (x = 4, top right plot), and in the conformal window (x = 4.5, bottom plot). For small m q the dependence is weak, but when m q /Λ UV 1 the decay constant vanishes very fast with increasing m q . This signals the decoupling of the pion mode in regime C, and is consistent with the findings of Appendix B (see Eq. (B.45)). Actually, all low-lying meson states are expected to decouple for potentials I.

Scaling of the S-parameter
We have also analyzed the S-parameter in the limit m q → 0. Numerical results are shown in the conformal window (for x = 4.5) in Fig. 19, where the red dots are our data and the lines are given by the functions where the parameters β 1,2 , as well as S( ), were fitted to the data. Here it is under-  stood that which is a finite number whereas the S-parameter vanishes at zero quark mass: S(0) = 0. Thus there indeed is a discontinuity at m q = 0. For both potentials β 2 0.08 fits the data very well. In order to understand the power law in (7.23), it is useful to first discuss in more detail how the discontinuity at m q = 0 arises. The mechanism is the same which was studied in detail in the case of m q = 0 and x → x c in section 6 and Appendix I of [23]. Here instead x ≥ x c and m q → 0. In both cases the RG flow of the coupling approaches the fixed point (λ = λ * ) but misses it finally (due to the finite tachyon). For such flows, it is useful to divide the background to the UV and IR sections, having λ < λ * and λ > λ * , respectively. Considering flows which get closer and closer to the fixed point, the S-parameter can be computed more and more precisely in terms of the IR section. The IR part takes a fixed shape in this limit, explaining the finite value of the S-parameter. For exactly zero quark mass the IR section of the background becomes disconnected from the UV section and is therefore not present in the physical vacuum solution, which now ends at the IRFP. This is reflected in the vanishing value of the S-parameter. By using similar arguments, we can also sketch how the power law in the mass dependence arises. It is understood to be the leading "perturbation" of the IR background due to the fact that the fixed point was not reached exactly. In the conformal window, the flow towards the fixed point is given by where δ is related to the dimension of the TrF 2 operator at the fixed point. It can be computed as the derivative of the holographic beta function at the fixed point (when the tachyon is set to zero exactly) [14]. One finds that where V i are the coefficients of the Taylor expansion of V eff at λ = λ * : Notice that V 2 < 0. Flow toward the fixed point ends when r ∼ 1/Λ IR . The difference of λ with respect to the fixed point value when this happens is given by where (A.14) was used to obtain the last expression. This difference controls the deviation of the IR section of the background from its limiting shape as m q → 0 and correspondingly the deviation of the S-parameter from the limiting value S(0 + ). Indeed by using the explicit expressions for the potentials one obtains δ ∆ * 0.0780 , (x = 4.5) , (7.29) a value in a very good agreement with the fit from above, shown in Fig. 19. Finally notice that the difference λ * − λ IR not only controls the corrections to the S-parameter at small mass, but also to other quantities that can be defined in terms of the IR section of the background. Examples are decay constants and meson masses in IR units, which are therefore expected to have qualitatively similar m q dependence to the S-parameter at small masses 20 .
The mass dependence of the S-parameter can be analyzed similarly in the walking regime (x → x c − ). The above calculation is approximately valid in regime B, i.e, when the quark mass controls the amount of walking. One can use (3.9) together with (7.28) to obtain Therefore one should effectively take ∆ * → 2 in (7.28), meaning that the power is continuous over the conformal transition at x = x c , as ∆ * = 2 at the transition. For even smaller m q , i.e., in regime A, the mass term of the tachyon can be treated as a linear perturbation to the whole background, and therefore the mass dependence of the S-parameter is linear. Both the scaling of (7.28) and the linear dependence can be seen in the numerical results in the walking regime, see Fig. 20.

Convergence of the spectral representation
Finally we will analyze the spectral representation of the S-parameter in order to understand better why it increases with x in the QCD and walking regimes. By inserting (7.4) and (7.5) in the definition for S, one obtains (7.31) This series may, however, be ill defined. For example, the decay constant approach asymptotically constant values whereas the masses have linear trajectories m 2 n ∼ n for potentials I (see Appendices E and F in [23]). One can check that (7.31) is convergent thanks to the asymptotic cancellation of the vector and axial terms, but it is not absolutely convergent, and therefore the result may depend on the ordering of the terms. The definition of (7.31), where the states are ordered in terms of their excitation numbers n, would work for potentials with linear trajectories if the slopes of the vector and axial spectra are the same. This is not 21 the case for potentials I, and consequently (7.31) is incorrect. The solution to this issue is simple: the contributions should be ordered by the meson masses rather than by the excitation numbers.
It is straightforward to verify numerically that (7.31) converges towards the Sparameter if the terms are ordered according to the masses. It turns out to be convenient to define a mass dependent cutoff, A smooth cutoff function was chosen instead of a step function because it improves convergence drastically. The convergence also means that the value of the S-parameter is determined through the dynamics in the deep IR, because the same holds for the masses and decay constants (the argument can be made precise in the walking regime, see Appendix I in [23]). The convergence of the regulated series (7.32) towards the S-parameter is demonstrated in Fig. (21). The resolution of the cutoff function was fixed to the mass difference of the two lowest vector states: δm = m 1 . The speed of the convergence is best visible from the right hand plot, which shows as a function of the cutoff. The convergence is exponential for all values of x, but becomes significantly slower as x increases and one moves from the QCD regime to the conformal window. The slowness of the convergence is not due to changes in the spectra. To show this, the filled circles and boxes marking the masses of the lowest five vector and axial mesons, respectively, were added in each plot. It is seen that the spectrum changes relatively little with x in unit of Λ IR . Finally let us try to extract the reason for the increase of S with increasing x in the region of low values of x from the plots of Fig. 21. Notice that the curves for x = 1, 3.5 and 4 essentially overlap at small values of m cut in the left hand plot. The curves for x = 1 and x = 3.5 deviate from that of x = 4 as m cut increases, and after deviating rapidly saturate to the final value of S. Therefore it appears that contributions to the S-parameter are roughly mass-independent up to a saturation scale, which increases with x. In order to have a good estimate for the S-parameter, a growing number of terms need to be included in the sum with increasing x, whereas the individual terms in the sum are of roughly constant size. Therefore the increase of S with x in the QCD and walking regimes can be seen to be due to slower convergence of the sum. In the conformal window, i.e., for the curve with x = 4.5, something different happens. The convergence of the sum is even slower, but the contributions at fixed m cut are suppressed, resulting in the decrease of the S-parameter with increasing x.

Finite temperature phase diagram
The finite temperature phase diagram has been studied in detail for IHQCD in [73] and for V-QCD at zero quark mass and at small values of the quark mass in [18]. Here this study is extended to large values of m q as well as very high values of x. The code for constructing solutions at finite temperature is available at [74]. Some extra tricks are necessary in order to obtain reliable results in the BZ region and at very large m q (see Appendix G).
First recall the generic structure of the (x, T ) phase diagram [18], which was already reviewed in Sec. 2. At zero quark mass, there is a first order deconfinement transition in the QCD and walking regimes, but there is also the possibility (depending on the choice of potentials and the value of x) of a chiral symmetry restoration at a separate second order transition. In the conformal window, there is a continuous phase transition at zero temperature. When a finite quark mass is turned on, chiral symmetry is always broken, and the second order chiral transition will become a fast crossover when the quark mass is small, and completely disappear at larger quark masses. The system is in a tachyonic thermal gas (TG) phase at small temperatures for all 0 < x < x BZ . As the system is heated, there is a first order deconfinement transition to the high temperature phase, which is implemented through a transition from TG phase to the black hole (BH) phase in holography.
The existence of the deconfinement transition requires an order parameter. While QCD at finite N c and N f has no order parameter related to deconfinement, the pressure acts as an effective order parameter at large N c . This is clear in the 't Hooft limit, where the number of degrees of freedom (and consequently the pressure) is O (N 0 c ) in the low temperate phase and O (N 2 c ) in the high temperature phase. In the Veneziano limit the number of degrees of freedom is of the same order in both phases, but the phase transition may still be identified as a discontinuity of the pressure. In fact, the pressure of the model is still exactly zero in the TG phase, because our approach does not capture the contributions corresponding to loops of pions (as well as mesons with higher masses) in this phase. Including these contributions in the model would affect the critical temperature, and potentially even alter the order of the transition [21].

Scaling laws at finite temperature
The critical temperature has nontrivial dependence on m q , which can be analyzed analytically. Notice that the temperature brings in an additional energy scale with respect to the zero temperature solutions. Another difference is that the definition of the standard reference scale Λ IR cannot be extended to the BH solutions in a natural manner (because the geometry now ends at a horizon rather than an IR singularity). Therefore it is understood that Λ IR is defined below through the TG solution (or equivalently through the zero temperature solution at the same values of m q and x).
Let us first discuss the mass dependence of the critical temperature, which can be inferred by using the results from [18] and from Sec. 3. The temperature of the black hole can be related to the metric through the formula where A h and r h are the values of the scale factor and the bulk coordinate at the horizon. In V-QCD models T (r h ) has a nontrivial minimum (for tachyonic BHs) and the transition takes place at the scale of the minimum. Indeed, the entropy density decreases monotonically (and fast) with r h , as suggested by the UV and IR (zero temperature) expansions of A(r) and as can be verified numerically. Further, the geometry of the BH solution approaches smoothly the TG solution as r h → ∞, and therefore p TG = lim r h →∞ p BH . Integrating p BH (r h ) = s(r h )T (r h ), a node p BH = 0, and consequently a first order phase transition, is found near the minimum of T (r h ).
In conclusion, one should locate the minimum of T (r h ) in order to determine the scaling of T c . When the geometry is close to AdS, (8.1) implies that T ∼ 1/r h . This result holds both in the UV asymptotic region and when there is an approximate IRFP, i.e., walking. When the quark mass is large, there is also an approximately AdS region where the flavors are already decoupled but the dilaton λ is still small. In summary, For r h 1/Λ IR the temperature increases with r h as seen by studying the IR expansions (see [18]). Therefore, the minimum of T (r h ) takes place at r h ∼ 1/Λ IR . By continuity, (8.3) for all m q > 0 and 0 < x < x BZ . The scaling results in units of Λ UV immediately follow by using the results from Sec. 3: and the dependence on m q is expected to be a linear perturbation. The dependence of T c on x at m q = 0 [18] is in qualitative agreement with analysis based on field theory (see, e.g., [75]).
• In regime B, • In regime C, In addition to the phase transition, also various crossovers can be identified as the maxima of the interaction measure As it turns out, such crossovers reflect the different regions of the zero temperature geometry. This can be understood by approximating the horizon as a sharp cutoff added on the zero temperature background. Substituting an AdS metric in the formulas (8.1) and (8.2), one finds that the interaction measure vanishes. First, there is the crossover which marks the transition from the quasi-conformal or walking phase (with approximate IRFP) to the asymptotic UV phase [18,76]. Such a crossover is found whenever there is walking, i.e., in regime B, and the x → x c edge of regime A. The UV asymptotics is valid for r 1/Λ UV , and the flow from the UV fixed point to the IRFP is characterized by Λ UV . Consequently, using (8.1) and (8.3), the crossover temperature is expected to be T co,qc ∼ Λ UV (8.10) independently of m q . Second, there is a crossover at large quark mass, corresponding to the transition from the region where the quarks are decoupled to the UV asymptotic region. The decoupling of the quarks takes place at r ∼ 1/m q as pointed out in Sec. 3. Consequently, the crossover temperature is given in terms of the quark mass: When T m q the quarks are effectively decoupled, and therefore the thermodynamics is the same as for pure YM (that is, the x → 0 limit of V-QCD). Notice that Λ UV , however, is defined in terms of the UV asymptotics, i.e., effectively at infinite energy, and different from that of YM even as m q → ∞: for YM, Λ UV ∼ Λ IR but in the limit of large m q these scales are related through (3.14) instead. Consequently, in order for the thermodynamics to smoothly approach YM thermodynamics as m q → ∞, dimensional quantities should be expressed in units of Λ IR or T c rather than in units of Λ UV .

Numerical results
Let us then illustrate the dependence of the various critical temperatures on x and m q numerically. The basic features of the phase diagram at small quark mass are demonstrated in the top-left plot of Fig. 22. The first order "deconfinement" transition temperatures for potentials I at zero and at tiny (10 −6 ) quark mass are shown as functions of x on the logarithmic scale. The first order transitions are shown as blue curves. The curves overlap at small x, but as the conformal transition is approached, the curves become separated. The lower curve (which overlaps with the red dashed curve and is therefore not well visible) is the transition temperature at m q = 0 which goes to zero with Miransky scaling as x → x c . For x > x c there is no transition when m q = 0. When a tiny quark mass is turned on the transition (upper blue curve) is present for all values of x. the critical temperature T c decreases with x inside the conformal window 22 . The dashed thick red curve is the second order chiral restoration transition, which also shows Miransky scaling as x → x c and is absent for x > x c . This transition only exists for m q = 0 in the walking regime.
The dependence of T c on m q is also demonstrated numerically for V-QCD in Fig. 22. The top-right, bottom-left, and bottom-right plots are in the running regime (x = 1), walking regime (x = 4), and in the conformal window (x = 4.5), respectively, and T c is given by the blue curve in each plot. The power laws are in agreement with the above formulas. For example, in the regime C with large m q we find that the exponent of (8.8) is 1 − b 0 /b YM 0 = 2x/11 0.182 at x = 1 and 0.727 at x = 4, which 22 In [18] also a region where T c increased with x was seen at high x 4.5 (see Fig. 27 there), but this effect turns out to be due to the UV cutoff being too low in the numerics. is consistent with the plots. The crossover between the quasiconformal and UV regions at T = T co,qc is seen as the horizontal lines in the bottom row of Fig. 22. The ratio T co,qc /Λ UV is constant as expected, but the value of the constant deviates significantly from its expected value, i.e., one. This happens because Λ UV deviates from the scale of the UV RG flow at large values of x, as was explained above in Sec. 4. The crossover due to the decoupling of the quarks at large m q at the temperature T ∼ m q is best visible in the top-right plot of Fig. 22.
In addition, in large part of the parameter space there is also a separate maximum of the interaction measure in connection to the first order transition. This kind of maxima have also been included as thin red curves in Fig. 22, and can be found close to the blue curves denoting the transition temperatures.
Recall that in regime C the mass gap of the mesons does not have the expected behavior m gap 2m q in V-QCD (if Sen-like tachyon potential is assumed). Nevertheless, the crossover due to the decoupling of the quarks takes place correctly at T co,mq ∼ m q . This is seen as a consequence of the decoupling of the mesons having an unphysically low mass, which was discussed above in Sec. 4.
Finally let us study in some details of the thermodynamic functions at very small and very large m q . We plot the (normalized) energy density, pressure, and interaction measure for potentials I with x = 4 in Fig. 23 (left). The solid blue, dashed red, and dotted magenta curves have m q /Λ UV = 0, 10 −12 , and 10 −10 , respectively. The value of x = 4 was chosen such that the model shows the separate second order chiral restoration transition, which appears as a kink in the energy density and the interaction measure of the functions at zero quark mass. The vertical thin dashed black lines mark the locations of the phase transitions. As a tiny quark mass (m q /Λ UV = 10 −12 ) is turned on, the second order transition turns into a crossover. Indeed the dashed curves follow closely the solid curves except for very close to the kink, where the curves have a smooth behavior instead. At m q /Λ UV = 10 −10 a much larger deviation is seen already.
The approach to YM theory can be seen by studying the thermodynamics at very large quark mass. Thermodynamic functions at sizeable m q are shown for the same choice of potentials and compared to the results for the limit of YM theory (x → 0) in Fig. 23. As the quark mass increases, the quarks are decoupled and the thermodynamics near T = T c is expected to converge to that of YM theory. Indeed, the (dashed red) curves for m q /Λ UV = 10 10 lie essentially on top of the (solid blue) curves of the YM thermodynamics 23 . For smaller quark mass (m q /Λ UV = 10 8 , dotted magenta curves) a much large deviation is seen already.
where * and m τ, * are the AdS radius and tachyon mass at the IRFP, and ν = Im∆ * . In the deep IR, i.e., for r 1/Λ τ , the tachyon and the dilaton have calculable IR asymptotics which depend on the details of the action [23]. These asymptotics shall not be needed here.

A.2 Scaling results
The dependence of the various scales on the quark mass can now be found by requiring continuity (and the continuity of the derivative) of the tachyon solution. More precisely, one should require that the dominant and subdominant tachyon solutions are both continuous, but in the cases studied here this is equal to requiring that continuity of the derivatives.

A.2.1 QCD regime
The QCD regime is defined by 0 < x < x c and x c − x 1. Let us first consider the case of small mass, m q /Λ UV 1. In this case there is (almost) spontaneous chiral symmetry breaking, which in V-QCD means that the (vev term of the) tachyon grows large at some value of the radial coordinate (energy scale) and triggers a nontrivial IR geometry [14]. The fields obey their IR asymptotics exactly when the tachyon is large. Therefore both the UV and IR scales are comparable to the scale of the tachyon, Λ UV ∼ Λ IR ∼ Λ τ , and σ ∼ Λ 3 UV . In the dynamic AdS/QCD models [33] similar results are obtained by introducing an IR cutoff where the tachyon grows large. From (A.3) one indeed sees that the mass term of the tachyon is suppressed with respect to the vev term for r 1/Λ UV which implies that m q can be treated as small perturbation.
When the quark mass is large, m q /Λ UV 1, the tachyon grows large at small r ∼ 1/m q and therefore Λ τ ∼ m q . The RG flow of coupling is determined by the QCD beta function for r 1/Λ τ and by the YM beta function for r 1/Λ τ as the growing tachyon decouples the flavor sector. Explicitly 24 , < 1, requiring continuity leads to the counterintuitive result that the scale of the IR expansions is larger than that of the UV expansions: This is however consistent with the fact that decoupling the flavor degrees of freedom makes λ to run faster than it does at small m q .

A.2.2 Walking regime
The walking regime is defined by x < x c and x c − x 1. For very small quark mass it is expected that the scaling is the same as for m q = 0 and that Λ τ ∼ Λ IR -the IR geometry is linked to the growth of the tachyon as in the QCD regime. Requiring continuity 25 with generic IR boundary conditions at r ∼ 1/Λ IR gives C w ∼ Λ 2 IR /Λ 2 UV in (A.4). Continuity of (A.3) and (A.4) at r ∼ 1/Λ UV further leads to σ ∼ Λ 2 IR Λ UV , and from (A.3) it is seen that the quark mass term is small if where we used (3.5). When the amount of walking is controlled by the quark mass. We still require that the IR geometry is exactly where the tachyon is large so that Λ τ ∼ Λ IR and C w ∼ Λ 2 IR /Λ 2 UV , but now continuity of (A.3) and (A.4) at r ∼ 1/Λ UV yields and there is no Miransky scaling. When m q /Λ UV 1, the coupling never flows close to the fixed point and the walking behavior is therefore absent. One finds the same results as in the QCD regime above.

A.2.3 Conformal window
Conformal window is the regime with x c < x < x BZ . Let us first assume that we are not in the perturbative BZ regime, x BZ − x 1. When the quark mass is small, m q /Λ UV 1, one again requires that Λ τ ∼ Λ IR and that there is walking. Continuity at r ∼ 1/Λ IR gives for the coefficients of (A.5). Matching the leading and subleading terms of (A.3) and (A.5) at r ∼ 1/Λ UV one further obtains (A.14) When m q /Λ UV 1, the flow does not become close to the fixed point, and the results are the same as in the QCD and walking regimes.

A.2.4 BZ regime
Let us then discuss the BZ regime (x < x BZ and x BZ − x 1). One might expect that Λ τ ∼ m q independently of the value of m q since ∆ * approaches one in this limit. We have, however, defined m q asymptotically in the UV, and the UV RG flow of the quark mass is singular in the BZ limit. Therefore m q and Λ τ are not simply related as x → x BZ -see the top-left plot of Fig. 3, where Λ τ increases with x in the BZ region instead of approaching m q . Therefore, we will use the scale Λ τ in the analysis below.
When the quark mass is small, the coupling flows toward the BZ fixed point when r 1/Λ τ , and like in YM when r 1/Λ τ . Since λ * is small, one can use (A.8) to describe the YM flow. Requiring the flow to start at λ * when r ∼ 1/Λ τ gives the exponential scaling In the opposite limit, i.e., large quark mass, the reasoning leading to (A.9) applies, if one replaces m q by Λ τ in the formulas. There is, however, also a subtlety in the definition of Λ UV . In [14], Λ UV was defined essentially as the scale where b 0 λ becomes O (1). However in the BZ regime the coupling reaches the fixed point well before reaching the value 1/b 0 ∼ 1/(x BZ − x). Consequently, Λ UV is exponentially suppressed with respect to the true characteristic scale of the RG flow in the UV. Let us instead denote by Λ UV the scale defined as [52] (roughly corresponding to the scale where λ/b 0 becomes O (1)): This formula is the RG flow given by the two-loop BZ beta function, and it is reproduced in V-QCD up to correction suppressed by x BZ − x (for the whole flow when τ = 0 and in the UV for all backgrounds), whereas the geometry is AdS, A − log(r/ ), up to highly suppressed O ((x BZ − x) 4 ) corrections. Therefore (A.9) is a better estimate after the replacements Λ UV → Λ UV and m q → Λ τ . In order to take into account the extremely slow RG flow in the BZ region, the conditions for the validity of (A.15) and (A.9) should actually be written as respectively. Here b 1 is the NLO coefficient of the QCD beta function and b 2 0 /b 1 ∼ (x BZ − x) 2 . Notice that the slowness of the RG flow is already visible at not so large values of x (as demonstrated by the plots in the text, see, e.g., Fig. 17), because we use parameters m q and Λ UV defined asymptotically in the UV. Indeed the equations (7.28), and (A.9) depend on the ratio m q /Λ UV through combinations of the type (m q /Λ UV ) O(x BZ −x) .

B. Schrödinger potentials and mass scales
Let us then analyse the behavior of mass gaps and mass splittings (separation of the lowest bound state masses) in V-QCD. For the flavor nonsinglet fluctuations this is rather straightforward, as the fluctuation equations can be transformed into the Schrödinger form. The singlet fluctuations are more involved, because there is in general nontrivial mixing between the meson and glueball states. Only in the probe limit x → 0 and in the limit of large quark mass m q → ∞ the glueballs and mesons are decoupled. We will restrict here to the nonsinglet states and study the scalar singlet states only numerically in Sec. 4.3.3.
For every flavor nonsinglet sector (vectors, axial vectors, pseudoscalars, and scalars) the fluctuation equation can be written as where V S (u) is the Schrödinger potential, and m 2 is the mass of the fluctuation. This form is obtained after a coordinate transformation defined by and u(r = 0) = 0. The Schrödinger potential is given by where Ξ(u) and H(u) are different functions for each sector and can be explicitly expressed in terms of the potentials (see Appendix A in [23]). For example, in the vector sector we find that and for the axial vectors Ξ A = Ξ V but The vector/axial decay constants are given by (see Appendix F in [23]) Notice that by using (B.3) the Schrödinger equation (B.1) may be written as Integrating this over u and inserting in (B.6), we obtain (B.8)

B.1 Small quark mass
Let us then analyze the dependence of the bound state masses on the quark mass in the three regimes of Fig. 2. In regime A, the quark mass is a small perturbation and the results from [23] can be used directly. The potential V S can be computed both when r 1/Λ UV and r 1/Λ IR by using the asymptotic expansions of the background. The results in the UV are independent of the potentials of the V-QCD action: where v UV = −1/4 for the pseudoscalars and 3/4 for other sectors, and we used the fact that u r in the UV region.
In the IR the asymptotics of the tachyon, and consequently the coordinate dependence of the Schrödinger potential, strongly on the choice for the potential functions in the V-QCD action. All regular choices considered in Appendices D and E of [23] lead to a confining potential, which grows as a function of u in the IR regime.
If Λ UV ∼ Λ IR , as is the case at small x and m q , we immediately notice that the Schrödinger potential has its bottom for all sectors expect for pseudoscalars at r ∼ 1/Λ UV ∼ 1/Λ IR . The value of the Schrödinger potential is O (Λ 2 IR ) near the bottom, as can be estimated from (B.9) above by requiring continuity at u ∼ 1/Λ IR , and therefore the mass gap is ∼ Λ IR . Similarly it can be seen that the mass splittings and vector/axial decay constants are O (Λ IR ). This is rather evident as only one energy scale enters the definitions above, both the tachyon and λ are O (1) when r ∼ 1/Λ IR , and no cancellations are expected in the formulas.
The pseudoscalar sector is special because of the pion modes which obey the GOR relation, as shown in Sec. 5.3. The excited pseudoscalar states appear at mass scale Λ IR , as can be seen numerically.
In order to study regime B and the remaining part of regime A (near x = x c where walking is seen) we need to check the Schrödinger potential for 1/Λ UV r 1/Λ IR , i.e., in the near-conformal region. The tachyon remains small also here and r u. The Schrödinger potentials have been derived in [23] and read where v w = 3/4 for vectors and axials, and depends on the anomalous dimension ∆ * at the (approximate) fixed point for scalars and pseudoscalars. For scalars v w < −1/4 when x < x c , whereas it is positive for pseudoscalars. Therefore the coefficient in the scalar potential is critical [28], and might potentially lead to an instability or a light state (see Sec. 5.3 of [23]). Numerically it is seen, however, that this is not the case and the spectrum of scalars is not qualitatively different from that of the vectors, for example.
We see from (B.10) that V S has similar dependence when the coupling constant walks as in the UV region in (B.9), for the vector and axial states. Therefore the bottom of the Schrödinger potential is at u ∼ 1/Λ IR . Similar arguments as above show that the mass gap, mass splitting, and decay constants are O (Λ IR ). As we pointed out above, for the scalars and pseudoscalars more careful or numerical analysis is needed. The numerical result is that the masses are similarly O (Λ IR ), with the exception of pion masses in regime A, which obey the GOR relation. In regime B, pions also have masses O (Λ IR ).
In conclusion, mass gaps, mass splittings, and decays constants are O (Λ IR ) at small quark mass, with the sole exception of the pions in regime A.

B.2 Large quark mass
Next we shall study the dependence on the (flavor nonsinglet) meson mass spectrum on the quark mass in the regime C (m q /Λ UV 1). Recall that in this limit some features of the bound state can be analyzed starting from field theory, because the low end of the spectrum becomes nonrelativistic. The expected mass gap is roughly equal to 2m q , and the states can be studied by using the Schrödinger equation. For QCD states at large mass, one expects to have two main contributions in the Schrödinger potential. First, at very short distances one has Coulomb potential as perturbative gluon exchange dominates. Second, there is a nonperturbative confining potential which is expected to be linear in the distance and arise from a flux tube between the quarks. While it is difficult to derive such a linear potential from first principles, it is consistent with the observed quarkonium spectra, lattice simulations, and also found in holographic calculations.
In the limit m q → ∞, the lowest states are therefore governed by the Coulomb potential, and one finds typical Hydrogen-like spectrum with negative binding energies. For slightly higher states, the linear potential dominates. Assuming precisely linear confining potential ∼ Λ 2 IR r, where Λ IR is the scale of glueballs in QCD, it is straightforward to solve the Schrödinger equation. The masses of the states obey roughly the scaling law For charmonium and bottomonium the observed states are in the region where both the Coulomb exchange and nonperturbative effects are important (so that (B.11), for example, is not a good approximation) and have positive binding energies. The behavior at large quark mass in V-QCD is somewhat dependent on the choice of the potentials. As pointed out above in (A.9), the scales m q and Λ IR become separated as m q → ∞. For r 1/m q one obtains the usual UV asymptotic solution, and for r 1/Λ IR the background obeys the IR asymptotics, leading to the confining Schrödinger potential. In the middle, however, the behavior of the Schrödinger potential is nontrivial. To compute it, we first need to solve the tachyon from its EoM.

B.2.1 The tachyon and the mass gap
The tachyon background EoM may be written as [14,23] When 1/m q r 1/Λ IR , the tachyon has already grown large, and decoupled the quarks from the glue. The dilaton is still small, and its evolution is governed by the YM RG flow, as the quarks are decoupled. More precisely, in this regime (B.13) where 0 = (x = 0) is the UV AdS radius for YM, and b YM 0 is the leading coefficient of the YM β-function. The appearance of the IR scale Λ IR in these expressions, which are the UV expansions for YM theory, may be surprising. In fact, it would perhaps be more appropriate to denote the scale of the expansions by a new quantity Λ YM UV (which is not the same as Λ UV , the scale of the UV expansions for r 1/m q ). But YM has only one energy scale, which is therefore the only scale in the model smaller than m q thanks to the decoupling of the quarks. That is, Λ YM UV ∼ Λ IR , and for simplicity we have already used this in (B.13).
Let us then insert the expansions in (B.12) in order to solve for the tachyon. We keep the leading behavior of all coefficients in 1/ log(rΛ IR ). In particular, since λ /A = O (log(rΛ IR ) −2 ), the terms involving λ will be dropped. As we shall see below, this works at least for r close to the lower end of the scaling region. We also expect that the tachyon grows large, and therefore 1 e −2A κ(τ ) 2 . The last term on the first line of (B.12) and the term on the second line dominate. Consequently the tachyon satisfies where we used the expansions (B.13), and κ 0 = κ(λ = 0). Let us study the asymptotics where we require v p > 1 in order to ensure that the tachyon grows fast enough 26 . Then the "asymptotic" solution of (B.14) reads 27 For v p > 2 there are no regular solutions. Let us then compute the Schrödinger potentials. First, the Schrödinger coordinate is given by The behavior of Ξ is dominated by its dependence on the tachyon in all sectors (flavor nonsinglet vectors, axials, scalars, and pseudoscalars): Since log Ξ increases fast enough with τ , the first term in (B.3) can be approximated by Combining these we obtain This result can be expressed in terms of the Schrödinger coordinate by using (B.16), (B.17) and (B.19), but the above form turns out to be more useful. The second term of (B.3) is negligible for the vectors and scalars, but important for the axials (and pseudoscalars), for which we find where w 0 = w(λ = 0). 26 The analysis can be extended to 0 < v p ≤ 1 where some assumptions made above fail and a more careful analysis is needed. The tachyon still grows only logarithmically and the conclusions are similar to the case 1 < v p < 2. Let us first assume that the results (B.16) and (B.17) hold in the whole regime 1/m q r 1/Λ IR . We will later discuss when this is not the case. For 1 < v p < 2 we then obtain that up to logarithmic corrections, V S ∼ 1/r 2 . The potential decreases with r, and reaches its bottom at r ∼ 1/Λ IR , where V S ∼ Λ 2 IR . Therefore the meson mass gaps are characterized by Λ IR , and only logarithmically enhanced with increasing m q .
For v p = 2, the estimates (B.22) and (B.23) match up to the multiplicative coefficient. The Schrödinger potentials behave as If ξ < 1, the result decreases with r. The bottom of the potential is reached at r ∼ 1/Λ IR , which leads to the mass gap If ξ > 1 the result increases with r. The bottom of the potential is therefore at r ∼ 1/m q , and m gap ∼ m q . (B.26) In the marginal case ξ = 1 the Schrödinger potential is flat and O m 2 q in the whole regime. The mass gap is therefore also m gap ∼ m q . (B.27) In conclusion, the mas gap has the power law expected for nonrelativistic states when ξ ≥ 1. However, as we shall see below, for ξ > 1 the mass splitting is O (m q ), i.e., larger than that of nonrelativistic bound states. Therefore only the marginal case ξ = 1 may potentially reproduce both realistic mass gap and splitting.
Let us then discuss logarithmic corrections to the critical case v p = 2. That is, we assume as τ → ∞. In this case, the tachyon solution is given by 30) and no regular solution is found when v > 1. The leading term of the Schrödinger potential is now By similar analysis as above, the bottom of the Schrödinger potentials lies at r ∼ 1/Λ IR and the mass gap is O (Λ IR ) when v < 0, up to corrections which grow slower than any power of m q as m q → ∞. When 0 < v ≤ 1, the bottom of the the potential is at r ∼ 1/m q , and the mass gap is O (m q ).
Recall also that, as we pointed out above, in some cases the above analysis is not valid in the whole regime 1/m q r 1/Λ IR . This can happen if the terms involving λ (r) in (B.12) start to dominate as r approaches the value 1/Λ IR . The growth of these terms requires that the factors in the square brackets depend on the tachyon, because λ is known to obey the RG flow of YM in this regime, and remains smallactually r ∼ 1/Λ IR is exactly the region, where λ finally reaches values O (1). Such a dependence on the tachyon may arise from the logarithmic derivative ∂ ∂λ log V f . The most natural Ansatz which leads to this is an exponential where the crucial point is that the factor a(λ) depends on λ. For the sake of generality, we however assume that ∂ ∂λ log V f ∼ −a 1 τ dp (B.33) as τ → ∞, where d p may differ from v p . From (B.12) we see that the terms which were neglected above become important when If d p > 0 and for any growing tachyon solution this condition is saturated within the range 1/m q r 1/Λ IR , because at r ∼ 1/Λ IR the left hand side is O (1) and the tachyon is already 1 (and for r ∼ 1/m q , the left hand side is sizeable whereas the tachyon is, by definition of m q , O (1)).
For the sake of concreteness, let us consider the case v p = 2 so that the tachyon has the r-dependence of (B.16) when the additional terms ∝ λ are still small. Then (B.35) is saturated for That is, the tachyon behaves as r ξ only up to r ∼ r c , which is only enhanced with respect to r ∼ 1/m q by a logarithmic term. For r c r 1/Λ IR (i.e., for almost the whole interval 1/m q r 1/Λ IR ) we need to solve the tachyon EoM with different dominant terms. Assuming that the right hand side in (B.34) dominates over the left hand side, we obtain which becomes for the current case This equation is solved by so that the tachyon only increases logarithmically for r c r 1/Λ IR . Inserting this expression in (B.22), and recalling that r c is close to 1/m q , we see that the bottom of the Schrödinger potential is found at r ∼ 1/Λ IR , and the mass gap is again O (Λ IR ) up to logarithmic corrections.
In conclusion, the phenomenologically interesting large mass gap (with power law dependence on m q ) is obtained only when v p = 2, 0 ≤ v ≤ 1, and when ∂ ∂λ log V f is suppressed (e.g., d p ≤ 0). The last requirement means that for the critical exponential asymptotics log V f ∼ −aτ 2 , the factor a cannot 28 depend on λ.
Interestingly we notice that potentials I, the construction of which was motivated by using completely independent arguments in [14,23], do produce a large mass gap. In this case we have where a 0 is indeed independent of λ. The value of a 0 was fixed in the UV (where the tachyon is small) to reproduce the UV dimension of the chiral condensate and quark mass so that The value of ξ is therefore slightly below the critical value ξ = 1 (because also 2 0 / 2 < 1). This is not surprising since holographic models are not expected to work perfectly in the UV region, where the coupling λ is small, so fixing the parameters by UV arguments may lead to some tension in the model. Recall, for example, that the UV asymptotics of the pressure and the UV asymptotics of the correlators for the energy momentum tensor give slightly different numbers for the normalization of the glue action [77]. For potentials II, we used instead a(λ) with nontrivial dependence on λ, and consequently the meson mass gap will be small.
Finally, let us point out that in all cases classified above, H A of (B.5) contributes at leading order to the Schrödinger potentials near their bottoms. This is seen from the estimates (B.22) and (B.23) when the bottom is at r ∼ 1/Λ IR , and by inserting the UV expansions in (B.5) if the bottom is at r ∼ 1/m q . Consequently, axial vectors and pseudoscalar mesons have larger mass gaps than vectors and scalars: the ratio of, say, the axial and the vector mass gap approaches a finite number which is larger than one in the limit of large quark mass 29 .

B.2.2 Mass splittings
Let us then discuss the mass splittings between the lowest meson states for the various cases described above. When the Schrödinger potential V S (u) has a clear and regular minimum, the splittings can be estimated by computing the second derivative V S (u 0 ) at the minimum u = u 0 . Most of the potentials discussed above fall into two categories, where the Schrödinger potential has a clear minimum either at r ∼ 1/m q or at r ∼ 1/Λ IR .
When the minimum is at r ∼ 1/m q , we found that the mass gap was O (m q ). From the definition (B.2) of the Schrödinger coordinate we see that u ∼ r and the scale of the derivatives of all relevant functions is given by d/dr ∼ m q . Therefore m q is the only scale that enters the analysis, and the splitting is also O (m q ), i.e., much larger than expected for nonrelativistic bound states. Notice that this also includes some of the asymptotics with critical v p = 2, i.e., those with 0 < v ≤ 1, and those with v = 0 and ξ > 1.
When the minimum is at r ∼ 1/Λ IR , the computation is more involved. We shall only discuss the case v p = 2 and v = 0, for which τ ∼ (rm q ) ξ and the mass gap is large, Recall that the minimum is at r ∼ 1/Λ IR when 0 < ξ < 1. At the minimum, all fields (τ , λ, and A) have logarithmic or power-law dependence on the coordinate r so the scale of the r-derivatives is Λ IR . The relevant quantities are, however, the u-derivatives of the potential, which can be estimated by using the chain rule and du/dr ∼ τ ∼ (m q /Λ IR ) ξ . We find that where n > 2. The extent of the Schrödinger wave functions around the minimum is determined for the lowest fluctuation modes by V S (u 0 ), and roughly given by ∆u = u − u 0 ∼ 1/Λ IR . Outside this region, the wave functions vanish very fast. The higher order derivatives in (B.43) vanish as m q → ∞, and therefore higher order terms in the Taylor expansion of the Schrödinger potential are suppressed (for u − u 0 ∼ 1/Λ IR ) and the potential takes the Harmonic oscillator form. In particular, the masses are given by where n = 0, 1, 2, . . .. From here we see that the mass splitting is suppressed as In the remaining case (v p = 2, v = 0, and ξ = 1) the bottom of the Schrödinger potential, as obtained from the leading tachyon solution, is flat for 1/m q r 1/Λ IR , suggesting a relatively small level splitting. There are, however, subleading logarithmic corrections to the solution due to the YM RG flow, which we have not computed. They will be important, and are expected to cause a minimum either at r ∼ 1/m q or at r ∼ 1/Λ IR . Then the mass splittings of the very lowest states will be as in one of the cases discussed above (with ξ = 1 if the minimum is at r ∼ 1/Λ IR ), up to logarithmic corrections, unless the subleading corrections cancel miraculously.

B.2.3 Decay constants
Finally let us discuss the mass dependence of the (vector/axial) decay constants, given in (B.8). The Schrödinger wave function vanishes very rapidly whenever m 2 n < V S (u), which limits the integral to the "classically allowed" region.
First we notice that when the Schrödinger has its bottom at r ∼ 1/Λ IR , the decay constants of the lowest states are very small, because Ξ in the above integral formula contains the factor V f which is exponentially suppressed because the tachyon is also large near r ∼ 1/Λ IR . For the interesting case of v p = 2, v = 0, and 0 < ξ < 1 we find that One can also show that the pion decay constant has similar dependence on m q . Recall that H V = 0 for the vectors and H A is given in (B.5). The low-lying states are therefore asymptotically decoupled. Sizeable decay constants are only found for highly excited states. More precisely, the suppression factor disappears when the classically allowed region extends to r ∼ 1/m q where the tachyon is no longer small. Since V S ∼ 1/u 2 in the UV, this requires m n ∼ m q (whereas the lowest states had m n ∼ Λ IR (m q /Λ IR ) ξ ). Therefore the coupled states appear at the scale where mesons are expected in QCD. By using (B.8), it is possible to show that for such states the decay constants are When the Schrödinger potential has its bottom at r ∼ 1/m q , it is easy to see from the above expressions that for the lowest states, because m q is the only scale which enters the formulas.
C. Analysis of the free energy C.1 Chiral condensate as the derivative of free energy The chiral condensate (including the sum over flavor) can be defined as 30 where E is the (zero temperature) energy density of QCD, and the subscript E denotes that Euclidean signature was used (so that the sign of the action is opposite with respect to the Minkowski signature, which was used in the main text). By computing the renormalized on-shell action using the identity (C.1) one can find the relation between the chiral condensate and the coefficient σ in the vev term of the tachyon. This is, however, slightly complicated in V-QCD. First, since there is full backreaction: changing the value of the quark mass will affect the geometry, possibly leading to nontrivial contributions in the m q derivative of (C.1). Second, the counterterms needed to regularize the on-shell action depend on the quark mass, and may also contribute in the derivative. As it turns out, these issues can be fully solved in the limit of zero quark mass.
Below we will first demonstrate how the free energy, and the backreaction in particular, can be analyzed by using the EoMs for the background. This will be compared to the direct computation of the regularized on-shell action done in Appendix C.3. In the calculations below the fields shall be decomposed as asymptotically in the UV (where r → 0). Here the terms with subscript zero (one) are the source (vev) terms. The terms A τ and λ τ and are sourced by the leading quadratic terms in the tachyon in the UV (∝ τ 2 0 ) and are therefore ∝ m 2 q . The term τ q is the leading nonlinear term of the tachyon. It is ∝ m 3 q and computed in Appendix C.2. It is shown to be subleading with respect to the vev term τ 1 and therefore irrelevant in the calculations below. It is also argued that the same is true for similar terms due to quartic and higher order tachyon perturbation in the expansions of A and λ. That is, the terms sourced by τ 0 τ q or τ 4 0 are subleading with respect to the vev terms A 1 and λ 1 . The UV expansions of the various terms will be given in Appendix C.3.

C.1.1 Chiral condensate at vanishing quark mass
Let us start by computing the chiral condensate at zero quark mass. This is the simplest case because the square of the source term of the tachyon does not contribute in the on-shell action, and as we shall see, there is no issue with extracting value of σ (the proportionality coefficient of the vev term) from the UV asymptotics.
Consider a perturbation of a generic background (around m q = 0) which keeps Λ UV (but not the quark mass) fixed. The variation of the background solution in the UV is found by using the UV expansions from Appendix C.3: The chiral condensate may then be computed by studying the variation of the on-shell action as seen from (C.1). The computation is analogous to the holographic derivation of the first law of thermodynamics at finite temperature and chemical potential: the differential of the free energy equals the variation of the action, and expressing this in terms of the various UV and IR boundary terms gives the terms in the desired expression. In order to be as precise as possible, we will formulate the computation in terms of a conserved infinitesimal current, which can also be easily used in the computations later. It can be found as follows.
The on-shell Lagrangian can be expressed as a total derivative [14]: But any leading variation of a generic Lagrangian around its on-shell value is a total derivative as well, given formally by 31 where the sum is over all fields in the Lagrangian. Requiring equality of the generic expression and the variation of (C.7), one identifies the following infinitesimal conserved "current": Indeed it is straightforward to check that J = 0 by using the EoMs (and their variations).
We will then require that the variation of the background is regular in the IR. By using the IR expansions of the background [14] in the expression (C.9) one sees that J vanishes in the IR (r → ∞) so it must vanish everywhere for regular variations: (C.10) In the UV the above expansions are inserted, leading to The expression for the free energy density in terms of G, σ, and m q can be extracted from the finite temperature computation in Appendix C.3 (see the expression for the free energy (C.63)). At small quark mass and at T = 0 (so that also C = 0) one obtains the expression 32 Inserting here (C.11), Therefore the chiral condensate is given by The coefficient in this equation can be fixed by using the asymptotics of the scalarscalar correlator (see Appendix C in [23]), which leads to In the case of the gravitational action there is a complication because R contains second derivatives. It is, however, well known that the second derivatives can be isolated in another total derivative term, related to the Gibbons-Hawking boundary term. 32 We omitted contributions from the reference background used to regulate this expression, but as we shall argue below, this does not affect the result at m q = 0.

C.1.2 Chiral condensate at finite quark mass
Let us then try to generalize the above computation to finite m q . Recall that there is an issue in the definition of the vevs G and σ. In principle, they are well defined as the coefficients of the vev terms A 1 and τ 1 , respectively. However such definitions are useless in practice, because the source terms A 0 and τ 0 cannot be solved perturbatively to high enough orders in order to separate them from the vev terms in the UV. In practice only differences of the vevs can be computed. Therefore we will define the values of G and σ with respect to some reference solutions. Let us first discuss how we define a reference solution for all values of m q .
First a (IR regular) solution with exactly zero tachyon (and therefore zero quark mass) is picked. This solution is chosen to have, by definition, G = 0 (and trivially σ = 0). Therefore A τ and A 1 in (C.2) are zero, and the solution for the scale factor A defines the source term A 0 . Let us call this solution 1. Then solution 2 is chosen which has finite m q , defines the nonnormalizable term τ 0 of the tachyon in the UV, and has σ = 0 by definition. Further one requires that G = 0 also for solution 2. Since A 0 was already defined, this choice also defines A τ .
For 0 < x < x c , the solution 2 can be identified with the solution that was used to define the vanishing of σ in the plots of Fig. 9. In principle any solution can be picked, but as we argued in Sec. 5, a choice which avoids fine tuning is to pick the solution 2 is the crossover between regimes A and B when 0 < x < x c and a solution with tiny quark mass in the conformal window.
We have argued that "nonlinear" terms of the tachyon are suppressed with respect to the vev terms. Similarly, the terms sourced by quartic tachyons in the UV expansions for A and λ are suppressed with respect to the vev terms. Therefore fixing the values of G and σ for solutions 1 and 2 is enough to define the vevs for all backgrounds. It is also possible construct a reference background which has vanishing G and σ for arbitrary m q , as a (generally not IR regular) combination of the two IR regular solutions 1 and 2 constructed above by appropriately scaling the constructed A τ and τ 0 to have the desired quark mass.
The above construction implies that the vevs can be written as where the hatted quantities are the exact coefficients of the vev terms, and the quantities without hat are given by the above subtraction procedure. The coefficients k i and G 1 are independent of m q .
Let us then proceed with the calculation. Consider again a perturbation of a background with δm q = 0, but now at generic value of m q . The UV expansion of the perturbation has additional terms, in particular where we chose to use the hatted vevs. The variation of b 0 λ includes a term ∝ m q dm q which turns out to only contribute at subleading orders in the UV, and the term (compare to (C.4)) It is again useful to study the current (C.9). As above, J = 0 because it vanishes in the IR limit. In the UV one obtains In particular, we expect that no terms ∝ m q δm q arise in the UV becauseσ andĜ are exactly the coefficients of the vev terms. This identity may be rewritten as where the expression in the square brackets has similar structure to the regularized vacuum energy given above in (C.12). The natural expectation is that indeed but in principle the expression could also contain additional terms ∝ m 2 q (or constants) which would cancel in the regularization.
We see that where extra terms ∝ m 2 q in the definition of E would effectively change the value of k τ which remains unknown in any case 33 . Therefore the chiral condensate is given by

C.1.3 Free energy on the Efimov spiral
Let us then discuss how the free energy is computed for the configurations of Fig. 9, which are found when 0 < x < x c . Let us assume for simplicity that we are able to define the vevs such that the constant k τ of (C.23) vanishes. In this case the chiral condensate is simply given by the derivative of the free energy density with respect to the quark mass: Inserting the asymptotic result of (5.14) in (C.24) and integrating, one can readily find the free energy for asymptotically small quark mass: where E 0 is the free energy of the solution having m q = 0 = σ, and u = log Λ UV /Λ IR . In general, free energy differences are given by the area between the spiral and the horizontal axis (see, for example, Fig. 12). Interestingly, the result for both zero m q and zero σ simplifies to where we used the handedness of the spiral in (5.16): sin(φ m − φ σ ) > 0. Therefore these solutions are dominant over the solution with zero tachyon. By using the UV boundary conditions, this identity can be rearranged to read n − 1 n g n c n σ σ n = −c σ σδm q − nmax n=2 1 n c n σ σ n δg n .
(C.31) Therefore the quantity in the square brackets is identified as the new free energy (over N f N c ). It matches with the free energy obtained by renormalizing the action (C.12) up to the last term involving the couplings g n . This extra term is proportional to where W is as in (6.2). Therefore the term agrees with that found in [61]. The condensates are given by where the normalization factors were read from (6.1). Therefore we find agreement with the large N factorization of the expectation values. The free energy on the Efimov spiral asymptotically close to the origin can be computed as above. That is, by inserting the spiral equations (6.8) to (C.33) and integrating, we obtain Notice that the term involving g n in (C.35) vanishes for n = 2, which is consistent with nonzero g 2 amounting to a redefinition of the parameters of the Efimov spiral. The higher order terms are suppressed as u grows because σ ∼ exp(−2u). Therefore we conclude as above that the solutions with m q = 0 but σ = 0 dominate over the solution with σ = 0, and that the solution with the smallest u (for u 1 so that the parametrization of the spiral is accurate) has the lowest free energy.

C.2 Contributions ∝ m 3 q in the on-shell action
The on-shell action may involve terms ∝ m 3 q which arise due to the nonlinear nature of the tachyon EoM. We will now check how these terms behave. We will need the UV asymptotics of the source and the vev terms of the tachyon which are given in equations (C.48) and in (C.54) in Appendix C.3. Using results from [14], the tachyon EoM can be written in the UV as where κ 0 = κ(0). All contributions suppressed by 1/ log(rΛ) 2 or more are hidden in the coefficients of the terms linear in the tachyon, as well as all contributions suppressed by 1/ log(rΛ) in the nonlinear terms (including the complete terms ∝ τ 2 τ and ∝ τ 2 (τ ) 3 ). The "nonlinear" tachyon solution is then found by writing where the UV expansion of τ 0 is given in (C.48) and τ q is the term which needs to be solved. Using the fact the τ 0 solves the linear equation, one finds that One finds that (dropping the terms corresponding to the source and vev terms, which arise as solutions to the homogeneous equation) , (C. 39) which agrees with the result of [46] when ρ = 0. Notice that the nonlinear term is subleading to the vev term in (C.54) provided that ρ > 1 4 , (C. 40) which is satisfied for QCD in the Veneziano limit as where b 0 and γ 0 are the leading coefficients of the beta function and the anomalous dimension of the quark mass, respectively. Therefore one can conclude that only linear terms of the tachyon EoM are relevant in the computation of the on-shell action 34 . The conclusion is different from that of [46] where the quark mass did not run. Notice, however, that these terms are only logarithmically suppressed and would be important for large values of m q if a finite value of the UV cutoff was used. 34 Notice that even though the nonlinearities in the tachyon asymptotics do not appear directly in the computation of the on-shell action, nonlinear terms of the EoM are important in general because they affect the values of the vevs such as the chiral condensate through IR boundary conditions.

C.3 Regularized on-shell action
The free energy is given by the on-shell value of the action. It is a UV divergent quantity and needs to be renormalized 35 . We work here at finite temperature, and subtract the UV divergences by considering the difference with respect to a reference (zero temperature) background with the same values for the sources. We could also consider the difference between two different zero temperature backgrounds.
It is not difficult to show that the Lagrangian is a total derivative even in the presence of the tachyon (as already shown at zero temperature in [14]). After a straightforward computation, and taking into account the Gibbons-Hawking term, one finds that When the quark mass is finite, the tachyon contributes at O (r 2 ) (terms ∝ m 2 q ) and at O (r 4 ) (terms propto m q σ) in the UV expansions. These contributions need to be taken into account in the holographic renormalization procedure. The calculation of [73] at N f = 0 used a reference (thermal gas) solution to subtract the divergences. We will use the same technique here. Denoting the fields and other variables of the reference solution by tildes, the following conditions need to be fulfilled at the UV cutoff:β where the possibility that the cutoffs of the two solutions are different was included, = , as in [73]. In the absence of the tachyon this was convenient since the UV scales Λ = Λ UV of the two solutions could be chosen to be the same,Λ = Λ. The fact that the present system has two scalars suggest that it is better to choose herẽ = and satisfy the last two conditions by varying the sources Λ and m q . However, it turns out to be convenient to still require thatΛ = Λ and vary the cutoff (and m q ) instead. This is so because the tachyon contributions are independent on whether one chooses to vary Λ or : The variation if suppressed by O ( 4 ), i.e.,˜ / = 1 + O ( 4 ), so that the effect on the tachyon contributions will be down by O ( 6 ) and therefore negligible. Keeping Λ fixed one can maintain very close contact to the computation of [73].
Turning on the quark mass will modify the UV expansions of A and λ. at this point it is useful to write down the complete UV expansion including all relevant terms. One can decompose, as r → 0, 35 For detailed analysis of the holographic renormalization of dilaton gravity, see [60,78] The source terms have the expansions [14] A 0 = − log r + log + 4 9 log rΛ and ρ = γ 0 /b 0 with γ 0 being the leading coefficient of the anomalous dimension of the quark mass. The O (r 2 ) terms A τ and λ τ were also included here which are sourced by the O (r 2 ) tachyon perturbation 36 . As it turns out, λ τ is suppressed by logarithms of log rΛ with respect to A τ so that it will not enter the calculation and its expansion will not be needed. A τ is given by where W 0 = V f 0 (0), κ 0 = κ(0), and as above The vev terms have the expansions (compare to [73]) Let us then go on with the renormalization procedure. In the expressions below only the differences of the vevs G and σ with respect to the values of the reference solutions will appear. Therefore, without loss of generality, one can set the vevs of the reference solution to zero. As motivated above, one can chooseΛ = Λ. Then the relation between˜ and is fixed 37 by the first condition in (C.44): There are also O r 4 terms which arise due to the perturbation from the quartic terms in the tachyon. As we argued above, these contributions are subleading with respect to those arising from the vev terms. 37 It can be checked that the variation of m q , which would enter through the term λ τ is subleading.
(C. 56) It is then straightforward to calculate the renormalized pressure. Inserting the expression of the (unrenormalized) on-shell action from Eq. (C.42) one finds that −βV 3 p = lim

D. Gell-Mann-Oakes-Renner relation
In this Appendix we check the GOR relation explicitly. The starting point is the normalization integral (E.13) for the pion wave function. When m q is small, it is dominated at small u r. Also the wave functionψ P is constant in the UV up to corrections O (r 2 ) (see Appendix G in [23]). Therefore integral one needs to compute is (the UV contribution to) (D.1) 38 Recall that the vevs here should be interpreted as their differences with respect to the values of the reference solutions, because we set the vevs of the reference solution to zero.
Let us first try to compute the integral by using the UV expansions of the various fields, given in Appendix C.3. We obtain Choosing the cutoff such that r cut m q /σ (but so that it is not too large so that our approximation are still valid) we obtain (D.11) Our approximations will break down at r cut ∼ 1/Λ UV , where the r-dependence of the tachyon changes qualitatively. In order to complete the calculation, one should estimate the contributions to the normalization integral (E.13) for r 1/Λ UV . But in this regime the dependence on the tachyon is regular and one can just do a Taylor expansion at m q = 0. The cutoff dependence should cancel against the UV contribution given in (D.11). As the dust clears, we can effectively set r cut ∼ 1/Λ UV , in (D.11), obtaining The normalization condition (E.13) becomes for small quark mass with corrections suppressed by m q Λ 2 UV /σ. For the pion, n = π, this is the GOR relation. It can be checked that the proportionality factor is correct for our definitions of f π and qq . For higher pseudoscalar mesons this relation implies that f 2 n vanish as m q → 0, which is consistent with the vanishing of Π L in this limit.

E. Fluctuation equations, f π , and the S-parameter
The radial wave functions for the flavor nonsinglet fluctuations satisfy the following equations [23]: In addition it is convenient to define the pseudoscalar wave function which satisfies the single equation (E.9) At nonzero quark mass the pion decay constant may be defined in terms of the pole of Π L in (7.16) at q 2 = −m 2 π . As the computation is similar also for higher modes, a generic pseudoscalar fluctuation can be considered. To compute the decay constants one needs to study the fluctuations for small values of ψ P = k L ψ P,n + δq 2ψ P + O (δq 2 ) 2 , (E.12) where k L is a normalization constant which will be fixed below. When δq 2 = 0 the wave functions are normalizable. One can therefore choose the fields ψ I,n to satisfy the usual normalization condition withψ P,n defined as in (E.7) for the normalizable mode. At finite but small δq 2 we impose the standard normalization condition in the UV: 1 = ψ L (0) k L δq 2ψ L (0) ; 0 = ψ P (0) ∝ψ P (0) , (E.14) where we used the fact that the normalizable wave functions vanish in the UV. Expanding the fluctuation equations at small δq 2 gives H A ∂ uψP = m 2 n ∂ uψL − ∂ u ψ L,n . (E. 16) By using these equations and the fluctuation equations for ψ I,n one finds the identity ∂ u m 2 n C V ψ L ∂ u ψ L,n − ψ L,n ∂ uψL + H A C V ψ P ∂ u ψ P,n − ψ P,n ∂ uψP = H A C V ψ P,n (ψ P,n − ψ L,n ) (E.17) Integrating this over u, using the boundary conditions from above, and further using the fluctuation equations to simplify the result one obtains Inserting (E.14), integrating partially, and using again fluctuation equations this relation simplifies to By using this result and the definitions from (E.11), Π L from (7.16) can be written as for small δq 2 . The decay constants are given in terms of the residue at δq 2 = 0: From (E.20) we see that for q 2 ∼ m 2 π , as m q → 0, Π L (q 2 ) f 2 π m 2 π q 2 (q 2 + m 2 π ) = f 2 π q 2 − f 2 π q 2 + m 2 π (E.22) so that S 0 (defined as the zero momentum value of q 2 Π L (q 2 )) approaches f 2 π . The same factor for the axial-axial correlator does not have a pion node, so that simply 23) in the same scaling limit. Consequently, the axial-axial correlator has the correct structure at small momentum: q 2 η µν − q µ q ν Π A (q 2 ) + q µ q ν Π L (q 2 ) f 2 π η µν − q µ q ν q 2 + m 2 π . (E.24) That is, the longitudinal term arises from the coupling of the pion to the axial current.
Let us then derive a few integral representations which are useful when computing observables numerically. First we analyze (E.25) where we inserted (7.16), and used the fact that r u near the boundary. Recalling the normalization ψ V (u = ) = 1 = ψ V (u = ) we obtain where fluctuation equations (E.1) were used at the last step. Taking → 0, it follows that D(q 2 ) = 1 4 A rather similar formula can be derived for the S-parameter. Let us denote Then the S-parameter (7.18) can be written as Derivating the fluctuation equations with respect to q 2 results in where it is understood that ψ V,A are evaluated at q 2 = 0. The combination appearing in the S-parameter can then be written as where we dropped higher order terms in in the first step (noticing thatψ V,A vanish fast in the UV) and used the fluctuation equations (E.1) as well as equations (E. 30) in the second step. Taking → 0, we obtain the final result Here one could also insert that ψ V = 1 when q 2 = 0.

F. Free field (one-loop) computation of the S-parameter
The free field result for the vector-vector and axial-axial correlators is given by Tr [γ µ (γ 5 )(/ k + m)γ ν (γ 5 )(/ k − / q + m)] where the γ 5 's are only present in the axial-axial correlator. The loop is divergent but the contribution to the S-parameter will be finite. Doing the integral with dimensional regularization one obtains where the value at q 2 = 0 was subtracted in order to remove the logarithmic divergence.
For the S-parameter one obtains the well-known result whereas the corresponding finite difference reads G. Details on numerics The numerical results in this paper were computed for sets of potentials termed "potentials I" and "potentials II". They are exactly the sets given in [23], but we repeat their definitions here for completeness: • Both Potentials I & II. Here the most of the coefficients are fixed by matching to perturbative QCD. Exceptions include the normalizations factors V 0 , which fixes the UV AdS radius, and W 0 , which remains as a free parameter. We also set (x = 0) = 1, and choose the parameter λ 0 , which only affects the higher order coefficients of the UV expansions, such that the higher order coefficients have approximately the same relative size as with standard scheme choices in perturbative QCD. Explicitly, the coefficients satisfy Two qualitatively different choices for W 0 are possible: either constant W 0 , which satisfies 0 < W 0 < 24/11 , (G. 5) or W 0 fixed such that the pressure agrees with the Stefan-Boltzmann (SB) result at high temperatures [18] (without the need to include x dependence in the normalization of the action). The latter option is given explicitly (when (x = 0) = 1) by (Stefan-Boltzmann) , (G.6) so that the AdS radius is (x) = 3 1 + 7 4 x . (G.7) In this article we have always chosen W 0 = 3/11 for potentials I, and the SB normalized W 0 for potentials II. For the coupling of the gauge fields w(λ) which is required for the computation of the vector correlators, we used w(λ) = κ(λ) for potentials I and w(λ) = 1 for potentials II. The numerical result were mostly computed as detailed in [14,18,23]: coupled ordinary differential equations were solved by shooting from the IR, either starting from near the IR singularity (at zero temperature) or near the horizon (at finite temperature). In order to obtain accurate and reliable results, some tricks had to be used in various cases, in particular near the BZ point x = x BZ and at large quark masses.
In general, the difficulties of the numerics in V-QCD arise from two sources: the IR singularity and the logarithmic corrections to the UV asymptotics. The latter easily lead to sizeable numerical errors when one tries to extract the values of the sources or vevs near the boundary. To improve the accuracy of the numerical analysis, we did the following: • The scale factor A was used as the coordinate instead of r. This makes it easier to analyze the UV asymptotics, because the mapping from r to A is logarithmic and expands the details near the boundary.
• Near the IR singularity (at zero temperature), where possible, EoMs were implemented such that they do not contain the extremely small factors ∝ exp(−aτ 2 ). In many cases, this could be done by writing all expression in terms of V log (λ, τ ) = log V f (λ, τ ) = −aτ 2 + log V f 0 (λ) rather than in terms of V f .
• Near the BZ region and at very large m q , the zero temperature backgrounds were constructed by shooting from the IR toward the UV in four steps. Very close to the singularity, where the tachyon is large and completely decouples the flavors, YM solution was used for the geometry and the dilaton, whereas the tachyon was solved by using a simplified EoM with only the dominant terms on top of the YM background. In the next step, the complete tachyon EoM together with the decoupled YM flow for the geometry was used. In the third step, all fields were solved from the full coupled EoMs. In the final step, the tachyon was again solved separately since it is decoupled from the other fields near the boundary. Actuallyτ = e A τ was used as the field, because it decreases much slower than τ , and the flow can then be tracked closer to the boundary.
The first two steps were necessary because the tachyon EoM becomes stiff in the IR, in particular for large m q , so that the tachyon takes larger values in the IR than otherwise. The last step was useful in particular at large x, where the RG flow of the fields is slower, and one needs to solve them closer to the boundary in order to control the leading logarithmic behavior of the various terms.
• For the finite temperature backgrounds in the BZ region and at large m q , the background was solved in two steps, which were analogous to the two last steps of the zero temperature construction. That is, the full system was used in the IR, and the tachyon was treated as a decoupled field near the UV.
• The S-parameter was computed by using the integral formula (7.21) rather than UV asymptotics of the fluctuation wave functions.
• The nonsinglet meson wave functions were computed by rewriting the fluctuation equations into a system of two coupled first order equations (rather than a second order equation). Again the equations were written in a form which does not contain explicitly the tiny factors exp(−aτ 2 ).