Holographic QCD in the Veneziano limit and neutron stars

We use the holographic V-QCD models to analyse the physics of dense QCD and neutron stars. Accommodating lattice results for thermodynamics of QCD enables us to make generic predictions for the Equation of State (EoS) of the quark matter phase in the cold and dense regime. We demonstrate that the resulting pressure in V-QCD matches well with a family of neutron-star-matter EoSs that interpolate between state-of-the-art theoretical results for low and high density QCD. After implementing the astrophysical constraints, i.e., the largest known neutron star mass and the recent LIGO/Virgo results for the tidal deformability, we analyse the phase transition between the baryonic and quark matter phases. We find that the baryon density $n_B$ at the transition is at least 2.9 times the nuclear saturation density $n_s$. The transition is of strongly first order at low and intermediate densities, i.e., for $n_B/n_s \lesssim 7.5$.


Introduction
Last year witnessed a remarkable achievement that marked the beginning era of multimessenger astrophysics: the observation of a neutron star collision by the LIGO and Virgo collaborations. The aftermath of analysing the gravitational wave data and the associated electromagnetic kilonova observations still continues. A particularly thrilling fact is that already this first observation of a neutron star merger sets bounds on the properties of the underlying theory, i.e., quantum chromodynamics (QCD) at high densities and low temperatures. In the future more neutron star mergers are to be observed, which will then dictate how matter will behave under extreme conditions as quantified in the Equation of State (EoS).
Connecting these observations with the theory is, however, hard because QCD cannot be solved in the relevant regime. The cores of cold neutron stars are extremely dense so that they cannot be described in terms of weakly interacting baryons, but not dense enough for perturbative QCD (pQCD) to be applicable. Lattice calculation are not available in this regime either, and consequently it is not even known what the relevant microscopic degrees of freedom are at the center of a typical neutron star.
The gauge/gravity duality may provide us with a way forward. The applications of the gauge/gravity duality to quark-gluon plasma has been a very active field [1,2]. Very recently these techniques have also been applied in the context of neutron stars [3,4] and stiff phases of matter [5,6] which are likely to be needed to achieve large enough stars [7]. In [4], a holographic model dual to N = 4 super-Yang-Mills and quenched supersymmetric matter was studied, and it was observed that neutron stars can only accommodate very small amounts of quark matter.
In this article, we aim at progress in predicting the EoS of dense QCD matter by using more realistic holographic models, which break conformal symmetry, include confinement, and are able to model the dynamics of real QCD with good precision. The class of models that we will adopt is called V-QCD, which are effective, bottomup holographic models for QCD [8]. The letter "V" in the acronym stands for the Veneziano limit: we take N f , N c → ∞ while keeping x f = N f /N c fixed. This means that unlike most of the holographic models for QCD considered in the literature, it contains a fully dynamical quark sector on top of the gluons. The gluon sector is given by improved holographic QCD [9][10][11], i.e., five dimensional Einstein-dilaton gravity with a properly chosen dilaton potential, whereas the quark sector is given by tachyonic DBI (Dirac-Born-Infeld) actions [12][13][14][15]. Naturally, taking the backreaction of the quarks into account is particularly important at high densities.
The description of dynamical quarks comes with a price to pay: these models are also clearly more complicated than typical holographic models of QCD. Being effective models, the actions contain several a priori unknown potentials which need to be chosen by comparing the model predictions to QCD physics. With a correct choice of action the models have known qualitative features of QCD such as asymptotic freedom, linear confinement, and discrete spectrum of hadrons and glueballs [8][9][10]16,17]. The gluon sector has also been shown to agree well with lattice thermodynamics for pure Yang-Mills theory and with lattice results for the glueball masses [18].
In this article, we take the first steps of comparing the full, backreacted V-QCD models directly to lattice QCD data. Specifically we choose a set of models which fit well the three-flavor lattice EoS at µ = 0 and its first nontrivial cumulant around µ = 0. Remarkably, we show that a very good fit of the lattice data is possible by using relatively simple Ansätze for the potentials. We demonstrate that already after this fit, the EoS at low temperatures and high densities is realistic, and constrained enough to make nontrivial predictions.
We then use the holographic method to model the EoSs in the deconfined quark matter phase, and a set of interpolated EoSs (see, e.g., [19]), which are constrained by matching with effective theories for nuclear matter at low densities and with pQCD at high densities, to model the equation state in the baryonic phase. This construction results in a family of all possible EoSs which are consistent with the holographic model and our knowledge of QCD. Further demanding that neutron stars (obtained by solving the TOV equations) satisfy the known tidal deformability constraints and are able to support two-solar-mass stars enables us to make generic predictions. In particular, we consider the possibility of finding quark matter deep in the cores of neutron stars. We find the answer to be negative. Even though this means that the state of the neutron star is determined by the interpolated EoSs rather than holography, we can still make nontrivial predictions on the location and nature of the phase transition between the hadronic and quark matter phases. Most importantly, we can set relatively tight constraints on the latent heat at the transition, which is an important parameter in the neutron star mergers [20]. The basic picture which arises is the following: the nuclear (holographic) matter has a stiff (soft) equation of state, and the latent heat at the transition between the phases is sizeable.
The article is organized as follows. In the following section we will review the holographic model. In Section 3 we will determine the forms of the potentials in the V-QCD action by fitting them to available lattice QCD data. Section 4 discusses the thermodynamics of the holographic model at finite density and vanishing temperature. We also explain how we determine that the resulting pressures are physically viable by comparing them to a family of polytropic equations of state that are thermodynamically consistent and have the correct asymptotic behavior both at high and low density. In Section 5 we construct hybrid equations of state where the baryonic EoS at low density, given in terms of the polytropes, has a single transition to the quark matter phase, which is modeled by our holographic framework. The EoSs are then confronted with the astrophysical constraints, and predictions for the phase transition and neutron star masses are derived. Finally, Section 6 both contains the discussion and outlines some thoughts about future outgrowths of our work. Details of the fitting of the potentials are relegated in Appendix A.

Holographic model
In this work we will be using a class of holographic models for QCD (V-QCD) [8]. These models are effective bottom-up models containing both gluons and dynamical quarks which are fully backreacted to the glue in the Veneziano limit: The glue sector is described in terms of improved holographic QCD (IHQCD), i.e. five dimensional Einstein-dilaton gravity with properly tuned dilaton potential [9,10]. The quark sector is described by a tachyonic DBI + CS (Chern-Simons) action which is associated to a pair of space filling D4-D4 branes [14,15].
The dictionary is as follows: • The dilaton λ = e φ is dual to the TrF 2 operator and therefore sources the 't Hooft coupling g 2 N c of the Yang-Mills theory.
• The tachyon τ is dual to theqq operator and sources the quark mass. In this article, we will be considering chirally symmetric backgrounds so that the tachyon vanishes.
• The (vectorial) gauge field A µ is dual to the current J µ =qγ µ q. The temporal component therefore sources the quark chemical potential.
Here we already suppressed the flavor structure -our solutions will be flavor independent since all quark masses vanish.
The terms in the action of the model which are relevant for the current study are where are the actions for the glue and quark sectors, respectively. For the tachyon potential, we take the Ansatz V f (λ, τ ) = V f 0 (λ)e −τ 2 . We are left with four potential functions V g (λ), V f 0 (λ), w(λ), and κ(λ) which need to be determined by qualitative and quantitative comparison to QCD data. In this article we will concentrate on chirally symmetric backgrounds describing the deconfined phase of QCD so that τ = 0. Consequently, the results actually do not depend on κ(λ) directly. The Ansatz for the metric can be written as The thermodynamics in this model has been discussed at zero µ in [21] and at nonzero µ in [22], by using potentials that reproduce various features of QCD at qualitative level. It was shown that depending on the precise choice of potentials, there can be a single first order order phase transition, or separate phase chiral and confinement-deconfinement phase transitions. 1 In the latter case, the critical temperature of the (second order) chiral transition was found to be larger than that of the (first order) deconfinement transition in the regime of low µ. In other words, an intermediate deconfined but chirally broken phase was found. At high enough µ, only the deconfined chirally symmetric quark matter phase was present for all potentials. A similar phase structure has been found also in models based on a D3/D7 brane system [23,24].
In the next section, we will compare the predictions of the V-QCD model to lattice data at zero chemical potential. This comparison strongly favors potentials producing phase diagrams with no intermediate (deconfined but chirally broken) phase. In the absence of the intermediate phase, the model will only have nontrivial thermodynamics in the deconfined phase. 2 We will concentrate on analysing the results in this phase in the rest of the article.

Comparison to lattice data
We then proceed to the determination of the potentials. The main idea is as follows. For each of the potentials, there are various qualitative constraints which restrain the dependence of the potentials on λ in the UV (λ → 0) and/or in the IR (λ → ∞). We choose Ansätze which satisfy these qualitative constraints, and also include extra parameters which can then be fitted to lattice Yang-Mills and QCD data at intermediate values of the coupling. We proceed by first fitting the potential V g of the glue action to lattice data at large N c . Then, keeping V g fixed, we fit the remaining potentials to QCD data at N c = 3 (because there is not much data available at higher N c ). Details will be given below and in Appendix A. We will argue that using data at N c = 3 is a reasonable approach in Sec. 6.
In this work we will only consider lattice data for thermodynamics at µ = 0. Comparison to other data such as meson and glueball spectra will be considered in future work [26].
Earlier related works considered fitting QCD lattice data at µ = 0 using a somewhat simpler holographic model based on Einstein-Maxwell-dilaton actions and studying the continuation of the equations of state to the (µ, T )-plane [27][28][29][30]. These arti-1 Notice that while ordinary QCD does not have a well defined order parameter for the deconfinement transition, there is such a parameter at large N c : the pressure in the confined (deconfined) phase scales as N 0 c (N 2 c ). 2 The pressure in the confined phase is dominated by string loop contributions which have been modeled in a simple approximation in [25]. cles concentrate on the physics near the critical point, whereas in the current article we will be interested in the physics at low temperatures. The main difference in the current model to the earlier work is the presence of the first order confinementdeconfinement transition even at small chemical potentials. We will discuss how the transition is treated below.
There are also numerous other approaches which aim at modeling the thermodynamics of QCD in the deconfined phase, also taking into account constraints from lattice data. These include extensions of pQCD to relatively low temperatures and chemical potentials using quasiparticle models based on hard thermal loop computations [31,32] and analyses based on Nambu-Jona-Lasinio model [33]. We are planning to carry out a comparison of our results for the thermodynamics to other approaches in future work.

Glue sector
First, the dilaton potential V g is chosen after a careful comparison to Yang-Mills data [9,10,18,34,35] at x f = 0 = N f . We choose here a potential [25] which has been shown to work well with backreacted flavors. The Ansatz reads This Ansatz has been chosen to reproduce qualitative features of QCD such as confinement, magnetic charge screening, and linear glueball trajectories (determining IR asymptotics) as well as asymptotic freedom (determining UV asymptotics). Here V 1 and V 2 are fixed by requiring the UV RG flow of the 't Hooft coupling to be the same as in QCD up to two-loop order. The fit parameters in the potential are then V IR and λ 0 , but also M and the energy scale of the UV RG flow Λ UV need to be determined through the fit (see Appendix A for details).
The result of the fitting to lattice data extrapolated to N c → ∞ [36] is shown in Fig. 1. Notice that the overall normalization of the pressure determines the value of M , the critical temperature fixes Λ UV whereas the parameters V IR and λ 0 affect the shape of the curve. The fit is stiff: producing results which differ qualitatively from lattice data would require input functions which have highly nontrivial λ dependence.

Flavor sector
The flavor sector of the model can be fitted to QCD data (now including flavors) following a rather similar strategy as for the glue sector. However since there is no lattice data available at large numbers of N c and N f , we choose to compare the model to lattice data with 2 + 1 flavors. Figure 1: Comparison of the thermodynamics of the glue sector of V-QCD to Yang-Mills lattice data [36]. The solid lines and error bars represent the extrapolation of lattice data to N c = ∞ and the dashed curves are the results for the holographic model.
Our Ansatz for the flavor potentials reads The asymptotics of the potentials were chosen to reproduce known features of QCD such as linear radial meson trajectories, correct behavior at large quark mass, compatibility with flavor anomaly structure, and UV properties of the flavored correlators [8,17,37] (see also [38,39]). We will also consider potentials with W 0 = 0 and in this case we choose a slightly modified Ansatz for w(λ), see Appendix A.
The parameters W 1 , W 2 , κ 0 , and κ 1 were determined by the UV dimension of thē qq operator and the RG flow of the coupling and the quark mass in the UV. The remaining additional parameters are W 0 , W IR , w 0 , w 1 ,w 0 , w s ,κ 0 , andκ 1 .
These parameters are then fitted to lattice data at µ = 0. The fit is carried out in stages: first we consider the interaction measure and pressure of QCD, with lattice data with 2+1 dynamical flavors from [40]. This fit is independent of the shape of w(λ), and κ(λ) affects only the location of the phase transition of V-QCD from the confined to deconfined phase. The most relevant parameters are therefore W 0 and W IR . There is an important complication with respect to the fit to Yang-Mills data: QCD at physical quark masses has a crossover, whereas the holographic model has a first order phase transition. The analysis of [25] suggests that this is due to the neglect of stringy loop contributions which are important in the confined phase. We follow the notion of the reference, and fix the phase transition temperature to around 120 MeV rather than near the peak of the susceptibilities in the lattice data near 155 MeV. The location of the transition is then expected to shift to higher temperatures and become weaker due to the aforementioned loop contributions. This strategy produces much better overall fit than fixing the transition of the holographic model to lie at, say, 155 MeV.
We choose to work at x f = 1, roughly corresponding to 2 + 1 flavors at N c = 3. We choose samples with three values for the critical temperature (110, 120, and 130 MeV) and four values for the parameter W 0 (0, 1, 2.5, and 5.886 with the last value giving the "automatic" normalization to Stefan-Boltzmann thermodynamics at high temperatures [21]). The parameter W IR is fitted to the functional form of the interaction measure. The parametersκ 0 andκ 1 are chosen to set the location of the phase transition at the desired point. Also M and Λ UV are refitted reflecting their possible dependence on x f . We stress that the fit is again stiff and it would not be possible to find good fits to arbitrary data without introducing complicated structure in V f 0 . In particular, after tuning V IR such that the shape of the normalized interaction measure is roughly correct, the bump in the interaction measure is automatically reproduced by the dynamics of the model rather than put in by the fit Ansatz.
The comparison to the lattice data for basic thermodynamic functions is shown in Fig. 2 (top row). The various curves correspond to different values of the transition temperature and W 0 from the samples given above. The thick curves (for T > T c with T c = 155 MeV) represent the fit and we also show the result of the model as the dashed lines for the deconfined phase below T c : as explained above, in this regime the deconfined phase should be replaced by the confined phase with added loop contributions. We plan to return to this issue in future work.
The shape of w(λ) can then be fitted to the first nontrivial cumulant of the pressure at µ = 0, i.e., the baryon number susceptibility χ B = d 2 p/dµ 2 . In our model, this can be computed at µ = 0 through [22] where boundary was set to be at r = 0 and r h is the location of the horizon.
For each of the V f 0 (λ) and κ(λ) determined above we carry out two additional fits of w(λ) to χ B with lattice data from [41]: one rough fit with three parameters, setting w 1 = 0, and a more precise four parameter fit including also w 1 . Notice that w 0 only affects the normalization of χ B . The results are shown in Fig. 2 (bottom row). The blue curves are the three parameter fits and green curves are four parameter fits. Comparison to higher order cumulants is left for future work.

Holographic thermodynamics and cold QCD matter
Having discussed the thermodynamics of deconfined holographic quark matter at zero chemical potential and finite temperature, we then discuss the predictions of the holographic model at zero temperature and at high baryon density. It is important to check that the result roughly agrees with what is already known about the QCD equation of state at zero temperature. As we will show in this section, the agreement with both nuclear matter EoS (at low densities) and pQCD (at high densities) is remarkably good, given that the only input data of the model are lattice results at µ = 0 and qualitative constraints.
In order to carry out the comparison, we construct a family of equations of state which interpolate between the known results from effective field theory at low densities and pQCD at high densities [19,42,43].

Construction of the interpolated equations of state
We are interested in quiescent neutron stars, so we focus on electrically neutral strongly interacting matter at beta equilibrium. 3 State-of-the-art nuclear physics methods are only reliable up to around the nuclear saturation density n s ≈ 0.16fm −3 , whereas the densities reached at the cores of neutron stars, and where the possible quark matter phase is commonly expected to reside, are few multiples of n s . This means that extrapolating even the well-established low density EoS to the regime of interest comes with huge uncertainty.
We will not choose a particular model but instead consider a whole family of EoSs, closely following [19] (for earlier pioneering work, see [42,43]). The key idea is that the thermodynamic functions can be arbitrary as long as the following conditions are met: • The pressure and the quark chemical potential should be continuous throughout the star and the thermodynamics has to be consistent (i.e., ∂ µq p = n q > 0, ∂ 2 µq p > 0).
• The speed of sound cannot exceed unity anywhere inside the star.
• At low density the EoS should be consistent with experimentally bound phenomenology and conform with chiral effective theory (CET).
• At asymptotically high density the EoS should agree with perturbative QCD.
We first describe how to select the hadronic phase EoS. At low densities we use the EoS obtained from the chiral effective theory following [44,45] (denoted below by HLPS), which is reliable until nuclear saturation density n s . Beyond this the systematic errors will begin to grow significantly, and at 1.1n s they are 24%. We choose this density as the starting point of the interpolation.
On the other asymptotic regime, i.e., at very high densities the EoS we employ for polytropes is the perturbative QCD one [46]. There too, there is theoretical uncertainty coming from the renormalization scale dependence. At quark chemical potential values µ q ≈ 0.87GeV, the uncertainty, similarly to CET, is 24%. This point is then chosen as the end point of the interpolation.
The family of EoS curves that we will consider presents all the possible interpolations between the low and high density limits. This is achieved by the following procedure. We begin by dividing the chemical potential range [µ q (1.1n s ), 0.87GeV] into several intervals. Then on each interval we consider a polytropic EoS of the form p = κn γ which are smoothly joined at the endpoints either to adjacent polytropic Ansatz or on the asymptotic regions. On each of the intervals, the parameters κ and γ can take arbitrary positive values as long as the requirements listed above are met. We will use polytropic Ansätze with four intervals. These "quadrutropes" then act as our proxies for the hadronic matter EoSs and one could continue along the lines of [19] and implement further constraints from astrophysics to figure out what is the family of allowed EoSs that the dense matter could follow.

Comparison between the holographic and interpolated equations of state
The vantage point of holographic models is that the quark matter is strongly interacting in the cores of neutron stars, i.e., at densities which are not asymptotically high. Indeed, this is the regime where we believe V-QCD to be applicable. But in the regime of asymptotically high densities, the predictions of V-QCD may not be reliable. Therefore we do not necessarily expect the thermodynamics of V-QCD to follow the Stefan-Boltzmann (SB) behavior in the extreme densities, where the quarks are expected to become free. It has been shown, however, that for the choice of action which holographically reproduces the logarithmic UV RG flow of the coupling (and which we also use in this work) it is possible to reproduce the SB law both at infinite T and infinite µ [21,22]. In the current article, we have chosen not to match the normalization of the pressure to the asymptotic SB law. This was done because this matching is in slight tension with the fit to lattice data, and we have preferred the quality of the fit over the asymptotic agreement with pQCD. Consequently, the various sets of potentials constructed above will not precisely agree with pQCD at high densities. This should be kept in mind when interpreting the comparison between the interpolated and holographic EoSs which we discuss next.
We present the numerical results for thermodynamics that stem from extrapolating the holographic EoS down to finite chemical potentials and to zero temperature. We depict the pressure of quark matter phase in this regime in Fig. 3 (left). The different curves represent the pressures from all 16 potentials constructed in Appendix A. The four bunches of curves correspond to different choices for the parameter W 0 (see Sec. 3 and Appendix A). The gray band is then generated by parametrizing the low density EoS using quadrutropes as discussed in Sec. 4.1. The parameters were chosen from a set such that we allowed a first order phase transitions to occur anywhere. The bands thus represent all the physically viable EoSs.
Remarkably, we see that the pressures from V-QCD for all choices of the potentials are in the same ballpark as extrapolated nuclear matter pressure at chemical potentials right above the vacuum-nuclear matter transition. All of the holographic curves also intersect the HLPS pressures in this regime. Moreover, the holographic curves, except for potentials 1-3, remain inside the gray band up to relatively high chemical potentials and pressures, where the band is already given by pQCD prediction. The agreement with pQCD is therefore good even if we did not require precise agreement between pQCD and V-QCD when choosing the potentials in Sec. 3. Notice also that none of the V-QCD pressures are positive for µ 370 MeV which is well above the vacuum-nuclear matter transition. Therefore stable quark matter at low densities, and in particular quark stars, are strongly disfavored by the holographic model.
We also show the pressure as a function of energy density in Fig. 3 (right) with the same notation as in the left hand plot. Note that in this plane the quark matter pressures do not intersect any of the HLPS curves. Therefore combining the HLPS equation of states with V-QCD necessarily leads to first order phase transitions -we will discuss this further in the next section. All the quark matter EoSs eventually match up on the pQCD region, depicted as blue band on the right asymptote.
It is noteworthy that all our holographic EoSs are inside the band (apart from the regime of low pressures where they cannot represent the dominant phase as we shall demonstrate in the next section) in the pressure-energy density plane of in Fig. 3 (right) even though they deviate 4 from the band in the left hand plot. This is possible because the mapping from one figure to another is nontrivial: in particular since = µ dp/dµ − p, simple rescalings µ → Cµ move the curves on the left hand plot but leave the right hand plot unchanged. 5 Notice also that as we pointed out above, the parameter that affects the EoS most is W 0 which creates the four families of curves in the left hand plot. At fixed W 0 , the various choices done when fitting the lattice data are clearly a subleading effect, as can be seen from the small spread of each family of curves. The effect of varying W 0 is however similar to rescaling µ: in the right hand plot most of the effect is gone and the families are no longer distinct.
Also the pressures for potentials 1a-3b, which have W 0 = 0, stay inside the interpolated gray band in the pressure-energy density plane. However, we exclude them because of their behavior on the pressure-chemical potential plane. They fail to follow 6 the allowed band already at relative low µ. Therefore the existence of an equation of state which would follow the prediction from V-QCD with these potentials at intermediate µ, i.e., the regime relevant for neutron stars, but still would conform to pQCD thermodynamics at extreme densities looks highly questionable.
We therefore conclude that only the EoSs for the potentials 4-9 are thermodynamically consistent and in the correct ballpark in the sense of conforming to pQCD bound at extreme densities. However, it is immediate that some of these potentials are in tension with astrophysical constraints if a particular model from the hadronic side is chosen. This, together with finite temperature physics associated with neutron star mergers will be discussed in [20]. In the following we are content with discussing the EoS at zero temperature and with drawing lessons of physics interest for quiescent neutron stars.

Hybrid equations of state
In this section we will discuss our results in matching the holographic equations of state to what is known at low densities or constrained from astrophysics. That is, we take the approach [4] suggested by Fig. 3 (left): we use a family of interpolated QCD EoSs to describe the low density baryonic phase, and a family of EoSs determined by V-QCD to describe the high density phase, matching the pressures in the middle where the curves intersect on the p − µ-plane. We aim at making generic predictions and asking if there are some universal results that can be inferred from holographic methods.

Setup
We compute the pressure as a function of chemical potential using our holographic framework, choosing one of the potential sets 4-9, and contrast it with the pressure coming from a given polytropic construction. Unlike in Fig. 3, we use polytropic interpolations without a first order transition, as matching between the polytrope and the holographic pressure will necessarily lead to such a transition. If there is a single phase transition to the holographic quark matter phase at some quark chemical potential ∈ [µ q (1.1n s ), 0.87GeV], the constructed EoS is accepted. In particular we exclude hybrid EoSs with multiple crossings between the deconfined quark matter and the hadronic matter phase. We then check the combinations of all holographic EoSs, using each of the potentials 4-9, with all the 86k interpolated polytrope EoSs, and survey the characteristic numbers for all the acceptable EoS potentials constructed this way.
To summarize, the full pressure as a function of quark chemical potential is then a piecewise hybrid curve: CET up to 1.1n s , then polytropes up until a critical chemical potential µ c , where there is a first order phase transition to a holographic quark matter phase. Notice, we do not consider those polytropes that would intersect the holographic pressures at even higher chemical potentials than 0.87GeV. We expect the pQCD approach to be valid in this regime and, moreover, at such high chemical potentials the corresponding deconfinement phase transition would occur at extremely high baryon densities ( 20n s ) and would presumably leave no imprint for observers.
The last step then is to ask if such an EoS, as constructed above, will be physically plausible, given the astrophysical constraints? At current stage there are two such constraints. The first criterion for the EoS to pass is the ability to support large enough neutron stars. This is tackled as follows. One solves the full Einstein equations for spherical configuration of a self-gravitating perfect fluid with a given energymomentum tensor. These equations are called the Tolman-Oppenheimer-Volkov (TOV) equations and they are ordinary coupled differential equations. The solutions of these equations one obtains the mass-radius (M − R) relationship from where one can infer if a star of a given maximal mass is obtained (and stable). The gravitational mass of the most massive neutron star known is at least 1.97M [47,48]. Another criterion for the viability of a given EoS comes from the gravitational wave observations of merging neutron stars [49]. An updated analysis of [50] states that a neutron star of mass 1.4M should have a tidal deformability in the range 70 ≤ Λ(1.4M ) ≤ 580.
We therefore solve the TOV equations for the all constructed hybrid EoS as well as compute the tidal deformability for a neutron star of mass 1.4M . If we do not find the EoS to produce a neutron star of at least of mass 1.97M or to satisfy the LIGO/Virgo constraint for the tidal deformability, we will discard the corresponding EoS.

Results
The bands spanned by the final hybrid EoSs (which also pass the astrophysical constraints) are shown in Fig. 4. These bands are shown as blue and compared to the gray bands of the interpolated EoSs with a phase transition, considered in Sec. 4. The holographic model restricts the band in the intermediate region between the baryonic phase and the region governed by pQCD, and in particular removes the EoSs with lowest pressures on the p − µ plane. On the p − plane, the excluded gray re- gions below the allowed band are due to the constraint on the maximum mass of the neutron star (low energy densities) and matching with the holographic model (high energy densities). The excluded regime above the allowed band is mostly due to the constraint on tidal deformabilities.
It is generally believed that the baryonic phase cannot exist at densities which exceed the nuclear saturation density by more than roughly by a factor of ten as the baryons would need to be extremely tightly packed. Therefore, we have also shown the red bands where we added the additional constraint that n B /n s < 10 for the hadron phase at the transition. This further reduces the bands, in particular at larger pressures, as the transitions with very high values of µ are excluded. Notice that the blue bands also include the red bands even in those regions of the plots where this does not seem obvious.
One of our most important results is shown in Fig. 5 (left), where we surveyed all the produced EoSs and calculated the latent heats at the critical point. Typical latent heats are sizeable, of the order 1 GeV/fm 3 . The phase transition between the hadronic and quark matter phases at low and intermediate densities, i.e., up to n B /n s 7.5, is strongly first order. As can be verified by studying the solutions of the TOV equations for any of the EoSs in this regime, cold neutron stars cannot sustain quark matter deep in their cores due to the great energy barrier. That is, the slope of the mass-radius curve is positive whenever the center of the star is in the quark matter phase.
In the regime of high densities, i.e., for n B /n s 7.5, we can only set a high upper limit for the latent heat. This is not surprising because as n B grows, the transition is pushed towards higher chemical potentials and the pQCD regime, where Figure 5: Left: We plot the latent heat of the deconfinement transition versus the transition density for all the hybrid EoSs (corresponding to potentials 4-9 and all the quadrutropes). The blue points correspond to those EoSs that survived our criteria of sustaining massive enough neutron stars together with having tidal deformabilities in the accepted range, while the orange points fail either criterion. For n B /n s 7.5, the latent heats are bounded both from below and above. Right: The maximum mass of the neutron star versus the transition density. The plot shows that neutron stars cannot be heavier than 2.6M . The horizontal dot-dashed line shows the minimal allowed value of the maximum mass.
the spread of the interpolated EoSs grows and the holographic model becomes less predictive. Also weakly first order and second order transitions are possible in this region, suggesting that these EoSs can support neutron stars with stable quark matter cores. We have, however, checked the mass radius curves for selected EoSs in the range 7.5 n B /n s 10, and quark matter cores are unstable by a clear margin for all of them. 7 In Fig. 5 (right) we have collected the data for maximum masses of the neutron stars as inferred from the critical pressures of the transitions. We find that neutron stars cannot exceed 2.6 solar masses. There is an additional uncertainty in this plot which does not affect any of our other results: When constructing the polytropes, following [19], we chose as the starting point either the softest or the stiffest HLPS EoS at low densities. The maximum mass of the neutron star also depends on the structure of its crust and outer core which are determined by the low density end of the EoS. Therefore adding EoSs with low density tails between the soft and stiff HLPS EoSs would cause some additional spread of the maximum mass. Including this effect would slightly change the shape of the region with blue dots, mainly at low transition densities.
We show the transition density as a function of the critical quark chemical potential in Fig. 6 (left). For n B /n s < 10, the critical chemical potentials are found to be between 460 MeV and 680 MeV. The spiky structure of the dotted regions at low densities is due to our choices of potentials in the holographic model, in particular the fact that we used potentials with three different values for the parameter W 0 . Letting W 0 vary continuously would smooth out the regions and fill the spaces between the spikes. Fig. 6 (left) shows the central pressure of the star as a function of the transition density.

Discussion
In this paper we took the first steps in using lattice QCD results at zero chemical potential to constrain the various scalar potentials in the V-QCD action. We demonstrated that a very good fit to the QCD pressure and its first cumulant at zero µ was possible, with errors comparable to those of the lattice data. We then projected the resulting thermodynamics in the quark matter phase to the region relevant for quiescent neutron stars, i.e., to vanishingly small temperatures and finite quark chemical potentials.
It was a priori not at all clear how well this approach would work given that this region is well beyond the range of applicability of, e.g., simple minded analytic continuations of the lattice results to finite µ. Despite this we found that after the fit to lattice data, the predictions of V-QCD were already tightly constrained in the region of interest. The dominant leftover freedom was seen to be the choice of the parameter W 0 which controls the relative normalization of the dilaton potential in the flavor term, i.e., in the tachyonic DBI action. Moreover, the resulting pressures from V-QCD for all choices of W 0 compared well with extrapolations of effective field theory results for baryonic matter, and as we demonstrated in Sec. 5, could be used to construct "hybrid" equations of state which meet available constraints. This is a highly nontrivial consistency check of the model. We then continued by analysing the results from the constructed hybrid equations of state. Once we carefully took into account the constraints coming both from the gravitational wave observations and the existence of two-solar-mass neutron stars from Shapiro delay measurements, we found that the holographic model is predictive. We found strong support to the fact that the latent heat of the transition from the hadronic to a three-flavor quark matter phase is well bounded from below at moderate baryon number densities. This then implies that the maximum mass of stable, cold, neutron stars directly follows from the location of the deconfinement phase transition.
It is important to bear in mind that our holographic analysis contained some uncontrolled approximations. Naturally, there is uncertainty from finite N c effects. For the glue sector, we were able to use directly lattice data extrapolated to N c = ∞. For the case of full QCD, the best available data is however at N c = 3, and we used this data to constrain the holographic model. In the Veneziano limit, corrections at finite N c and N f are challenging to analyse on both sides of the gauge/gravity duality. To the contrary, in the case of the thermodynamics of the pure Yang Mills theory the corrections due to finite N c have been analysed on the lattice in [36] and found to be small: already the results at N c = 3 serve as a good model for the limit N c → ∞. Therefore we expect that using N c = 3 lattice data for the thermodynamics of QCD in the Veneziano limit is also reasonable.
As another approximation, we kept the quark masses of all flavors zero on the holographic side while fitting to lattice data with 2+1 quarks at their physical masses. While the perturbation due to light quark masses is expected to be insignificant compared to other uncertainties in our approach, this is not the case for the strange quark mass. Since our holographic model is effective, the error due to setting the quark masses to zero. That is, the effects due to the finite strange quark mass may be captured by changes in the various potential functions as determined by the fit to the lattice data (and the same holds also for the effects due to finite N c ). While we have no way to estimate the size of the errors due to this procedure in the current setup, the fact that we obtain a reasonable model for the thermodynamics at all µ and T suggests that we capture at least some of the quark mass effects without including the masses explicitly.
The above issues motivate ways to make the holographic model increasingly realistic. Addressing the 1/N c corrections would be technically demanding as it requires studying stringy corrections on the holographic side. Therefore the first step would be to add flavor dependent quark masses, and in particular to include the effect due to the strange quark mass. This would generalize the setting with flavor independent quark masses in V-QCD which has already been studied in [8,37,51].
There are also several other possible extensions and related future projects. An ongoing work will carry out more detailed comparison of the model to QCD data [26], including meson and glueball spectra, decay constants, lattice results for higher order cumulants of the pressure, etc. It is likely that this comparison will further narrow the family of allowed holographic EoSs in the regime of low temperatures.
One should keep in mind that "exotic" color superconducting and color/flavor locking phase are expected for QCD at low temperatures and high chemical potentials. Such phases should also be present in the holographic model, as is suggested by the observation that the entropy does not vanish as T → 0 in the relevant region [22].
We have discarded some of the holographic EoSs in the current work because they did not compare well with the results from pQCD at high chemical potentials. This constraint could be made more systematic, e.g., by requiring the existence of smooth polytropic interpolation between the holographic result at intermediate µ and pQCD at high µ.
Moreover, we used poorly constrained interpolated EoSs to model the baryonic phase in the current article. It may be possible to develop a reliable holographic description for this phase also. Unfortunately, holographic baryon physics is technically demanding as baryons are usually realized as soliton configurations of the gravitational dual. An ongoing work discuss baryons in V-QCD in simple approximations [52]. Studying magnetic effects would also be interesting, and relevant in particular to study the physics of magnetars. As a first step in this direction one could analyse the case of constant external magnetic field, which can be relatively easily added in the current holographic setup.
where the first term 8 dominates near the UV (λ → 0) and the second term which represents nonperturbative contributions, dominates in the IR (λ → ∞). Potentials will tend to constants in the UV, and the powers α and β, which control the IR asymptotics, are determined by requiring the model to produce qualitatively different features of QCD [8,10,17,37]. Interestingly, the values of α found by comparing to QCD also typically match the first expectations from string theory.
In the Ansatz (11), many of the UV coefficients are determined by matching the UV dimensions of the operators and their UV RG flow with those found in perturbative QCD, and the IR coefficients remain as fit parameters which need to be compared to lattice QCD data or to experimental data for QCD. The general idea is that near the UV, at weak coupling, holographic predictions may not be reliable but we match our model there with pQCD in order to choose the best possible "boundary conditions" for the more interesting IR physics. In the IR, we can with the functions by comparing to some basic observables and then make predictions for more complicated observables. Relatively simple interpolations for the potentials between their expected UV and IR asymptotics produce a nice fit to the lattice QCD data.
In addition to the expansion parameters in the potentials there are two other parameters: the 5 dimensional Planck mass M and the dynamically generated energy scale Λ UV of the background solutions. The latter may be defined in terms of the UV expansion of the dilaton: where b i are the coefficients of the QCD beta function in the Veneziano limit, Moreover we restrict to zero quark mass in this article, so there is no source for the tachyon field.

A.1 Glue sector
For the dilaton potential we choose where requiring the holographic RG flow of (12) to match with perturbative Yang-Mills theory fixes [9] In order to fit the data we choose a value for λ 0 (which should be close to the characteristic value 4π 2 in perturbation theory) and fit the remaining parameters V IR , M , and Λ UV to lattice data. We find that the following values [25] lead to a good fit, shown in Fig. 1: where T c is the critical temperature of the lattice data. The value of M 3 which reproduces exactly the Stefan-Boltzmann law for the pressure at high temperatures is 1/45π 2 , so our fitted pressure will overshoot this limiting behavior by a factor of 1.3.

A.2 Flavor sector
Let us first discuss the fitting of V f 0 and κ to QCD lattice pressure with 2+1 dynamical flavors at zero chemical potential. As we are considering the chirally symmetric phase, the background solutions are independent of the tachyon kinetic coupling κ(λ). Apart from the parameters of the glue sector, the entropy at µ = 0 therefore only depends on the function V f 0 (λ). In addition, since the tachyonic (thermal gas) vacuum is used as a reference when computing the pressure, it contains a constant which also depends on κ(λ).
Here requiring the UV dimension of the quark mass and theqq operator to match with the asymptotically free flavored theory sets for x f = N f /N c = 1. Moreover requiring the holographic RG flow of (12) to match with that of full QCD with dynamical flavors at x f = 1 gives and requiring the RG flow of the quark mass to agree with the leading order anomalous dimension from pQCD sets κ 1 = 1 3π 2 , (W 0 = 0) ; κ 1 = 11 24π 2 , (W 0 = 0) . As discussed in the main text, we choose values for the parameter W 0 from the set {0, 1, 2.5, 5.886} where the last value leads to high temperature pressure at µ = 0 having the correct flavor dependence (Stefan-Boltzmann law) with x f -independent M . For the temperature of the first order phase transition in the holographic model we pick values from the set {110, 120, 130} MeV, therefore fixing the physical energy scale. We keep λ 0 fixed to the value preferred by the Yang-Mills fit. For each pair of these reference values we then fit the parameters W IR , again the normalization M , and eitherκ 0 orκ 1 as follows.
The decoupling of the tachyon field corresponds to the annihilation of the D4 and D4 branes in the IR (at zero temperature) is required for correct physics such as reasonable meson spectra and correct anomaly structure [15]. For our choice of potentials, this sets a constraintκ 0 > 1.295 [8,17]. Even though no issues appear in the thermodynamic analysis of this article even if this bound is violated, we choose to satisfy it. Already values ofκ 0 very close to the lower bound lead to clear changes in the meson spectrum -therefore we requireκ 0 > 1.5. We fit the data takingκ 1 = 0 if this is possible withκ 0 > 1.5. If this is not the case, we setκ 0 = 1.5 and fitκ 1 to the data. Notice that κ(λ) only affects the thermodynamics through the constant reference pressure of the confined phase, i.e., the thermal gas background which is independent of T and µ, so that fitting more than one parameter of κ(λ) to data does not make sense.
The fit parameters are given in Table 1. The fit was carried out by changing the values of the parameters "by hand" and comparing the result visually to the lattice results, see Fig. 2 (top). For the nonzero values of W 0 we excluded the largest reference value of T c because the quality of these fits was not high enough. The combination 45π 2 M 3 3 /(1 + 7/4) in the last column comes from the normalization with respect to the Stefan-Boltzmann pressure: the listed numerical values indicate how large the pressure is at asymptotically high T compared to the Stefan-Boltzmann (i.e., pQCD)  Table 2: Fit to the baryon number susceptibility at µ = 0: values of the parameters in the function w(λ).
In the Ansatz for W 0 = 0, we could leave out the denominator factor in the term ∝ w 1 since the IR term will be dominant as λ → ∞ even without it, but we choose to keep it for the two Ansätze to be as similar as possible.