Layers of deformed instantons in holographic baryonic matter

We discuss homogeneous baryonic matter in the decompactified limit of the Sakai-Sugimoto model, improving existing approximations based on flat-space instantons. We allow for an anisotropic deformation of the instantons in the holographic and spatial directions and for a density-dependent distribution of arbitrarily many instanton layers in the bulk. Within our approximation, the baryon onset turns out to be a second-order phase transition, at odds with nature, and there is no transition to quark matter at high densities, at odds with expectations from QCD. This changes when we impose certain constraints on the shape of single instantons, motivated by known features of holographic baryons in the vacuum. Then, a first-order baryon onset and chiral restoration at high density are possible, and at sufficiently large densities two instanton layers are formed dynamically. Our results are a further step towards describing realistic, strongly interacting matter over a large density regime within a single model, desirable for studies of compact stars.


Introduction
Cold and dense matter in the interior of compact stars is strongly interacting, governed by Quantum Chromodynamics (QCD). Its phases and properties are poorly known because it is much denser than ordinary nuclei on earth, but not asymptotically dense and thus not quantitatively accessible with weak-coupling methods. The gauge/string duality [1][2][3] is a powerful tool to study strongly interacting matter and has proven to be very useful to get insight into hot QCD matter at low baryon densities produced in heavy-ion collisions [4]. It is thus natural to ask whether we can use it to learn something about cold and dense matter too [5][6][7]. Dense QCD is expected to have a very rich phase structure, including color-superconducting quark matter [8], and at present there are no holographic approaches that can be expected to predict reliably any details of this phase structure. Here we are asking a more modest question, which nevertheless may turn out to be valuable for the study of compact stars. We are asking whether the Sakai-Sugimoto model [9][10][11], a certain realization of the gauge/string duality that comes as close to QCD as currently possible, can be used to understand at least the gross thermodynamic properties of dense nuclear and quark matter and possibly the transition between them, ignoring all complications such as Cooper pairing of nucleons or quarks. The main point of our current effort, started in Ref. [12], is to first find a feasible approximation that gets the basic properties of dense matter right, and then, in future studies, to apply this approximation to the physics of compact stars. In particular, we are interested in the onset of nuclear matter, which must show a discontinuity in the baryon density due to the finite binding energy, and in the transition to quark matter, which is expected to happen at high densities and which is needed to investigate hybrid stars, i.e., compact stars containing quark matter in the core, surrounded by nuclear matter. Although the Sakai-Sugimoto model is a top-down approach, our study should not be understood as a first-principle calculation because we apply various approximations and simplifications. We rather aim at a model description of dense matter, which has some advantages over many of the field-theoretical models used in the same context: we employ a genuine strong-coupling formalism, we can account for nuclear and quark matter in a single model, and the model has very few parameters (3 in the version we consider: the 't Hooft coupling λ, the Kaluza-Klein mass M KK , and the asymptotic separation of the D8-and D8-branes L).
Baryons in the Sakai-Sugimoto model are introduced as D4-branes wrapped around the 4-sphere of the background geometry, following the general concept of baryons in the gauge/string duality [13,14]. Here, this is equivalent to gauge field configurations with nonzero topological charge on the connected flavor branes of the model [10], and various properties of baryons in the vacuum have been studied within this approach [10,11,[15][16][17]. Baryonic matter at nonzero density and temperature was first considered in a pointlike approximation of the instantons on the flavor branes [18], and this approach was improved and complemented by a number of studies [12,[19][20][21][22][23][24][25][26][27][28][29][30][31]. (For studies of baryonic matter using a different holographic approach, based on a D3-D7 setup, see for instance Refs. [32,33].) The idea of the present paper is to improve the instanton gas approach, introduced in Ref. [26] and further developed in Ref. [12]. More specifically, it is known that away from the λ = ∞ limit the Sakai-Sugimoto instantons are anisotropic in the sense that they break the SO(4) symmetry of rotations in the space of the holographic coordinate and the three spatial dimensions [28,29]. We shall account for this anisotropy by introducing a "deformation parameter" into the standard flat-space instanton solution. Furthermore, it has been argued that the repulsion between the instantons makes them spread out in the holographic direction [19], which is realized for instance in crystalline structures in the confined phase of the model [24,31]. We introduce this repulsive effect by allowing for an arbitrary number of instanton layers in the bulk and determine this number and the distance between the layers dynamically as a function of the baryon chemical potential.
Unless backreactions of the flavor branes on the background geometry are taken into account, cold matter in the Sakai-Sugimoto model does not deconfine, which is in accordance with expectations from large-N c QCD [34]. In order to allow for a transition between nuclear and quark matter at low temperatures, we work in the "decompactified" limit of the model: if the separation of the flavor branes L is sufficiently small, the deconfined geometry of the model has a chirally symmetric and a chirally broken phase. Therefore, we are able to include the transition from nuclear matter to quark matter without the complications of backreacting flavor branes. Varying L from its maximum value, as used in the original works [10,11], to very small values is best understood as changing the dual field theory: the limit of maximal L, i.e., an antipodal separation of the flavor branes in the space of the compactified extra dimension of the model, is related to large-N c QCD; the limit of very small L, on the other hand, corresponds to a field theory comparable to a Nambu-Jona Lasinio (NJL) model [35][36][37], and its rich phase structure in the deconfined geometry is possibly closer to real-world QCD than the antipodal limit, at least with respect to the chiral phase transition.
The paper is organized as follows. In Sec. 2 we explain our ansatz and derive the free energy and its stationarity equations. This is done by first discussing the Dirac-Born-Infeld action in Sec. 2.1, including a very general form of the non-abelian field strengths in Sec. 2.1.1, our specific ansatz for the anisotropic instantons in Secs. 2.1.2, the approximations for our many-instanton system in Sec. 2.1.3, and the symmetrized trace prescription in Sec. 2.1.4. Then, in Sec. 2.2, we add the Chern-Simons contribution to obtain the full Lagrangian, and in Sec. 2.3 we explain how we solve the system, including the minimization of the free energy. Sec. 3 is devoted to the numerical results and is divided into two subsections: in Sec. 3.1 we minimize the free energy with respect to all parameters of the ansatz, while in Sec. 3.2 we impose certain constraints on the shape of the single instantons, increasing the number of free parameters of our model to 5. We give our conclusions in Sec. 4.

Setup
The general setup follows numerous other works in the Sakai-Sugimoto model, and for all details and foundations of the model we refer the reader to the original works [10,11] or reviews [38][39][40]; the notation we are using is consistent with Ref. [12]. Our starting point is the action for the gauge fields on the flavor branes, which consists of a Dirac-Born-Infeld (DBI) and a Chern-Simons (CS) part, (2.1) We now discuss these two contributions separately.

Dirac-Born-Infeld action
The DBI part is given by Here, the integral is taken over imaginary time τ with the temperature T , over position space X = (X 1 , X 2 , X 3 ), and over the holographic coordinate U ∈ [U c , ∞]. In this section, we discuss the chirally broken geometry, where the D8-and D8-branes are connected, with U c being the location of the tip of the connected branes. The model contains a compactified direction X 4 , whose radius is expressed in terms of the inverse Kaluza-Klein mass M KK , Figure 1. Illustration of the three phases whose free energies are compared in this paper. It shows the cylinder-shaped subspace of the deconfined geometry, spanned by the compact extra dimension X 4 (with radius M −1 KK ) and the holographic coordinate U , and the D8-and D8-branes, which can either be connected (chiral symmetry spontaneously broken, left and middle figure) or disconnected (chiral symmetry restored, right figure). They are asymptotically separated by a distance L, and in the chirally broken phases their (density-dependent) embedding, including the location of the tip U c , has to be determined dynamically (we assume no backreaction on the background geometry). Baryon number in the chirally broken phase is introduced through instantons on the flavor branes, here symbolized by two circles. Our ansatz allows for an arbitrary number of instanton layers N z in the bulk (see Sec. 2.1.3), but we shall find that more than two are never energetically preferred. We assume that L π/M KK , which is the "decompactified" limit of the model, where the deconfined geometry extends down to arbitrarily small temperatures (apart from using this fact, which allows us to set the temperature to zero, U T = 0, we never make any assumptions about L and M KK in our calculation). For the discussion of the instantons we sometimes switch to an alternative holographic coordinate Z along the connected branes, as indicated in the middle figure.
X 4 ≡ X 4 + 2π/M KK , and the embedding of the flavor branes in the background geometry is given by X 4 (U ). This function has to be determined dynamically and is subject to the boundary condition X 4 (U → ∞) = ±L/2, where L is the asymptotic separation of the D8and D8-branes. The chirally broken geometry accommodates the baryonic phase, discussed in this section, and the mesonic phase, which is well known and whose free energy we shall later need and simply quote from the literature. In the chirally restored phase the flavor branes are straight, X 4 (U ) = ±L/2, and disconnected, and again we will quote the corresponding free energy later and use it in the energy comparison with the baryonic phase. The geometry of the three different phases is shown schematically in Fig. 1. The dilaton field is e Φ = g s (U/R) 3/4 , where R is the curvature radius and g s the string coupling. Moreover, α = 2 s with the string length s , T 8 = 1/[(2π) 8 9 s ] is the D8-brane tension, V 4 = 8π 2 /3 is the volume of the 4-sphere, and "str" denotes the symmetrized trace (we shall discuss below the procedure that we follow to evaluate this trace). We work in the deconfined geometry, whose induced metric on the flavor branes g is given by where i = 1, 2, 3, dΩ 2 4 is the metric of the 4-sphere, and we have abbreviated where U T is related to temperature T and curvature radius R via In our final results we shall restrict ourselves to T = 0. Strictly speaking, the preferred geometry at zero temperature is the confined one. However, in the decompactified limit L π/M KK , the critical temperature for deconfinement is much smaller than the critical temperature for chiral restoration (at zero chemical potential). One may think of letting M KK → 0 while keeping L fixed; this renders the region of the confined geometry in the phase diagram arbitrarily small and justifies our zero-temperature approximation. In the decompactified limit, the structure of the phase diagram (without baryonic matter) is similar to what is obtained in an NJL model, where there is no confinement either. Two differences to NJL are that our formalism allows for a well-defined way to implement baryons (which are rarely included in NJL studies, although it is possible [41]) and that in the NJL model it is easy to include nonzero current quark masses (which is difficult in the Sakai-Sugimoto model, although it is possible [42][43][44][45]). As a consequence, in our calculation, the chiral phase transition is always a phase transition in the strict sense (in fact, it turns out to be always a first-order phase transition), while in the NJL model with quark masses (and in nature) this transition is allowed to be a continuous crossover because chiral symmetry is not an exact symmetry.
We work with two flavors, N f = 2, and express the U (2) field strengths in terms of the gauge fields A µ , µ = 0, 1, 2, 3, U , in which we separate the abelian from the non-abelian part, with the Pauli matrices σ a , normalized such that [σ a , σ b ] = 2i abc σ c . Consequently, with In our ansatz the only nonzero abelian field strength isF 0U , where the quark chemical potential will be included as the boundary value ofÂ 0 , and the baryons are implemented through the non-abelian field strengths F iU , F ij . All other field strengths are set to zero.
The trace over the square root in the DBI action is not uniquely defined in the nonabelian case, and thus we need to follow a certain prescription. Here we follow Tseytlin [22,46]: we first compute the determinant over space-time indices as if the field strengths were numbers, This expression factorizes if (F ij F kU ijk ) 2 = 2F 2 iU F 2 ij (which shall be fulfilled by our ansatz), The prescription then requires us to expand the square root over this determinant to all orders in α , apply the symmetrized trace for each term separately, and then resum the resulting infinite series. Within our ansatz, this resummation can be done analytically and leads to a relatively simple analytic form for the DBI action, see Sec. 2.1.4. Nevertheless, the numerical evaluation turns out to be more difficult compared to the simpler prescription that uses the standard (unsymmetrized) trace [12,19]. Thus, after showing numerically in Sec. 3.1 that the results of the two prescriptions do not differ much we shall resort to the unsymmetrized prescription in Sec. 3.2.

General form of non-abelian field strengths
The ansatz for the non-abelian part is best discussed in the new holographic coordinate Z, defined as 1 such that Z = 0 corresponds to the tip of the connected flavor branes, and Z = ±∞ to the holographic boundary on the D8-and D8 branes, see Fig. 1. The most general ansatz 1 Later we shall come back to using the coordinate U in many equations. This is partly to connect to previous literature, and partly because of convenience. Neither of the two coordinates turns out to be overly superior when it comes to compactness in notation or convenience in the calculation.
that is SO(3) symmetric in the spatial directions for an instanton located at X = 0 can be written as [17,28,29] where X = | X|, and a Z , a X , φ 1 , φ 2 are all functions of X and Z. From this ansatz we compute the non-abelian field strengths needed in Eq. (2.8) (summation over a = 1, 2, 3) whereˆ X ≡ X/X, and 14) The field strengths squared are obviously linear combinations of the products σ a σ b . Except for the non-diagonal structureˆ X · σ in F ij F kZ ijk , there are only diagonal terms, σ 2 a = 3, X 2 a σ 2 a = X 2 . We have written the Pauli matrices explicitly in the results because this is needed for the discussion of the symmetrized trace, see Sec. 2.1.4.

Anisotropic instantons
We now specify our ansatz for the gauge fields. To this end, it is convenient to work with dimensionless coordinates defined as (2. 15) In these coordinates, our ansatz for a single instanton placed at (2.17) This corresponds to the Belavin-Polyakov-Schwarz-Tyupkin (BPST) instanton [47], where the z coordinate is rescaled with respect to x by a factor γ. In appendix A we derive the single-instanton solution in the deconfined geometry, which does show a nontrivial (and temperature dependent) γ, see Eq. (A.7). Here we will treat γ as a free parameter, accounting for the "deformation" of the instanton. Such a deformation has also been observed in the full solution of a single baryon in the vacuum [28,29]. The instanton width in the spatial direction is ρ/γ, while the width in the holographic direction is ρ (for convenience, we shall often simply refer to ρ as the instanton width). For a given ρ, the deformation parameter thus has the effect of stretching the instanton along the holographic direction z (large γ) or along the radial spatial direction x (small γ). A single instanton becomes elongated along x and wider in both x and z for values of the 't Hooft coupling λ away from infinity, which was shown in a full numerical evaluation of the equations of motion, based on the most general ansatz for the gauge fields (2.12), see Ref. [28]. (Already from the SO(4) symmetric case we know that only the finiteness of λ prevents an instanton in the vacuum from being pointlike [15].) Translated to our parametrization, we thus expect a smaller γ and a larger ρ for finite λ than for λ = ∞. Our approximation, based on the ansatz (2. 16), and extended to a many-instanton system below, is too simplistic to allow for a dependence on λ apart from a trivial rescaling (this is in contrast to the "homogeneous ansatz" [12,19], which is not based on any instanton solution). Therefore, besides computing ρ and γ dynamically in Sec. 3.1, we shall impose certain constraints on ρ and γ in Sec. 3.2, with the idea of capturing some of the λ < ∞ physics, which seems to be crucial to obtain more realistic results, already for baryons in the vacuum [17,28,29]. With the ansatz (2.16), the field strengths become particularly simple. The nondiagonal term and all terms proportional to X a σ a vanish, 18) and the remaining terms become proportional to the same function of x and z, The field strengths now fulfill the relation (F ij F kZ ijk ) 2 = 2F 2 iZ F 2 ij , and thus we may use Eq. (2.9). This leads to the DBI action (2.20) where we have replacedÂ 0 → iÂ 0 since we work in Euclidean space, have introduced the dimensionless quantitieŝ and have denoted the derivative with respect to u by a prime. Also, we have abbreviated N ≡ 2T 8 V 4 R 5 (M KK R) 7 /g s and and used the relation λ 2 s = 2R 3 M KK .

Spatial average and instanton layers
Next, we go from a single instanton to a many-instanton system. We do so on the level of the field strengths squared. We place i n many instantons at the points z n in the bulk, n = 0, . . . , N z − 1 (N z ≥ 1), and distribute them at the points x in , i = 1, . . . , i n , in position space. The total number of instantons is N I ≡ i 0 +. . .+i Nz−1 . One can think of an instanton lattice sitting at each point z n in the bulk. In this general notation, the lattice structure is allowed to be different at different points in the bulk. However, we shall drastically simplify this general picture in our calculation such that the lattice structure in position space becomes irrelevant: we shall average the field strengths squared over position space before solving the equations of motion [12,26], and as a consequence it does not matter at which points x in the instantons sit. The many-instanton system within our approximation is thus obtained by replacing where, in the second step, we have assumed that the same number of instantons sits at every z n and denoted this number by N x , such that the total number of instantons is now N I = N z N x , and where we have defined the normalized instanton profiles It is important that the deformation parameter γ has not dropped out, although we have averaged over position space. We can thus later compute the deformation of the instantons Figure 2. Instanton distribution along the holographic direction z. The number of layers N z and their extension z 0 will be determined dynamically. The solid lines are the separate terms in the sum of Eq. (2.24b). For N z → ∞ at fixed z 0 , the sum over all instanton layers approaches the function D ∞ (z) (2.31), here shown as a dashed line. Each instanton layer in the holographic coordinate accommodates i n instantons in position space, and we assume i 0 = . . . = i Nz−1 ≡ N x , the total instanton number thus being N I = N x N z . even though our simplified equations of motion only involve the holographic coordinate z, and not x.
For the instanton distribution in the bulk we employ the following ansatz. We assume the layers of N x instantons to be separated equidistantly by a distance ∆z from each other, and to be centered at the points such that they are spread over a symmetric interval of length 2z 0 around the tip of the connected flavor branes z = 0, and z 0 = (N z − 1)∆z/2. This is illustrated in Fig. 2. This ansatz, where the instanton layers all have the same shape and are separated by the same distance, allows for a continuous transition between N z = 1 and N z = 2, but all other transitions -if they occur -will necessarily be discontinuous. Inserting Eq. (2.23) into Eqs. (2.22) yields where we have introduced the dimensionless instanton density (per flavor, hence the division by N f = 2) and defined ).] For N z = 1 we recover the function q(u) from Ref. [12], where the instanton repulsion was neglected, We shall treat z 0 and N z as dynamical parameters with respect to which we minimize the free energy. We include the possibility of infinitely many instanton layers, i.e., a smoothly smeared instanton distribution. In this limit, letting N z → ∞ while keeping z 0 fixed, we can approximate the sum in Eq. (2.24b) by an integral and we obtain which is also shown in Fig. 2.

Symmetrized trace
As explained above, we need to decide on a certain prescription to evaluate the non-abelian DBI action. We expand the square root and take the symmetrized trace in each term separately [46]. It is known from string theory that this prescription is accurate up to O(F 4 ) [48]. For the structure we have in Eq. (2.20), this yields [22] str (1 + ϕσ a σ a )(1 + ψσ a σ a ) The equations of motion in this prescription as well as the stationarity equations for the free energy are worked out appendix B. If we instead take the standard trace in this series, we obtain Since ϕ, ψ ∝ F 2 , this result differs from Eq. (2.32) starting from terms of order F 4 . We thus expect different results for large densities. Below we shall present a comparison of the two prescriptions, showing that there is indeed a difference. However, this difference turns out to be small and the results are qualitatively the same, see Fig. 4. Therefore, we shall mostly (in the equations in the main text and in all results except for Fig. 4) use the unsymmetrized prescription, which leads to significantly simpler equations, resulting in much faster numerics.

Chern-Simons action and full Lagrangian
Within our ansatz, the CS action is Having computed the field strengths and having introduced convenient dimensionless quantities, we can easily compute this contribution. We obtain F ij F kZ ijk from Eqs. ( Putting this together with the DBI action in the unsymmetrized prescription, this yields the action with the Lagrangian with g 1 and g 2 from Eqs. (2.27). This Lagrangian has exactly the same form as the one used in Ref. [12], see Eq. (30) in that reference. The extensions of the present approach are hidden in the functions g 1 , g 2 : we reproduce the functions g 1 , g 2 of Eq. (31) in Ref. [12] by considering only one instanton layer, N z = 1, z 0 = 0, and by choosing the instanton deformation to be γ = 3u 3/2 c /2 (this specific value was chosen by transferring the BPST result of the confined geometry to the deconfined geometry) 2 .

Minimizing the free energy
The equations of motion forâ 0 and x 4 , obtained from the Lagrangian (2.37), are, in integrated form, where k is an integration constant, and The equations of motion (2.38) can easily be solved forâ 0 and x 4 algebraically. The resulting expressions can then be inserted into the (dimensionless) free energy density of the baryonic phase, where the Lagrangian is given in Eq. (2.37) and where, in the second line, we have employed partial integration in the CS term and used the boundary conditionsâ 0 (∞) = µ, x 4 (∞) = /2, with the dimensionless chemical potential µ and the dimensionless asymptotic separation of the flavor branes = M KK L. (The complete, dimensionful free energy density is obtained by multiplying Ω baryon with N N f .) We note that the asymptotic behavior of a 0 and x 4 is given by This confirms that n I is the (dimensionless) baryon density, which is also given by the derivative of the free energy with respect to the chemical potential, This equation seems like an obvious thermodynamic relation, but there are some subtleties in the Sakai-Sugimoto model if baryon number is (partially) created through a magnetic field [49,50]. In that case, a modified Chern-Simons term has been used to ensure the relation (2.43) [23,[49][50][51][52]. Here, no such modification is necessary.
To be precise about the meaning of our dimensionless quantities, we notice that µ is a dimensionless quark chemical potential, while n I is a dimensionless baryon number density. The physical quark chemical potential is related to µ by the factor introduced in the definition of the dimensionless abelian gauge field in Eq. (2.21), Inserting this relation into Eq. (2.43) and using that the dimensionful free energy density is N N f Ω baryon , we read off the physical quark number density, i.e., With Eq. (2.28) we conclude that the dimensionful baryon number density (= quark number density divided by N c ) is exactly the instanton density N I /V . The free energy (2.41) is a function of the parameters k, n I , u c , ρ, γ, z 0 , N z . They are independent of each other except for the obvious condition that the separation of instanton layers vanishes, z 0 = 0, if and only if there is exactly one instanton layer, N z = 1. We shall discuss the following two approaches and present their results in Secs. 3.1 and 3.2, respectively.
(i) Minimize Ω baryon with respect to all seven parameters.
(ii) Impose the following constraints on the parameters that determine the shape of a single instanton, i.e., the instanton width ρ and instanton deformation γ, and fix ρ 0 , γ 0 . Then minimize Ω baryon with respect to the remaining five parameters k, n I , u c , z 0 , N z .
Approach (i) requires no further motivation, it yields the ground state that the system chooses to be in within the given approximation. The idea behind approach (ii) is as follows. We do not know how our many-instanton ansatz is related to the full solution of the problem. Therefore, we have to take the result of the straightforward minimization of scenario (i) with some care: the minimum in our restricted parameter space might be very different from the minimum in the full functional space. However, as mentioned below Eq. (2.17), we do know some features of the full solution of single instantons in the vacuum, in particular we know that the width and the deformation change as a function of λ away from the λ = ∞ limit. In order to consider a many-instanton system, we have given up some complexity, in particular we do not expect our ansatz to reproduce these important λ < ∞ features of single instantons. Therefore, we choose to impose external constraints on width and deformation and scan through the resulting parameter space. We might simply have "rigidly" fixed ρ and γ. However, we do expect these quantities to change with density. Therefore, we have chosen a particular scaling with (the density-dependent) u c . This "natural" scaling is chosen such that u c drops out of but one minimization equations, as we shall see below. (The factor 3/2 in the scaling relation for γ is chosen such that γ 0 = 1 corresponds to the calculation done in Ref. [12].) The use of the constraints on ρ and γ is justified a posteriori by the observation that only in approach (ii) we do find a layered structure of the instantons in the bulk, whose existence is suggested from other, complementary, approximations in the literature.
In both approaches (i) and (ii), the parameters are determined by setting the various derivatives of the free energy to zero. The parameter N z is special because it is an integer and thus we cannot simply take the derivative of the free energy with respect to N z . Instead, we will solve the below equations for various values of N z , including N z = ∞. It turns out that there is a clear tendency in the behavior of the free energy as a function of N z , and thus this procedure is sufficient to determine the preferred N z . Since we have written the free energy in the same form as in Ref. [12], we can skip the details of the derivation of the stationarity equations and directly quote the results. Anyway, only the derivative with respect to u c is not completely straightforward, see Sec. III of Ref. [12] or appendix B of the present paper, where we go into some details in the context of the symmetrized trace prescription. The resulting equations (in the order: minimization with respect to k, n I , ρ, γ, z 0 , u c ) are where we have used the abbreviation and where, in Eq. (2.47f), we have defined with c 1 giving the behavior of x 4 close to u c , Notice that the minimization with respect to k (2.47a) is nothing but the condition that the asymptotic separation of the flavor branes be . We thus have to solve 6 coupled equations for k, n I , u c , ρ, γ, z 0 in approach (i) and 4 coupled equations for k, n I , u c , z 0 in approach (ii) (and do so for various values of N z ). However, in both cases, two equations decouple. First, we observe that the only explicit appearance of µ is in Eq. (2.47b). Therefore, rather than fixing µ we can fix n I and determine the corresponding µ with the help of Eq. (2.47b) after we have solved the other equations. (This is also advantageous because µ is always a single-valued function of n I , while n I can become a multi-valued function of µ.) Second, we can rescale all quantities with appropriate powers of u c and introduce the new integration variable u/u c . One can show that this eliminates u c from all equations except for Eq. (2.47a). Hence, Eq. (2.47a) also decouples, we can solve the remaining equations for the rescaled quantities, then compute u c from Eq. (2.47a) and then undo the rescaling with the help of the resulting u c . A similar rescaling of all quantities with the externally given parameter eliminates from all equations. As a consequence, all results scale with in a simple way; different 's do not lead to qualitatively different results.
Since the rescaling with u c , in particular together with our two approaches (i) and (ii), may be somewhat confusing, let us explain this in more detail. In deriving Eq. (2.47f) we have first applied Eq. (2.46) and then taken the derivative with respect to u c . For scenario (i) this is not very crucial because minimizing with respect to u c , ρ, γ is equivalent to minimizing with respect to u c , ρ 0 , γ 0 . [To see this, consider the function Ω baryon = Ω baryon (u c , ρ, γ) and take the derivatives with respect to u c , ρ, γ on the one hand and, via the chain rule, with respect to u c , ρ 0 , γ 0 on the other hand. The apparent additional terms created in the latter procedure are all zero because the derivatives with respect to ρ 0 and γ 0 are required to vanish.] In scenario (ii) it is crucial to correctly capture the dependence on u c within ρ and γ, because we do not minimize with respect to ρ 0 and γ 0 . In both scenarios, the eventual rescaling with u c (where ρ 0 and γ 0 by construction do not scale anymore with u c ) is then merely a convenient trick to simplify the numerical evaluation. As a check, we have also evaluated the equations without this eventual rescaling and found the same result.
Once the minimum of the free energy is found within our ansatz for baryonic matter, we need to compare the value of Ω baryon at that minimum with the free energies of the mesonic phase (= chirally broken phase without nuclear matter) and the quark matter phase (= chirally restored phase), see Fig. 1. We shall restrict ourselves to zero temperature, although the equations derived above for the baryonic phase provide the full temperature dependence. At zero temperature, the free energy of all three phases at the stationary point can be written in the very compact form where Λ is an ultraviolet cutoff, replacing the boundary u = ∞, and where n I and k remain to be determined numerically in the baryonic phase and have simple analytic forms in the mesonic and quark matter phases, see below. All free energies show the same constant divergence which becomes irrelevant when we compare them which each other. In the baryonic phase, the form (2.51) is derived as follows. We start from Eq. (2.41), rescale c γ, and introduce the new integration variable u/u c . We do not rescale the externally given quantities µ and . Neither would we rescale T , but we have already set T = 0 and thus f T = 1. Then, we extremize the resulting expression with respect to u c , taking into account the u c dependence in the upper boundary of the integral, which now is Λ/u c . The condition that the derivative of Ω baryon with respect to u c vanishes, gives exactly Eq. (2.51). As an aside, this procedure also gives an alternative form of the minimization with respect to u c (2.47f). For the mesonic and quark matter phases we use the well-known results (for instance from appendix B of Ref. [12]) and note that at T = 0 they have the form (2.51) with quark matter: (2.52b) In particular, the free energy of the quark matter phase does not depend on , while the free energy of the mesonic phase does not depend on µ.

Fully dynamical instanton width and deformation
In this section, we evaluate and discuss approach (i), i.e., we minimize the free energy with respect to all free parameters, including the instanton width ρ and the deformation γ. For the minimization with respect to z 0 (2.47e) we observe that for small z 0 This implies that z 0 = 0 always solves Eq. (2.47e), and thus one solution of our system is always found by solving the remaining equations with z 0 = 0. In those equations, then, the number of instanton layers in the bulk N z only appears implicitly in n I ∝ N x N z /V , which is determined dynamically anyway. We thus do not have to compute the free energies for various values of N z separately to find the ground state. This has to be done only for the solutions with z 0 > 0. To search for such solutions let us first assume that there is a continuous transition from z 0 = 0 to z 0 > 0 as a function of µ. The critical chemical potential for this transition can be found by dividing Eq. ; without loss of generality, we have set N z = ∞ for this plot.] A zero of this function would give a critical chemical potential at which a solution with more than one layer, z 0 > 0, starts to exist. The numerical result shows the absence of such a critical chemical potential (and suggests that a layered structure is approached asymptotically at µ = ∞). The plot does not exclude the possibility of a discontinuous transition to z 0 > 0, but we have not found such a transition.
numerically by computing the right-hand side of Eq. (2.47e), with n I , ρ, γ and k from the z 0 = 0 solution inserted. We have plotted the result as a function of µ in Fig. 3. A zero of the plotted curve would correspond to a critical chemical potential for the onset of a second instanton layer. We see that there is no zero. Interestingly, multiple layers seem to be "postponed" to infinitely large densities because the plotted function approaches zero asymptotically for µ → ∞. This argument does not exclude that there is a discontinuous transition to a phase with multiple instanton layers. In a numerical search we have not found any solution z 0 > 0, but a rigorous proof for that absence is difficult. We discuss the only solution we have found, N z = 1, z 0 = 0, in the following, and come back to solutions that show instanton repulsion in Sec. 3.2, where we work with approach (ii), in which case we do find solutions with z 0 > 0, both in a continuous and a discontinuous transition, depending on the values of ρ 0 and γ 0 .
The results for N z = 1, z 0 = 0 are shown in Fig. 4, which leads to the following observations 3 .
3 All quantities in this and all following figures are rescaled with appropriate powers of , i.e., µ stands for 2 µ, nI for 5 nI etc. If we wish to assign physical units to the plot, we have to choose values for the three parameters of the model, say MKK = 950 MeV, λ = 16, and L = 0.3 π/MKK. Then, using the relations (2.44) and (2.45) for the dimensionful chemical potential and density, the maximum baryon chemical potential in the two upper panels is 3.3 GeV (about 3.6 times the chemical potential of the realworld baryon onset), while the baryon density at that point is 5.3 fm −3 (about 35 times real-world nuclear saturation density). Here we are only interested in the qualitative properties of our approximations -which, in Fig. 4, are not in agreement with real-world baryonic matter -and thus these numbers do not mean much. • The transition from the mesonic to the baryonic phase is second order, as can be seen from the continuity of the baryon density in the upper left panel. This implies the absence of a binding energy, in contradiction to real-world nuclear matter. This result is qualitatively the same as for the pointlike approximation of baryons (shown as dashed lines). At the baryon onset, our solution approaches that of the pointlike approximation. Qualitatively, this is exactly the same observation that had already been made in Ref. [12], where the instanton deformation was not determined dynamically and where only one instanton layer was taken into account. We thus conclude that allowing for a dynamical instanton deformation and multiple instanton layers is not sufficient to turn the unphysical second order baryon onset into a physical first order transition.
• The instanton deformation γ decreases just after the onset and then increases for large baryon densities, for very large µ we find γ ∝ µ 2 . Hence at large densities the instanton becomes elongated along the holographic direction. More specifically, the instanton width in the holographic direction ρ increases monotonically with density, with asymptotic behavior ρ ∝ µ 3/2 , while the instanton width in the spatial direction ρ/γ behaves non-monotonically, increasing for small µ and decreasing like µ −1/2 for very large µ. At the baryon onset, the width of the instanton is zero in all directions.
Since the density just above the onset is infinitesimally small, our approximation thus predicts a pointlike baryon in the vacuum. We know that holographic baryons do acquire a width if corrections of finite λ are taken into account. This indicates that our present approximation is too simplistic to yield realistic isolated baryons.
• There is no chiral restoration at large chemical potentials. This can be seen from the lower right panel, where the ratio of the baryonic pressure over the pressure of the chirally restored phase (quark matter) is shown (P = −Ω). Chiral restoration would occur if that ratio were to decrease below 1. While in the pointlike approximation this ratio approaches 1 for µ → ∞ (which can be shown analytically [23]), the ratio appears to saturate at a much larger value, P baryon /P quark ∼ 2.4, for our extended, and deformed, instantons. Again, it is instructive to compare this result to that of Ref. [12], where the deformation was fixed. In that case, chiral restoration did occur.
Here, we allow the system to settle at a lower free energy by adjusting its instanton deformation. As a consequence, the transition to quark matter has disappeared.
• In the lower panels we compare the results for the two different prescriptions for the non-abelian DBI action. We see that they do differ for large chemical potentials, the free energy from the symmetrized prescription is somewhat larger (smaller ratio P baryon /P quark ), however the difference is small and not relevant for our main conclusions. Had we plotted both results in the upper panels, the curves would have been indistinguishable by naked eye. Since the unsymmetrized prescription is much simpler, we shall in the following section only work with it and discard the symmetrized prescription.

Constraints on instanton shape
We now turn to approach (ii), where we impose the constraints (2.46) on ρ and γ. One of the crucial differences to the previous section is that now we do find solutions for all N z . We thus have to compare the free energies of all phases with different numbers of instanton layers in the regime where their solutions coexist. Let us start with a specific choice of parameters, ρ 0 = 2.5, γ 0 = 4. The results for the free energies and corresponding densities are shown in the left column of Fig. 5. [The ratio of the free energy densities of two phases is obviously the same as the ratio of pressures because the minus signs in the free energies simply cancel. However, since a ratio larger than 1 (and a larger pressure, but a lower free energy) is favorable, it is somewhat more natural to label the vertical axis with the ratio of pressures, as in Fig. 4, although in the text we continue to speak of the free energy of the phases.] The first important result is   Figure 5. Pressures and corresponding densities for an instanton width ρ 0 = 2.5 and two different instanton deformations, γ 0 = 4 (left column) and γ 0 = 6.2 (right column). Upper row: comparison of pressures of baryonic matter P Nz , with N z = 1, 2, 3, 4, 5, ∞ instanton layers in the bulk (solid lines) and the mesonic phase P meson (dashed lines). In both panels, there is a first order phase transition from the mesonic phase to baryonic matter. In the left panel, this transition is to the N z = 1 phase, and the solutions for N z > 1 only start to exist at a larger chemical potential. In the right panel, the solutions N z > 1 start to exist below the baryon onset, and the transition is to the N z = 2 phase. There are energetically disfavored branches whose pressure we have not shown (for N z = 1 in the left panel and for all N z in the right panel). Lower row: corresponding baryon densities for N z = 1 and N z = 2, showing all solutions, including the energetically disfavored branches. The dashed vertical lines indicate the baryon onset. The curves for N z = 3, 4, 5, ∞ are not shown since they would be difficult to distinguish from the N z = 2 curve on the given scale. that the solution for N z = 1 is multi-valued in a certain regime of chemical potentials and as a consequence there is a first-order baryon onset, in contrast to the result of the previous section. In the left column, the solutions for all N z > 1 are single-valued and start to exist at the same point. As soon as they exist, their free energy is lower (and their density larger) than that of the N z = 1 solution. We find that the free energy of the N z = 2 solution is lowest, and Ω Nz=2 < Ω Nz=3 < . . . < Ω Nz=∞ < Ω Nz=1 . Therefore, there is a transition from the baryonic phase with a single instanton layer to a baryonic phase where two instanton layers separate in the bulk. This transition is smooth. The right column of Fig. 5 shows the result for the same ρ 0 , but a larger deformation parameter, γ 0 = 6.2. Now, the solutions for all N z are multivalued, and there is a first-order baryon onset directly to the phase with N z = 2. In both cases, we have thus, by imposing the constraint (2.46), arrived at a first-order phase transition to baryonic matter, as expected from real-world nuclear matter.
In Fig. 6 we show the details of the solution obtained for the parameters from the left column of Fig. 5. In the upper row, the instanton profile in the holographic direction D(z), see Eq. (2.24b), is plotted for three different chemical potentials. The chemical potentials chosen here are all above the onset of the N z = 2 solution, i.e., the (red) thick line that represents the solution with lowest free energy has always two maxima, symmetrically placed around z = 0. These maxima move apart with increasing chemical potential. The spreading of the instantons in the bulk with increasing density was observed previously in the Sakai-Sugimoto model and related models in various different approximations. Firstly, it was suggested in the confined phase of the Sakai-Sugimoto model within a simple approximation unrelated to any single-instanton solution [19]. A similar observation, taking into account a crystalline structure in both spatial and holographic coordinates, led to the term "baryonic popcorn" [24,31], referring to a successively increasing number of instanton layers in the holographic direction with increasing density -in contrast to our approximation, where at most two layers are favored. The observation of baryonic popcorn was confirmed in a full numerical calculation within a simpler, 2+1 dimensional model [29,53]. These results in the literature suggest that the occurrence of multiple instanton layers, or, more generally, the spreading of the instantons away from the tip of the connected flavor branes, at large baryon density is a general feature, and it is intriguing that our simple approximation for homogeneous baryonic matter, based on the flat-space BPST instanton solution, shows the same feature, if we enforce the constraints (2.46).
In the second row of Fig. 6 we show the same instanton profiles, but now in the two-dimensional space of holographic and spatial directions. We recall that even though we have averaged over position space before solving the equations of motion, we can still go back to the instanton profile D(x, z) from Eq. (2.24a) and ask how this profile looks for different chemical potentials. Although we have fixed ρ 0 and γ 0 , the instanton width ρ and deformation γ remain nontrivial functions of the chemical potential due to their dependence on u c , which is determined dynamically. The figure shows that the instantons not only develop a second layer at large µ, but also become elongated in the holographic direction: with increasing density, the instantons get wider in the holographic direction because ρ ∝ u c increases (already obvious from the first row of the figure), and narrower in the spatial direction because ρ/γ ∝ u −1/2 c decreases. The third row of Fig. 6 shows the embedding of the flavor branes in the background geometry for the same chemical potentials as the first two rows. These plots illustrate the instanton layers in the subspace spanned by u and x 4 . They show in particular that, due to the flatness of the brane profile just above u c , a small distance between the instanton layers in the z (or u) coordinate can result in a large distance along the brane profile, if the instantons are close to u c . We also see that the shape of the profile does not change qualitatively with density (apart from moving up), even though the instantons move from the flat part of the profile up to the almost vertical segments of the branes. It does not change much with the number of instanton layers either, which can be seen from the comparison with the N z = ∞ embedding (we have checked that the embeddings for other values of N z  Figure 6. Instanton profiles (first two rows) and corresponding embedding functions of the flavor branes (third row) for ρ 0 = 2.5, γ 0 = 4 and three different chemical potentials, µ = 1 (left column), µ = 3 (middle column), µ = 6 (right column). First row: profiles D(z) in the holographic direction z. For chemical potentials above a certain critical potential (µ 0.43 for the given parameters, see left column of Fig. 5) solutions for all N z exist, but the energetically preferred solution is N z = 2 [thick (red) line]. We have also plotted the energetically disfavored solutions N z = 1, 3, ∞ [thin (black) lines]. Second row: energetically preferred profiles D(x, z) in the holographic and radial directions, showing that the instanton gets elongated along the z direction for large densities. Note that the vertical axis in the first row and the color scale in the second row is adjusted for each panel (whereas the z and x intervals are fixed). Third row: the solid (red) curve is the embedding of the preferred N z = 2 solution, with the dots marking the centers of the two instanton layers (this is the calculated version of the cartoon in Fig. 1). The dashed (black) curve, which is barely distinguishable from the solid curve, is the embedding for the (energetically disfavored) N z = ∞ solution. are also barely distinguishable from the shown curves). In other words, the flavor branes do not seem to care much about how many instanton layers they carry and where they sit.
As we have seen in Fig. 5, different choices of ρ 0 and γ 0 can lead to qualitatively different behaviors of the system. In principle, we can now scan the entire two-dimensional parameter space, and determine the phases and phase transitions for all chemical potentials. Keeping T = 0, this would result in a three-dimensional ρ 0 -γ 0 -µ phase diagram. We present two two-dimensional slices of this phase diagram in Fig. 7, where we have fixed ρ 0 and computed all phase transition lines in the γ 0 -µ plane. In the right panel, ρ 0 = 2.5, i.e., the results of Figs. 5 and 6 are obtained along two vertical lines in that panel. The left panel has been calculated with a smaller value, ρ 0 = 1.5. In appendix C we explain how we have computed the various phase transition lines and critical points.
Before we come to the observations, let us add a remark regarding the connection of the phase diagrams to the results of the previous section. The minimization carried out in Sec. 3.1 determines a trajectory of the system through the three-dimensional space spanned by µ, ρ 0 , γ 0 . This trajectory intersects each of the slices shown in Fig. 7 in a point. We find that for both slices this point lies in the region where N z = 2 (for the left slice this can be read off of the upper right panel of Fig. 4). One might thus naively conclude that we have found a stationary point of the free energy with N z > 1, even though we have argued in Sec. 3.1 that such a point does not exist. But this conclusion is not correct: there may very well be a minimum at N z = 2 under the constraint of a fixed pair (ρ 0 , γ 0 ), but no minimum for the same µ with N z = 2 if we search for the minimum in the entire parameter space, including ρ 0 and γ 0 .
The main observations of Fig. 7 are as follows. We see that the baryon onset can be of first order, as in Fig. 5, but also of second order, as in Sec. 3.1. It appears that smaller values of ρ 0 can produce a second order onset. This makes sense because by decreasing ρ 0 the instanton width becomes smaller, i.e., we approach the pointlike limit, and we know that the pointlike approximation predicts a second order onset. Small values of ρ 0 (like high densities and large values of γ 0 ) also seem to prefer instanton spreading. This suggest that, if we were to extend the pointlike approximation by allowing for instanton repulsion, we would presumably find the degenerate, delta-peaked instantons move up in the holographic direction with increasing chemical potential. This calculation might be of some interest because it would be the simplest system in which the instanton repulsion could be observed, possibly allowing for some analytic results, at least in certain limits such as large densities. On the other hand, our present results show that the pointlike approximation is not a good approximation for large densities and thus we do not include this calculation here. In the present scenario, the instanton width ρ = ρ 0 u c is always nonzero because u c never goes to zero. Thus, one way of thinking about the results is that in approach (i) of Sec. 3.1, the system chooses to have pointlike baryons at infinitesimally small densities, resulting in a second order onset, and here, in approach (ii), we forbid pointlike baryons via the external constraint ρ 0 > 0 and in accordance with expectations from finite-λ corrections, and thus we are able to see a first order onset.
For sufficiently small values of γ 0 , i.e., stretching the instanton in the spatial direction compared to the holographic direction, there is no baryonic phase at all: the vacuum is directly superseded by quark matter. For slightly larger values of γ 0 the baryon onset is followed by a transition to chirally restored matter (upon increasing µ at fixed ρ 0 and γ 0 ). The numerical results suggest that there is no chiral transition at all for sufficiently large γ 0 . However, the numerics become difficult at very large µ, and thus we cannot say this with certainty. In any case, except for a very narrow regime, the chiral phase transition occurs -if at all -at much larger values of the chemical potential than the baryon onset (note the logarithmic µ scale in Fig. 7).
It is interesting to note that the overall structure of the phase diagrams is very similar to the one obtained with the "homogeneous ansatz", with γ 0 replaced by the 't Hooft coupling λ, see Fig. 7 of Ref. [12] (by comparing the equations it is obvious that the instanton deformation γ in the present instantonic ansatz plays a very similar role to λ in the homogeneous ansatz). One difference is that the chiral phase transition line bends in the other direction: within the homogeneous ansatz the baryon onset is never followed by a chiral phase transition if λ is held fixed. In both cases, µ-dependent parameters [here ρ 0 (µ), γ 0 (µ), there λ(µ)] would allow for an equation of state that shows a firstorder baryon onset and a transition to quark matter at moderately large densities. From a purely phenomenological point of view, it might be tempting to search for such a suitable µ dependence: one might add ρ 0 and γ 0 to the three free parameters of the model and fit them to known properties of nuclear matter at the saturation density, or to the (poorly known) critical chemical potential of the chiral phase transition. However, fitting them in a density-dependent way would necessarily include some arbitrariness or, at best, some extrapolation to large densities. And, from a theoretical point of view, there is no reason for using the phase diagrams of Fig. 7 and moving through them with externally given functions ρ 0 (µ), γ 0 (µ). We do know how ρ and γ "want" to behave as a function of µ in the given approximation. This was discussed in Sec. 3.1 and has led to unphysical results, a second-order baryon onset and no chiral restoration. Therefore, the results of the present section should be understood as a step towards a better understanding of the instanton approach to baryonic matter, its relation to the homogeneous ansatz and, in future work, towards further improvement of the approximation, rather than a straightforward recipe for constructing a strong-coupling equation of state for dense matter inside compact stars.

Summary and outlook
We have investigated homogeneous baryonic matter at zero temperature in the decompactified limit of the Sakai-Sugimoto model. Our main point was to improve existing approximations based on flat-space instantons and to ask whether these improvements bring us closer to real-world nuclear matter, in particular whether they give rise to a first-order baryon onset and a chiral transition to quark matter at high densities. Motivated by results of holographic baryons in the vacuum, we have introduced a deformation parameter into the instanton ansatz that allows for an anisotropy in the space of holographic and spatial directions. While in the spatial direction we have employed an averaging procedure, accounting for homogeneous matter in a very simple way, we have introduced instanton repulsion in the bulk: the instantons are allowed to spread out in the holographic direction in the form of a number of instanton layers, this number and the distance between the layers being determined dynamically.
We have found that if we minimize the free energy with respect to the width and the deformation of the instantons (and with respect to the various other parameters of our ansatz) that (1) there is a (unphysical) second-order baryon onset, and at the onset the instantons are pointlike, reproducing the approximation of degenerate, delta-peaked instantons, (2) baryons (unphysically) refuse to go away at large densities, i.e., there is no chiral restoration, (3) at large densities, the instantons tend to get elongated along the holographic direction, and (4) the instantons prefer to sit all at the same point in the bulk, i.e., the number of instanton layers remains 1 for all densities.
Besides this most straightforward approach we have also worked with external constraints on the width and the deformation of single instantons -being aware that our simple approximation cannot capture the shape of the full solution -and studied their effect on the many-instanton system. More precisely, we have constrained the width ρ to be proportional to the (density-dependent) location of the tip of the connected flavor branes u c and the deformation parameter γ to be proportional to u 3/2 c , and treated the proportionality constants ρ 0 and γ 0 as free parameters. Interestingly, with these constraints, which in particular result in non-pointlike instantons at all densities, we do find that at sufficiently large densities the instantons are divided into two layers. With increasing density these layers move up in the holographic direction, away from the tip of the connected flavor branes. All higher numbers of layers turn out to be energetically disfavored. In particular, we have included the possibility of infinitely many layers, corresponding to a smeared distribution along the holographic direction. Our results regarding the instanton layers is interesting in view of various other approximations and approaches in the same model, and complete solutions in simplified models, which show similar effects [19,24,29,31,53,54]. We have also found that in the presence of the constraints on width and deformation, we do see a first-order baryon onset and chiral restoration at large densities, as expected from QCD: for sufficiently large ρ 0 , the baryon onset is first order, and for a small range of values of γ 0 , baryonic matter is superseded by chirally restored quark matter, with the critical chemical potential for chiral restoration being very sensitive on γ 0 .
In conclusion, we have shown that the present instanton approximation -if evaluated at its stationary point -is too simplistic to show multiple instanton layers (which are suggested by other approximations of the model), and it does neither show a first-order baryon onset nor chiral restoration (which occur in the real world). We have shown that by imposing external constraints on the shape of the single instantons, these features do appear.
For a short discussion of future perspectives we first recall that this work has been, to a large extent, motivated by a phenomenological question: can we come up with a strongcoupling model description of dense nuclear and quark matter within a single model? In principle, we could use our results to find such a model description. However, to fulfill even the most fundamental requirements of real-world matter, this would require a theoretically ill-motivated, density-dependent choice of ρ 0 and γ 0 , which renders any resulting predictions questionable. Therefore, the results should mainly be considered as a further theoretical step towards such a model. We have embedded our ansatz into a more general setup, which allows for systematic improvements. For instance, it would be very useful to extend our ansatz to one that incorporates a nontrivial dependence on the 't Hooft coupling λ, even though a systematic treatment of finite-λ effects would require string corrections, which is very difficult. It would also be interesting to gain a deeper understanding of the relation between our instanton approach and the results of the "homogeneous ansatz" [12,19] and of the relation of both to the full solution. For the latter it is useful to retreat to simpler models or simplifications of the Sakai-Sugimoto model, where the full solution is available and can be compared to the various approximations [29,53]. Other promising extensions, within the present ansatz or one that is further improved, are to include nonzero temperatures (all our equations contain the full temperature dependence, but we have restricted ourselves to zero temperature in the numerical results), to include an isospin chemical potential, the possibility of a chiral density wave, or to add an external magnetic field.
for the embedding function for the flavor branes. As a consequence, we are able to solve the equations of motion for all gauge fields analytically, including the dependence on the position space coordinates. The analogous calculation in the confined geometry can be found in the literature, for maximally separated flavor branes [15] and with a general, not necessarily maximal, separation [16]. The present calculation in the deconfined geometry in particular yields a nontrivial temperature dependence for the instanton width and the instanton deformation.
In the main part, we have ignored any dependence of the abelian gauge fieldÂ 0 on X, and thus we first have to reinstateF 0i into Eq. (2.8) [or Eq. (2.9), which only differs from Eq. (2.8) by terms of order F 4 , which shall be neglected in this appendix] and insert the result into the DBI action (2.2). Then, we derive the YM action by expanding in the field strengths up to order F 2 , which results in the action with λ 0 ≡ λ/(4π), is a purely geometric term, not depending on any gauge fields, and where Implicitly, we have introduced the dimensionless gauge fields a i = A i M KK , a u = A U R(M KK R) 2 , and, as in the main text,â 0 =Â 0 λ 0 M KK . In the absence of instantons, the equation of motion for x 4 (u) is given solely by S 0 . With the boundary condition x 4 (u c ) = ∞ we find .

(A.4)
This is the leading order result for small instantons. (There are subleading contributions which we ignore, i.e., we work without backreactions of the instanton on the embedding of the flavor branes.) For small instanton widths, the integrands in S YM and S CS are nonzero only in a small vicinity around z = x = 0. This renders the abelian termsF 2 0i ,F 2 0z of higher order than the non-abelian ones F 2 ij , F 2 iz . Anticipating the eventual solution, one can do this systematically by rescaling x → x/ √ λ, z → z/ √ λ, and the gauge fields accordingly, and applying a systematic expansion for large λ. To keep the notation simple, we do not introduce rescaled quantities, but keep this expansion in mind, which yields leading and subleading contributions to the energy of the instanton that are eventually of first and zeroth order in λ, To compute these contributions, we first insert (A.4) into the YM action and change the integration variable from u to z. Then, S YM is obtained by an expansion around z = 0 and dropping the abelian field strengths and all higher order terms O(z 2 ), Consequently, to leading order, the equations of motion for the non-abelian gauge fields yield the flat-space BPST solutions, with field strengths equivalent to the ansatz (2.16) in the main text. In order to compute the abelian field strengths, we go to subleading order, The equation of motion forâ 0 becomes [recall that we work in Euclidean space, wherê where the first term in the square brackets on the right-hand side comes from S Minimizing E with respect to ρ yields This solution justifies our expansion a posteriori since we now confirm that S [to be compared to the Lagrangian from the unsymmetrized prescription (2.37)], where, for the sake of a compact notation in this appendix, we have abbreviated This completes the set of equations needed to find the ground state for the case of the symmetrized trace prescription. The equations are very similar to the simpler ones of the unsymmetrized case. The main difference are the cubic equations whose solutions enter the integrands of all stationarity equations, making the numerical evaluation more involved.

C Calculation of phase transition lines and critical endpoints
In this appendix we explain our calculation of the various phase transition lines for the phase diagrams in Fig. 7. It is obviously very tedious to compute the ground state on a grid in the ρ 0 -γ 0 -µ space and deduce the phase transition lines from this calculation. More efficiently, after getting an idea of the overall structure of the phase diagram, one may proceed as follows.
• Chiral phase transition. This is the first-order phase transition between the chirally symmetric phase and the baryonic phase, i.e., the phase transition line is defined by Ω baryon = Ω quark . We solve this equation simultaneously with Eqs. (2.47b), (2.47e), and (2.47f) for the variables k, n I , z 0 , and µ [all rescaled with appropriate powers of u c , such that Eq. (2.47a) decouples]. The baryonic phase at this transition can have either z 0 = 0 (N z = 1) or z 0 > 0 (N z = 2). The transition with z 0 = 0 can be computed separately without the minimization with respect to z 0 (2.47e), and the intersection of the two phase transition lines defines a critical point, see upper left corner of both panels in Fig. 7.
• Onset of second instanton layer. This is the transition within the baryonic phase which separates z 0 = 0 from z 0 > 0 (z 0 > 0 only appears for N z = 2, no higher number of instanton layers is preferred). This transition can be either second or first order. The second order transition line is determined as explained in the context of Fig. 3: after dividing the minimization with respect to z 0 (2.47e) by z 0 , we take the limit z 0 → 0, i.e., we approach the transition form the z 0 > 0 side. The resulting equation is solved simultaneously with Eq. (2.47f) in the limit z 0 → 0 for the variables n I and k, which are then used to compute the critical chemical potential with the help of Eq. (2.47b). The calculation of the first-order onset of the second layer is a little more complicated because we have to compare the free energies of the phases with one and two layers. We need to solve the following 7 equations simultaneously: Eqs.
(2.47b) and (2.47f) for z 0 = 0, Eqs. (2.47b), (2.47e), and (2.47f) for z 0 > 0, plus the two conditions that the chemical potentials and the free energies of both phases are identical. The corresponding 7 variables are k, n I , µ for both phases and z 0 for the z 0 > 0 phase (as always we work with quantities rescaled by u c , so, more precisely, the two chemical potentials used as variables in our calculation areμ 1 = µ 1 /u c,1 , µ 2 = µ 2 /u c,2 , and at the phase transition µ 1 = µ 2 , while in generalμ 1 =μ 2 ). The critical point that separates the second-order from the first-order transition line can be found by asking at which γ 0 the curve z 0 (µ) becomes multivalued, which is equivalent to asking at which γ 0 this curve bends to the left at z 0 → 0.
• Baryon onset. This is the phase transition between the mesonic and baryonic phases. Therefore, it is defined by Ω baryon = Ω meson , and this equation has to be solved simultaneously with Eqs. (2.47b), (2.47e), and (2.47f) for the variables k, n I , z 0 , and µ. This is completely analogous to the chiral phase transition. Again, we need to compute the two cases z 0 = 0 and z 0 > 0. The former transition can be either first order or second order, and the second order transition can be found by solving the single equation (2.47f) for k with an infinitesimally small n I → 0.