Towards a holographic quark-hadron continuity

We study dense nuclear and quark matter within a single microscopic approach, namely the holographic Sakai-Sugimoto model. Nuclear matter is described via instantons in the bulk, and we show that instanton interactions are crucial for a continuous connection of chirally broken and chirally symmetric phases. The continuous path from nuclear to quark matter includes metastable and unstable stationary points of the potential, while the actual chiral phase transition remains of first order, as in earlier approximations. We show that the model parameters can be chosen to reproduce low-density properties of nuclear matter and observe a non-monotonic behavior of the speed of sound as a function of the baryon chemical potential, as suggested by constraints from QCD and astrophysics.


Contents
1 Introduction

Context and purpose
If nuclear matter is compressed to sufficiently large densities it turns into quark matter, which is weakly interacting at asymptotically large densities. The nature of this transition is unknown because first-principle calculations within Quantum Chromodynamics (QCD) are currently inaccessible for matter at low temperature and large baryon density due to the strong-coupling nature of the problem and the so-called sign problem of lattice gauge theory [1]. In a simple picture, extreme compression makes the nucleons overlap, and at some point it is the quarks, no longer the nucleons, which are the relevant degrees of freedom. Since there is no strict order parameter for confinement and chiral symmetry breaking, it appears that there is no qualitative difference between low-density nuclear matter and ultra-dense quark matter, such that this transition is allowed to be continuous. This scenario is realized at zero baryon density, where we do know from first principles that there is a crossover from the hadronic phase to the quark-gluon plasma. At large baryon densities, the situation is more complicated due to the presence of Cooper pairing. Since asymptotically dense and sufficiently cold matter is in the color-flavor locked (CFL) phase [2,3], the actual question is whether one can go continuously from nuclear matter to CFL rather than to unpaired quark matter [4][5][6][7][8][9][10]. Here we will ignore Cooper pairing and address the question of a possible continuity between unpaired, isospin-symmetric nuclear matter and unpaired quark matter, as for instance envisioned within a percolation picture [11,12].
Combining nuclear matter and quark matter in the same calculation on a microscopic level is very challenging, even if the rigor of a first-principle calculation is given up. However, besides the theoretical motivation just discussed, such studies are highly desired in the context of neutron stars, where for instance mass and radius, but also gravitational wave signals from neutron star mergers, depend on thermodynamic properties of the matter in the core of the star and on the possible existence of a first-order phase transition between nuclear matter and quark matter [13][14][15][16]. Therefore, various studies exist which combine two distinct models, one with nucleons as degrees of freedom, one with quarks as degree of freedom. They are then either glued together at a first-order phase transition, or a smooth interpolation is introduced between the two phases [17][18][19]. In either case, these approaches have no predictive power for the nature and the location of the quark-hadron phase transition. Other studies attempt to find a sufficiently general parameterization of the equation of state without relying on a microscopic theory, using microscopic and/or astrophysical input as constraints for the parameterization [13,20,21]. Another strategy is to use a single model that has either quark degrees of freedom or nucleonic degrees of freedom, but not both. In that case, potential chiral and deconfinement phase transitions are determined consistently within the given model, but at least one side of the phase transition is then only a very rough approximation to the real world [22][23][24][25]. An example is the Nambu-Jona-Lasinio (NJL) model [26][27][28], which usually does not contain baryons, although it is conceivable to include them [29][30][31]. Our present approach, based on the gauge/gravity duality, is at best a distorted version of QCD, in some sense similar to an NJL model, which however does have a well-defined concept of baryons in a chirally broken phase and of chirally symmetric quark matter. Our main results are qualitative and thus at this point we do not attempt to make quantitative predictions for the physics of neutron stars. We do make a step in this direction, however, by fitting the parameters of our holographic model to properties of nuclear matter at saturation and by calculating the speed of sound for all densities. The results indicate that -together with future improvements -our approach can be a useful alternative to the models currently employed for dense matter in neutron stars. There exist other recent studies that have connected results from the gauge/gravity duality to neutron star physics [32][33][34]. However, these studies have combined a holographic equation of state for quark matter with a traditional approach for nuclear matter and thus fall in the above category of combining two distinct theories.

Method
We employ the gauge/gravity correspondence ("holography") which allows for a non-perturbative calculation of strong-coupling physics via the classical gravity approximation of certain string theories [35][36][37]. Specifically, we work with the model developed by Sakai and Sugimoto [38,39], who introduced fundamental chiral fermions through N f D8 and D8 branes ("flavor branes") into Witten's original model for low-energy QCD based on a background of N c D4 branes [40], see Refs. [41][42][43] for reviews. Here, N f and N c are the numbers of flavors and colors, respectively, and the model is usually evaluated in the limit N f N c , which is a necessary condition to neglect effects of the flavor branes on the background geometry. We shall work in the so-called decompactified limit of the model, where the background is given by the deconfined geometry. In this limit, the gluon dynamics are decoupled and the dual field theory is expected to resemble the NJL model in certain aspects [44][45][46]. The price we have to pay for employing this limit is that the rigorous connection to large-N c QCD is lost -which in any case only exists in the limit of small 't Hooft coupling, which is inaccessible within the gravity approximation. But, in the decompactified limit, the chiral phase transition depends nontrivially on the baryon chemical potential even if changes of the background geometry are neglected. This results in a richer phase structure which may well be closer to real-world N c = 3 QCD, at least with regard to the chiral phase transition.
Baryons are introduced following the general recipe of the gauge/gravity correspondence, i.e., through branes that wrap around an internal subspace of the ten-dimensional geometry and that have N c string endpoints. In the non-supersymmetric Sakai-Sugimoto model, this picture becomes equivalent to instanton solutions of the gauge theory on the flavor branes [47]. Here, "instanton" refers to an object localized in position space and the holographic direction, i.e., they are similar to ordinary Yang-Mills instantons with the time direction replaced by the holographic coordinate. The instantons are introduced in the chirally broken phase, where the D8 and D8 branes are connected, and they carry topological baryon number due to the presence of a Chern-Simons term in the action. This is our holographic version of "nuclear matter". In the chirally restored phase, the flavor branes are disconnected and baryon number is carried by quarks, represented by strings attached with one end at one of the N f flavor branes (either D8 or D8, depending on the chirality of the quarks) and the other end at one of the N c D4 branes that are responsible for the curved background geometry. This is our holographic version of "quark matter". One of our main results is the observation that the two embeddings of the flavor branes (connected and disconnected) can be continuously transformed into each other in the presence of instantons.
Throughout the paper we shall work in the chiral limit. Since the D8 and D8 branes intercept the D4 branes, quarks are massless in the chirally restored phase. For discussions about nonzero quark masses in the Sakai-Sugimoto model see Refs. [48][49][50][51][52]. In the massless limit, chiral symmetry is exact, which suggests a true phase transition between (chirally broken) nuclear matter and (chirally restored) quark matter, either of first or second order. It thus appears contradictory to address the question of the continuity between nuclear and quark matter in the chiral limit. However, our present goal is more modest: we seek to connect chirally broken and chirally restored phases in the "order parameter space", along stationary, but not necessarily stable, points of the free energy. This is one step towards an effective potential defined in the entire "order parameter space", which depends on chemical potential and temperature and where nuclear and quark matter correspond to local minima that can be continuously connected whenever they exist. Whether the model allows for a quark-hadron continuity between stable phases can be addressed more comprehensively only after including quark masses.
The main difficulty of our approach is the treatment of the baryonic phase. Single baryons in the Sakai-Sugimoto model have been studied extensively in the literature, from analytical flat-space approximations to purely numerical studies [38,39,47,[53][54][55]. Here we are interested in a many-baryon system, which requires various approximations. We follow the approach developed in Refs. [56][57][58][59], which is based on the pioneering work [60], where pointlike instantons were considered. Our work is a direct improvement of Ref. [58], whose approach we extend by accounting for the interactions of the instantons in the bulk. This is done by using the exact two-instanton solution in flat space, which is a special case of the general Atiyah-Drinfeld-Hitchin-Manin (ADHM) construction [61], and which has been discussed previously in the context of the Sakai-Sugimoto model to study the nucleonnucleon interaction [62,63]. As we shall see, the instanton interactions are crucial for our main observation. There are various other, in some aspects complementary, approximations to many-baryon phases in the Sakai-Sugimoto model, see for instance Refs. [64][65][66][67][68][69][70]. One of them is based on a homogeneous ansatz for the gauge fields in the bulk [57,71,72], which is expected to yield a better approximation for large densities, but which is less transparent from a physical point of view because it is not built from single instantons.
Our paper is organized as follows. The main point of Sec. 2 is to introduce our ansatz for the non-abelian gauge fields on the flavor branes, based on the single instanton and twoinstanton solutions in flat space. This is discussed in Sec. 2.2 with the technical background deferred to appendix A. In Secs. 2.3 and 2.4 we solve the equations of motion and set up the calculation of the stationary points of the free energy. Sec. 3 contains our main results. In particular, we discuss the continuity between the chirally broken and chirally symmetric phases in Sec. 3.2, and we compute the speed of sound in Sec. 3.3. We give our conclusions in Sec. 4.

Action
The general setup of the model in the given limit has been used in various previous works. The study most relevant to the present paper is Ref. [58], whose notation we shall employ. Here we will not go into details regarding the construction of the model itself and refer the reader to the reviews and original works quoted in the introduction and, more specifically, to Refs. [57,58]. The most important difference to Ref. [58] is our ansatz for the non-abelian field strength, which we will explain in detail.
Our starting point is the action for the U(N f ) gauge fields on the D8 and D8 branes in the deconfined geometry. In our treatment of baryons we shall restrict ourselves to N f = 2, i.e., we shall not include hyperons, and thus the non-abelian field strengths will be approximated with the help of SU(2) instantons. For consistency we shall use N f = 2 also for the quark matter phase, although for our purpose it would be trivial to add a third massless flavor in the chirally symmetric phase. We shall briefly comment on the threeflavor case in Sec. 3.1. The action has a Dirac-Born-Infeld (DBI) and a Chern-Simons (CS) contribution, with the 't Hooft coupling λ. For convenience, we have started with a mostly dimensionless formulation, the relation to the corresponding dimensionful quantities can be found in Refs. [57,58]. The only dimensionful quantities in the present formulation are in the prefactor of the action, the temperature T and the Kaluza-Klein mass M KK , which is the inverse radius of the compactified extra dimension with coordinate x 4 . In the dimensionless coordinates used here, x 4 ≡ x 4 + 2π. The integration over imaginary time has already been performed, yielding the prefactor 1/T . The remaining integrals are taken over position space R 3 and the holographic coordinate u ∈ [u c , ∞], where u c is the location of the tip of the connected flavor branes 1 and the holographic boundary is located at u = ∞. The embedding of the flavor 1 The notation "uc" originates from Ref. [60], where pointlike instantons were used as an approximation, which cause a cusp at the tip of the connected flavor branes, hence the subscript c. We keep using this notation even though the embedding is smooth in the presence of instantons with nonzero width. branes, given by the function x 4 (u) and by u c itself will later be determined dynamically, with the boundary condition x 4 (∞)−x 4 (u c ) = /2, where is the dimensionless asymptotic separation of the D8 and D8 branes, its dimensionful version given by L = /M KK . We have denoted derivatives with respect to u by a prime. The DBI action has been written in terms of the symmetrized trace "str", to be taken over the 2-dimensional flavor space. For simplicity, we shall later work with the ordinary trace. This is a tremendous simplification, and it has been shown for a very similar calculation that the use of the symmetrized trace prescription does not make a large quantitative difference in the results [58]. The expression in the square root is obtained from det(g + 2πα F), where g is the induced metric on the flavor branes, F is the field strength tensor, and α = 2 s with the string length s (which is absorbed in the definition of the dimensionless field strengths). For the explicit form of the metric see for instance Ref. [58]. It contains the temperature-dependent function with the dimensionless temperature t = T /M KK . The field strengths are decomposed according to U(2) ∼ = U(1) × SU (2). Within our ansatz, the only nonzero abelian field strength isF u0 = iâ 0 with the temporal component of the abelian part of the gauge field a 0 (u) (the factor i is needed in the imaginary time formalism). The quark chemical potential enters the action as a boundary condition for this abelian gauge field,â 0 (∞) = µ. The only nonzero non-abelian field strengths are F ij and F iu , i = 1, 2, 3. They are responsible for the nonzero baryon number through the topological charge provided by the CS contribution, and their form plays the main role in our calculation. We do not attempt to solve the full equations of motion for all gauge fields plus the embedding function in x and u. (In the action (2.2) we have already neglected the spatial derivatives ofâ 0 and x 4 .) We shall rather insert an ansatz for the square of the non-abelian field strengths, which depends on x and u, and then average over position space such that the spatial integration in the action becomes trivial. Then, we will proceed without further approximation and solve the equations of motion forâ 0 (u) and x 4 (u) fully dynamically.
We will exclusively work with the deconfined geometry, whose metric has been used to derive the DBI action (2.2a). Since we are interested in the chiral phase transition, we need to work in a regime where the critical temperature of the chiral phase transition is larger than that of the deconfinement transition. Since T deconf c = M KK /(2π) (for all µ unless backreactions are taken into account) and T chiral c 0.15/L (at µ = 0) [73], this requires /π 0.31 ( = π is the antipodal limit used in the original works by Sakai and Sugimoto). We shall also assume that the deconfined geometry is preferred down to arbitrarily small temperatures, which is justified in the decompactified limit, where the radius of the x 4 circle is taken to infinity (at fixed L) such that M KK and thus T deconf c goes to zero. In this limit, the chiral phase transition is entirely determined by the dynamics on the flavor branes, not by the background geometry.

Ansatz for non-abelian field strengths
With our ansatz for the non-abelian field strengths we attempt to capture as many effects as possible from a many-instanton system and at the same time try to keep the expressions as simple as possible to make feasible the full evaluation of the equations of motion for the abelian gauge field and the embedding function of the flavor branes. We will construct our ansatz from the 1-instanton and the 2-instanton solutions in flat space, which we now discuss.

Single instanton
For large, but not infinite, 't Hooft coupling, the leading-order expression for a single instanton has been discussed for the confined geometry in Ref. [47] for maximally separated flavor branes (= on opposite ends of the x 4 circle, = π) and for non-maximal separation in Ref. [53]. Here we need the analogous result for the deconfined geometry from Ref. [58], which reads where the superscript "(1)" is added to indicate the single-instanton solution, where x 2 = x 2 1 + x 2 2 + x 2 3 , and where z ∈ [−∞, ∞] is the holographic coordinate along the connected flavor branes, defined by such that z = 0 corresponds to the tip of the connected branes, u = u c . The shape of the instanton is given by where, as Eq. (2.5) shows, ρ can be interpreted as the width in the holographic direction, while ρ/γ corresponds to the spatial width. Therefore, γ is a "deformation parameter", which characterizes the deviation of the instanton from SO(4) symmetry. This deviation has also been observed in a purely numerical evaluation of a single instanton [55]. In the large-λ expansion, ρ 0 = 6 2π/ √ 5 10.06, γ 0 = 2 2/3 1.633 [58], and we have abbreviated The solution (2.5) is the flat-space Belavin-Polyakov-Schwarz-Tyupkin (BPST) instanton [74], which is exact for a single instanton in the limit of large λ because the width goes to zero for infinite λ and thus the effect of curved space is negligible. For our many-baryon system this is no longer true, and in principle we would have to construct a (numerical) solution in curved space. For simplicity, however, we shall keep the form (2.5) as an ansatz. In particular, we shall keep the dependence of ρ and γ on u c from Eq. (2.7), but treat ρ 0 and γ 0 as free parameters, which can be adjusted to capture effects that are known from the numerical solution of a single instanton in curved space and can also be adjusted to reproduce known properties of realistic nuclear matter 2 . Since u c will be determined dynamically for each µ and T , the instanton width in the spatial and holographic directions will depend on density and temperature. Self-duality of the instanton solution implies Note that all field strengths squared are proportional to the 2×2 unit matrix in flavor space. At this point, γ appears as a mere rescaling of the holographic coordinate z. However, when we insert our ansatz into the DBI action γ becomes a nontrivial parameter that cannot be eliminated by rescaling z. The single-instanton solution has topological charge 1, as it should be, where, for later convenience, we have denoted 3 (2.11)

Two instantons
The second ingredient for our ansatz is the flat-space two-instanton solution. If the separation of the two instantons is sufficiently large, this solution is well approximated by two single instantons. Let us for now consider SO(4) symmetric instantons, denote the width of both instantons by ρ, and let us place one instanton at the origin of the four dimensional ( x, z) space and the other one separated by a distance δ in the x 1 direction. Then, the exact flat-space solution that we will use is parameterized by ρ and δ and reads The parameters ρ and δ lose their interpretation of width and separation as δ is decreased and the instantons start to overlap. The derivation of Eq. (2.12) is explained in appendix A. It is based on the general ADHM construction [61] and its application to the Sakai-Sugimoto model [62,63]. For simplicity, we have assumed the instantons to have the same orientation, termed "defensive skyrmions" in Ref. [62]. 2 This approach was already suggested and explored in Ref. [58]. However, in this reference, ρ ∝ uc, γ ∝ u 3/2 c was used, i.e., the dependence of ρ on uc was chosen differently for simplicity. Here we work with the "correct" dependence suggested by the lowest-order, single-instanton solution, which leads to a somewhat more difficult calculation, but which turns out to be crucial (together with taking into account instanton interactions) for our main observations. 3 In Ref. [58], q0 was denoted by D.
This renders the subsequent calculation much simpler, but we should keep in mind that in principle the relative orientation of the instantons in the many-instanton system has to be determined by minimizing the free energy. In this sense, the free energy we calculate has to be understood as an upper limit, likely to be reduced by choosing a nontrivial relative orientation. We introduce the deformation parameter γ by the same rescaling as in the single-instanton solution (2.5): we replace ρ → ρ/γ, z → z/γ and divide F iz F (2) iz by γ 2 . Then, we can apply the self-duality (2.9) for the two-instanton solution, and we verify that the topological charge is 2, Here we have introduced the function with the polynomials The particular way (2.14) of writing the function q int is not only motivated by compactness, but also turns out to be useful for integrating and differentiating q int . We have also introduced the overlap parameter which is the distance between the instantons normalized by twice their spatial width, such that, if they were spheres with radius ρ/γ, they would overlap for d < 1.

Many-instanton approximation
The ADHM construction in principle provides exact many-instanton solutions, which however become increasingly complicated and are much too unwieldy for our purpose. Elements of these general solutions have been used in the Sakai-Sugimoto model in Refs. [75,76], see also a brief discussion in Ref. [70]. Here we construct a simple approximation for a manyinstanton system, taking into account only 2-body interactions. To this end, we first define an interaction energy density by subtracting the two individual single-instanton contributions from the full two-instanton solution, Via the self-duality (2.9), analogous relations hold for the other two relevant quadratic forms F ij F ij and F ij F kz ijk . The two instantons are centered at the points ( x 1 , z 1 ) and ( x 2 , z 2 ). Each term in Eq. (2.18) also depends on x and z, which we have not indicated explicitly in the arguments. The field strength squared for a system of N I instantons is then approximately constructed as Let us first comment on the holographic location of the instantons. Previous studies have shown that at low baryon density the instantons are located at z = 0, the tip of the connected flavor branes, while at large density additional instanton layers at |z| > 0 are expected to appear [58,67,72,77]. Allowing for multiple layers makes our calculation somewhat more complicated, both the numerical evaluation and the following equations. Therefore, and since our emphasis is on the new interaction terms, we discuss all modifications that arise due to multiple layers in appendix B. With the help of this appendix, we have checked that allowing for more layers does not change our main conclusion regarding the continuous connection between nuclear and quark matter. Here we continue with a single layer, i.e., we center all instantons around the point z = 0 in the holographic direction. Consequently, z n = z m = 0 in Eq. (2.19), and we shall omit the arguments referring to the holographic coordinate from now on. Next, we put the instantons on a regular lattice in position space and apply the nearest-neighbor approximation. We denote the number of nearest neighbors by p and the lattice vectors that connect each instanton with its nearest neighbors by δ m , with the nearest neighbor distance δ, i.e., | δ 1 | = . . . = | δ p | ≡ δ. Then, Eq.
In the second line we have employed a further approximation by averaging the field strengths squared over position space, denoted by angular brackets, and used Eq. (2.18). This has allowed us to shift every single term of the sum to the origin, x n → 0, with δ connecting an instanton in the origin to one if its nearest neighbors, say along the x 1 direction. We have used the spatial integrals for the one-instanton and two-instanton solutions from Eqs. (2.10) and (2.13), respectively. Moreover, we have introduced the dimensionless instanton density (per flavor) where v is the dimensionless three-volume of the system, related to the dimensionful volume is of course related to the instanton density. This relation depends on the lattice structure. We write the volume of a Wigner-Seitz cell as v such that r is the inverse Wigner-Seitz volume in units of the nearest-neighbor distance δ. Consequently, with Eq. (2.17) we can express d as The parameters p and r carry all information about the lattice structure in our approximation. For a cubic, a body centered cubic, and a face centered cubic crystal we have (p, r) = (6, 1), (8, 3 √ 3/4), (12, √ 2), respectively. By setting the number of nearest neighbors p to zero the result becomes independent of r and we recover the non-interacting approximation of Ref. [58]. Replacing the field strengths squared by their spatial average is an enormous simplification since now the equations of motion will be ordinary differential equations, only depending on the holographic coordinate u. A more refined approximation would be to average not over the squared field strengths, but over the entire square root of the DBI action (2.2a) before solving the equations of motion. Although still resulting in ordinary differential equations, this procedure would be much more complicated since there is no simple analytic form for the spatial integral over the DBI Lagrangian, even if the above analytic ansatz for the non-abelian field strengths is used.

Solving the equations of motion
We now insert our ansatz (2.20) together with the corresponding expressions for F 2 ij , F ij F kz ijk into the action (2.2). As already mentioned, instead of the symmetrized trace we take the standard trace for simplicity. Since all field strengths squared are diagonal in flavor space, this simply gives an overall factor N f , and we can write with the Lagrangian where (Here we are switching back and forth between the coordinates u and z, whose relation is given in Eq. (2.6).) We have brought the Lagrangian into the same form as in Ref. [58], with all effects of the instanton interaction absorbed in the function q(u). Therefore, for the solutions of the equations of motionâ 0 (u) and x 4 (u), we can simply follow this reference.
With the integration constant k we obtain where we have abbreviated and with the polynomials and S 1 , S 2 , a, b defined in Eq. (2.16). At first sight, finding an analytic expression for the integral in Eq. (2.31b) seems hopeless, due to the complicated form of q int . As explained in appendix C, however, an analytical form can be found from the observation that taking successive derivatives does not change the structure of the function, only the exponents in the denominator and the explicit form of the polynomials that multiply S 1 and S 2 . The solution (2.28) includes the case of a trivial embedding function, x 4 = k = 0, which corresponds to the chirally symmetric quark matter phase. It also includes the mesonic phase, which has vanishing baryon number, n I = 0, and thus a constant fieldâ 0 (u) = µ. In the baryonic phase both k and n I are nonvanishing. For all cases, the asymptotic behavior at the holographic boundary u → ∞ can be written aŝ This expansion follows from Eq. (2.28) and Q = 1 + O(u −6 ), ζ = 1 + O(u −5 ) and confirms the interpretation of n I as the baryon number density.

Stationary points of the free energy
According to the on-shell action (2.24) the free energy is N N f Ω with Here, the Lagrangian is obtained by re-inserting the solutions (2.28) into Eq. (2.25). As written, the free energy is divergent. We can easily deal with this divergence by introducing a cutoff for the upper boundary of the u integral. This introduces a cutoff dependent term that turns out to be independent of chemical potential and temperature. We can therefore simply subtract this vacuum term since it does not affect the physics. The free energy depends on various parameters. Firstly, there are the model parameters λ, M KK , and . Secondly, there are the quantities k, n I , u c , ρ, γ, which, in principle, have to be determined dynamically. However, as shown in Ref. [58], this leads to pointlike instantons at the baryon onset (which, unphysically, is of second order), while we know that at finite λ the instantons should have a nonzero width. We have checked that this observation remains unchanged in the present approach and thus we shall use the form of ρ and γ from Eq. (2.7) and treat ρ 0 and γ 0 as parameters, increasing the number of model parameters from 3 to 5 [58]. This can be understood as an attempt to capture nontrivial effects, for instance known from other approximations and numerical studies, for which our original ansatz is too simplistic. It can also be understood as a phenomenological approach, since it gives us additional freedom to reproduce properties of real-world nuclear matter. Within this approach, it remains to find the stationary points of Ω with respect to k, n I , and u c . The three stationarity equations can be written more explicitly as where the derivatives with respect to d are taken at fixed ρ and vice versa. The first two equations, the derivatives with respect to k and n I , are straightforwardly derived. The derivative with respect to u c is more complicated, and we explain the derivation in appendix D. We have denoted the dimensionless entropy density by s, Note that s is computed solely from the derivative with respect to the explicit dependence on the dimensionless temperature t because the implicit dependence through k, n I , u c vanishes at the stationary point. We see from Eq. (2.36a) that stationarity with respect to k is equivalent to the boundary condition x 4 (∞) − x 4 (u c ) = /2, given by the fixed asymptotic separation of the flavor branes. Therefore, Eq. (2.36a) should be viewed as a constraint, it is not a minimization of the free energy. This constraint restricts the parameter space to a two-dimensional surface in (k, n I , u c ) space. Local minima of the potential are points that fulfill Eqs. (2.36) and in whose vicinity the potential increases in all directions on this two-dimensional surface (and not necessarily in all directions in the three-dimensional (k, n I , u c ) space).

Choice of parameters
The stationarity equations (2.36) have to be solved simultaneously for k, n I , u c for given chemical potential µ and temperature t. In the remainder of the paper we shall set t = 0. Moreover, we shall work with the cubic lattice, p = 6, r = 1. We have checked that other lattice structures slightly change our results quantitatively, but not qualitatively. It remains to fix the model parameters λ, M KK , , ρ 0 , γ 0 . Here we do not attempt to study this parameter space systematically, as it was done in Ref. [58] within a simpler version of the model.
We shall rather work with a single parameter set, which we determine by requiring the model to reproduce the basic properties of nuclear matter at saturation (a first step in this direction was made in Ref. [59]). This is, firstly, a useful preparation for future extensions and applications of this model. For instance, if the equation of state is calculated and used to describe dense matter inside neutron stars, it is important to anchor the approach to known low-density properties. And, secondly, the fact that it is possible to reproduce properties of real-world nuclear matter can be viewed as a phenomenological validation of the model. We will reproduce 4 properties with 5 free parameters, and thus the fact that the matching works might not seem very impressive. But, due to the complicated, highly nonlinear structure of the stationarity equations it is a priori not obvious that the equations set by the constraints do have a solution.
The physical constraints we impose are as follows. We require the vacuum mass of the nucleon to be m N = 939 MeV. The remaining conditions are properties of infinite, isospin-symmetric, zero-temperature nuclear matter at saturation: the binding energy E B = −16.3 MeV, the saturation density n 0 = 0.153 fm −3 , and the incompressibility, for which only a certain range is known, K (200 − 300) MeV. Relating these quantities to our model gives the following four conditions, The right-hand sides depend on the parameters of the model, through the dimensionless quantities µ 0 , m 0 , n 0 I , κ, which are functions of ρ 0 , γ 0 , λ, and . Here, µ 0 is the quark chemical potential at the baryon onset, and n 0 I is the baryon density just above the onset. To relate µ 0 and n 0 I to their dimensionful counterparts, we have used Eqs. (2.44) and (2.45) of Ref. [58]. Moreover, m 0 is the constituent quark mass in the vacuum, such that N c m 0 is the vacuum mass of the baryon. In our approach, it can be computed from the chemical potential in the limit of vanishing baryon density, n I → 0. Since in this limit the instantons are infinitely far away from each other we can obviously use the non-interacting limit to compute the vacuum mass. We find that the solution of Eqs. (2.36) for n I → 0 with µ approaching a nonzero value is (for zero temperature) We mention already here that there is a second solution for n I → 0, where µ does go to zero as well. This solution will play crucial role later, see Eq. (3.1). Finally, κ is the dimensionless version of the incompressibility at saturation, This derivative with respect to the baryon density is taken at fixed temperature, but with µ, k, u c being functions of n I . In contrast, the same derivative, used above for minimizing the free energy, was taken at fixed temperature and fixed µ, k, u c . We compute the derivative in (2.40) purely numerically.
The 4 conditions (2.38) are used to fix ρ 0 , γ 0 , λ/ , and /M KK = L. This matching is not unique because K is not known exactly and because it is conceivable that different, disconnected regions in parameter space satisfy the physical constraints. It is beyond the scope of this paper to present a systematic analysis of the multi-dimensional parameter space. The parameter set that will be used throughout the rest of the paper is It can be checked numerically that these values satisfy Eqs. It is then easy to check that these values satisfy the physical conditions at saturation; in particular K 220 MeV, which is within the experimentally known range [78,79], possibly somewhat smaller [80]. The predicted value for L from Eq. (2.41) can be used to compute the critical temperature for the chiral transition at zero chemical potential. Equating the free energies of the mesonic phase and the quark matter phase -for instance using appendix B of Ref. [57] -yields a first-order phase transition with T chiral c 0.15384/L, from which we obtain T chiral c 156 MeV. This is very close to the transition temperature for the chiral crossover in real-world QCD. We should keep in mind, however, that we work in the chiral limit with two quark flavors. Since lattice simulations are not available in this limit, the exact value and even the order of the transition is not rigorously known. It is nevertheless worth pointing out that from fixing our parameters to cold nuclear matter, we obtain a prediction for the finite-temperature phase transition which can in principle be compared to QCD.
Besides the absolute value of L let us also briefly comment on the relation between and M KK in Eq. (2.41). This relation suggests that by choosing M KK arbitrarily any value of can be obtained. However, is constrained because the asymptotic separation of the flavor branes can obviously not be larger than half of the circumference of the x 4 circle, and hence < π is an absolute limit. As discussed at the end of Sec. 2.1, a more restrictive condition comes from the observation that there is a chirally broken phase in the deconfined geometry only if /π < 0.30768. Interestingly, at this upper limit, we have λ 14.6 and M KK 980 MeV, which is close to the original fits in the confined geometry, which were obtained using completely different constraints based on the rho meson mass and the pion decay constant [38]. The decompactified limit, which allows us to apply the deconfined geometry to arbitrarily small temperatures, implies π. It is a matter of interpretation whether this stronger constraint should be obeyed since one might interpret the decompactified limit as a separate model, allowing for extrapolations of /π beyond the originally allowed regime. Nevertheless, it is reassuring that π does not contradict our physical constraints.

Solution of the stationarity equations
Having fixed all parameters, we can solve the stationarity equations and compute the free energy. This has to be done numerically in general, but in one instance we have found an analytical expansion. This is the limit of small µ, where the solution to Eqs. (2.36) is where we have abbreviated With Eq. (2.17), this solution yields in particular d ∝ µ −1/2 , i.e., the ratio of the distance over the width of the instantons becomes infinitely large as µ goes to zero. As a consequence, the instantons are essentially non-interacting in this limit, and Eq. (3.1) does not depend on the interaction terms 4 . We show the full numerical solution in Fig. 1. In this figure, we have plotted n I (upper panels) and u c (lower panels) 5 . The third variable k is not plotted because it shows a qualitatively identical behavior as u c , although there is no simple analytical relation between them. Each quantity is shown in a linear plot (left) and a double-logarithmic plot (right) to make all important features of the solution visible. In the left plots, solid lines corresponds to the stable solution, dashed lines indicate metastable and unstable branches of the solution. Fig. 1 also includes the result from the mesonic phase, which has zero baryon number, n I = 0 and where u c is constant in µ and assumes the value (2.39a). We see that there is a first-order baryon onset, which occurs by construction due to our matching of the parameters to real-world nuclear matter at saturation. The double-logarithmic plots show that the two branches in the left plots are continuously connected. As a check, we have also included the analytical result (3.1) to confirm that it is in exact agreement with the full result at (very) small µ. Moreover, we see that at small, but not too small, chemical potentials, n I and u c become three-valued in the baryonic phase. This multivaluedness is not visible in the linear plots. We do not have any particular interpretation for this behavior. It does not play any role for the ground state since in this regime the mesonic phase has a smaller free energy than any of the three baryonic solutions. We shall thus not further discuss this feature, although it is prominently visible in the double-logarithmic plots.
In the left panel of Fig. 2 we show the pressure P = −Ω for the mesonic phase (straight red line) and the baryonic phase (solid black curve), both divided by the pressure of the 4 The solution (3.1) does not exist if ρ ∝ uc instead of ρ ∝ u 3/4 c , which explains why this solution was not found in Ref. [58]. 5 Here and in all following figures, the axes labels have to be understood to include appropriate powers of , i.e., nI stands for nI 5 , µ for µ 2 , uc for uc 2 , and Ω for Ω 7 ; the overlap parameter d does not scale with .
where the subscripts indicate the shape of the flavor branes of the two phases. The favored state for a given µ is the one with the largest value of P/P || . We see that for small µ the mesonic phase (the vacuum) is preferred, until at µ 0 0.253/ 2 there is a first-order phase transition to the baryonic phase. This first-order phase transition is weak in the sense that the chemical potential at the transition is not much smaller than the vacuum mass of the nucleon, as required from the matching to QCD. Therefore, in the double-logarithmic plot, the multivaluedness of the baryonic pressure around this first-order phase transition is not visible. At µ 42.7/ 2 there is a first-order phase transition to the chirally symmetric phase.
In physical units, this corresponds to a critical chemical potential of about 160 GeV or, using the value for the dimensionless baryon density at the transition n I 3780/ 5 , to more than 10 5 times saturation density. This is an enormously large value, for instance compared to the chiral transition suggested by phenomenological models, but also because at these ultra-high densities we are approaching the regime where perturbative QCD is expected to be valid. If two-layer or multi-layer instanton solutions are taken into account, which decrease the free energy of the baryonic phase [58], the chiral phase transition will be even further shifted to larger densities. On the other hand, it has been shown, at least without instanton interactions [58], that the chiral phase transition can be arbitrarily close to the baryon onset for certain choices of the parameters γ 0 and ρ 0 , which then, however, lead to unrealistically large binding energies of nuclear matter [59]. In any case, we should take the large value of the critical chemical potential with care for various reasons. Firstly, we have neglected backreactions on the background geometry, which can be expected to become large at large baryon densities. Even if our approximation is taken seriously, the large value may be a consequence of the large-N c limit that is inherent in our calculation and/or the lack of asymptotic freedom of the model. We should also recall that for consistency we work with two flavors also in the chirally symmetric phase. In reality, however, strange quark matter is the ground state at high densities. (Even more quark flavors appear if the chemical potential is sufficiently large, but, eventually having in mind neutron star applications, let us ignore them.) In our approximation of massless quarks, an additional quark flavor in the chirally symmetric phase simply changes the prefactor of the free energy from N f = 2 to N f = 3. Since the free energy is negative, this reduces the free energy, i.e., makes strange quark matter more favorable, as it should be, and we find a critical chemical potential of about 30 GeV, still large, but reduced by about a factor 5 compared to two-flavor quark matter. Since here we are mostly interested in qualitative results, we leave a more systematic discussion of the chiral phase transition for the future.
In the right panel of Fig. 2 we show the instanton overlap parameter d. The (black) solid curve, bounded by the two dots, is the stable baryonic segment, from the baryon onset (large d) to the chiral transition (smaller d), while the dashed segments are unstable or metastable. The result shows that the instantons "refuse" to overlap strongly. Even at the highest densities we have d > 1, because, as the instantons are squeezed, they become smaller. As a consequence, the non-interacting solution is a very good approximation for large parts of the curve. This solution, obtained in our approximation by setting the number of nearest neighbors to zero, p = 0, is shown by the two thin (red) curves. The right branch of the non-interacting result ends at large µ, we have not found any solutions beyond this point. Since the numerics get increasingly difficult at these large chemical potentials we cannot exclude that this branch continues beyond the point where we have stopped. On the other hand, we were able to continue the left branch to much smaller values of µ than shown in the plot. Therefore, if instanton interactions are neglected these two branches appear to be disconnected. Only in the presence of interactions all baryonic solutions we have found are continuously connected. And, they are also continuously connected to the chirally symmetric phase, as we shall discuss now.

Connecting nuclear matter with quark matter
Our results show a continuity between nuclear and quark matter in the following sense: the curve [k(µ), u c (µ), n I (µ)] describing the stationary points of the free energy is continuous and is a path that can be traced from the baryonic phase into the quark matter phase (along which µ is non-monotonic). This continuous path manifests itself also in the free energy, which is shown schematically in Fig. 3. We have chosen a schematic representation for convenience, which is a distorted, but not disrupted, version of the path obtained from the actual calculation. We have also indicated the geometries of the flavor branes in the two-dimensional (u, x 4 ) subspace of the model, together with a "legend" that explains these geometries. If we follow the state with the lowest free energy, the transition from nuclear to quark matter is not continuous. For instance n I , which is the (negative) derivative of Ω with respect to µ, as well as u c are discontinuous. Therefore, the continuous path just described passes through metastable and unstable states. At zero chemical potential, the baryonic phase connects continuously to the chirally symmetric phase. With the help of the analytical solutions (3.1) this statement is exact and does not rely on a numerical calculation: the small-µ expansions show for instance that u c indeed vanishes as µ goes to zero. We can investigate this small-µ regime in more detail for instance by studying the behavior of the instantons. In Fig. 4 we show the instanton profiles at certain locations of the continuous path.
This figure shows the instanton profiles via cuts in the (z, x 1 ) and (x 1 , x 2 ) subspaces in the first two columns. The profiles are given by the field strengths squared from the first line of Eq. (2.20), and in each row the color scale is adjusted to the maximal value of the profile. For simplicity, we have calculated the profiles from the non-interacting limit, which makes almost no difference since, as we have seen and as this figure further illustrates, the instantons never overlap significantly. The figure also shows the corresponding embedding of the flavor branes (third column) and indicates the value of µ and d for each row (fourth column). The first three rows correspond to stable solutions, the first row being located at the baryon onset and the third row at the chiral phase transition. As we move towards the chirally symmetric phase the instantons become infinitesimally thin in the holographic direction, while they spread out in the spatial direction and become infinitely wide 6 . This is a consequence of the scaling of the holographic width ρ ∝ u 3/4 c and the spatial width ρ/γ ∝ u −3/4 c . Since u c goes to zero according to Eq. (3.1) we have ρ ∝ µ 3/2 and ρ/γ ∝ µ −3/2 . Note that nonzero-temperature effects will slightly change this picture. The geometry dictates that u c ≥ u T and thus u c > 0 for T > 0. Hence the instantons will always retain a nonzero width in the holographic direction and a finite spatial width within the present 6 From a top-down string-theoretical point of view one might argue that the resulting large derivatives of the field strengths require derivative corrections of higher order in α to the DBI action. Here we do not include such terms for simplicity.  ansatz. We leave a more thorough study of nonzero temperatures to the future.
Having emphasized the continuous connection between nuclear and quark matter, let us discuss and interpret this result. It is striking that the connection relies on the nuclear and quark matter branches to "meet" at zero chemical potential. The reason why the quark matter solution exists all the way down to µ = 0 is that we work in the chiral limit and thus quarks can be placed into the system at infinitesimally small chemical potential. Of course, at small chemical potential the quark matter phase is energetically disfavored compared to the mesonic phase (in real-world QCD one would expect an infinite energy cost due to confinement). Nucleons, on the other hand, do have a mass and are placed into the system at a nonzero chemical potential. Therefore, the uppermost branch in Fig.  3 should not be confused with ordinary nuclear matter where the density is reduced as we approach the origin. Ordinary low-density nuclear matter sits on the lower branch and is superseded by the mesonic phase as the chemical potential is decreased. We thus emphasize that we are not simply decreasing the density of two phases until all matter is gone and then claim a continuous connection between them. We rather observe that all branches, from the mesonic phase through the nuclear matter phase up to the quark matter phase are continuously connected if we include all solutions of the equations of motion and all stationary points of the free energy, which are not necessarily stable. In fact, the uppermost branch in Fig. 3 can be expected to be unstable because it has the largest free energy of all solutions we have found. It should thus be a local maximum, not a local minimum of the free energy. By plotting the free energy for a fixed µ as a function of n I and u c , with k determined such that the constraint (2.36a) is fulfilled, we have confirmed this expectation. Despite being unstable, the existence of this branch is interesting for the following reason. Ultimately, we are interested in a single continuous potential, defined in the entire "order parameter space". If there is a first-order phase transition between nuclear and quark matter, one should be able to use this potential to connect the two phases continuously. The knowledge of the full potential is for instance needed for the calculation of the surface tension of an interface between nuclear and quark matter [23,25]. Our approximation does not yield such a potential. It is easy to check that there are regions in parameter space where the potential becomes complex. However, compared to previous approximations within the same model (and compared to the vast majority of field-theoretical models) we have made a step forward by at least finding a path in parameter space along stationary points of the potential that connects nuclear and quark matter.
Since we are working with massless quarks, the transition from the chirally broken to the chirally restored branch at zero chemical potential can be expected to involve a discontinuity in second derivatives of the free energy, just like an ordinary second-order phase transition. We shall see in the next section that this is indeed the case, by computing the speed of sound. In other words, since chiral symmetry is an exact symmetry, instantons are topologically protected and a transition from a state with nonzero instanton number to a state without instantons must involve some kind of discontinuity. These arguments suggest an obvious improvement for the future, namely to include nonzero quark masses, which breaks chiral symmetry explicitly. In this case, firstly, the quark and nuclear branches may connect without reaching back to zero chemical potential (since massive quarks require a finite energy and thus the quark matter branch cannot start at zero chemical potential). And, secondly, the transition between the chirally broken and (approximately) symmetric phases is allowed to be smooth, including the second derivatives of the free energy. It is then even conceivable that the multivalued curve in Fig. 3 "straightens out" at large densities, such that an actual quark-hadron continuity occurs in the phase diagram. Therefore, and despite the various approximations we have made and despite the various differences to real-world QCD, we believe that our approach can be very useful in understanding the transition between nuclear and quark matter from a physical point of view.

Speed of sound
The general form of the speed of sound is where the derivative of the pressure P with respect to the energy density = −P + µn + sT is taken at fixed entropy per particle s/n, and where the derivatives of entropy density s and baryon number density n with respect to the baryon chemical potential µ are taken at fixed temperature T and vice versa. We derive Eq. (3.4) in appendix E, see also Refs. [81,82] for similar derivations in the context of superfluids and appendix A of Ref. [83]. As above, we restrict ourselves to zero temperature. In this case the speed of sound in the baryonic phase and the chirally restored phase becomes which is a useful relation to keep in mind for the following results. Since a phase can only be thermodynamically stable if its density increases monotonically with the corresponding chemical potential, this relation also shows that a negative speed of sound squared indicates an instability (assuming that n and µ themselves are both positive). We compute the speed of sound numerically for the baryonic solution discussed in the previous subsection, i.e., with the parameters from Eq. (2.41). The results for the mesonic phase and the quark matter phase are derived in appendix F and read mesonic: chirally symmetric: The result for the mesonic phase away from the zero-temperature limit can only be computed numerically, but is not needed here. The general result for the chirally symmetric phase is written in terms of the baryon density n I , which depends on µ and T in a nonanalytical way. (Although there are no instantons in the chirally symmetric phase we have kept denoting the dimensionless baryon density by n I for consistency.) In Fig. 5 we show the speed of sound squared as a function of chemical potential in the baryonic phase and, in the regimes where they are favored, in the mesonic and quark matter phases. Close to the baryon onset, barely visible on the plot, c 2 s is negative. The solution is unstable for densities smaller than about 66% of the saturation density, then becomes metastable and finally stable at saturation. Besides this instance, there are three other turning points of the curve n I (µ) close to which its derivative is negative. In contrast to the turning point close to the baryon onset, in these three instances the curve n I (µ) goes through a point where its derivative is zero such that, according to Eq. (3.5), the speed of sound diverges. Apart from small vicinities of the turning points the speed of sound is real for the entire nuclear matter solution. We should keep in mind, however, that c 2 s > 0 is a necessary, but not sufficient, stability criterion. As we know from the previous subsection, certain segments of the thin (black) curve in Fig. 5 correspond to maxima of the free energy and thus are unstable.
At very small µ, i.e., in the regime of the continuous transition from the connected to the disconnected flavor branes, the speed of sound squared approaches 1/6. This can be seen from the analytical result (3.1), which shows n I ∝ µ 6 , from which c 2 s = 1/6 follows with the help of Eq. (3.5). This value is different from the speed of sound in the chirally restored phase, which is c 2 s = 2/5, as shown for large chemical potentials in the figure. Consequently, as we go from the chirally broken to the chirally restored phase at zero chemical potential, the speed of sound is discontinuous. As explained above, it is expected that second derivatives of the free energy show discontinuities due to the different symmetries of the phases. Since the speed of sound (3.5) contains a second derivative, its discontinuity at zero chemical potential is no surprise. Interestingly, 1/6 is the value of the speed of sound in the chirally restored phase if we set µ = 0 and then take the limit T → 0 (note the different values of c s at T = µ = 0 depending on the order in which the limits of the function in Eq. (3.6b) are taken).
Let us now comment on the stable (red) branch. At the baryon onset, the speed of sound jumps from its mesonic value down to c 2 s 0.025, then increases to a maximum of about c 2 s 0.47, then decreases and finally jumps down to the constant value of the quark matter phase at the chiral transition. In QCD, we know that for asymptotically large µ the speed of sound assumes the conformal value c 2 s = 1/3 and perturbative corrections decrease this value. We also know that nuclear matter at saturation has a much smaller -non-relativistic -speed of sound, close to the value we have found here. It is unknown how the low-density and high-density regimes are connected. It has been pointed out that the simplest scenario, connecting the two regimes with a monotonic function, creates tension with astrophysical observations [84]. More precisely, a sufficiently stiff equation of state (= a sufficiently large speed of sound) is required to obtain observed neutron star masses of twice the solar mass [85,86], possibly larger [87]. Hence, in the intermediate density regime the speed of sound most likely has to exceed the conformal value. There have been suggestions in the literature that the conformal limit may be an absolute upper bound, but counterexamples have been pointed out within the gauge/string duality [88,89]. Values larger than c 2 s > 1/3 are routinely reached in phenomenological models of nuclear matter or in extrapolations of low-density effective theories, but in these calculations the speed of sound typically behaves monotonically and cannot be traced into the chirally symmetric phase. Therefore, phenomenologically motivated interpolations between the low-density and high-density QCD results have been studied and the consequences for masses and radii of neutron stars have been tested [21]. It is intriguing that the nuclear matter of our holographic calculation does show a non-monotonic behavior, as required from putting together constraints from QCD and neutron stars. Of course, our result is different from QCD at asymptotically large µ. The reason is that the mass scale M KK plays a role even in this asymptotic regime, allowing the speed of sound to be different from 1/3. This might be interpreted as a non-perturbative strong-coupling effect, having in mind that we need to stop trusting the model for asymptotically large densities, where QCD exhibits asymptotic freedom. We also mention again that the chiral transition occurs at an extremely large density, a fact easy to be overlooked due to the use of the logarithmic scale in Fig. 5, and, as explained above, that one should treat the quantitative interpretation of the model with some care, especially at large densities. Nevertheless, our approach allows for a consistent microscopic calculation of nuclear and quark matter with a qualitative prediction for the speed of sound that is very interesting in view of recent astrophysical constraints.

Summary and outlook
We have shown that chirally broken and chirally symmetric phases in the holographic Sakai-Sugimoto model can be continuously connected in the presence of instantons on the flavor branes. Geometrically, this continuity is realized by transforming the curved, connected flavor branes for left-and right-handed fermions into straight, disconnected branes. The transformation appears dynamically by solving the equations of motion and computing the embedding of the flavor branes in the deconfined geometry of the model and in the presence of an interacting many-instanton system. We have approximated the instanton interactions with the help of the exact flat-space two-instanton solution. We have found that, as the connected flavor branes become straight and approach the disconnected embedding, the instantons become infinitesimally thin in the holographic direction, but spread out to become infinitely wide in the spatial direction. Nevertheless, they barely overlap since in this limit also their density goes to zero. Most previous related studies had included instantons only in the confined geometry or did not include interactions between the instantons, i.e., to the best of our knowledge this continuity has been observed for the first time in the present paper.
Our observation is closely related to the transition between nuclear matter (= instanton system with connected flavor branes) and chirally symmetric quark matter (= disconnected flavor branes) in dense QCD. While our continuity does not directly connect stable nuclear matter at the lowest densities with stable ultra-dense quark matter, it does so on a path that goes through metastable and unstable stationary points of the potential. In particular, the actual phase transition between nuclear and quark matter turns out to be of first order. Since we have worked in the chiral limit, it is obvious that the actual transition cannot be continuous, since chiral symmetry is exact. It would thus be very interesting to include quark masses into our approach and see whether and how the continuity is affected. This can be done with worldsheet instantons [52], which have been used to study the effect of quark masses on the hadron spectrum [51,[90][91][92], or with a tachyonic effective action [49,50].
Besides this main theoretical result we have also pointed out phenomenological applications of our approach. We have shown that it is possible to fit the parameters of the model to reproduce the basic properties of isospin-symmetric nuclear matter at saturation. In doing so, it was crucial to introduce two additional parameters that characterize the instanton shape, increasing the number of model parameters from 3 to 5. We have worked with a single parameter set, and future studies are required for a more systematic investigation of the parameter space. For instance, we have seen that for the chosen parameters the chiral transition occurs at an extremely large baryon density -probably unrealistically large, at least compared to predictions of other models. It would be interesting to see whether different parameter sets can be found that are also in agreement with nuclear physics and yield different values for the chiral transition. We have also calculated the speed of sound to further connect our results to the phenomenology of dense matter. Our result shows a non-monotonic behavior of the speed of sound in nuclear matter that reaches a maximum larger than that of the quark matter phase, before it jumps to the quark matter value at the chiral transition. Interestingly, this non-monotonic behavior is suggested by putting together constraints from QCD and astrophysical data, most notably the largest observed neutron star mass of about two solar masses. It is therefore a natural extension of our work to compute the equation of state and the resulting maximum mass of a star that can be reached with it.
Let us conclude with mentioning some further extensions that may improve our calculation. Obvious technical improvements concern our many-instanton approximation. We have constructed the interaction terms solely from 2-instanton contributions, and those, in turn, are based on the flat-space solutions. Moreover, we have simply averaged over position space before solving the equations of motion and the embedding of the flavor branes. Also, we have kept the background geometry fixed although the gauge fields on the flavor branes become very large at large densities and backreactions may be non-negligible. Improvements on all these points are difficult, but since we are working within a top-down approach it is in principle known how to proceed. More direct extensions are for instance the study of two-or multi-layer solutions, which we have briefly discussed but not included into our numerical evaluation; isospin-breaking terms, which are needed for a more realistic study of neutron star matter; and nonzero-temperature effects, which we have included in almost all our equations, but not in the numerical results.

Acknowledgments
We would like to thank N. Evans, A. Rebhan, and S. Reddy for valuable comments. A.S. is supported by the Science & Technology Facilities Council (STFC) in the form of an Ernest Rutherford Fellowship. KBF thanks the University of Southampton for their hospitality during his sabbatical.

A Instantons from the ADHM construction
In this appendix we derive the two-instanton solution from the ADHM construction [61]. This is in large parts a recapitulation of derivations from Refs. [62,63], adapted to our purposes. In particular, we derive Eq. (2.12) in the main part.
According to the ADHM construction, the gauge fields for an Sp(n) k-instanton solution, where Sp(n) is the unitary symplectic group, are where M = 1, 2, 3, z, and U is a (n + k) × n matrix with quaternionic entries, which satisfies and Here, ∆ is a quaternionic (n + k) × k matrix, such that the quaternionic k × k matrix is symmetric and real, i.e., it commutes with quaternions. The relevant case for us will be n = 1 since Sp(1) ∼ = SU (2). We write the quaternions as with M = 1, 2, 3, 4, q M ∈ R, and the basis vectors e i = iσ i for i = 1, 2, 3 and e 4 = 1, where σ i are the Pauli matrices, and the 4-direction corresponds to the holographic z coordinate. We denote q † = q MēM withē i = −iσ i for i = 1, 2, 3 andē 4 = 1. Below we shall need the scalar product and cross product of two quaternions, defined as where we have used {σ i , σ j } = 2δ ij , [σ i , σ j ] = 2i ijk σ k . Note in particular that the scalar product is proportional to the unit matrix, while the cross product is spanned by the other three basis vectors of the quaternionic space. It is important to keep in mind that the 2 × 2 space introduced by the representation of the quaternions occurs for every gauge group, it has nothing to do with the SU(2) gauge group relevant for our case. We will represent the four-dimensional vector x M = ( x, z) as a quaternion as The matrix ∆ can be written as ∆ = a + b ⊗ x with quaternionic (n + k) × k matrices a, b that are constant in x, and where in b ⊗ x each element of b is multiplied by x from the right. One can choose b † = (0 k×n , −1 k ), while a remains generic and encodes the widths and locations of the instantons (more precisely, a depends on parameters that correspond to widths and locations of the instantons if the instantons are sufficiently far away from each other). Thus, with a quaternionic k × k matrix X and a quaternionic n × k matrix Y we can write Using Eq. (A.1) and where η ijk = ijk , η ij4 = −η i4j = δ ij , η i44 = 0. Consequently, the ADHM construction ensures the self-duality of the field strengths, A useful relation for the trace over the field strengths is [95] Tr[F 2 M N ] = − 2 ln det L , (A.11) where = ∂ 2 M . Below we shall compute F 2 M N explicitly for the 2-instanton solution and thus we do not really need to employ this relation. Nevertheless, calculating the trace of the field strength squared via L can be used as a check for our explicit form.

A.1 Single instanton
As a warm up, it is useful to start with a single instanton. In this case, k = 1 and thus X and Y are quaternions. Anticipating the notation used for k = 2 below let use denote y = Y . We determine U by writing U † = (α † , β † ) with quaternions α, β (in general, α is a quaternionic n × n matrix and β is a quaternionic k × n matrix). Then, Eq. (A.3) reads After multiplying this equation by y from the left we obtain α, which we insert into Eq. (A.3) to obtain |β| 2 . Denotingβ = β/|β| we thus find where ξ 2 = |x − X| 2 and ρ 2 = |y| 2 , and where we have used y † = ρ 2 y −1 . Inserting this into the gauge field (A.1), we compute where f ≡ ξ 2 /(ξ 2 + ρ 2 ), g ≡ (x − X)/ξ. This is the well-known BPST solution. The field strengths can easily be computed from Eq. (A.8). In the single-instanton case L = ρ 2 + ξ 2 is a scalar and taking its inverse is trivial. We find With η aM N η bM N = 4δ ai δ bi we obtain

A.2 Two instantons
For the two-instanton solution we write with quaternions y 1 , y 2 , X 1 , X 2 , w, from which we compute L according to Eq. (A.7), The second line is obtained as follows. From the requirement that L be symmetric, we find where we set c = 0 [63]. This yields where we have defined ρ 1 ≡ |y 1 |, ρ 2 ≡ |y 2 |. Finally, we have rewritten where the right-hand side is manifestly real because it is composed of scalar products of quaternions.
We can now compute the field strengths analogously to the case of the single instanton. We first need to compute U † = (α † , β † ). For the two-instanton solution, α remains a quaternion, while β becomes a 2 × 1 quaternionic matrix, which we write as β † = (β † 1 , β † 2 ) with quaternions β 1 , β 2 . Then we can express Eq. (A.3) solely through quaternions, We manipulate the resulting two quaternionic equations as follows. We multiply the first one from the left with y 1 , and the second one with y 2 , and we defineŷ i = y i /ρ i for i = 1, 2.
Then, the first equation yields and subtracting the second from the first equation gives where we have defined with the definition Q ≡ŷ From the condition 1 = U † U = |α| 2 + |β 1 | 2 + |β 2 | 2 we find and thus we obtain the final result for U , where we have denotedβ ≡β 2 for brevity. We can insert this result into Eq. (A.8) to compute the field strengths. Denoting the elements of L −1 by L −1 ij ∈ R, with L −1 12 = L −1 21 , we obtain (A.31) After some algebra we find Since all terms are scalar products of quaternions, F 2 M N is proportional to the unit matrix (while F M N isn't), as for the single-instanton case. In contrast to the single instanton, there is a nontrivial dependence on the orientations of the instantonsŷ 1 ,ŷ 2 ∈ SU(2). Since only the relative rotation between the two instantons matters, we parameterize [62] y 1 = 1 ,ŷ 2 = e i θ· σ = cos θ + iˆ θ · σ sin θ , (A.33) which yieldsŷ 1 ·ŷ 2 = cos θ ,ŷ 1 ×ŷ 2 = iˆ θ · σ sin θ . (A.34) An SU(2) transformationŷ can be formulated in terms of a rotation R ∈ SO(3) as follows. Start from an SU(2) matrix G and its transformationG =ŷ † Gŷ, insert G = G a σ a and G =G a σ a , multiply by σ b , and take the trace to obtaiñ Thus, in our case, the rotation corresponding toŷ 2 is Tr[σ a e −i θ· σ σ b e i θ· σ ] = δ ab (cos 2 θ − sin 2 θ) − 2 abcθc cos θ sin θ + 2θ aθb sin 2 θ . (A.36) We now write X 1/2 = X ± D and set X = 0, i.e., the origin of the four-dimensional space is chosen to be the center point of the line that connects the two instantons. In this case, X † 1 X 2 = −|D| 2 , and thus (X † 1 X 2 ) · (y 2 × y 1 ) = 0 since any cross product is orthogonal to the unit matrix direction, as remarked below (A.6), and from Eq. (A.19) we recover Eq. (21) in Ref. [62]. We choose the separation to be in the x 1 direction and to have length δ, i.e., D = iδσ 1 /2. This yields such that L from Eq. (A.19) is given by the entries We can now invert L and compute the field strength squared (A.32) explicitly. As a check, this can also be computed with the help of Eq. (A.11). Here we restrict ourselves to "defensive skyrmions" [62], where θ = 0, i.e., R = 1. Also setting ρ 1 = ρ 2 ≡ ρ for simplicity, we obtain Eq. (2.12) in the main text.

B Multiple instanton layers
If we allow for multiple instanton layers, Eq. (2.19) is generalized to where we have written the intra-layer and cross-layer interaction terms separately. We have denoted the number of layers by N z and their location by z n , and in this general notation each layer is allowed to contain a different number of instantons, denoted by i n . From now on we assume all layers to have the same number of instantons, i 0 = . . . = i Nz−1 ≡ N x , such that N I = N x N z is the total instanton number. Moreover, we neglect the cross-layer interactions. Including them does not pose a conceptual problem, we leave this improvement for future work. As in the main part, we employ the nearest-neighbor approximation for the intra-layer interactions and take the spatial average to obtain which is the generalization of Eq. (2.20). Now, since the instanton number in the spatial lattice is different from the total instanton number, Eq. (2.22) is modified to v/N x = δ 3 /r, and thus Eq. (2.17) becomes d = γ/ρ[6π 4 rN z /(λ 2 n I )] 1/3 . Following Ref. [58], we may use the following ansatz for the location of the layers, with n = 0, . . . , N z − 1, and where the distance between the layers ∆z has to be determined dynamically. This ansatz is completely general for one and two layers (apart from the obvious assumption that the two layers are symmetric with respect to z → −z). For N z > 2 it assumes equidistant layers, which is a simplification in order to limit the number of free parameters. It has been shown that within this ansatz and neglecting instanton interactions three or more layers are never energetically preferred [58]. Although this might be a consequence of the simplistic ansatz, the fact that the same observation was made in a completely different approximation [72] suggests that this result might hold more generally. The Lagrangian (2.25) has the same form as in the single-layer case, with the generalized functions where we have used that q 0 (z) and q int (d, z) are symmetric under z → −z.
Finally, the multi-layer scenario within the ansatz (B.3) requires an additional equation for the minimization of the free energy with respect to ∆z, i.e., we have to add the following equation to the set of equations (2.36), All results in Sec. 3 in the main part concern the single-layer case. With the help of the results of this appendix we did check that including multiple layers does not change our main conclusions: the two-layer solution is preferred for large baryon densities, and thus the critical chemical potential for the chiral phase transition is affected, but this does not disrupt the continuity between the chirally broken and chirally symmetric geometries.

C Integrating the interaction term
Here we derive the result (2.31b), i.e., we compute the z integral over the function q int defined in Eq. (2.14). The integration can be done with the help of the following observations. With the definitions from Eq. (2.16) we computė where the dot denotes derivative with respect to z. Next, for any α, β we have We also find Consequently, combining these relations yields The significance of this relation is that the square root is "absorbed" into S 1 and S 2 , such that only polynomials appear as prefactors of S 1 and S 2 . Putting all this together yields Therefore, for the function with polynomials H 1 , H 2 , we can write the derivative aṡ with the polynomials Now, by identifyingQ int with q int , with n = 3 and m = 1, we find an analytic expression for the integral over q int : the left-hand sides of Eqs. (C.8) are given by the polynomials from Eqs. (2.15). Then, making a polynomial ansatz for H 1 , H 2 , we have reduced the problem to a system of linear equations for the coefficients of these polynomials, which can easily be solved. The result is Eq. (2.31b) in the main text. From this result one can straightforwardly check thatQ int = 2q int (the factor 2 appears due to the requirement that q(u) be normalized to one in the domain u ∈ [u c , ∞]).

D Derivatives of the free energy
Here we explain the derivation of Eqs. (2.36) in some detail since in particular taking the derivative with respect to u c is somewhat complicated. We start by rewriting the free energy ( c η(u → u c ), which is divergent. This divergence is eventually canceled by another term, such that there is nothing wrong with this evaluation, which has in fact been used in Refs. [57,58]. However, to avoid dealing with divergent terms we may evaluate the derivative with respect to u c in the following equivalent, but more elegant, way. We introduce rescaled quantities by k =ku 4 c , n I =n I u 5/2 c , ρ =ρu c , such that Ω = Ω(u c ,ku 4 c ,n I u 5/2 c ,ρu c , γ, d) , (D.9) and by the chain rule the derivative with respect to u c at fixedk,n I ,ρ, γ 0 is relativistic hydrodynamics. We start from the general forms of the current j µ and the stress-energy tensor T µν j µ = nu µ , T µν = ( + P )u µ u ν − g µν P , (E.1) with the number density n, the pressure P , the energy density , the metric tensor g µν = diag(1, −1, −1, −1), and the four-velocity u µ = γ(1, u), where γ = (1 − u 2 ) −1/2 is the usual Lorentz factor with the modulus of the three-velocity u = | u|. The sound modes are obtained from the hydrodynamic conservation equations Using dP = ndµ + sdT and + P = µn + sT , where s is the entropy density, we express all derivatives in terms of derivatives of the independent variables four-velocity u µ , chemical potential µ, and temperature T , (By contracting the second equation with u ν one obtains the equation for the entropy current, which is conserved in the present non-dissipative case. But for our purpose this reformulation is not necessary.) We now introduce fluctuations in the three-velocity as well as in the chemical potential and temperature, u( x, t) = δ u e i(ωt− k· x) , µ( x, t) = µ + δµ e i(ωt− k· x) , T ( x, t) = T + δT e i(ωt− k· x) , where µ and T are equilibrium values (we do not impose an external velocity u on the fluid). Linearizing in the fluctuations yields ∂ µ j µ ie i(ωt− k· x) ω ∂n ∂µ δµ + ∂s ∂µ δT − n k · δ u = 0 , (E.4a) ∂ µ T µ0 ie i(ωt− k· x) ωT ∂n ∂T δµ + ∂s ∂T δT − sT k · δ u + µ∂ µ j µ = 0 , (E.4b) where w = µn + sT is the enthalpy density. We contract the last equation with k i to arrive at three equations for the three fluctuations δµ, δT , k · δ u.  where, in the second step, we have used P = P [µ(n, s), T (n, s)]. All partial derivatives with respect to n are taken at fixed s and vice versa. Next, we change variables from (n, s) to ( , s/n), where s/n is the entropy per particle. We use d = T ds + µdn and d(s/n) = ds/n − sdn/n 2 to find where the derivative is taken at fixed entropy per particle s/n.

F Speed of sound in mesonic and chirally symmetric phases
This appendix gives a brief derivation of the speed of sound in the mesonic phase and the chirally symmetric phase. The relevant free energies are taken from the literature, most conveniently from appendix B of Ref. [57], which uses the same notation as the present paper.

F.1 Mesonic phase
In the mesonic phase, the location of the tip of the connected flavor branes u c is determined from . (F.1) In general, this has to be solved numerically. For small temperatures t we find 2 u c = u (0) c + u (1) c ( t) 6 , (F.2)