Building a realistic neutron star from holography

We employ the recently improved description of dense baryonic matter within the Witten-Sakai-Sugimoto model to construct neutron stars. In contrast to previous holographic approaches, the presence of an isospin asymmetry allows us to implement beta equilibrium and electric charge neutrality. As a consequence, we are able to model the crust of the star within the same formalism and compute the location of the crust-core transition dynamically. After showing that a simple pointlike approximation for the baryons fails to satisfy astrophysical constraints, we demonstrate that our improved description does account for neutron stars that meet the current experimental constraints for mass, radius, and tidal deformability. However, we also point out tensions in the parameter fit and large-$N_c$ artifacts and discuss how to potentially resolve them in the future.


I. INTRODUCTION
Dense matter in the interior of a neutron star is notoriously difficult to understand from first principles. The strong coupling nature of the problem -neutron star matter is dense, but not asymptotically dense -makes the AdS/CFT correspondence ("holography") [1,2] a viable theoretical tool. While the holographic principle allows for a rigorous strong-coupling calculation, studies are currently and for the foreseeable future constrained to string models whose dual is more or less different from the relevant underlying field theory, Quantum Chromodynamics (QCD). Moreover, at least in the most accessible approximations, the results are strictly valid only in the strong-coupling limit and for a large number of colors N c of the dual field theory, while N c = 3 in QCD. Nevertheless, for instance in the context of the quark-gluon plasma in heavy-ion collisions [3] or deep inelastic scattering in the so-called small-x regime [4,5], considerable qualitative, and in some cases quantitative, progress has been made with the help of holography, in conjunction with more traditional field-theoretical approaches.
Recently, holographic models have also increasingly been used to understand the regime of large baryon densities and the properties of compact stars. Several studies use holography for a description of quark matter, and combine it with "ordinary" nuclear matter in the outer layers of the star, for instance within a D3/D7 approach [6,7], which can be improved by implementing a running coupling [8,9]. "Bottom-up" models, inspired by but not rigorously based on a string theory, have also been employed, for example the so-called V-QCD model. It is valid in the Veneziano limit, which alleviates the large-N c artifacts, and has mostly been used for quark matter, * Electronic address: nicolas.kovensky@ipht.fr † Electronic address: a.d.poole@soton.ac.uk ‡ Electronic address: a.schmitt@soton.ac.uk but also allows for a description of nuclear matter and the high-density deconfinement phase transition. In combination with known low-density approaches such as chiral effective field theory it has been used to model compact stars, including simulations of neutron star mergers [10][11][12][13][14][15]. For other holographic approaches to compact stars see for instance Refs. [16,17].
In this paper, we employ the Witten-Sakai-Sugimoto model [2,18,19], a "top-down" approach based on type-IIA string theory. We go beyond the previous holographic studies of neutron stars in the sense that our approach allows us to implement equilibrium with respect to the electroweak interaction and electric charge neutrality, two essential conditions for realistic neutron stars that are routinely implemented in field-theoretical approaches. We can therefore construct a mixed phase of our holographic baryonic matter with a gas of leptons, which models the crust of the star. Such a unified approach is highly desired even beyond the realm of holography [20][21][22][23][24][25] and enables us to determine the location of the crust-core boundary fully dynamically.
Baryons in the Witten-Sakai-Sugimoto model correspond to instantons of the gauge theory in the bulk [26,27], and dense baryonic matter has been described in approximations of various degree of sophistication. The pioneering works used a pointlike approximation [28] and -in a complementary approach -a homogeneous ansatz for the gauge fields [29], later refined in Ref. [30]. Improvements to the pointlike approach included a multilayered structure in the holographic direction [31] and the construction of quarkyonic matter [32]. Non-pointlike instantonic matter was considered in Refs. [33,34], including interactions of the instantons in the bulk [35].
Here we shall, firstly, consider the pointlike approximation. It was shown already in Ref. [36] that this approximation cannot account for realistic tidal deformabilities and masses at the same time. We confirm this conclusion and show that the energetically more favorable 2layered solution [31,32] leads to similar results and does not produce realistic neutron stars either. Secondly, and mainly, we then use the homogeneous ansatz of Ref. [29] in the improved version of [37], where an isospin asymmetry for two-flavor baryonic matter was introduced. This improvement provides the setup to compute the thermodynamic properties of baryonic matter for arbitrary baryon and isospin chemical potentials. Thus, by adding a lepton gas we can construct (globally or locally) neutral dense matter in beta equilibrium. The calculation will be performed in the confined geometry of the model with antipodally separated flavor branes, which is the original version of the model introduced by Sakai and Sugimoto and which only has two free parameters, the 't Hooft coupling and the Kaluza-Klein scale. As a consequence of this simple setting, our calculations will be independent of temperature. Also, we shall not include strangeness, quark matter, quarkyonic matter or any form of Cooper pairing of the baryons.
With the equation of state computed from the holographic model and with the help of the Tolman-Oppenheimer-Volkoff (TOV) equations we then compute masses, radii, and tidal deformabilities of the resulting stars and discuss their properties, most notably those of the holographic crust and the comparison of our results with the latest experimental constraints from gravitational waves [38,39], the mass measurement of the heaviest known neutron star [40], and the estimates for neutron star radii from the NICER mission [41][42][43][44].
Our paper is organized as follows. Sec. II is devoted to the pointlike approximation, including the holographic setup in Sec. II A and the evaluation and discussion of the results in Sec. II D. Secs. II B and II C contain the TOV equations and collect the relevant astrophysical constraints, which are also needed for Sec. III. Apart from these subsections, readers only interested in our more realistic approach may go directly to Sec. III, where we discuss the neutron stars constructed from the homogeneous ansatz for baryonic matter. We briefly summarize the holographic setup in Sec. III A and discuss the addition of leptons and the construction of the crust in Secs. III B and III C. The numerical results for the thermodynamics are presented in Sec. III D and for the properties of the resulting stars in Sec. III E. A summary and an outlook are given in Sec. IV.

II. POINTLIKE BARYONS
We start with the discussion of the simplest approximation of holographic baryonic matter, where the instantons in the bulk are approximated by delta peaks and assumed to be non-interacting [28] (this does not mean that the baryons on the field theory side are non-interacting). In this approximation, the number of flavors N f simply appears as a prefactor in the action, and thus the baryons can be thought of as objects composed of N c quarks of a single flavor. Moreover, the baryon onset turns out to be of second order. In other words, matter at a nonzero baryon density always has nonzero positive pressure and thus cannot coexist with the vacuum. Therefore, baryonic matter in this approximation is obviously different from ordinary nuclear matter, which is made of neutrons and protons, and which, in the isospin symmetric case, does have a nonzero saturation density with zero pressure. Despite these shortcomings it is a useful first step to start with this approximation. It allows us to connect our results to the previous literature, and it will be instructive to contrast our more realistic calculation in Sec. III with the results from pointlike baryons.

A. Holographic calculation
The holographic setup needed for this section can be found in previous works [28,32,33,45]. We shall therefore simply collect and only briefly explain the main equations needed to compute the equation of state, without going through the derivations or the details of the model.
In this section, we work with the deconfined background geometry and allow for the D8-and D8-branes ("flavor branes") to have an arbitrary asymptotic separation L, such that together with the 't Hooft coupling λ and the Kaluza-Klein mass M KK the model has 3 free parameters. (We shall later set N f = 2 and N c = 3.) The flavor branes connect in the bulk at the point u = u c , where u is the coordinate of the holographic direction. This is interpreted as spontaneous breaking of chiral symmetry since the gauge theory on the flavor branes corresponds to the global flavor symmetry of the field theory on the boundary. We shall work with zero current quark masses throughout the paper for simplicity (a nonzero mass can be included following Refs. [32,46]). In the following we may restrict our calculation without loss of generality to one half of the connected flavor branes, The action, composed of Dirac-Born-Infeld and Chern-Simons terms, including the source term for the baryons, is where prime denotes derivative with respect to u, V is the spatial volume, T is the temperature, and with λ 0 = λ/(4π), and where is the blackening factor of the background metric, with the location of the horizon u T , related to temperature via 4πT /M KK = 3u T . The action is a functional of the abelian part of the gauge fieldâ 0 (u) and the embedding function x 4 (u) of the flavor branes. The embedding gives the geometric shape of the branes in the x 4 -u subspace of the 10-dimensional background, where the x 4 direction is compactified with a radius given by M −1 KK . Pointlike baryons are placed at the point u = u B , which, in principle, has to be determined dynamically. Due to the symmetry of the two halves of the branes, u B > u c corresponds to a configuration of two layers of baryons. We shall also include the simplest case of a single baryon layer, where the baryons are forced to sit at the tip of the connected branes, u B = u c , as in the original work [28]. We have denoted the baryon number density, which is proportional to the number density of the delta-like instantons, byn B . Following Refs. [32][33][34][35]46], the action is formulated in terms of dimensionless quantities (the only dimensionful quantities in Eq. (1) are V , T , and N ), and the physical, dimensionful baryon number density n B is The equations of motion forâ 0 and x 4 can be solved algebraically forâ 0 and x 4 , where k is an integration constant, and The on-shell action (times T /V ) is identified with the free energy density with the dimensionless version Here we have used that stationarity of the on-shell action with respect ton B yieldŝ and we have subtracted the infinite, medium-independent vacuum contribution 2 7 Λ 7/2 −p 0 , where Λ is an ultraviolet cutoff (which is sent to infinity after the subtraction). The finite contribution has been included in the subtraction to normalize the vacuum pressure to 0, i.e.,Ω = 0 for u T =n B = 0. Here, = M KK L is the dimensionless version of the asymptotic separation of the flavor branes, and we have used that in the vacuum k = u 4 c and The dimensionless (quark) chemical potential is introduced as the boundary value of the abelian part of the gauge field,μ B =â 0 (∞), and is related to its dimensionful (baryon) counterpart by Using Eq. (9), we can thus write the boundary conditions for x 4 andâ 0 as For the minimization with respect to u c one has to treat the 1-layer and 2-layer cases separately. In both cases, this condition can be solved explicitly for k, where η = 1 (η = 0) for 1 layer (2 layers) 1 . In the 2-layer case, the additional condition of minimizing the free energy density with respect to u B yields where we have abbreviated One now proceeds by solving Eqs. (13) with the help of Eqs. (14) and (15). The asymptotic separation can be eliminated by working with the variablesñ B =n B 5 , µ B =μ B 2 ,ũ c = u c 2 ,ũ B = u B 2 , which yields the free energy densityΩ =Ω 7 . Since our choice of the dimensionless variables already absorbs all other model parameters, Eqs. (13) only have to be solved once, without fixing λ, M KK , and L, which only become relevant in the translation to physical results. Having in mind the "decompactified limit" of the model 1, where the confined phase is ignored [32,45], we shall restrict ourselves to zero temperature. In this case, the second-order baryon onset occurs atμ with u c from Eq. (11). This is true for both 1-layer and 2-layer configurations. In the 2-layer case the baryons sit at u B = u c at the onset and move away from this point towards the ultraviolet as the chemical potential (and thus the baryon density) is increased. The free energies of the two configurations are identical only at the onset; as soon as the chemical potential increases, the 2-layer solution is energetically preferred. Multi-layer solutions are expected to take over at large densities [31,32], and they should be taken into account in principle. However, we shall demonstrate that the 1-layer and 2-layer solutions do not differ much in their effect on compact star properties. And, since we shall see that realistic constraints are far from being met within the present approximation, we do not expect multi-layer solutions to change this conclusion. Also, although the pointlike approximation suggests that multi-layer solutions become preferred at large densities, calculations allowing for finite-width instantons indicate that more than 2 layers are never preferred [30,34]. We therefore ignore configurations with more than two baryon layers.
To obtain the equation of state, the solutions of Eqs. (13) are inserted back into the free energy density (8), which yields the pressureP = −Ω and the zerotemperature energy densitȳ =Ω +μ BnB .
Dimensionless pressure and energy density are both translated to their dimensionful counterparts P and with the factor given in Eq. (7). Finally, we will also need the speed of sound Here, the derivative with respect to is in general taken at fixed entropy per particle, and the right-hand side is valid at zero temperature. We show the results for the 1layer and 2-layer configurations in Fig. 1. In both cases the speed of sound is non-monotonic as a function of density and approaches the value c 2 s = 2/5 for asymptotically large densities (in an asymptotically free theory such as QCD the result approaches 1/3). Interestingly, in the 2-layer configuration the speed of sound approaches the asymptotic value much faster. We have indicated the value of the central density within the most massive compact star -to be computed in the following subsections -for both cases, which shows that the non-monotonicity plays no role in the interior of any star.

B. TOV equations and tidal deformability
We now combine the thermodynamics with gravity to construct stars in hydrostatic equilibrium. To this end, we employ the well-known TOV equations [47][48][49], restricting ourselves to the static, spherically symmetric case, where G = 6.709 × 10 −39 GeV −2 is the gravitational constant. The third equation (20c) is added to the TOV equations (20a) and (20b) to compute the tidal deformability [50][51][52][53]. It contains the function y(r), which is related to the metric perturbation from tidal deformations. The equations are solved for mass M (r), pressure P (r), and y(r) as functions of the radial coordinate r after employing an equation of state (P ) and using the boundary conditions in the center of the star M (0) = 0, P (0) = P c , y(0) = 2, where P c is the central pressure. Varying P c yields all possible stars for a given equation of state. Then, P (R) = 0 (the vacuum pressure being set to zero) defines the radius of the star R, and the total gravitational mass of the star is M (R), which we shall simply denote by M in the following. The tidal deformability is then computed from the compactness of the star, and the so-called tidal Love number [54] where y R ≡ y(R). We shall perform the calculation with dimensionless quantitieŝ where we choose the scales r 0 , M 0 , 0 to obey The benefit of this choice is that the dimensionless version of Eqs. (20) does not explicitly depend on r 0 , M 0 , 0 , which only become relevant to transform back to the dimensionful results. Since Eq. (25) leaves one freedom for the three scales, we will choose 0 to our convenience according to the holographic calculation and then determine M 0 and r 0 from that choice.

C. Astrophysical constraints
Here and in the rest of the paper we shall have in mind the following constraints from astrophysical data: • Each mass-radius curve must allow for a mass of about 2.1 solar masses according to the heaviest currently known neutron star [40]; if the compact object from the merger GW190814 with a black hole is a neutron star, the lower limit for the maximal mass is even 2.5 solar masses or higher [39].
• The tidal deformability for a roughly 1.4 solar mass star was constrained by the merger of two neutron stars GW170817 to be about 70 Λ 1.4 580 [38].
• Data from the NICER collaboration was used to estimate the radius of the above mentioned heaviest neutron star to be about (11.4 − 13.7) km [43] and ( [42].

D. Absence of realistic stars with pointlike baryons
In the case of pointlike baryons, we compute the free energy densityΩ =Ω 7 , withΩ defined in Eq. (7). Therefore, the obvious choice for the scale 0 is where we have introduced the energy scale With Eq. (25) we then find where M = 1.988 × 10 30 kg denotes the mass of the sun. We observe that mass, radius, and tidal deformability of the star depend solely on the particular combination of the model parameters in Eq. (27).
Solving the coupled equations (20) in their dimensionless version numerically for different central pressures yields the mass-radius curves shown in Fig. 2. We see that 1-layer and 2-layer phases yield very similar results, with the 2-layer solution allowing for a slightly smaller maximal mass. In the plot we have extended the curves to radii smaller than that of the maximum mass star, although this part of the curve corresponds to unstable stars with respect to radial oscillations [55,56]. As the central pressure is decreased, the radius of the stars asymptotes to a finite value of aboutR 23.9. To understand this we observe (numerically) that for small densities the equation of state behaves like P ∝ 2 . In this regime the Lane-Emden equation is applicable, obtained in the Newtonian limit and with a polytropic equation of state from the TOV equations, see for instance Ref. [57]. From this equation one can check analytically that a quadratic equation of state results in a fixed radius independent of the mass. Figure 2 also includes the tidal deformability as a function of compactness. For both plots no choice of the model parameters was necessary. To connect our results to the astrophysical data we need to discuss their dependence on the energy scale K.
As a first choice, let us fit our model parameters to basic QCD vacuum properties. For instance we may use the pion decay constant f π 93 MeV and interpret the value of the chemical potential at the baryon onset as the vacuum mass of the nucleon m N 939 MeV. Using the expression for the pion decay constant for the deconfined geometry of the Witten-Sakai-Sugimoto model [32,58], reinstating the dimensionful factor for the chemical potential (12) and using Eq. (11), this yields from which we obtain K 0.842 GeV. Here we have set N c = 3 (extrapolating down from large N c ) and N f = 2. We find a maximal mass M max 3.0 M (3.2 M ) for 2 layers (1 layer), a tidal deformability of a 1.4 solar mass star Λ 1. 4 6.6 × 10 4 (8.2 × 10 4 ) and a corresponding radius R 1. 4 30 km (31 km). The results for Λ 1.4 and R 1.4 are far beyond the constraints from astrophysical data, and thus we do not even come close to reproducing vacuum and neutron star properties at the same time within the pointlike approximation.
We may lower our expectations and ask if there is any parameter choice that at least fulfills the astrophysical constraints, ignoring m N and f π . Since the translation of the dimensionless results into physical units only requires the choice of the single scale K, it is easy to explore the entire parameter space. Each K gives a maximal mass M max and a tidal deformability Λ 1.4 . By varying K we can thus construct the curves shown in Fig. 3 for 1-layer and 2-layer configurations. The shaded regions show the constraints from the data, and we see that our curve does not enter the area were both constraints are fulfilled simultaneously. For instance, the largest possi-  ble mass compatible with the realistic band for the tidal deformability is about 1.62 M (1.67 M ) for 2 layers (1 layer), in accordance with the conclusion for 1 layer in Ref. [36]. Since the curves capture all possibilities, we can rigorously conclude that the present approximation, while being useful as a first attempt due to its simplicity, cannot produce realistic neutron stars. We therefore move on to a different approximation of holographic baryons which shall turn out to be more relevant for astrophysical applications.

III. BETA-EQUILIBRATED, ELECTRICALLY NEUTRAL, HOLOGRAPHIC NEUTRON STARS
We have seen that the pointlike approximation of the previous section can only be a first step at best towards a holographic description of dense matter inside neutron stars. More sophisticated approaches in the Witten-Sakai-Sugimoto model based on instantons with a nonzero, dynamically determined, width do exist in the literature. However, so far none of these instantonbased approaches have included isospin, except for studies of single baryons [26]. Since an isospin asymmetry is a crucial ingredient for neutron star matter we will use a recently developed approach which is not explicitly built on instanton solutions but which does account for isospin [37]. This approach makes use of an ansatz for the nonabelian part of the gauge fields that is homogeneous in position space [29,30,33].
Here, we shall only use the simplest version of the for-malism developed in Ref. [37]. Namely, we work in the confined geometry, with antipodal separation of the flavor branes, and employ the Yang-Mills approximation for the action. Moreover, we ignore the possibility of a pion condensate, which has been shown to coexist with baryonic matter in large parts of the phase diagram in the presence of an isospin chemical potential [37]. This is an interesting observation also for applications to compact stars since it is still an open question whether pions form a condensate in neutron star matter [59][60][61]. However, in the approach of Ref. [37] the holographic pion condensate cannot easily be separated from the baryonic contribution, and thus it is not straightforward to couple the system to electromagnetism, i.e., assign the correct electric charges to the nucleons and pions separately. Therefore (and since Ref. [37] only provides the necessary equations in the limit of massless pions), we leave holographic pion condensation in neutron stars for future studies. Also, our restriction to the confined geometry of the Witten-Sakai-Sugimoto model is done for simplicity. The deconfined geometry -used in the previous section for the pointlike baryons -has the advantage to be easily generalizable to nonzero temperatures and to allow for chirally symmetric phases. It may thus be used in principle to study hybrid stars with a quark matter or quarkyonic core. (However, the transition densities to chirally symmetric matter appear to be very large, perhaps too large for neutron star interiors, at least in the approximations currently available in the literature.) For our purposes, the deconfined geometry leads to numerically more challenging equations, which make a systematic study very cumbersome. We therefore work with the simpler equations from the confined background, where in particular the antipodal separation of the flavor branes simplifies the problem because their embedding does not have to be calculated dynamically.
Moreover, in this section we do not include the possibility of more than one baryon layer in the holographic direction. As for the pointlike case of Sec. II, this possibility should, strictly speaking, also be taken into account in the approach of the homogeneous ansatz. It was indeed studied in Ref. [30], albeit within a slightly different homogeneous ansatz, made on the level of the field strengths, not the gauge fields. A generalization of this ansatz to nonzero isospin density, however, has not yet been performed. We recall from Sec. II that within the pointlike approximation the results for 2 layers are very similar to the 1-layer case. Therefore, although a full calculation has to ultimately show the quantitative effect of additional baryon layers, we do not expect them to alter our main conclusions significantly.
Before collecting the relevant equations, let us point out one important shortcoming of the approach used here. Since it does not employ any quantization in the bulk, which can account for finite-N c effects and results in discrete states for neutron and proton [26], the baryon spectrum in the approximation used here is continuous in the isospin direction [37]. As a consequence, symmetric nuclear matter can be thought of as being composed of predominantly isospin-symmetric baryons, and an isospin asymmetry in the many-body system requires the population of heavier, isospin-asymmetric states. This is in contrast to the real world, where isospin-symmetric nuclear matter is already made of baryons with nonzero isospin number, neutrons and protons, and an isospin asymmetry can be created by changing their population, without creating new baryons with different masses. An important consequence is that our holographic nuclear matter has an unrealistically large symmetry energy [37] -of the order of the baryon mass rather than about 30 MeV -and we shall see that as a consequence our neutron star matter will have a very large proton fraction, much closer to symmetric matter than to pure neutron matter. Nevertheless, we shall see that qualitatively, and in many aspects also quantitatively, our approach yields realistic properties of neutron stars.

A. Holographic setup
In this subsection, we collect and briefly comment on the relevant equations of the holographic setup. More details and derivations can be found in Ref. [37]. The action in the Yang-Mills approximation is with the abbreviations As in Sec. II, prime denotes derivative with respect to u. The structure of the action is similar to the one used for the pointlike baryons in the deconfined geometry (1). Besides the absence of the square root due to the Yang-Mills expansion, the different background geometry leads to a different, temperature-independent, metric factor where u KK is the location of the tip of the cigar-shaped u-x 4 subspace, i.e., u ∈ [u KK , ∞], again working on one half of the connected flavor branes. In our units, u KK = 4/9. The action depends on the temporal component of the abelian U (1) gauge fieldâ 0 with boundary conditionμ B =â 0 (∞) and the temporal component of the non-abelian SU (2) part a a 0 σ a , where σ a are the Pauli matrices. The isospin chemical potential acts as a boundary value for a 0 ≡ a 3 0 ,μ I = a 0 (∞), and it is consistent to set a 1 0 = a 2 0 = 0, at least in the chiral limit employed here. The spatial components of the non-abelian part are written as a i = −λ 0 hσ i /2, with a single function h(u) that generates nonzero baryon number through a discontinuity at u = u KK . We work with the same dimensionless units as in Sec. II, such that the physical chemical potentials are obtained via Eq. (12).
The equation of motion forâ 0 can easily be integrated to giveâ where The dimensionless baryon densityn B and the isospin densityn I (given below) are related to their dimensionful versions via Eq. (4). The equations of motion for h and a 0 are These equations have to be solved numerically for a 0 (u) and h(u) with the boundary conditions already given above, together with h(∞) = 0 and a 0 (u KK ) = 0. From the solutions we then compute baryon chemical potential and isospin number densitȳ where h (1) is numerically extracted from the behavior close to the tip of the connected branes, Inserting the solutions back into the action yields the (dimensionless) free energy densitȳ The divergent part of the action is absent from the beginning due to the use of the Yang-Mills approximation, and the vacuum pressure is automatically normalized to zero,Ω B (n B =n I = 0) = 0, i.e., no further vacuum subtraction is necessary. We have added a subscript B to the free energy density to indicate that this is the baryonic part, to which the leptonic part is added now.

B. Adding leptons
In order to account for realistic neutron star matter we need to add leptons, which will serve to neutralize the system. The total (dimensionless) free energy density is thus the sum of baryonic and leptonic contributions, whereΩ =Ω(m e , µ e ) +Ω(m µ , µ µ ) (38) is the free energy density of non-interacting electrons and muons with masses m e = 511 keV, m µ = 106 MeV, with the zero-temperature Fermi gas expression Here we have divided by the factor N f N to use the same dimensionless form as dictated by the holographic setup. The corresponding (dimensionful) lepton density is with We shall require equilibrium with respect to the electroweak processes of beta decay and electron capture, n → p + e +ν e , e + p → n + ν e and the purely leptonic processes e → µ +ν µ + ν e , µ → e +ν e + ν µ . At zero temperature this implies the following conditions for the chemical potentials ("beta equilibrium"), where the neutron and proton chemical potentials are given in terms of baryon and isospin chemical potentials by We shall assume that neutrinos leave the system once they are created, i.e., their mean free path is of the order of or larger than the size of the star, which is a good approximation except for the very early stages in the life of the star or for mergers. Therefore we set µ ν = 0, and beta equilibrium yields In writing down the neutron and proton chemical potentials (43) we have interpreted the two isospin components in the space spanned by 1 and σ 3 as neutron and proton contributions. This seems natural, but we should keep in mind that our formalism does not exhibit actual neutron and proton states, as discussed at the beginning of this section. Within this identification, the proton and neutron densities are given by n B = n n + n p , n I = n n − n p .
Then, assigning the electric charges 0 and +1 to the neutron and proton components, local charge neutrality n p = n e + n µ can be written as where the dimensionless lepton densityn is related to n also by the factor given in Eq. (4). So far, the entire setup is taken -with some notational adjustments for our purposes -from Ref. [37]. For our application to compact stars we also need to compute the energy density, = Ω + µ n n n + µ p n p + µ e n e + µ µ n µ = Ω + µ n (n n + n p ) where, in the second line, beta equilibrium and charge neutrality have been used, and in the third line we have introduced the dimensionless energy densitȳ and the corresponding dimensionful factor With the help of the equations from the previous subsection, we can compute the energy density and the resulting equation of state: for instance, for a givenμ I we solve Eqs. (34), compute the correspondingμ B ,n B ,Ω B via Eqs. (35), (36), and the leptonic part from Eqs. (38) and (44).

C. Constructing a holographic crust
The construction described so far accounts for homogeneous nuclear matter in the interior of a neutron star. The outer layers are, however, expected to have a crystalline structure, which requires a different equation of state. Many previous holographic studies employed crust equations of state (plus a low-density layer of the core) from "traditional" methods and added the holographic part for high densities. This is a sensible way to proceed because the underlying low-density microscopic physics are well understood. Therefore, at least for the outer part of the crust, one might argue that holography is not needed. Nevertheless, here we construct the entire star from holography, which has the advantage of using a single microscopic approach for the entire star. As a consequence, the crust-core transition -where uncertainties in the more traditional approaches already become sizable -can be determined dynamically.
To this end, we construct the crust as follows. Firstly, we observe that our beta-equilibrated, neutral matter shows a first-order transition from the vacuum to nuclear matter. This allows us to construct a mixed phase in the vicinity of this phase transition, where nuclear matter (plus leptons) is spatially separated from a lepton gas. The construction we use here is often employed at the transition between nuclear matter and quark matter in the core, see for instance Refs. [56,[62][63][64]. The resulting structure of the crust is somewhat simplistic. Most notably, our construction is, at best, an approximation for the outer crust -clusters of nucleons immersed in an electron gas -while we shall not attempt to construct an inner crust, where the nucleon clusters coexist with a pure neutron phase, see Ref. [65] for a review of the physics of the neutron star crust. Another simplification we will make is the use of a single geometric structure, namely spherical bubbles of nuclear matter immersed in the lepton gas. Again, this is realistic for the outer crust, but in the inner crust more complicated structures are expected, in particular in the vicinity of the crust-core transition ("nuclear pasta" [23,[66][67][68]). It is possible within the given holographic model to make further refinements along these directions in future studies. However, given the crude approximation our nuclear matter represents to begin with, it seems questionable to refine the details of the inner crust and include different geometric structures before more fundamental improvements have been made.
Let us now turn to the calculation of the equation of state of our crustal mixed phase. We denote the volume fraction occupied by the leptonic phase with χ ∈ [0, 1], such that the volume fraction of the baryonic phase (nucleons plus leptons) is 1 − χ. Then, the conditions for the mixed phase read The first equation is the condition that the pressureand thus the free energy density -of the baryonic phase, Ω B +Ω , be identical to that of the leptonic phase,Ω , where, by construction, both phases must have the same electron chemical potential and thus the sameμ I . Consequently, the pressure in both constituents of the mixed phase is identical to the leptonic pressure, and the pressure of the baryonic component must be zero for any baryon density and any proton fraction that occurs in the crust. The second equation (50b) is the global neutrality condition (allowing each phase on its own to be electrically charged). In order to solve these coupled equations for χ andn B , we first fix a value forμ I . Then, with the help of (35b) and (36) we writen I andΩ B as results of a numerical routine that involves solving the differential equations (34), which in turn also depend onn B . Inserting this routine into Eqs. (50) then allows us to solve them for χ andn B . Afterwards, we can computeμ B from Eq. (35a), the pressureP mix = −Ω , and the energy densitȳ Repeating this procedure for many values ofμ I , we find the equation of state. So far our equations know nothing about the particular geometric structure of the mixed phase. This structure becomes important if surface and Coulomb effects are taken into account. We shall do so in the Wigner-Seitz approximation, where the shape of the unit cell is chosen according to the geometric structure of the mixed phase, e.g., spherical for bubbles. For a given volume fraction χ, the size of the Wigner-Seitz cell is then determined by the competition between surface tension (preferably large unit cells) and electrostatic Coulomb energy (preferably small unit cells). In principle, these effects can be included fully dynamically if the interface profiles are calculated, which also includes screening effects automatically [64,69]. Here we proceed with a simple approximation that assumes the interfaces to be sharp surfaces, with spatially uniform charge density in either phase, as often used for the quark-hadron mixed phase [56]. A brief summary and derivation of this approximation can be found in Ref. [64]; here we simply quote the relevant results.
The cost in free energy from surface and Coulomb effects in Heaviside-Lorentz units is with the charge densities of the baryonic and leptonic phases, respectively, where d is the co-dimension of the geometric structure, d = 1 for slabs, d = 2 for rods, and d = 3 for bubbles. We have chosen χ in Eq. (52) such that it corresponds to the volume fraction of the phase in the outer region of the Wigner-Seitz cell, i.e., if we are interested in bubbles of the baryonic phase immersed in the lepton gas, χ is the volume fraction of the lepton gas, as in Eq. (50). Moreover, Σ in Eq. (52) is the surface tension, which, in the given step-like approximation of the interface profiles, is simply an external parameter. For simplicity we shall assume Σ to be independent of baryon and isospin density. As a benchmark, it is useful to keep in mind that the empirical value for symmetric nuclear matter at saturation is about Σ 1 MeV/fm 2 [70,71] and that the surface tension in the inner crust is expected to become smaller with increasing neutron excess [72]. Due to our simple construction of the crust and since we shall not adopt a density-dependent surface tension, it makes sense for us to vary Σ around its empirical value and check the dependence of our results on this variation.
For the practical calculation we write the dimensionless free energy cost as with the numerical constant A 197.327. This expression can be used to insert the holographic results from Eqs. (50). Then, the equation of state including surface and Coulomb effects is obtained from the pressureP mix,sur =P mix − ∆Ω and the energy densitȳ mix,sur =¯ mix + ∆Ω.

D. Thermodynamics and equation of state
We now evaluate the equations of the previous subsections numerically and discuss the thermodynamic properties of our system. The coupling with gravity and the resulting properties of the holographic stars will be discussed separately in Sec. III E.
For an explicit evaluation we need to choose values for λ and M KK . In the present approach -and in contrast to the pointlike approximation -the equations of motion (34) already depend on λ explicitly; there is no rescaling of the variables that eliminates λ. The equations of motion do not depend explicitly on M KK , but this scale enters when we add the lepton gas (due to the lepton masses) and when we compute the surface effects in the mixed phase (due to the surface tension). In the present section we shall work with λ = 10 and M KK = 949 MeV. This particular choice is not very crucial for now because the main purpose of this section is to point out qualitative properties and establish our construction of the crust. In Sec. III E we shall discuss several parameter sets and also present a systematic study in the λ-M KK plane. To put our parameter choice into context, we have listed three parameter fits in Table I, two of which are taken from the literature and one that is obtained by computing the saturation properties of symmetric nuclear matter within the present setup. This table shows that already without any astrophysical constraints the current version of the model is too simplistic to account for correct QCD vacuum properties and properties of nuclear matter at the same time. We will further discuss this tension in the parameter space in Sec. III E.
In the upper panel of Fig. 4 we present various free energy densities as a function of the neutron chemical po-  Table I: Parameter fits in the confined geometry with maximally separated flavor branes, using pion decay constant fπ and rho meson mass mρ (original work by Sakai and Sugimoto [18,19]), including the QCD string tension σ for a fit with large-Nc lattice data (used in a study of glueball decay rates [73]), and using saturation density n0 = 0.153 fm −3 and binding energy EB = −16 MeV within the homogeneous ansatz for baryonic matter employed in this paper.
tential (Eqs. (12) and (49) have been used to obtain the physical units for the given parameter set). The black, two-valued curve is the result for the pure nuclear matter phase, i.e., homogeneous, beta-equilibrated, locally charge neutral matter, including a metastable/unstable branch with positive free energy density. The other black line is the vacuum, Ω = 0. If we were to ignore the other curves, there would be a first-order phase transition from the vacuum to nuclear matter. The blue curve, with endpoints indicated by the dots, is the free energy density for the mixed phase without surface and Coulomb effects. Where it exists it has lower free energy than the pure phases and thus it replaces the first-order transition with two transitions where the baryon density is continuous -vacuum/mixed phase and mixed phase/nuclear matter. The low-density end of the mixed phase consists of a lepton gas whose volume fraction χ approaches 1 and whose density approaches 0 and a baryonic phase whose volume fraction 1 − χ approaches 0 with n B remaining nonzero (such that the spatially averaged density n B = (1 − χ)n B goes to 0). The high-density end of the mixed phase can be found by solving Eqs. (50) at χ = 0 forn B andμ I . The blue curve represents the unphysical form of the mixed phase because surface and Coulomb effects are not included. The result of these effects is shown by the red curve, which is computed with the help of the energy cost (55), where we have set d = 3 -considering spherical nuclear matter bubbles immersed in a lepton gasand Σ = 1 MeV/fm 2 . Since the red curve is barely distinguishable from the blue curve in the upper panel, we plot the difference of all free energy densities relative to the blue curve in the lower panel. As expected, the added energy cost decreases the interval in µ n where the mixed phase is preferred, but it still survives. We also see that the transition between the vacuum and the mixed phase as well as between the mixed phase and homogeneous baryonic matter is now of first order again (the slope of the free energy curves changes at the transition points). Anticipating our application of these results to compact stars, the low-density transition point corresponds to the surface of the star, while the high-density transition point corresponds to the crust-core transition inside the star. We have indicated these transitions by vertical dashed lines. The corresponding equation of state is shown in the upper panel of Fig. 5. The fact that the realistic mixed phase (red) has nonzero energy density at zero pressure confirms that there is a (very weak) first-order transition between the vacuum and our crust. Also, the (small) jump in at the crust-core transition clearly shows that this transition is of first order. The speed of sound squared and the adiabatic index are shown in the lower panel. The adiabatic index has been used as a criterion for the distinction between nuclear and quark matter [74]. Since weakly interacting quark matter has a smaller γ than low-density nuclear matter, the value γ = 1.75 was, somewhat arbitrarily, chosen as a value for the transition in Ref. [74], where a set of parameterized equations of states are used and thus there is no microscopic criterion for this potentially smooth transition. We find that our holographic nuclear matter persists down to the conformal value γ = 1 and even in the center of the star we find values as small as γ 1.4.
In the low-density regime, we see that both c 2 s and γ show a curious structure within the mixed phase (with a close look a corresponding cusp-like structure can be seen in the equation of state). This is due to the onset of muons, i.e., at that point the electron chemical potential becomes larger than the muon mass. Translated to the stars we shall discuss in the next subsection, this means that muons do appear in the inner part of the crust within our approximation, and not only in the core. The reason 0.97 n0 in the uniform phase (close to but not exactly at the intercept of black and blue curves). As in Fig. 5, the star corresponds to the center of the most massive star, where we find nB 3.45 n0. of this surprisingly "early" appearance of muons is the large proton fraction x p = n p /n B , which can be traced back to the large symmetry energy of our holographic matter, as discussed at the beginning of Sec. III. Since, within our holographic approach, it is energetically very costly to move far away from symmetric nuclear matter, the system decides to keep the proton fraction high, which results in a large electron chemical potential and thus in an early onset of muons. In Sec. III E we shall, for comparison, also consider the case where muons are omitted altogether and we find that our main conclusions are not very sensitive to their presence.
It is thus instructive to compute the proton fraction directly. We show the result as a function of baryon density in Fig. 6. Here we have normalized the baryon density by the saturation density n 0 of symmetric nuclear matter. Within the given parameters, we findn 0 0.089, which corresponds to n 0 0.21 fm −3 (as discussed above, the present parameters are not chosen to reproduce the physical saturation density exactly). While the overall magnitude of the proton fraction is larger than expected in real-world neutron stars, it is worth noticing the qualitative behavior: as expected from a realistic crust, we start off with nearly symmetric matter at very low (spatially averaged) densities and produce more neutron-rich matter as we move towards the core. Within the core, the proton fraction increases before it starts to decrease again at ultra-high densities. This decrease is relevant for the most massive stars, as we have indicated in the figure.
E. Holographic neutron stars

Mass and radius
We now insert the equation of state and the speed of sound into the TOV equations and the equation needed for the tidal deformability (20). The natural choice of the dimensionless pressure and energy density from Eq. (24) isP =P ,ˆ =¯ , such that the scale for the energy density 0 in Eqs. (24) and (25) is given by Eq. (49). (This scale differs from the one used in the pointlike approximation by the factor 7 ; in the present approach does not appear explicitly because the asymptotic separation of the flavor branes is fixed to be antipodal from the beginning.) From this scale and Eq. (25) we determine the scales for mass and radius For now we continue with the parameters of the previous subsection, λ = 10, M KK = 949 MeV. Several massradius curves with these parameters are shown in Fig. 7.
In this and the following plots, segments of the curves representing stars unstable with respect to radial oscillations have been omitted, each curve terminates at its maximal mass. Let us first focus on the black curves, which are all for beta-equilibrated, charge neutral matter, but which are obtained with different descriptions for the crust. All curves give almost the same maximal mass, and this maximal mass is consistent with astrophysical constraints if we assume that the 2.5 solar mass object in the merger GW190814 is not a neutron star (see Sec. II C for the list of astrophysical constraints we discuss here and in the following). If the crust is completely ignored, the mass-radius curve bends back to the origin, resulting in a relatively small radius for, say, a 1.4 solar mass star. This behavior is easy to explain since without the crust there is a first-order transition between our holographic nuclear matter and the vacuum. In other words, nuclear matter at a certain nonzero density can coexist with the vacuum, and that's what happens at the surface of the thus constructed star. This allows for arbitrarily small and light chunks of matter, where gravity does not play any role, close to the origin of the M -R plane, and the curve assumes the typical shape of self-bound stars often encountered for quark stars in simple models (if the crust of the quark star is neglected) [53].
The other extreme, giving rise to much larger radii, is shown by the curve labeled by Σ = 0, meaning we have included the crust but neglected surface and Coulomb effects (blue curves in Figs. 4 -6). This results in unphysically large crusts and overall radii. As we have seen in the previous subsection, surface and Coulomb effects reduce the regime where the mixed phase is favored (red curves in Figs. 4 -6), and as a consequence give rise to smaller, more realistic crusts. Examples are the two mass-radius curves in between the two extremes. While the surface tension Σ = 1 MeV/fm 2 is close to the empirical value for symmetric nuclear matter, we have also included a curve for Σ = 2 MeV/fm 2 to illustrate the effect of a variation in the surface tension, which for simplicity we keep constant throughout the crust for each star of a given curve. We already see that our holographic construction not only reproduces reasonable maximal masses but also, in the most realistic version of the crust, is able to reproduce realistic radii. We will come back to this observation below in a more systematic survey of our parameter space. Strictly speaking, the curves with realistic crust also bend back to the origin for extremely small central pressures. The reason is that even in the presence of the crust (if surface and Coulomb effects are taken into account) there is a (weak) first-order transition to the vacuum, see Figs. 4 and 5. However, this segment of the mass-radius curves is irrelevant for our purposes and thus not included in the figure.
As discussed above, our holographic nuclear matter is plagued by the large-N c artifact of a very large proton fraction. To get an idea of the size of this effect, we have added two more mass-radius curves in Fig. 7 where we fix the proton fraction by hand. We show the result for pure neutron matter (red) -defined in our holographic calculation by n B = n I -and for isospin-symmetric nuclear matter n I = 0 (blue). Both cases are artificial in the sense that the resulting stars are not electrically neutral and not in beta equilibrium. We see that, as expected, our isospin-asymmetric stars with neutrality and beta equilibrium are better approximated by isospin-symmetric matter than by pure neutron matter. The plot suggests that improvements of the holographic model to account for quantized isospin states will move our mass-radius curve more towards the pure neutron matter curve. This would imply that our current approximation somewhat overestimates the maximal mass. It is also worth noticing that the mass-radius curve for pure neutron matter shows the same qualitative behavior of reaching back to the origin as the other cases without crust. The reason is that within our holographic approach even pure neutron matter has a nonzero saturation density at which the pressure is zero. This is in contradiction with results from chiral effective theory, which show that pure neutron matter has positive nonzero pressure for all densities [75]. Given the large-N c artifacts of our isospin spectrum, this discrepancy is not surprising.
In Fig. 8 we only keep the most realistic version of our model, i.e., including a crust with surface tension Σ = 1 MeV/fm 2 , but now we study the dependence of our results on the model parameters λ and M KK . For this discussion we recall the results from Table I, which defines a parameter window where we expect the model to approximately reproduce properties of the QCD vacuum and QCD matter. The choice of the previous plots λ = 10, M KK = 949 MeV lies within this regime, and in Fig. 8 we include the mass-radius curves for three additional parameter pairs chosen such that all curves show more or less realistic maximal masses (one of the curves gives a maximal mass slightly below 2 solar masses). The step size of the variation in λ and M KK is chosen such that the variation in either variable has a similar effect. For instance, starting from the parameter set used in the previous plots (blue curve), the maximal mass is decreased by roughly the same amount by either increasing λ by 1 or increasing M KK by about 50 MeV.  Fig. 8. The dots at the endpoints (hardly distinguishable) mark the most massive stars. The dashed curve is a parameterization constructed to approximate a large number of equations of state taken from Ref. [76]. Lower panel: Ratio of the crust thickness ∆R over the radius R of the star as a function of compactness, for the same four parameter sets, again with dots indicating the 1.4 solar mass stars and the most massive stars. Again, the dashed line shows a simple parameterization, taken from Ref. [77].

Tidal deformability and crust thickness
We may also check the corresponding tidal deformability Λ against the known constraints. This is done in the upper panel of Fig. 9, where we plot Λ as a function of the compactness c. Here we see that for all four parameter sets Λ 1.4 lies in the band given by estimates deduced from the merger GW170817. Given that they also produce realistic masses, this is in stark contrast to the pointlike approximation of Sec. II. Remarkably, all curves Λ(c) are indistinguishable on the given scale (while the value for Λ 1.4 does depend visibly on the parameters). This ap-parent universal behavior is not obvious since here, in contrast to the pointlike approximation, the equation of state even in its dimensionless version explicitly depends on λ (via the holographic equations of motion) and on M KK (via the lepton masses and the surface tension). We have checked that for 't Hooft couplings of about λ = 50 variations of the curve are visible on the chosen scale, although they are still small. It is known that different equations of state show an approximately universal behavior if certain observables are set in relation. We have, for comparison to our results, included a simple parameterization of the function Λ(c) (dashed curve) taken from Eq. (78) in Ref. [76] that is believed to approximate a wide class of equations of states. (See also Ref. [78] for a very similar parameterization.) A similar universal behavior is seen in the lower panel of Fig. 9, which shows ∆R/R, where ∆R = R−R cc is the thickness of the crust, R cc being the radial location of the crust-core transition. Again, the values for M = 1.4 M are marked by dots. The absolute values for crust thicknesses and radii for the four cases in the order (blue, red, black, green) are ∆R 1.4 (2.48, 2.06, 2.06, 1.69) km and R 1.4 (12.9, 11.8, 11.7, 10.7) km. These crust thicknesses are somewhat larger than the usually assumed values centered around 1 km, which, however, are subject to significant uncertainties [21,25,79]. For the most massive stars in each case, also marked by dots, we find the absolute values ∆R = (0.88, 0.82, 0.83, 0.77) km. Again we show a simple parameterization, taken from Eqs. (B6) and (B8) of Ref. [77], that is believed to be a good approximation independent of the exact behavior of the equation of state. We see that although the qualitative behavior is the same, our results for ∆R/R are significantly larger than those given by this parameterization.
In Fig. 10, in order to further put our results into context, we compare the four equations of state used for Figs. 8 and 9 with the "allowed" band of Ref. [80] 2 . At low densities, n B 1.1 n 0 , this band contains the results for homogeneous, beta-equilibrated, charge neutral nuclear matter, obtained from a phenomenological extrapolation of pure neutron matter from chiral effective field theory [75,81], and shown to behave similarly to a widely used equation of state for the crust [82,83]. At extremely high densities, n B 40 n 0 , the band represents the predictions from perturbative QCD [84,85]. The intermediate regime is constructed from an interpolation using piecewise polytropic equations of state obeying thermodynamic consistency, causality, and the astrophysical constraints M max > 2 M and Λ 1.4 < 580. We see that our crust deviates from the low-density band, while the intermediate and high-density parts are mostly within the allowed region. The curve that slightly violates the band at intermediate densities corresponds pre-  Figure 10: Equations of state corresponding to the four parameter sets of Fig. 8 and 9, compared with the "allowed" band of Ref. [80], defined by low-density chiral effective theory, high-density perturbative QCD, and polytropic interpolations between them, constrained by astrophysical observations. The dots mark the 1.4 solar mass stars and the most massive stars for each equation of state.
cisely to the parameter set that does not produce stars above two solar masses, i.e., this behavior is in agreement with the findings of Ref. [80]. All our curves continue within the allowed region even beyond the highest densities in the most massive stars, before they deviate from the band at ultra-high densities (our curves stop where the numerical evaluation becomes problematic). Since our holographic model does not account for asymptotic freedom, it is no surprise that the weak-coupling results are not reproduced. More importantly, we conclude that the asymptotic constraints do not invalidate our results: it is conceivable that our holographic equation of state approximates real-world QCD up to (and even somewhat past) the highest densities in neutron stars, and beyond that, where we do not intend to apply our results anyway, it is possible to connect them to the first-principle QCD calculation without violating any theoretical or astrophysical constraints.

Discussion of the parameter space
Finally, let us discuss our parameter space more systematically. To this end, we determine the regions in the λ-M KK plane where astrophysical constraints on mass and tidal deformability are obeyed, see Fig. 11. Before we discuss the results let us briefly explain how this figure was produced. For simplicity, here we have neglected the electron mass (which makes almost no difference for any observables we discuss) and ignored muons (which does make a small difference, but does not affect the main conclusions). With these simplifications, the holographic calculation itself, which is the most time-consuming part of our numerics, becomes independent of M KK , and we  Table I: fit to mρ and fπ (circle), to mρ and σ (diamond), and to the saturation properties of nuclear matter (square). In this plot, muons are neglected for simplicity.
only have to do it once for each λ. The energy cost from Coulomb and surface effects unavoidably induces a dependence on M KK . Adding this cost within our approximation to obtain the equation of state is trivial, but we have to solve the TOV equations -which is numerically less demanding than the holographic part of the calculation -on a suitably fine grid in the λ-M KK plane for each grid point separately. In this way we can, for each λ, determine the value M KK needed to obtain a given maximal mass and a given tidal deformability of a 1.4 solar mass star. As a result, we can identify the windows where the maximal mass is larger than 2.1 solar masses (shaded red) and where the tidal deformability is within the bounds 70 < Λ 1.4 < 580 (shaded blue). For completeness we have also added the contour for the maximal mass of 2.5 solar masses (red dashed curve), having in mind the possibility of an ultra-heavy neutron star in the merger event GW190814 (although a general upper bound for the maximal mass below that value was suggested [86]). We have also used three different values of the surface tension. This creates a small variation in the deformability and essentially no change in the maximal mass, see right panel for a zoom-in, where the (blue) lines for the three different values of Σ become distinguishable. The main observation of this figure is that there is a region in the parameter space where the constraints of the 2.1 solar masses and the tidal deformability can be met, as already seen for selected parameters in the previous plots. However, a 2.5 solar mass star is more difficult to reconcile with the constraint for Λ 1.4 , as the red dashed line barely enters the blue area.
We have indicated the three particular parameter choices of Table I in the plot. We see that none of them lies in the overlap region of the red and blue areas. For the interpretation of this result we should keep in mind that they already disagree with each other, i.e., had we looked for a parameter set within the given approximation that can fit "everything" it would not have been necessary to calculate properties of compact stars. On the other hand, we may use Table I to define a window in parameter space in which QCD properties of the vacuum and of nuclear matter at saturation can be reproduced simultaneously at least in an approximate way. It is reassuring that this window overlaps with the one defined by the astrophysical constraints in Fig. 11. We thus conclude that the Witten-Sakai-Sugimoto model -even in its simplest version and by fitting only two parameterscan account approximately for observables in regimes as diverse as the vacuum, nuclear matter at saturation, and basic properties of neutron stars. Moreover, the small number of parameters distinguishes our approach from many phenomenological models used for dense matter in neutron stars. This allows us also to ask whether we can make some parameter-independent observations, as we shall discuss next.
It is striking that the curves in Fig. 11 appear to be straight lines to a good approximation on the doubly logarithmic plot (one can check that the allowed region is roughly approximated by M KK √ λ ∼ (3.0 − 3.2) GeV). In particular, lines of constant maximal mass appear to be also lines of constant tidal deformability Λ 1.4 . To further quantify this observation, in Fig. 12 we show Λ 1.4 and the corresponding radius R 1.4 as a function of the 't Hooft coupling along the lines of constant maximal mass 2.1 M and 2.5 M , i.e., moving along the red solid and red dashed lines in Fig. 11. We see clearly that Λ 1.4 and R 1.4 are not exactly constant, but their variation is intriguingly small, at least in the shown range, omitting values of the 't Hooft coupling which are orders of magnitude larger or smaller than λ ∼ 10, which we know is the regime where the model reproduces properties of QCD. We can thus use these plots to make approximate parameter-independent predictions of the model. For instance, the model predicts that an equation of state whose most massive star has a mass of 2.1 solar masses will produce a 1.4 solar mass star with tidal deformability Λ 1.4 ∼ 230 − 350 and radius R 1.4 ∼ (11.4 − 13.5) km (upper and lower limits of the blue bands shown in Fig.  12). Interestingly, this radius prediction is in very close agreement with the range given in one of the two interpretations of the NICER data [41]. Since we know from astrophysical observations that 2.1 M is the lower limit for the maximal mass, the plot also demonstrates that the lower boundaries of the blue bands are absolute lower bounds predicted by our model. It is interesting to note that the lower bound for the tidal deformability, Λ 1.4 230, is identical to the lower bound found in the V-QCD approach [15]. The figure also illustrates once more the tension of producing stars above 2.5 solar masses: as the left panel shows, it is possible to produce them in a certain regime in λ, but with values for Λ 1.4 at the very upper end of the allowed window and only with sufficiently large, perhaps unrealistic, values of the surface tension.

IV. SUMMARY AND OUTLOOK
We have constructed neutron stars from the holographic Witten-Sakai-Sugimoto model. Our approach employs the same holographic formalism from the densest matter in the center of the star up to the low-density surface. This includes a construction of the crust, which had not been done before within holography. Our holographic nuclear matter is based on a spatially homogeneous ansatz for the non-abelian gauge fields in the bulk, which yields an approximation for a many-baryon system that is expected to work best at high densities. In particular, we have used a recent extension of this approach to isospin-asymmetric matter. After adding a non-interacting lepton gas to our holographic setup, this allows us to incorporate equilibrium with respect to the electroweak interaction and electric charge neutrality. The crust is then constructed in the Wigner-Seitz approximation as a mixed phase of the lepton gas and nuclear matter, assuming sharp, step-like interfaces between the two phases. The comparison of the free energies of pure and mixed phases is used to determine the crust-core transition dynamically.
After showing that a simple pointlike approximation of the holographic baryons is not in agreement with astrophysical data, we have computed mass-radius relations and tidal deformabilities of the neutron stars constructed from the homogeneous ansatz. Our calculation is performed in the confined geometry of the model with antipodal separation of the flavor branes. In this scenario, the model only has two parameters, the 't Hooft coupling λ and the Kaluza-Klein mass M KK , supplemented in our approach by the surface tension of nuclear matter, which we treat as an external parameter. We have shown that using a realistic surface tension the model parameters can be chosen to meet the known astrophysical constraints on maximal mass, tidal deformability and radius of the star.
We have systematically studied the parameter space to compute the astrophysically allowed window in the λ-M KK plane. In this simple version of the model, the corresponding parameters reproduce approximately, but not exactly, the vacuum properties used in the original works as a fit. This tension in fitting to different properties is already apparent if the saturation properties of nuclear matter are taken into account. We have also extracted parameter-independent predictions of the model. For instance, we have shown that any parameter set that produces a maximum mass of 2.1 solar masses or higher produces a radius for a 1.4 solar mass star larger than 11.4 km and a tidal deformability larger than 230. More, and more stringent, bounds can be obtained by a more exhaustive analysis, also making use of the NICER data for the radius of a 2.1-solar-mass star. This was done in our follow-up study [87].
In constructing the entire neutron star from a single model our work goes beyond previous holographic -and most traditional -approaches. Nevertheless, several improvements are desirable. Most notably, as discussed in the main part of the paper, we know that our approach shows large-N c artifacts in the isospin spectrum. While we can define neutron and proton number and assign electric charges to them, additional states in the continuous isospin spectrum unavoidably become relevant in the simple approximation used here. A very large symmetry energy and a resulting large proton fraction are manifestations of this artifact and an important step would be to improve the treatment of isospin, for instance by including the quantization in the bulk of the holographic setup [26], which however is difficult to generalize to dense matter. Moreover, our construction of the crust is simplistic in the sense that we did not include an inner crust, which contains a mixture of nearsymmetric nuclear matter and a pure neutron fluid. One can try to construct such an inner crust within the setup used here. Other extensions would be to compute the surface tension dynamically within the given model, or to include pion condensation, which may coexist with nuclear matter in the relevant part of the phase diagram within the present approximation [37]. It is also possible to extend our approach to include strangeness and check whether hyperons appear in the stars constructed from our holographic model. Also, one can straightforwardlybut with more numerical effort -extend our study to the deconfined geometry of the model. The advantage would be that a nontrivial temperature dependence can be included, relevant for the simulation of neutron star mergers, and that the transition to a chirally restored phase could be included. This would be interesting for a hybrid star containing quark matter (or quarkyonic matter [32]) in the core. However, in currently employed approximations the chiral phase transition occurs at relatively large densities, perhaps too large to be relevant for the interior of neutron stars. One could also consider computing transport properties of our neutron star matter. Finally, it would be interesting to see if our construction of the holographic crust is also applicable to similar holographic models.