{'fraction_non_alphanumeric': 0.03163332100065458, 'fraction_numerical': 0.01727523692973219, 'mean_word_length': 4.526936688950059, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unified weak/strong coupling framework for nuclear matter and neutron stars

Ab initio methods using weakly interacting nucleons give a good description of condensed nuclear matter up to densities comparable to the nuclear saturation density. At higher densities palpable strong interactions between overlapping nucleons become important; we propose that the interactions will continuously switch over to follow a holographic model in this region. In order to implement this, we construct hybrid equations of state (EoSs) where various models are used for low density nuclear matter, and the holographic V-QCD model is used for non-perturbative high density nuclear matter as well as for quark matter. We carefully examine all existing constraints from astrophysics of compact stars and discuss their implications for the hybrid EoSs. Thanks to the stiffness of the V-QCD EoS for nuclear matter, we obtain a large family of viable hybrid EoSs passing the constraints. We find that quark matter cores in neutron stars are unstable due to the strongly first order deconfinement transition, and predict bounds on the tidal deformability as well as on the radius of neutron stars. By relying on universal relations, we also constrain characteristic peak frequencies of gravitational waves produced in neutron star mergers.


Introduction
The era of multimessenger astrophysics has enabled theorists working on neutron stars to scrutinize their models against new available data. In particular, the announcement of the tidal deformability constraint by the LIGO/Virgo collaboration [1,2] has excluded numerous models, some of which were already proposed decades ago. Future observations of neutron star mergers and successfully combining the information coming in from the gravitational waves with the electromagnetic kilonova observations will further narrow down the possible models, with the prospect of nailing down the underlying equation of state (EoS). In addition, constraints to the EoS from spectroscopic measurements of accreting neutron stars, X-ray bursts, or high precision timing observations of rotating stars, are expected to become more accurate in the near future. Discerning the EoS alone, however, will not be fully satisfactory as it does not give direct information on the low energy degrees of freedom. To achieve this, progress in theoretical understanding of dense QCD is needed.
In principle, the theory governing the composition of neutron stars is well-known: it is the Standard Model of particle physics, in particular QCD in the regime of dense matter. However, first principles solutions to the properties of matter in this regime have escaped from theorists due to technical difficulties, due to the famous sign problem of lattice QCD [3] and due to perturbative QCD being invalidated by the sizable value of the gauge coupling [4]. In the current work we will seek another approach and solve a theory which is not QCD but has been phenomenologically tuned to mimic QCD as closely as possible. Our approach is to rely on holographic duality and map the equations of this QCD-like theory to a problem in classical general relativity which is much easier to tackle. The merit of our program is that we can easily provide, e.g., the EoS in the strong coupling regime of cold and dense matter where standard perturbative methods are not applicable.
The holographic approach has been surprisingly successful in describing the hot quarkgluon plasma [5][6][7]. It is hence irresistibly captivating to ponder if holography could also teach us useful lessons at finite baryon density. At face value, this is far from clear as the description of matter at nuclear densities, especially for nucleons, would greatly depend on N c . The limit where the dual gravity framework is solvable rests heavily on taking N c → ∞. Nevertheless, the original work [8] applied the famous N = 4 super Yang-Mills, with quenched fundamental flavors [9], in the context of cold and dense regime of QCD matter showing that the holographic approach continues to fare very well. The microscopic mechanisms of many complex phenomena occurring at finite density quark matter have since been addressed [10][11][12][13].
However, any top-down string theory construction, even if backreaction is taken into account, cannot possibly capture all the rich phenomena of QCD. An alternate approach is to relax the stringent rules stemming from string theory and to take a more phenomenological bottom-up approach: follow the rules of gauge/gravity duality as closely as possible but allow some freedom in model building. Such an effective holographic framework for dense and cold QCD matter has been studied recently in [14][15][16][17] based on models introduced in [18][19][20][21], i.e., improved holographic QCD and V-QCD. The idea is to first tune a relatively complex generic five dimensional gravity model to describe well generic, salient features of QCD and lattice data in the low density region where it is available. The model can then be used to compute precise predictions for dense matter. This approach was seen to work surprisingly well: First, a very precise fit to the lattice data for QCD thermodynamics was seen to be possible [14]. Second, the predicted quark matter EoS at finite density was seen to be tightly constrained by the fit, and in accordance with known constraints both at zero [14] and at finite [16] temperature. Third, including nuclear matter in the same framework by using a simple approximation scheme leads to a stiff EoS (i.e., a high speed of sound) for dense nuclear matter [15], which is what is needed for the EoS to pass known astrophysical constraints [22,23]. This progress has made it possible to construct examples of phenomenologically viable "hybrid" EoSs, 1 where traditional nuclear matter models are used at lowest densities and the holographic model is used for both high density nuclear and quark matter, matching the EoSs from different sources continuously in the middle [17]. Moreover, the mergers of neutron stars governed by such a hybrid EoS were also simulated [17]. The analysis of the produced gravitational wave spectrum therefore opened a new way of comparing predictions of holographic models to experimental data.
In this article we will continue this work by carrying out a detailed analysis of the hybrid EoSs. The goal is to explore all remaining freedom in this construction, arising from the uncertainties in the nuclear matter model as well as the model dependence of the holographic approach, while also taking into account the known constraints to the EoS coming from the LIGO/Virgo observations and neutron star mass measurements. Apart from constraints from the holographic approach and observations, by using a set of nuclear matter models with rather regular EoSs and because of the matching with the holographic model, we also effectively restrict to a family of EoSs with regular behavior: apart from features linked to changes in physics (e.g., phase transitions), the EoSs are continuous monotonic functions without surprising features. Our construction allows us to pin down the predictions of the approach not only for the EoS, but also for the properties of the nuclear to quark matter phase transition as well as for the masses, radii, and deformabilities of the neutron stars.
The rest of this paper is organized as follows. In Sec. 2 we will explain how to construct hybrid equations of state that will conform with known theoretical bounds. At low density we will consider several nuclear physics models that are compatible with chiral effective field theory, the unitary gas conjecture, and are otherwise not too drastic. On the high density side we will use the equations of state following the holographic approach. In Sec. 3 we will first discuss constraints on the equations of state from astrophysics. We will draw lessons of physics interest for neutron star properties and make predictions, e.g., on characteristic frequencies of gravitational waves produced in neutron star binary mergers for future observations. Finally, in Sec. 4, we will summarize our main results and discuss possible outgrowths of our work. The paper is supplemented with two appendices containing technical details and discussion on the adiabatic index.

Hybrid equations of state 2.1 Setup
The equations of state used in this work are constructed from two ingredients, namely a part that comes from a nuclear matter model, and a part that comes from a holographic model. Nuclear matter models are based on well understood physics, where the approximations used are well under control in the regime of small densities. These models can be extrapolated to higher densities, but the underlying assumptions, specifically that the interactions between the nucleons are weak enough to allow a perturbative treatment, eventually break down. In fact, there is a recent body of work [24][25][26][27] showing that non-perturbative effects can lead to drastic differences for the equation of state for both hadronic and quark matter degrees of freedom.
Holographic models, instead, incorporate the assumption that the coupling is strong, and this property can be used in a way such that a holographic description can complement a nuclear matter model. Notice that the breakdown of weakly coupled nuclear matter models takes place at relatively low densities: these models are reliable at best up to about two times the nuclear saturation density [28,29], a density much smaller than the potential transition density from nuclear to quark matter. For realistic holographic modeling of neutron stars it is therefore not enough to just include quark matter, but the dense regime of the nuclear matter phase should also be modeled.
Holographic baryons and nuclear matter have been studied extensively in the literature. It is well established that a single baryon is dual to a soliton configuration of bulk gauge fields, similar to the BPST instanton in Yang-Mills theory [30]. This description has been analyzed in detail in the Witten-Sakai-Sugimoto model [31][32][33][34]; the nuclear matter phase turns out to be an instanton crystal [35,36]. In this article, however, we will employ a simpler approach which avoids the need to consider inhomogeneous configurations and therefore solving partial differential equations. Specifically, we use a homogeneous approximation originally developed in the context of the Witten-Sakai-Sugimoto model [37][38][39][40]. This approximation is expected to work best in the regime of large density, and therefore the holographic homogeneous approach and traditional models of nuclear matter complement each other's regime of validity. By combining the two approaches in their respective regimes into one equation of state with a matching procedure, we can then obtain a hybrid equation of state.
Below, we will describe in more detail the model input, the precise matching procedure, and we will explore how large the set of equations of state is that one obtains in this way.
To determine how generic predictions from such a hybrid construction are, we perform this construction for a wide range of different nuclear matter models and holographic models. For the nuclear matter models, we take care to include both hard and soft choices, as well as several in between, so as to properly parametrize the uncertainty in these models. The models which we use, roughly from soft to stiff, are the following: Hebeler-Lattimer-Pethick-Schwenk (HLPS) soft variation [41], Akmal-Pandharipande-Ravenhall (APR) [42], Skyrme Lyon (SLy) [43,44], HLPS intermediate, IUF [45,46], and DD2 [47]. There are even stiffer EoSs available than DD2 in the range of densities where we will use them, but as we shall see, DD2 is already too stiff for our purposes: the hybrid constructions using it will fail to meet all constraints. 2 It is also important to note that these equations of state are well above the theoretical lower limit derived from the unitary gas conjecture [48]. 3 For the holographic model, we use V-QCD. This is a bottom-up holographic model, which is strongly inspired by string theory and top-down holographic models, but has been phenomenologically tuned to reproduce qualitative features of QCD. The model is in turn based on earlier models: improved holographic QCD (IHQCD) for the gluon sector [18,19], and tachyonic brane constructions for the flavor sector [51][52][53][54][55][56]. We work at zero quark mass. For a detailed exposition of the model, see Appendix A.
For a successful modeling of properties of cold QCD matter it is crucial to carry out a careful comparison to QCD data, in particular to the available lattice data for the thermodynamics of QCD at small densities. For the full V-QCD model this was done in [14] following similar work for IHQCD in [57,58]. A similar approach has been used to study, among other things, the critical point on the QCD phase diagram at finite temperature and density [59][60][61][62][63]. In V-QCD, the data fit was seen to produce a tightly constrained and viable prediction for the EoS of dense quark matter both at zero [14] and at finite [16] temperature. In the current article we will use V-QCD with three sets of potentials arising from the lattice fits of [14], see Appendix A for details.
As the final ingredient for the current setup, the implementation of nuclear matter in the homogeneous approach was studied for V-QCD in [15]. It was shown that this approach produces a physically reasonable phase diagram, and a stiff EoS for dense nuclear matter, i.e., an EoS with a high speed of sound, reaching values well above the value for conformal field theories (c 2 s = 1/3). Notice that the conformal value was long thought to set a ceiling, an obstruction for realistic holographic modeling of dense matter [64,65], but works [11,12] revived the interest in this context. High speed of sound makes it easier to alleviate the bounds to the EoS from the mass measurements of neutron stars [28]. In the current article we will demonstrate that the nuclear matter EoS is stiff independently of the values of those parameters of the holographic model which are left free by the fit to the lattice data.
We stress that, by choosing EoSs coming from specific models both for the low and high density regime and matching them in the middle, rather than using for example polytropic interpolations, we effectively restrict our approach to EoSs which are relatively smooth and regular (apart from the strong nuclear to quark matter phase transition predicted in the model). That is, we avoid EoSs where the speed of sound as a function of density, for example, 2 Notice that the DD2 EoS as such, without matching with holography, is actually too stiff to pass the astrophysical bounds. It could happen though that matching would make the EoS softer at high densities, so that the hybrid DD2+holography EoS would pass the bounds. We will see below that this is not the case. 3 The conjecture posits a lower bound on the energy at a given density. It follows from an idealized state of fermions interacting pairwise solely in the s-wave channel in the limit of infinite scattering length and zero range. Almost all deviations will make the interactions repulsive and hence increase the energy. It would be interesting to extend the recent holographic computations [49,50] to finite density and thus have a first principles value for the scattering length in our context. varies fast without obvious physical reason for the variation. In this sense our approach is similar to that of [23], which used families of analytic interpolating EoSs (including polytropes and piecewise continuous interpolations of the speed of sound), but regulated the EoSs by setting limits to the speed of sound as well as the parameters of the EoSs in consecutive intervals. Notice, however, that, as we shall see below, the use of the rather soft HLPS model (together with the stiff V-QCD predictions) is in slight tension with this idea: there will be a sizable jump in the speed of sound at the matching point when using this EoS at low density.
In summary, the constructed hybrid EoSs then include three regions: • For nuclear matter at low densities, up to densities equal to roughly 2 times the saturation density, we use the various weakly coupled models of nuclear matter.
• For dense nuclear matter, we use the V-QCD model with the homogeneous approach as established in [15].
• For quark matter we use the V-QCD EoSs constructed in [14].
We stress that the last two regions, as well as the nuclear to quark matter transition, are described by the same holographic model. Examples of such hybrid equations of state were already constructed in [17] and used as an input in the numerical simulation of neutron star mergers. In this article we explore how much freedom there is left in this construction, taking into account the uncertainties in the nuclear matter models, the parameter dependence of the holographic model, and astrophysical constraints.

Matching procedure
V-QCD contains several potential functions, which parametrize the remaining freedom in the model. In [14], these potentials have been quantitatively matched to available lattice data, resulting in several possible choices of potentials which are compatible with lattice data. Of these, we use potentials 5b, 7a, and 8b (see Appendix A). In this way, we also parametrize part of the uncertainty in the holographic model.
For each of these choices, the components need to be glued together to obtain a hybrid equation of state. Since the two parts are supposed to be two descriptions of the same matter, we would ideally want the density at which the matching is performed to be a crossover. This is not possible, however, so instead we demand that the transition is as smooth as possible, which is a second order phase transition. For a given density at which we want to perform the matching, we therefore require both the pressure and its first derivative with respect to the chemical potential to be identical for both the nuclear matter part and the holographic part. As these are two parameters, we need two parameters to tune to make this possible. These two parameters are: • The normalization of the pressure c b in the baryonic phase of the holographic model (see Appendix A). This maps to the normalization of the nuclear matter action of V-QCD, and letting it to be a free parameter amounts to taking the speed of sound as the input from the holographic model instead of the pressure itself.
• The parameter b as defined in [15] and in Appendix A. With the choice of potentials for V-QCD as discussed above, there is still this one degree of freedom left, which comes from the Chern-Simons part of the action of the holographic model. This parameter controls the location of the instanton in the bulk and its coupling to the tachyon field, which is dual to the chiral condensate in QCD.
In this way, we obtain for each combination of nuclear matter model and V-QCD potentials a family of equations of state parametrized by the matching density n tr . We choose the values of n tr to lie within the range 1.2n s ≤ n tr ≤ 2.6n s , where n s = 0.16fm −3 is the nuclear saturation density. The lower end of the interval is determined by the constraints: as we shall see only hybrids with n tr 1.4 are soft enough to pass the constraint coming from LIGO/Virgo. The upper end was chosen to lie above the maximal validity range of the nuclear models. For a discussion of the values of b and c b , determined by the matching, see Appendix A.
In summary the hybrid EoSs depend on the choices of the models and parameters as follows: • The low density region up to the nuclear saturation density n s , depends only on the choice of the nuclear matter model, HLPSs, APR, SLy, HLPSi, IUF, or DD2, from soft to stiff.
• The intermediate density region, n s n 2.5n s , depends mostly on the choice of n tr . Because the V-QCD EoSs are (up to few exceptions) stiffer than the nuclear matter models in this regime, increasing n tr means softer EoS.
• In the high density regime, n 2.5n s , the details and stiffness of the model depend mostly on the choice of the parameters of the holographic V-QCD model, i.e., the choice of the potentials 5b, 7a, or 8b (ordered here from soft to stiff). There is, however, also some dependence on the choice of the nuclear matter model and n tr because the parameters b and c b of the nuclear matter sector in the holographic model were determined through the matching procedure.

Results
We then discuss the results obtained by analyzing the hybrid EoSs. Apart from basic thermodynamic properties, we find the structure of nonrotating neutron stars for each of the EoSs by solving the Tolman-Oppenheimer-Volkov (TOV) equations. Some of the most important results are collected in Fig. 1, where the left hand plot depicts the pressure as a function of energy density in log-log scale and the right hand plot shows the mass-radius relations for  (1) and (2) span the light red band, and the excluded EoSs span the striped band. The light blue band represents the polytropic interpolations between the nuclear matter EoSs and pQCD satisfying the astrophysical constraints. Also presented are some examples of hybrid EoSs for the different nuclear matter models used, colored as in Fig. 2 and combined with potential 7a at n tr /n s = 1.9. (Right) The mass-radius cloud spanned by the EoSs following the same color coding as on the left. Also indicated are the 1σ contours for the NICER analysis on PSR J0030+0451 from Ref. [66] with the dashed blue line and from Ref. [67] with the dashed black line, along with the 1σ low-mass X-ray binary fits from the time-evolving X-ray burst spectra for 4U 1702-429 (dark cyan) from Ref. [68], 4U 1724-307 (magenta) and SAX J1810.8-2609 (light green), both from Ref. [69]. The examples of hybrid M − R curves presented correspond to the EoSs on the left.
nonrotating neutron stars. The light red bands are spanned by the viable hybrid EoS which pass known constraints. We will discuss the results shown in these plots in more detail below.

Implementing astrophysical constraints
In this subsection we will discuss which constraints from the astrophysical observations will be relevant in our analysis.
The most precise bounds on the EoS come from the mass measurements of neutron stars. Over the years, convincing evidence for massive neutron stars with M ∼ 2M has been accumulated from various astronomical observations. A compact star may exist in the LMXB (Low Mass X-ray Binary) 4U 1636-536 M = 2.0 ± 0.1M [70], the pulsar B1516+02B in the Globular Cluster M5 M = 2.08 ± 0.19M [71], as well as the millisecond pulsars J1614-2230 M = 1.97 ± 0.04M [72], J0348+0432 M = 2.01 ± 0.04M [73], and J0740+6620 M = 2.14 +0. 10 −0.09 M [74] accurately measured by using Shapiro delay. The existence of these massive stars sets a stringent lower bound to the maximum mass of neutron stars.
The gravitational wave observations by the LIGO/Virgo collaborations have also set con-straints on compact star properties. The (non-)observation of the squishiness of the star constrains the tidal deformability for the compact star of mass 1.4M in the range [2] 580 ≥ Λ(1.4M ) ≥ 70 (GW170817) . (1) This study assumed that both stars obey the same EoS, therefore effectively assuming that neither of the objects in the merger event GW170817 were black holes for which Λ be strictly vanishing (see, however, [75]). The upper bound essentially requires stars of mass 1.4M to be small, so that the tidal forces have very little effect: the matter is kept very tightly close to the center of the star. However, this is a merely a rule of thumb. The presence of a strong phase transition in the outer crust can lead to non-monotonic behavior of the radius R on Λ 1.4M [10].
Of the astrophysical observables relevant for the study of dense EoSs, reliable constraints on NS radii are more difficult to obtain than are the mass and tidal deformability estimates. There are, however, some interesting recent results for the isolated millisecond pulsar PSR J0030+0451 by the NICER collaboration. Both Riley et al. [67] and Miller et al. [66] have independently done Bayesian parameter estimation on the energy-dependent X-ray pulse waveform data with independently developed codes to determine the pulse waveforms and different modeling of the emitting spots on the stellar surface and the instrumental response. Both of these works find consistent estimates for M and R, with 1σ results for both studies shown in Fig There is an ongoing discussion on the atmospheric model used in the analysis to model the emitting hot spots [76], leading perhaps to slightly updated constraints in the future. The additional contours in Fig. 1 are Bayesian parameter fits applying state-of-the-art atmospheric models directly to observed time-evolving thermonuclear X-ray burst cooling tail spectra, obtained by the Rossi X-Ray Timing Explorer for low-mass X-ray binaries 4U 1702-429 (in dark cyan) from Ref. [68], 4U 1724-307 (in magenta) and SAX J1810.8-2609 (in light green), both from Ref. [69].
We then discuss how we take into account these astrophysical bounds in our setup. As we showed above, the most stringent bounds on the maximum mass of the neutron star come from the measurements of J0348+0432 and J0740+6620 [73,74]: the former has significantly smaller error bars, whereas the latter sets a higher lower bound at 1σ level. In the following we will not strictly follow either of these results but we adopt the following rule: we exclude any hybrid EoS if it does not support a star whose mass is at least We will also exclude all EoSs which do not satisfy the LIGO/Virgo bound (1). We will not impose any of the radius measurements directly, but will rather compare our results to them. As we will demonstrate below, the results constrained by maximum mass and tidal deformability are already consistent with the NICER analyses. The measurement of 4U 1702-429 appears more constraining, but one should recall that all the X-ray measurements may contain sizable systematic uncertainties due to the modeling of the neutron star atmosphere.
The effect of imposing the bounds (1) and (2) is shown in Fig. 2, where we plot the tidal deformability Λ(1.4M ) and the maximum mass M max as a function of the matching density n tr . The excluded zones are indicated in the plots by the gray regions, and the thick curves show the EoSs which pass both bounds. We also show the 1σ lower bounds from the mass measurements of J0348+0432 and J0740+6620, in addition to our mass cutoff, for comparison. Both figures exhibit a downward trend, as a higher matching density corresponds to a softer EoS, which supports smaller maximum masses and smaller tidal deformabilities. Thus, because we have an upper limit for Λ(1.4M ) and a lower limit for M max , and the curves in question are monotonic, we can bracket the values of n tr /n s from these astrophysical constraints. The LIGO/Virgo bound (left plot) excludes some of the stiffest EoSs, including all EoSs with n tr /n s < 1.4 and all hybrids with the DD2 nuclear EoSs. The effect of the maximum mass bound is less severe because it has the tendency of ruling out softer equations of state, and the V-QCD nuclear matter EoS is quite stiff. Some of the hybrids with HLPSs and APR nuclear matter models are, however, excluded for large values of n nt . From Fig. 2 we can also note that the non-excluded transition densities differ greatly between the different nuclear matter and potential combinations used. For example, HLPSs with potential 5b has a window of 1.4 < n tr /n s < 1.6, whereas the same nuclear matter EoS combined with potential 8b allows for 1.5 < n tr /n s < 2.1.

Physics lessons: thermodynamics and phase transitions
We then discuss the predictions of our construction to the equation of state and the nuclear to quark matter transition. First, in Fig. 1 (left), the light red band is spanned by all hybrid equations of state which satisfy the constraints of Eqs. (1) and (2). The light blue band is spanned by the quadrutropic interpolations [77] between low density nuclear matter and perturbative QCD results which satisfy the same constraints, 4 while the striped band is spanned by all the hybrid EoS of this article with 1.2n s ≤ n tr ≤ 2.6n s , also those which fail to satisfy the constraints. We also show examples of the hybrid EoSs as the colored curves, using V-QCD with potentials 7a and n tr = 1.9n s while varying the nuclear matter model.
In the leftmost section of the curves and the band, up to the matching density which is where the curves for the examples of EoSs have kinks (i.e., for 300 MeV/fm 3 ), the EoSs are determined by the nuclear matter models. In this regime our construction excludes the stiffest EoSs (having higher pressures). This happens because the V-QCD EoS at higher densities is stiff, which makes the effect of the LIGO/Virgo bound in (1) for the low density part of the EoSs more severe.
In the regime of dense nuclear matter up to the phase transition (which appears as a horizontal line on the curves) the EoS is given by the nuclear matter phase in V-QCD. Because this EoS is stiff and in part due to the maximum mass bound in (2), the softest EoSs are excluded in this region. In addition, the light red band of the hybrid EoSs is narrower than the light blue band of the quadrutrope interpolations, because our construction restricts to "regular" EoSs, excluding combinations of very soft and very stiff sections. As we pointed out above, our approach is therefore similar in spirit to that of [23], which used constraints on the speed of sound to regulate families of analytic interpolating EoSs. Our results in the nuclear matter regime are also similar to this reference. 5 Above the phase transition, the V-QCD quark matter EoS is very tightly constrained and leads to a narrow band. There is more spread though, if the thermodynamic potentials are plotted as a function of the chemical potential, as shown in [14].
In order to analyze the hybrid EoSs in more detail, we show the speed of sound squared c 2 s = dp/d as a function of the baryon number density in Fig. 3. In the top, bottom left, and bottom right plots we vary the nuclear matter model, the matching density n tr , and the potentials of the V-QCD model, respectively. As discussed in Sec. 2.2, varying the different parameters mostly affects the speed of sound in distinct regions: the nuclear matter model affects the low density region, n tr affects the intermediate density region, and the potential choice of V-QCD affects the high density region. The nuclear to quark matter transition density, however, depends on a combination of the parameters. We also note that for all hybrid equations of state, the speed of sound rises well above the conformal value c 2 s = 1/3 in the dense nuclear matter region, and even values c 2 s ≈ 0.6 can be reached with potentials 8b.
The matching between the nuclear matter and holographic model in effect leads to a second order transition, where the speed of sound jumps. Such a jump may be understood as a model for a region of rapid change of the speed of sound as we move from weakly coupled low density nuclear matter to the strongly coupled high density region. The jump is moderate for most of the hybrid models, and even practically absent for some hybrids with the HLPS intermediate EoS, but very drastic for hybrids with the HLPS soft EoS. An example of such an EoS is shown as the thin blue curve in the top row plot. Even though some of the EoSs with most extreme discontinuities (such as the one shown in the plot) are excluded by the maximum mass bound, we expect that the hybrids with the HLPS soft EoS are less likely to be good models of the nuclear matter EoS in this region than the other hybrids, and results based on them should be interpreted with caution.
Our results for the speed of sound have very similar features to the results obtained by using the functional renormalization group approach (see, e.g., [24,78,79] and references therein) in both the nuclear and quark matter phases [25][26][27]. This is particularly interesting as these studies use a somewhat similar approach to the current article: the use of a nonperturbative method in the regime where ab initio calculations cannot be trusted.
We also study the closely related quantity, the adiabatic index γ = d log p/d log = c 2 s /p in Appendix B. We note that it shows much weaker dependence on the parameters of V-QCD than the speed of sound. Near the phase transition we obtain that γ ≈ 1.5 universally for all EoSs, which is a very low value compared to typical predictions of nuclear matter models.
The plots in Fig. 4 show the effect of the constraints imposed on the nuclear matterquark matter first-order transition. In toto, the astrophysical constraints do not limit the QM transition parameters as strongly as they do the astrophysical observables. However, in specific cases, like for potential 7a with HLPSs, the constraint is quite stringent, limiting 4.7 < n b /n s < 6.8 for all allowed n tr (top left plot). 6 None of the constrained hybrid EoSs support a QM transition below 4 saturation densities, and for n tr /n s < 1.9, all constrained hybrid EoSs have n b /n s < 8.
The latent heat of the QM transition for each hybrid EoS is presented in the top right plot in Fig. 4. There we see that all hybrid EoSs produce a ∆ > 700 MeV/fm 3 , implying that the transition is strongly first order for all transition densities and EoSs. For n tr /n s < 2.2 the values of ∆ are bounded both from below and from above, so that 700 MeV/fm 3 < ∆ < 2000 MeV/fm 3 . This, combined with the solutions for the TOV equations, tells us that cold neutron stars cannot contain a QM core for these transition densities because of the energy barrier. This finding is consistent with the analysis of [14] which used V-QCD for quark matter and a model independent approach employing quadrutropes for nuclear matter. It is furthermore expected that the latent heats will decrease with a moderate rate with increasing temperature [16], leaving ample room for quasistable QM core generation in mergers.
Finally, the critical value of the chemical potential is shown in Fig. 4 (bottom). We observe that the critical value depends even more distinctly on the holographic model than, say, n b /n s . The dependence on n tr and on the nuclear matter model, which appears due to the matching procedure, is mostly a small correction. The results fall within the range 470 MeV µ c 680 MeV. The relatively large spread reflects a similar spread of the dependence of the thermodynamic potentials on the chemical potential in the quark matter phase [14] as the parameters of the holographic model are varied. Notice, however, that this dependence mostly cancels in the EoS, which is expressed as a relation between the potentials, = (p).
As a final remark, we notice that the upper bounds of all parameters shown in Fig. 4 are set by the hybrids using the soft variant of the HLPS EoSs for low density nuclear matter. As we have pointed out above, this means that there is a jump from a very soft to a somewhat stiff EoS at the matching density n tr , rendering the hybrid EoS potentially unreliable. The problem is enhanced by the fact that these EoSs (thick dashed blue curves) use the stiffest version of V-QCD with potentials 8b. Therefore the other curves should be regarded as more realistic predictions of our construction.

Physics lessons: neutron star properties
As mentioned above, given a hybrid EoS constructed in Section 2, we can solve the full Einstein equations, which for a spherical configuration of a self-gravitating perfect fluid are called the TOV equations. As a solution of these equations one obtains the mass-radius (M −R) relation for nonrotating neutron stars, examples of which are shown in Fig. 1 (right), along with a band spanned by all of the hybrid EoSs, NICER results for PSR J0030+0451, and a selection of results from measurements of cooling of X-ray bursts. We note that the hybrid EoSs favor relatively large neutron star radii, above 11 km for typical neutron star masses. Moreover, we note that the stars with lowest radii are obtained by using the HLPS soft model at low density, and as noted above, this means that such EoSs have a significant jump in the speed of sound at the matching density. The more regular EoSs (examples of which are shown in the plot) give neutron star radii around 12 or even close to 13 km. We will discuss this in more detail below. The radii being large is a consequence of the V-QCD nuclear matter EoS being stiff, and contrasts results obtained through the systematic use of effective field theory only [80] which gives R(1.4M ) < 11.9 km at 90% confidence level. Notice that the regular EoSs are however in good agreement with the radius measurement of 4U 1702-429 (dashed dark cyan curve in Fig. 1) which has the smallest errors from the X-ray results we have included, and are also consistent with the other direct measurements of radii.
In Fig. 5 we show the radii of the hybrid EoSs at neutron star mass M = 1.4M and M = 2.0M . From the left plot in Fig. 5 we notice, that the radius of a 1.4M NS for all hybrid EoSs is limited to 10.9 km R(1.4M ) 12.8 km . ( The lower limit on the radius of just under 11 km is achieved by HLPSs combined with the stiffest set of potentials 8b in the V-QCD model, meaning the the jump in the speed of sound at the matching density is highest among the hybrid EoSs we have constructed. Other hybrid EoSs produce minimal radii of around 11.5 km to 12.5 km, see 7 Table 1. Thus the constrained hybrid EoSs produce M − R curves that are bundled tightly within the limits of best current limitations on NS radii. It is noteworthy that the lower bound we obtain is larger than R(1.4M ) 10.3 km stemming from the unitary gas conjecture and by the hydrodynamical simulations of binary neutron stars that insist of not having a prompt collapse to a black hole [48,81,82].
For a given EoS and a star of a given mass and radius one can compute the second Love number k 2 and from that, the tidal deformability Λ = (2/3)k 2 [(c 2 /G)R/M ] 5 , which quantifies the effect of an external tidal field on the induced quadrupole moment of the star [83,84]. The tidal deformability Λ is one of the most important parameters describing the gravitational wave signal from neutron star mergers. Apart from characterizing the deviation from the point mass limit for the inspiral part of the signal, it can also be used to roughly predict some properties of the merger and postmerger parts (as we will discuss below).
We show the tidal deformability of a 1.4M NS and the maximum mass supported by the EoSs for various hybrid EoSs as functions of the matching density n tr /n s from NM to V-QCD in Fig. 2 (left). We also show the relation of the Λ parameter to the neutron star radius at three different neutron star masses, M = 1.4M , 1.8M , and 2.0M , in the top left, top right, and bottom plots of Fig. 6, respectively. We notice that the requirement of Eq. (2) also constrains the values of Λ. As an example of a constraint following from our construction, we obtain a lower bound for Λ at M = 1.4M : Here one should note that the lowest values are obtained (again) by the hybrid EoS using the HLPS soft model together with the stiffest (8b) version of V-QCD, so that there is a sizable jump in the stiffness of the EoS at the matching point. Using only more regular hybrid EoSs would push the lower bound close to 300 (see Table 1). Interestingly, these lower bounds are close to the value (∼ 300) obtained by analyzing the electromagnetic signal from the GW170817 event in [85]. Notice that for transition densities between 1.4 ≤ n tr /n s ≤ 2.2, even without constraints from the maximum mass, Λ(1.4M ) 150 for all hybrid EoSs. As is well known there is a strong correlation between the Λ parameter and the radius at fixed mass, as we can see from Fig. 6. Also, tidal deformability decreases fast with increasing mass.
As we commented above, one of the clearest outcomes of the our analysis for neutron stars is that quark matter cores are always unstable, i.e., dM/dR > 0 by a wide margin for solutions containing quark matter. Moreover, we do not find EoSs supporting additional stable branches, i.e., twin stars. The instability of the quark matter cores also limits the maximum mass of the neutron stars, shown in Fig. 2 (right), which marks roughly the point on the M − R curve where the star becomes unstable. For hybrid EoSs with potentials 5b and 7a, M max is set by the star where the central density reaches the phase transition density, so that the instability of the star around M = M max is driven by the phase transition. For potentials 8b, however, the maximum mass is reached within the nuclear matter phase, at densities slightly lower than the transition density. We also notice that the upper limit for Λ(1.4M ) limits the highest masses in the model. E.g., for HLPSs, we notice that for potential 5b the highest maximum mass of around 2.1 solar masses is achieved by n tr /n s ≈ 1.4 and for potential 8b the corresponding values are M max ≈ 2.4 and n tr /n s ≈ 1.5.
Features of the electromagnetic signal related to GW170817 suggest that the remnant collapsed into a black hole soon after the merger. It has been estimated that this sets an upper bound of around 2.2M to the maximal mass M max of nonrotating neutron star [86][87][88][89]. Since M max is sensitive to the parameters of the holographic model, such a bound would severely constrain these parameters. It would exclude most (but not all) hybrid EoSs using the stiffest version of V-QCD considered in the article defined using potentials 8b, see Fig. 2 (right). Almost all the the hybrids with potentials 5b and 7a would however pass this bound. Finally we also show in Table 1 how large a fraction of the neutron star is made out of holographic nuclear matter, i.e., has density above the matching density. This is given in terms of the minimal and maximal values for the ratio of the radius of the "holographic core" to the radius of the whole star. We see that typical numbers for this ratio are around 0.7 at M = 1.4M and around 0.8 or a bit above at M = 2.0M .

Frequencies of gravitational waves from neutron star mergers
In this subsection we discuss some of the key properties of the signal from mergers of neutron stars described by using the hybrid EoSs: the characteristic frequencies of the merger and postmerger parts of the signal. Such frequencies include the peak frequencies f 1 , f 2 , and f 3 of the power spectral density of the postmerger signal, as well as the value of the instantaneous frequency at the time of the merger f mrg [90,91]. Here we will restrict to the frequency of the most prominent peak f 2 , which is linked the rotation frequency of the hypermassive neutron star formed in the merger, and the merger frequency f mrg . Notice that f 2 is, however, absent if the binary promptly collapses into a black hole after the merger. Based on numerical simulations, it has been found that whether or not prompt collapse happens depends, to a good approximation, only on the parameter [92] κ T 2 ≡ 3 We follow here [93] and set the following threshold to single out the cases of prompt collapse to a BH: This limit excludes the mergers with the heaviest masses, including all mergers with average neutron star mass above ∼ 1.5M . For equal neutron star masses, M A = M B , we have κ T 2 = 3Λ/16, and the limit becomes Λ 370, i.e., a value well below the LIGO/Virgo bound of (1).
Comparison of static properties of neutron stars and the results of numerical merger simulations has shown that the frequencies can be estimated to a good accuracy by using universal relations involving the masses and tidal deformabilities of the individual stars [90,91,[93][94][95]. 8 We use here the relations for the frequencies f 2 and f mrg from [93], where they are fitted as rational functions of the variable where ν = M A M B /(M A + M B ) 2 is the symmetric mass ratio and c is also a fit parameter. The results for our hybrid EoSs are shown in Figs. 7 and 8. Fig. 7 shows the predictions for the frequencies f 2 and f mrg for equal mass binaries as a function of the mass of an individual star. Notation is as in Fig. 1: we show the bands spanned by all hybrid EoSs (striped band) and those passing the astrophysical constraints (1) and (2) (light red band) compared to the band for polytropic interpolations satisfying the same constraints, as well as examples of results for hybrid EoSs with n tr = 1.9n s . We also required the bound of (6) for f 2 which explains the sharp cutoff at large masses in the left plot. Notice that the bound contains a sizable uncertainty [92], which is not shown. We did not apply the cutoff for f mrg , which is also well-defined in the case of prompt collapse. The fits of [93] however used data with κ T 2 70, so that the results for large masses and frequencies in the right plot are based on extrapolation. We did not include the uncertainty of the universal relations in the plots -this uncertainty is less than 10% at 90% confidence level [93]. The simulations carried out in [17] typically found frequencies that were a bit smaller than those predicted by the universal relations. These deviations were also less than 10% in all cases. However this suggests that the bands in Fig. 7 (and also below in Fig. 8) systematically slightly overestimate the frequencies for the hybrid EoSs.  We note that the hybrid EoSs in general favor smaller frequencies than generic polytropic interpolations, because the hybrid EoSs are stiffer. This is in agreement with the results of [17] where mergers were simulated by using the hybrid EoS with the SLy model and V-QCD with potentials 7a, and it was found that the characteristic frequencies for the hybrid EoS were significantly lower than for the pure SLy EoS. Notice however that, compared to the polytropic results, our construction also excludes a bunch of stiff EoSs producing lower characteristic frequencies than the hybrid models. The excluded polytropic EoSs are very stiff at low densities but can still meet the LIGO/Virgo bound due to their softness at high densities, well above the saturation density. The hybrid construction does not allow for such softening: all hybrid EoSs are relatively stiff at high densities due to the input from holography.
Notice that at the average massM = 1.4M , corresponding roughly to the gravitational wave event GW170817, only a restricted set of hybrid EoSs passes the constraint (6), producing a frequencies f 2 in a narrow interval in Fig. 7 (left). That is, the binaries with softest hybrid EoSs (including typical EoSs using the soft variant of the HLPS models) are predicted to lead to prompt collapse, whereas stiffer hybrids lead to a (possibly short-lived) remnant neutron star. It has been estimated based on the observed electromagnetic signal from GW170817 that a hypermassive neutron star was formed, which collapsed into a black hole about 100 ms after the merger [100] (see also [101]). Consequently, this observation favors a hybrid EoS with one of the stiffer nuclear matter models, which (as we have argued above) are more regular than those obtained by using the soft HLPS model, and are also favored by the direct radius measurements using the X-ray channel.

Discussion
In this article we explored a framework where "traditional" models of nuclear matter were combined with predictions from a holographic model. In this approach, the low density (weakly coupled) EoS was given by a selection of well-established models of nuclear matter, whereas both the dense nuclear and quark matter EoS was given by the holographic V-QCD model. The "hybrid" EoSs for cold QCD matter were confronted with known constraints from the measurements of neutron star masses and radii, as well as gravitational wave and electromagnetic observations of the GW170817 neutron star binary merger. We found that all known astronomical bounds can be satisfied if the nuclear matter becomes strongly coupled at one-to-two nuclear saturation densities as described by holography.
The essential new ingredient in the EoSs studied in this article was the holographic modeling of nuclear matter, which was carried out by using a homogeneous approach in V-QCD [15,17]. We adopted a method where the EoS of the homogeneous holographic nuclear matter was continuously matched with the low density nuclear matter models. Effectively, this meant taking as an input intensive thermodynamics (e.g., the speed of sound) from the holographic model, while determining the normalization of the extensive thermodynamic potentials from the continuity conditions. We explored the model dependence (remaining after the fit to lattice data [14]) of the speed of sound, and showed that qualitative features were similar in all models (see Fig. 3): the speed of sound squared was seen to be an increasing function of the density with a roughly fixed slope, and would rise well above the conformal value 1/3. In other words, we demonstrated that the model has strong predictive power.
The shape of speed of sound curves in the dense nuclear matter region makes it easy to pass the well-known bounds coming from the Shapiro delay measurements of neutron star masses and the LIGO/Virgo event GW170817. More precisely, the LIGO/Virgo deformability bound for neutron stars having a relatively low mass of 1.4M requires the low density EoS to be soft, whereas the bound of the maximal neutron star mass being at least around 2.0M requires, on average, the EoS at a bit higher densities to be stiffer. Consequently, for a speed of sound increasing rapidly with increasing density, such as the speed of sound predicted by V-QCD, these bounds constrain the EoS only weakly. This leaves, however, constraints due to the functional form of the EoS and more generally from the "regularity" of the EoSs: using a selected model both in the low and high density regions guarantees that the EoS, as well as the speed of sound, is (except in some cases which we discuss below) a relatively smooth, regular, monotonic function apart from features induced by physical reasons such as the nuclear to quark matter transition or the transition from weakly to strongly interacting nuclear matter with increasing density.
The EoS as a whole depends on both the traditional nuclear matter models (at low density) and high the holographic models (at high density). After imposing the astrophysical bounds on the maximum neutron star mass and deformability at 1.4M , we found a family of viable hybrid EoSs, which spans the light red bands in Fig. 1. As we have pointed out, however, this family also contains some EoSs that are not so regular and smooth. In particular, many of the hybrids formed using the soft variation of the HLPS model have a sizable jump of the speed of sound at the matching point, whereas other models (such as SLy, the intermediate version of HLPS, and IUF) typically produce more regular EoSs. Interestingly, we noticed that such more regular EoSs are also slightly favored by X-ray measurements of neutron star radii (included in the left plot of Fig. 1) and the study of the electromagnetic signal related to GW170817 which appears to be consistent with a short-lived hypermassive neutron star remnant [100]. They have relatively large radii, around 12km and even above, see Table 1. Note also that the examples of hybrid EoSs constructed in [17] which were used in the merger simulations in the same article, belong to the latter class of regular EoSs favored by the data.
We also verified the expectation [14][15][16][17] that the deconfinement transition continues to be strongly first order, the neutron star cores are void of QM: the latent heat is over 700 MeV/fm 3 for all viable EoSs, making quark matter cores unstable. Because both dense nuclear and quark matter phases were modeled by V-QCD, the instability of quark matter cores (as well as other properties of the phase transition in Fig. 4) is mainly a prediction of the holographic model. Details in Fig. 4 depend on the low density model only weakly through the matching procedure.
The absence of quark matter in neutron star cores seems to be at odds with the recent work [23] which found strong evidence for the existence of sizable quark matter cores (see also [102]). This evidence was based on studying analytic interpolations for the EoSs between the low and high density regions, where nuclear matter models and perturbative QCD give reliable predictions for the EoS, respectively. One should note, however, that the latent heat of the nuclear to quark matter phase transition in V-QCD is well above the limiting value of 130 MeV/fm 3 , above which quark matter cores were indeed found to be unstable in [23].
Nevertheless, we find examples of hybrid EoSs which pass the astrophysical constraints but do not reach extremely high speeds of sound. The lowest max(c 2 s ) ≈ 0.41 of such EoSs is obtained from the hybrid EoS with the IUF model and V-QCD with potentials 5b at n tr ≈ 1.85n s . This EoS does not meet the regularity conditions of [23] though due to the jump in the speed of sound at the IUF to V-QCD matching point. However, in the hybrids with the intermediate HLPS model and V-QCD with potentials 5b the jump is absent for n tr ≈ 2.05n s , and the EoS still passes the bounds and has max(c 2 s ) ≈ 0.42. These numbers are considerably lower than those reported in [23]. The reason for this is apparently that V-QCD predicts somewhat lower values of the adiabatic index γ (see Appendix B) just below the deconfinement transition than what is classified as nuclear matter in this reference.
Apart from the properties of the phase transition, some other important observables such as the maximum mass of the neutron star in Fig. 2 are more sensitive to the physics of the core than the crust, and hence more to the holographic part of the hybrid EoSs than the nuclear model part. But there are also observables, such as the radii (Fig. 5) and deformabilities (Fig. 6), in particular at low mass, as well as the frequencies of Figs. 7 and 8, which are more sensitive to the low density nuclear matter model than the holographic model. However, interestingly the hybrid construction also adds information about the observables in the latter class: by determining with improved certainty how the models should be extrapolated to higher densities, we are able to study the effects of the astrophysical bounds on the low density region more reliably. In particular, the stiffness of the V-QCD EoS at higher densities means that the LIGO/Virgo deformability bound more severely restricts the stiff low density EoSs, than one might think based on the studies of polytropic interpolating EoSs, for example. As a concrete example we obtain the lower bound of 230 for the tidal deformability at 1.4M in (4), which is pushed up to about 300 if only more "regular" hybrid EoSs are considered (i.e., excluding those with the soft HLPS model and a huge jump in speed of sound at the matching point).
It is important to notice that since all results discussed in this article are functionals of the EoS only. Therefore, based on these results, it is difficult to point out "smoking gun" observables supporting the holographic (or hybrid holography+nuclear matter model) approach. Namely, it is always possible that some other model produces a very similar EoS based on a completely different approach, and therefore almost identical results for the neutron stars, which are therefore "masquerading" [103] as holographic (hybrid) neutron stars. Indeed, as we have pointed out, the functional renormalization group methods seem to lead to strikingly similar results to holography in both the dense nuclear and quark matter phases [25][26][27]. In principle it could be possible to differentiate between different setups by studying observables which cannot be directly derived from the EoS, such as transport properties. The key fact to recall from the holographic approach is that certain computations are relatively easy to perform, which are otherwise out of reach from "traditional" methods, and transport properties are exactly in this class of observables. In particular, all relevant transport coefficients for neutron star mergers in the (heated up) spatially homogeneous quark matter phase of V-QCD recently appeared [104].
A key ingredient in the analysis of this article was the implementation of nuclear matter in V-QCD, which was done by using a very simple approach following [15] where the nuclear matter is dual to a homogeneous bulk field. This approach depends, among other things, on the Chern-Simons term discussed in Appendix A which is not known for general backgrounds. One future project will explore the freedom in determining this term. Another direction would be to go beyond the homogeneous approximation, implement the baryons in terms of inhomogeneous solitonic gauge field configurations in the bulk. A first step would be to solve for the soliton dual to a single baryon, and further steps would consider interactions between the solitons.
Besides improvements in describing the nuclear matter phase, there are several other interesting issues that will require further study. In the current article we set the quark masses to zero, which presumably is a reasonable assumption at not-so-high densities. The generalization of the present work to non-zero flavor independent quark masses is straightforward [105] albeit somewhat tedious. A more challenging and ambitious improvement would be the inclusion of flavor dependence and in particular the strange quark mass. It may affect the comparison to lattice data carried out in [14], also leading to changes in the predictions of the holographic model at high density. Carefully studying this effect is an important future project.
When the chemical potential is increasing and one is creeping toward the deconfinement transition densities, one might also ask if some form of pairing between the quarks [106] could also manifest in V-QCD achievable through a mechanism recently revealed in [107,108] (see also [109][110][111][112][113] for earlier work in this direction). The details of the pairing mechanism is highly sensitive to the differences of the bare quark masses, especially those for the down and strange quarks.
Apart from the effects of the strange quark mass, there are other ways to refine the comparison of the holographic model with lattice and other data. The simple analysis of [14] omitted higher cumulants of the dependence of the lattice EoS on chemical potential, which are already known to a rather good precision [114,115]. The RHIC beam energy scan has also produced data which can be used to constrain the dependence on chemical potential [116]. One could also check if the available lattice data for magnetic effects [117] and conductivities [118] helps to reduce the uncertainties in the model parameters. Moreover, the physics in the low temperature confining phase should be carefully compared to experimental QCD data, in particular on particle spectra.
Finally, at very high densities, the Lorentz invariance of QCD could be spontaneously broken. Famous candidates include the spontaneously generated currents in a phase with a non-zero Kaon condensate [119,120], the spontaneously broken translational symmetry in the crystalline CFL phase [121][122][123], and in particular in the large-N c limit the chiral density wave [124,125]. We similarly expect that the Lorentz invariance might not be intact at very low temperatures in V-QCD, as suggested by several holographic flavor models [126][127][128][129][130], hence a detailed study of the global phase structure as well as of the transport coefficients of the ordered phases [131,132] is direly needed.
In this article we considered, apart from possible indirect effects due to finite temperature through the use of universal relations in the analysis of the frequencies of the gravitational waves in neutron star mergers, only the physics at zero temperature. In order to be able to directly describe the physics of neutron star mergers, where QCD matter is heated up to temperatures comparable to Λ QCD , it would be important to extend our description in all phases to include realistic finite temperature effects. We plan to report on this and also on other previously mentioned extremely interesting topics in future works.

A The holographic V-QCD model
The holographic model we use (V-QCD) [20] is obtained through a fusion of two building blocks: The first is improved holographic QCD (IHQCD) [18,19] which is a holographic model for pure Yang-Mills theory, loosely based on five-dimensional noncritical string theory. The second is a method for adding flavor and the physics of chiral symmetry breaking in such setups [51,53], which is realized through Sen-like tachyonic DBI actions for a pair of space filling D4-D4-branes (see also [55,56] for an analysis of this setup in a different geometry). The gluon and quark sectors are fully backreacted in the Veneziano limit, where one takes both the number of colors N c and the number of flavors to infinity keeping their ratio fixed.
The action of the model has two terms modeling the gluon and quark sectors: where, including terms which are relevant in the quark matter phase and the confined vacuum, Here the dilaton field λ is dual to the TrF 2 operator and sources the 't Hooft coupling on the boundary [18,19]. The "tachyon" field τ is dual to the quark bilinearqq and sources the quark mass on the boundary (which in this formulation is assumed to be flavor independent). It therefore controls the physics of chiral symmetry breaking [19,20,51,53]. In this work we will simply set the quark mass to zero. The quark chemical potential is realized by turning on a temporal component Φ for the gauge field such that F rt = −F tr = Φ (r) and µ = Φ(0), where r is the bulk coordinate.
In addition, we consider nuclear matter in a simple approximation scheme where it is dual to a homogeneous non-Abelian gauge field configuration. To be precise, we restrict to SU (N f = 2) and choose an Ansatz for the spatial components of the non-Abelian gauge fields [15,38,39]: where A L (A R ) are the gauge fields on the D4 (D4) -branes, and σ i are the Pauli matrices. This effectively gives rise to the following additional terms in the Lagrangian, see [15] for details: where the metric was taken to be and In particular, the Chern-Simons term S CS [53], through which the homogeneous Ansatz sources the Abelian gauge field Φ in the bulk. This term depends on the parameter b, which will be determined by matching the EoS with the nuclear matter models as explained in the main text. We also introduced a normalization parameter c b which is also determined by matching with the low density nuclear matter models -such a "correction factor" is needed for physically reasonable predictions for the thermodynamic potentials [15], which reflects the roughness of the homogeneous approximation. This means that we take the speed of sound (rather than the pressure) as a function of chemical potential as the input from holography in the nuclear matter phase. Notice that, as we demonstrate in the main text (see Fig. 3), the predictions for the speed of sound from the holographic model are robust in the nuclear matter phase, which supports this approach. Introducing such an extra parameter may also make sense because, as explained in [15], the nuclear matter was taken to have N f = 2 due to practical reasons whereas the normalization of the quark matter sector was determined using lattice data with 2 + 1 flavors.
The low temperature phases of the model at zero quark mass are then the following [15,21,133]: 1. Confining, chirally broken vacuum phase. This phase has a horizonless geometry ending in an IR singularity, and a nontrivial background configuration for the tachyon field, indicating that chiral symmetry is broken for the boundary theory. The gauge field Φ is constant and the baryon number density is zero.
2. Confining nuclear matter phase. This phase has the same geometry as the vacuum phase, and a nonzero tachyon, but in addition a homogeneous non-Abelian bulk gauge field, which models nuclear matter. The non-Abelian field acts as a bulk source for Φ due to the CS action (13), so that Φ has a nontrivial bulk profile and the baryon number density is therefore nonzero.
3. Deconfined, chirally symmetric quark matter phase. This phase has a geometry with a (planar) charged black hole, and the tachyon field is zero, indicating that chiral symmetry is intact on the boundary. The gauge field Φ is sourced only by the charge behind the horizon, so that the "instanton density" due to (13) is zero, but the baryon number density is still nonzero. In the zero temperature limit the IR geometry becomes AdS 2 × R 3 [21].
For the results in this article, the two latter phases with nonzero baryon number density are relevant.
The potentials V g , V f , κ, and w in the action (9), (10) are first constrained by requiring qualitative agreement with QCD physics [14, 18-20, 105, 134], which essentially fixes the UV (small λ) and IR (large λ) asymptotics of the potentials. The remaining degrees of freedom are then tuned to fit lattice data (at small density) [14,57,135].
We use three sets of potentials 5b, 7a, and 8b in this article, which were originally determined in [14] but require some tuning to be compatible with the nuclear matter setup of [15]. More precisely, [14] used lattice data in the chirally symmetric, high temperature phase where the tachyon is zero. Therefore the fit only directly constrained the potentials V g , V f , and w, whereas the tachyon kinetic coupling κ was tuned to set the critical temperature to the desired value. In [15] it was noted that a slight tuning of the potential κ (which keeps the critical temperature unchanged) was required for some potential sets to obtain physically reasonable nuclear matter configurations.
Here the following UV parameters were fixed to the values obtained by comparing to the UV RG flow of QCD [18,105]:  while the remaining parameters were tuned to match with the lattice data for thermodynamics at low density in the deconfined phase. We take the same values for the parameters of the gluon sector λ 0 = 8π 2 /3 , V IR = 2.05 , which were determined by comparing to lattice data for Yang-Mills thermodynamics at large N c [58], for all fits. The remaining fit parameters for the three potential sets are given in Table 2, where the AdS radius is given by The parameters of Table 2 were grouped in four different groups as follows: • The parameter W 0 was held fixed for each set of potentials, roughly corresponding to a flat direction in the fit.
• The IR parameter of the flavor potential W IR as well as the parameters of the w potential were fitted directly to the lattice data for the interaction measure ( − 3p)/T 4 and the baryon number susceptibility [136,137].
• The parameters of the κ potential were tuned to obtain the desired critical temperature and regular nuclear matter configurations [15]. The parameters are therefore slightly different from [14] for the potentials 5b and 7a.
• The global parameters Λ UV and M were fitted to the lattice data for the interaction measure.
Apart from these parameters, the nuclear matter sector [15] contains two additional parameters: b and the normalization of pressure c b , which are determined by matching the EoS with the various nuclear matter models as explained in the main text. We find that the value of b depends mostly on the choice of potentials. For potentials 5b, we find values within the range 9 < b < 10, for potentials 7a within the range 10 < b < 11, and for potentials 8b within the range 30 < b < 32. The values of 8b may appear unnaturally large; notice however that, as seen from the definition (13), the relevant parameter is actually the scale of the tachyon 1/ √ b. For potentials 8b this is 1/ √ b ≈ 0.18, i.e., still not very small. The normalization parameter c b is most sensitive to the value of n tr . The values which we find lie in the range 2.5 c b 5 for all models. The correction factor of the pressure is therefore clearly larger than the preferred value c b ∼ 1, but not by orders of magnitude. We stress that, as explained in [14], the fit to the lattice data is stiff: while the number of parameters seems large, the results depend on them only weakly. Actually the results are already determined at a qualitative level after fixing the asymptotics of the various potentials, and fitting the parameters, which control the potentials at intermediate values of λ, amounts to small tuning of the results. Moreover, we note that the lattice fit chooses such values of the parameters that all the potentials are regular, monotonic functions of λ. This in part explains why the model still gives very constrained predictions for the EoS.

B The adiabatic index
In this Appendix, we briefly discuss the adiabatic index γ = d log p/d log of the hybrid EoSs. The results are shown in Fig. 9. They should be compared to those for the speed of sound in Fig. 3 as the quantities are related by multiplication by p/ . Notice that the adiabatic index is roughly independent of the choice of potentials in V-QCD, whereas c 2 s shows moderate dependence. The roughness of the curves at low densities reflects the fact that our nuclear matter EoSs were extracted from discrete data which was somewhat sparse in some cases (but still dense enough to make reliable predictions). The structure near n = 0.5n s is due to the change from the crust EoS (which has small uncertainties) to the nuclear matter EoSs. Almost all models show a high peak of γ right above the matching density after which the index decreases fast, reaching values near 1.5 at the first order phase transition.