Dimension Changing Phase Transitions in Instanton Crystals

We investigate lattices of instantons and the dimension-changing transitions between them. Our ultimate goal is the 3D->4D transition, which is holographically dual to the phase transition between the baryonic and the quarkyonic phases of cold nuclear matter. However, in this paper (just as in [1]) we focus on lower dimensions -- the 1D lattice of instantons in a harmonic potential V M_2^2x_2^2+M_3^2x_2^2+M_4^2x_4^2 and the zigzag-shaped lattice as a first stage of the 1D->2D transition. We prove that in the low- and moderate-density regimes, interactions between the instantons are dominated by two-body forces. This drastically simplifies finding the ground state of the instantons' orientations, so we made a numeric scan of the whole orientation space instead of assuming any particular ansatz. We find that depending on the M_2/M_3/M_4 ratios, the ground state of instanton orientations can follow a wide variety of patterns. For the straight 1D lattices, we found orientations periodically running over elements of a Z_2, Klein, prismatic, or dihedral subgroup of the SU(2)/Z_2, as well as irrational but link-periodic patterns. For the zigzag-shaped lattices, we detected 4 distinct orientation phases -- the anti-ferromagnet, another abelian phase, and two non-abelian phases. Allowing the zigzag amplitude to vary as a function of increasing compression force, we obtained the phase diagram for the straight and zigzag-shaped lattices in the (force, M_3/M_4), (chemical potential, M_3/M_4), and (density, M_3/M_4) planes. Some of the transitions between these phases are second-order while others are first-order. Our techniques can be applied to other types of non-abelian crystals.


Introduction
In the ordinary N c = 3 QCD cold nuclear matter forms a quantum liquid, but for large N c it becomes a crystalline solid [2]. In many holographic models of QCD [3,4] baryons are represented by instantons of the U(N f ) gauge theory living on the flavor branes [5], so cold nuclear matter corresponds to a whole crystalline lattice of instantons [6]. The geometry and even the dimensionality of this lattice depend on the baryon density: at medium densities, the instantons form a 3D lattice in the x 4 = 0 hyperplane of the holographic 4D space (not counting the time); this corresponds to the baryonic phase of nuclear matter. At high densities, the baryons spread out into the x 4 dimension and form a 4D lattice; this is dual to the quark liquid phase of nuclear matter, or rather to the high-density quarkyonic phase [14] in which the quarks fill the Fermi sea, but the excitations near the Fermi surface are baryon-like bound states of N c quarks rather than free quarks by themselves. According to Rozali et al [7], the chemical potential of the quarks in this phase is related to the thickness of the lattice in the x 4 direction.
Most investigations of the low-temperature holographic nuclear matter focus on the extreme high-density limit where the x 4 dimension is thick and can be studied macroscopically (for example, see [9,8,7,10]). But our main interest is in the microscopic structures of the instanton lattices, and especially the phase transition between the 3D and 4D lattices.
In our previous paper [1] (with Dmitry Melnikov) we found that there is a whole sequence of such transitions: from a single 3D layer in the x 4 = 0 hyperplane, to two 3D layers in two parallel hyperplanes, to 3 layers, to 4 layers, etc., until the number of layers becomes too large to count individually. Alas, the 3D and 4D instanton lattices turned out to be too hard to analyze, so we resorted to oversimplified toy models. In the first toy model, we approximated the instantons as point-like charges repelling each other with 4D Coulomb forces; this approximation is very crude -it ignores the instantons' orientations, never mind the interference between the instantons -but it makes for a simple model of transitions between lattices of different dimensions. Indeed, in this model we found that increasing density makes a 3D lattice of instantons in x 4 = 0 hyperplane suddenly split into two sublattices in the x 4 = ±ǫ hyperplanes; pictorially, this looked like the popcorn suddenly jumping up in the popper, so we dubbed the dimension-changing transitions the popcorn transitions.
In our second toy model, we used actual instantons and the exact ADHM [11] construction for the multi-instanton system. However, we had confined the instanton centers to a 2D plane x 1 = x 2 = 0 by making the inverse 5D gauge coupling rise steeply in the x 1 and x 2 directions. We also had the 1/g 2 5 slowly rising in the x 3 direction, altogether Consequently, at low densities the instantons formed a 1D lattice in the x 4 dimension (which acted as the model's only flat space dimension), but for higher densities the instantons moved into the x 3 dimension (which acted as the holographic dimension of the model). In that model, we saw two phase transitions in response to increasing instanton density. The first transition -from a straight 1D chain of instantons to a zigzag-shaped chain -is second-order, so the zigzag amplitude behaves as ǫ ∝ √ ρ − ρ c . On both sides of this transition, the instantons' orientations form an anti-ferromagnetic pattern: two alternating orientations related by a 180 • rotation. At higher densities there is another transition, but this time it's first-order: the zigzag amplitude ǫ jumps discontinuously, and the instantons' orientations also change their pattern: the orientations of the nearest neighbors are now related by a rotation through angle φ < 110 • instead of 180 • .
Presumably, at still higher densities we would have seen popcorn transitions from the zigzag -which is a kind of two-layer lattice -to a 3-layer lattice, then to 4 layers, etc., etc. However, calculating the net energies of such multi-layer lattices turned out to be too hard, so we stopped at the zigzag.
Technically, the main difficulty in working with more complicated instanton lattices is setting up the ADHM construction [11] for the infinite number of instantons. Indeed, the ADHM construction for A instantons involves 4 non-commuting A × A matrices Γ µ mn whose off-diagonal elements follow from non-trivial constraint equations; for A → ∞ solving these equations becomes very hard. For a simple 1D lattice of instantons whose orientations are related by commuting SU(2) symmetries, the ADHM constraints were solved by Kraan and Baal [12], and we have re-derived and used their explicit formula for the instanton number density profile I(x) in [1]. For a zigzag-shaped instanton chain with link-periodic instanton orientations, we found an exact solution of the ADHM constraints in terms of Fourier transforms of some rather messy functions, but we could not use it to derive an exact formula for the I(x); instead, we had to use the small-instanton-size approximation. And it gets worse for more complicated instanton lattices and/or orientation patterns: For the simple 2D and 3D lattices and purely abelian orientation patters the ADHM constraints have known formal solutions 1 -which alas are way too complicated for any practical use -while for the non-abelian orientation patterns there are no known solutions at all. This is particularly unfortunate since the lowest-energy 3D configuration for low densities is likely to be an FCC lattice with non-commuting orientations of the instantons [13].
In this paper we follow a shortcut around the ADHM construction -and also around the more difficult part of the energy calculation. We assume from the beginning that the instanton density is not too high, so the distances between instantons are much larger than the instantons' radii. In §3 we show that for all such multi-instanton systems, interactions between the instantons are dominated by the two-body forces. The irreducibly-threebody forces, 4-body forces, etc. also exist, but they are suppressed by powers of the small (radius/distance) 2 ratio compared to the two-body forces. Consequently, we obtain a manageably simple approximate formula for the net energy of a multi-instanton system in terms of the instantons' positions and orientations.
Thanks to this formula, we do not need to assume (as we did in [1]) that in a lattice of instantons, their relative orientations have the same symmetries as their relative positions. Instead, we may use numerical methods to minimize the net energy with respect to all the orientations treated as independent parameters, without any symmetry assumptions. The algorithm we use in this paper is fairly simple: We set up a large lattice (200 sites is more than enough for 1 dimension) of fixed geometry, with a free SU(2) matrix y n at each site encoding the n th instanton's orientation (for N f = 2). We start with a completely random set of the initial y n and make them evolve along the steepest descent of the net energy function until they converge to a local minimum. Then we repeat the process many times with different random sets of initial y n . This gives us a short list of the local minima, from which we select the deepest -which is presumably the global minimum of the net energy function. This algorithm has convergence problems near the phase boundaries -where two or more local minima are almost degenerate -but it is good at finding all the orientation patterns that might be ground states for some combinations of background parameters and lattice geometries. And once we know all these patterns, we fit them all with common ansatz, calculate the net energy analytically as a function of the ansatz's parameters, and then we minimize it with a much higher precision to get the accurate phase diagram of the instanton lattice.
In this paper we put the instantons in a harmonic potential that rises at different rates along three out of four space coordinates; in terms of the 5D gauge coupling, Specifically, we shall see that for the 1D instanton lattice in the almost-holographic background M 2 , M 3 ≪ M 4 (where M 2 and M 3 help to line up the instanton centers along the x 1 axis but do not affect the two-body forces between the instantons), there is a large degenerate family of ground states for the instanton orientations y n ∈ SU(2)/Z 2 (for N f = 2).
This family includes many periodic patterns in which all the y n belong to a finite subgroup of the SU(2)/Z 2 : the anti-ferromagnetic pattern in which y odd n and y even n belong to 2 elements of a Z 2 subgroup, the period = 4 pattern in which the y n span the Klein group Z 2 × Z 2 , and two series of period = 2k patterns in which the orientations span a prismatic group Z k × Z 2 or a dihedral group D 2k . There are also many link-periodic patterns in which the relative orientations y † n y n+1 of neighboring instantons are periodic but the orientations y n themselves are not periodic. But the vast majority of the degenerate ground states are not periodic at all.
Increasing the M 2 and M 3 parameters lifts the degeneracy. For M 2 ≪ M 3 ∼ M 4 , the ground state of the 1D lattice is anti-ferromagnetic, while for M 2 ∼ M 3 ∼ M 4 the ground state is periodic with a longer period or link-periodic, depending on the M 2 /M 3 ratio. Specifically, for M 2 = M 3 the instanton's orientation in the ground state span the Klein group, for other rational M 2 2 /M 2 3 they span a dihedral group, while the irrational M 2 2 /M 2 3 favor the link-periodic patterns.
For the zigzag-shaped chains of instantons we need M 2 ≪ M 3 ∼ M 4 to make sure the transition from a straight chain to a zigzag happens when lattice spacing ≫ instanton radius -otherwise the two-body-force approximation would not work -but we may vary the M 3 /M 4 . In this setting, the ground state of instanton orientations may follow one of four different patterns, depending on the zigzag amplitude and on the M 3 /M 4 ratio. Besides the anti-ferromagnetic and the abelian φ < π patterns seen in [1], there are two non-abelian link-periodic pattern in which the y † n y n+1 "twists" for even and odd links of the lattice do not commute with each other. The non-abelian patterns appear only for M 3 < M 4 , that's why we did not see them in [1].
Focusing on the "popcorn transition" from a straight chain of instantons to a zigzagshaped chain (which is a two-layer lattice), we shall see that the thermodynamic order of this transition depends on the M 3 /M 4 ratio: the transition is first-order for M 3 < 0.725M 4 but second order for M 3 > 0.725M 4 . Likewise, when we further increase the instanton density -or rather, the 1D compression force which maintains it -the number and the orders of the subsequent phase transitions between different orientation phases of the zigzag-shaped chain also depend on the M 3 /M 4 ratio. This paper is organized as follows: For the reader's convenience, in §2 we review the basic features of holographic nuclear matter: effects of the large N c and large λ limits, the generalized Sakai-Sugimoto model of holographic QCD, realization of holographic baryons as instantons of the flavor gauge symmetry, basic features of the holographic nuclear forces, and our basic approach to the multi-baryon systems and the bulk nuclear matter. In §3, we argue that the interactions between multiple instantons include 2-body forces as well as irreducibly-3-body force, 4-body forces, etc., but when the distances between the baryons are much larger that the instantons' radii, the two-body forces dominate the net energy, while the 3-body,. . . , forces are suppressed by powers of (radius/distance) 2 ≪ 1. First we give a heuristic argument, and then we back it up with a technical proof for N f = 2; we also calculate the precise form of the net two-body force and its dependence on the two instantons' orientations.
In §4 we use the two-body-force approximation to study 1D lattices of instantons. We start with the M 2 , M 3 ≪ M 4 background and show that there is a very large degenerate family of ground states of the instantons' orientations. Then we move on to the M 2 ∼ M 3 ∼ M 4 backgrounds; this changes the precise form of the two-body force in a way that lifts the degeneracy between the orientation states. Consequently, the ground state of the lattice has a specific periodic or link-periodic pattern of orientations, depending on the M 2 /M 3 ratio.
In §5 we study zigzag-shaped instanton chains in the M 2 ≪ M 3 ∼ M 4 backgrounds. After a few preliminaries, we use numerical methods to find the lowest-energy state of instanton orientations without presuming any particular symmetries. We repeat the calculation for 50 different lattice geometries (ǫ/D = 0.01, . . . 0.99 by 0.02) and 50 values of the M 3 /M 4 ratio; the results are presented on the phase diagram 5 on page 63. Our numeric code is not very accurate and has convergence problems near the phase boundaries, but it is good at identifying the symmetries of the lowest-energy patterns. Using these symmetries, we construct a 3-parameter ansatz that includes all the preferred patterns, then calculate the net energy as analytic function of the parameters and minimize it. For the reader's convenience, the energy calculation and minimization is described in the Appendix, while the resulting phase diagram is presented in figure 6 (page 67) in the main body of §5.
In the last part of §5 we treat the zigzag amplitude and the lattice spacing as dynamic parameters, vary the compression force of the 1D instanton lattice, and study the popcorn transitions from the straight chain to the zigzag and between different orientations phases of the zigzag. We show that the net number of such transitions -1, 2, or 3 -and their orders depend on the M 3 /M 4 ratio. The overall phase diagram in the compression versus Finally, in section 6 we summarize our results and give a short list of open questions we hope to address in the near future.

Review of Holographic Nuclear Physics
In this section we review the basics of holographic nuclear physics that we shall need later in this paper. The holography here means gauge-gravity duality, and the semiclassical description of the gravity side of the duality requires the limits of large N c and large 't Hooft coupling λ = g 2 YM × N c [15]. In the first two subsections §2.1-2 we discuss the general consequences of these limits for the nuclear physics. In §2.3 we get more specific and review the Sakai-Sugimoto model of holographic QCD. There are many other models which differ in countless details -and none of them is exactly dual to the QCD, -but their basic features are rather similar, and the Sakai-Sugimoto model serves as a prototype. In §2.4 we introduce the baryons and explain their holographic duality to the instantons of the flavor gauge symmetry. We also review the interactions and the self-interactions of such instantons in the large λ limit and explain why they have small radii a ∝ λ −1/2 . Finally, in §2.5 we address the multi-instanton systems dual to nuclear matter and explain how the interactions between the baryons appear in at the first order in the 1/λ expansion. We also explain the technical difficulties in handling infinite crystals of baryons and preview a shortcut around those difficulties that we shall take in §3.
2.1 Nuclear matter in the large N c limit QCD perturbation theory in the large N c limit is dominated by planar gluon diagrams while contributions of the non-planar diagrams and of the quark loops are suppressed. At zero temperature and chemical potential below the mass of the baryon, the large N c QCD is a theory of free mesons and glueballs whose interactions are suppressed by powers of N c (except for interactions with the baryons). The masses of mesons and glueballs are O(Λ QCD ), and they scale as N 0 c . On the other hand, the baryons -which are made out of N c quarks -have masses of order N c × Λ QCD , so their relative abundance in thermal equilibrium is exponentially suppressed. The interaction energy between baryons also scales as N c .
The T − µ phase diagram of the large-N c QCD is believed to look like figure 1 below.
(For a review see [14] and references therein, also [22].) We assume the 't Hooft's N c → ∞   [14] limit in which the number of flavors remains finite. Consequently, the dynamics of the theory is dominated by the gluons, the quarks are sensitive to the gluonic background, but the backreaction from the quarks to the gluons is suppressed by N f /N c ≪ 1. Thus, at lower temperatures there is confinement but for increasing temperature there is a first order transition to the deconfined phase. The transition temperature T d is at the Λ QCD scale and it's almost independent of the quark chemical potential µ q = µ b /N c (as long as µ q does not grow with the N c ).
It is not clear whether for N c → ∞ the deconfining phase transition coincides with the chiral symmetry restoration for the light quarks. Several field theory arguments -for example [16,14] -suggest that for µ q = 0 the two transitions should happen at the same point. However, these arguments do not work for µ q > 0 [14], and there are other arguments for the existence of confined but chirally restored phases and/or deconfined phases where the chiral symmetry remains broken. Moreover, some holographic models -for example [4] -have deconfined but chirally broken phases even at zero chemical potentials.
For temperatures below the deconfining transition t d and baryon chemical potentials below the baryon mass -or equivalently for µ q m q ≡ M b /N c -the thermal state of the theory is a gas of glueballs and mesons with almost no baryons or antibaryons. But at µ q ≈ m q there is an abrupt phase transition to the bulk nuclear matter with finite baryon density. But unlike the ordinary nuclear matter which is in the quantum liquid state, the large N c nuclear matter is crystalline solid since the ratio of kinetic energy to potential energy decreases with N c . Indeed, the potential energy of baryon-baryon forces scales like N c ; more precisely [17], in the large N c limit the two-baryon potential becomes for some N c -independent profiles A C (r), A S (r), and A T (r) for the central, spin-spin, and tensor forces; their overall magnitudes are A ∼ Λ QCD for r ∼ 1/Λ QCD . Classically, this potential tries to organize the baryons into some kind of a crystal where the distances between neighboring baryons do not depend on the N c while the binding energy (per baryon) scales like N c × O(Λ QCD ). In quantum mechanics, the baryons in such a crystal behave like atoms in ordinary crystals -they oscillate in their potential wells with zero-point kinetic energies where d ∼ 1/Λ QCD is the N c -independent diameter of the potential well. Therefore, at zero temperature the ratio of kinetic energy to the potential energy scales like and becomes very small for large N c . At higher temperatures the kinetic energies of baryons become larger, K ∼ T , but in the confined phase we are limited to T < T d ∼ Λ QCD and hence K < O(Λ QCD ). Consequently, the kinetic to potential energy ratio scales like -which is larger than (2.3) but still becomes small in the large N c limit. Consequently, for large N c neither zero-point quantum motion nor thermal motion of baryons can melt the baryon crystal, so nuclear matter remains solid all the way to the deconfining temperature.
Before holography, the best models for the large-N c nuclear crystals were lattices of skyrmions. In this framework, Igor Klebanov [18] had found a curious phase transition from a lattice of whole skyrmions at low chemical potentials to a denser lattice of half-skyrmions at higher potentials. (According to [19], the whole-skyrmion lattice at low potentials has FCC geometry while the half-skyrmion lattice at higher potentials is simple cubic.) In the half-skyrmion-lattice phase, the order parameter for the chiral symmetry breaking vanishes after space averaging, so in QCD terms the transition is interpreted as chiral symmetry restoration at high µ q .
In QCD with N c = 3 and two massless flavors there is a similar chirally-symmetric phase at low T and high µ q . This phase is a quark liquid rather than a baryon liquid (the quarks are no longer confined to individual baryons) and there is a condensate of quark pairs making this liquid a color superconductor. But for large N c the situation is more complicated: there is no color superconductivity, and there is no deconfinement for T < T d . Instead, the dense cold nuclear matter forms a phase which combines the features of the baryonic and quark phases: the quarks fill up a Fermi sea, but the interactions near the Fermi surface are strong, so the excitations are not free quarks or holes but rather meson-like quark-hole pairs or baryon-like states of N c quarks. MacLerran et al have dubbed this phase "quarkyonic".
For µ q ≫ Λ QCD , the interior of the Fermi see is chirally symmetric, but near the Fermi surface the symmetry is broken by the chiral density waves [20]. (Although Son and Shuster [28] argue that such waves develop only for very large N c > 10 3 .) To be precise, the chiral density waves mix the chiral symmetry of the quarkyonic phase with the translational symmetry rather than simply break it. Averaging over space restores the chiral symmetry, just like it happens for the lattice of half-skyrmions.
On the other hand, for µ q just above m q = M b /N c , the baryonic crystal has a completely broken chiral symmetry. Thus, at some critical µ (c) q = m q × a few, there should be a chiral symmetry restoring phase transition from the baryonic crystal to a distinct quarkyonic phase.
In holography, this transition is believed to be dual to the "popcorn transition" from a 3D instanton lattice to a 4D lattice.

Effect of the large λ limit on the holographic nuclear matter
In holography, the semiclassical description of the gravity side of the gauge-gravity duality in terms of metric, fluxes, and branes requires the limits of large N c and also large 't Hooft coupling λ = N c g 2 YM . In the large λ limit, the baryons become very heavy: in units of the mesonic mass scale M ∼ Λ QCD , the baryon mass is M b ∼ λN c M. However, the interactions between the baryons do not grow with λ: even for two baryons right on top of each other, the repulsive potential between them is only V ∼ N c M ∼ M b /λ. At larger distances, the forces are even weaker since the hard-core radius of a holographic baryon shrinks with λ as R b ∼ M −1 λ −1/2 . Outside this radius, the repulsive potential decreases as 1/r 2 until r ∼ M −1 , at which point it becomes dependent on the meson mass spectrum of a specific holographic model: In some models, the potential becomes attractive for r M −1 while in others it remains repulsive at all distances. The overall picture is shown in figure 2: Since the nuclear forces are so weak in the holographic QCD, all transitions between different phases of the cold nuclear matter happen at chemical potentials µ b very close to the baryon mass: just below M b we have glueball/meson gas (or vacuum for T = 0), while just above M b we have dense quark matter. To see the baryonic matter phase (or any other intermediate phases) we need to zoom into the µ b ≈ M b region. Figure 3 (on the next page) illustrates this point: to see the baryonic-matter phase between the vacuum and the quark-matter phases on the plot of baryon density ρ as a function of the chemical potential µ, we need to zoom into the narrow range µ − M b = O(M b /λ). The figure also shows that the thermodynamic order of the phase transition between the vacuum (or the meson/glueball gas for T > 0) and the baryonic matter depends on the sign of the long-distances nuclear force. If the force becomes attractive at long distances, then bulk baryonic matter exists at zero external pressure and has µ = M b −binding energy. Consequently, the transition from the vacuum (or gas) to the nuclear matter (in the form of a baryonic crystal) is first-order as shown on figure 3(a).  Figure 3: Density as a function of the chemical potential in the large λ regime for attracting (a) and repelling (b) baryons. The nuclear matter phase is confined to a narrow window of the order ∆µ ∼ B B /λ. In the naive diagram the transition occurs directly from the no-baryon to a quark phase.
On the other hand, if the nuclear forces are repulsive at all distances, then the bulk nuclear matter does not exists except at positive external pressures, and its chemical potential must be µ > M b . Moreover, µ raises monotonically with the pressure and the density, so the transition from the vacuum to the bulk nuclear matter is second-order as shown on figure 3(b).
In this article we do not wish to focus on a particular model of holographic QCD, so we make no assumptions about the long-distance nuclear forces. Consequently, we cannot say anything specific about the transition from the vacuum to the baryonic phase. Instead, we focus on the transition from the baryonic phase to the quarkyonic phase, which correspond in the holographic picture to changing the lattice geometry of the instanton crystal, from a 3D lattice, through a sequence of intermediate steps, to a 4D lattice. Or rather, that's the goal of our program; in this article we shall focus on a simplified problem, namely the transitions between 1D and 2D instanton lattices.

Sakai-Sugimoto model as a prototype of holographic QCD
Some gauge theories -especially the N = 4 SYM -have exact holographic duals, where both sides of the duality follow as IR limits of the same string-theoretical construction, and all the extra degrees of freedom are superheavy. Alas, the Quantum Chromo Dynamics -or even the large N c limit of QCD -is not so lucky: it either does not have an exact holographic limit, or we have not found it yet. Instead, there is a large number of "holographic QCD" models, which are dual not to the QCD as such but rather to some QCD-like theories with a lot of extra junk that the real QCD does not have. Such models make for "spherical horse in a vacuum" approximations -good enough to get qualitative understanding, but not too accurate to make useful numerical predictions.
In a moment, we shall review the Sakai-Sugimoto model [3] of holographic QCD. That model should not be taken too seriously -its predictions for the hadronic spectrum are at best so-so, and for the nuclear forces are even worse, -but it serves as a type specimen for a large class of models. Here are some general features of this class: • The construction starts with N c coincident D-branes, which span the Minkowski space× a compact cycle. Typically, one uses D3 branes at a singularity, or D4-branes wrapping a circle, although some constructions use D5 branes on a 2-cycle, etc. For weak 't Hooft couplings λ = N c × g YM ∼ N c × g str ≪ 1, the open strings between the branes give rise to the gluons of the U(N c ) gauge theory (plus a lot of extra junk).
• For λ ≫ 1, the D-branes merge into a big fat black brane. All we see outside the horizon is a warped space-time geometry and the Ramond-Ramond fluxes induced by the conserved charges of the D-branes. Modes of the bulk-field excitations in this geometry are dual to the QCD glueballs. Or rather, some modes are dual to the glueballs while other models are dual to the extra junk that the real QCD does not have • To make the hadrons (mesons and baryons), add N f flavor D-branes, usually D5, D6, D7, or D8. At weak λ, the open strings connecting the color branes to the flavor branes give rise to the quarks and the anti-quarks.
• In the holographic limit N c → ∞ and λ → ∞, but the N f stays finite and the N f × g str remains weak. Consequently, the flavor branes remains D-branes (rather than merge into a black brane of their own), but now they live in a warped geometry with fluxes induced by the N c color branes.
• The open strings between the flavor branes give rise to N 2 f vector and scalar fields living on those branes. The 4D modes of these vector and scalar fields are dual to the QCD mesons.
• The YM instantons of the vector fields on the flavor branes are dual to the QCD baryons, see §2.4 for more details.
And now let us review the Sakai-Sugimoto model [3] in a little more detail. The color sector of the model is the Witten's model [23]: N c D4 branes spanning the Minkowski space×a circle; the circle has antiperiodic boundary conditions for the fermions, which breaks the N = 4 SUSY down to N = 0 * . In the field theory limit λ ≪ 1 =⇒ Λ QCD ≪ M KK the effective low-energy theory is pure U(N c ) Yang-Mills, but in the holographic limit λ ≫ 1 everything happens right at the Kaluza-Klein scale M KK = 1/circle radius, so the YM glueballs end up with similar O(M KK ) masses to a lot of non-YM junk.
On the gravity side of the duality, the D4 branes merge into a black brane which warps the 10D metric. Instead of flat 10D spacetime, we now have a warped product of R 3,1 Minkowski space, the S 4 sphere (originally surrounding the D4 branes), and a 2D cigar spanning the radial direction (⊥ to the branes) and the S 1 circle. The radial coordinate u runs from u Λ > 0 to infinity; at u Λ the S 1 circle shrinks to a point, hence the cigar. Altogether, we have warped metric, the 4-form flux, and the dilaton's gradient according to where x 4 was the coordinate along the S 1 circle and now is the polar angle on the cigar, (2.6) ℓ s = √ α ′ is the string length scale, and g s is the string coupling. The u Λ -the minimal value of the radial coordinate u at the tip of the cigar -is related to the original radius R of the S 1 circle as The same radius R also controls the 4D Yang-Mills coupling and hence the 't Hooft's coupling λ. Analytically continuing from λ ≪ 1 to λ ≫ 1, we have To add the flavor degrees of freedom to the model, Sakai and Sugimoto added N f D8 branes and N f anti-D8 branes. On the field-theory side of the holographic duality (i. e., for λ ≪ 1), they span all space coordinates except the x 4 coordinate along the S 1 circle. At the intersections of D8 and D4 branes, the open strings give rise to the massless quarks; likewise, at the intersections of the anti-D8 branes and the D4 branes we get massless antiquarks.
On the holographic side λ ≫ 1 of the duality, the exact solution for the flavor branes interacting with the warped metric and fluxes is not known, but for N f ≪ N c and N f × g s ≪ 1 we may use the probe approximation: The flavor branes seek the lowest-action configuration in the background metric (2.5), while their back-reaction upon the metric is neglected. Consequently, at low temperatures (below the deconfinement transition) 2 , the flavor branes span the Minkowski space × S 4 × a line on the cigar; the exact shape of this line follows from minimizing the branes' action, but its topology follows from the cigar itself: Since the D8 and the anti-D8 branes cannot reach all the way to the origin u = 0, they must reconnect to each other and form U-shaped configurations as shown on figure 4. The As shown on figure 4, the U-shaped profiles of reconnected branes depend on one parameter -the asymptotic distance L between the D8 and the anti-D8 branes along the S 1 circle for u → ∞. For L = πR the branes form the antipodal configuration in which the branes remain at opposite points on the circle for all u, all the way from u → ∞ down to u = u Λ where the branes reconnect; this is the original configuration of Sakai and Sugimoto. In the more general version of the model [4] we allow for the L < πR non-antipodal configurations. In such configurations, the distance between the branes in the x 4 direction depends on u -it becomes smaller for smaller u -and eventually the branes reconnect at u 0 before they reach the bottom of the cigar. The u 0 /u Λ ratio may be used to parametrize the non-antipodal configurations instead of the L/R.
The u 0 /u Λ or L/R parameter of the Sakai-Sugimoto model does not correspond to any adjustable parameters of the real-life QCD. Unfortunately, this parameter affects many physical properties of the model. For example, for (L/R) > 0.97 the deconfinement and the restoration of chiral symmetry happen at the same temperature, but for (L/R) < 0.97 they happen at different temperatures and the model has an intermediate deconfined but chirally broken phase [4]. Also, in the antipodal model the central nuclear forces are purely repulsive, while the non-antipodal models give rise to both repulsive and attractive nuclear forces [2] (but unfortunately the net force remains repulsive at all distances).
The low-energy dynamics of the flavor degrees of freedom living on the D8 branes is governed by the effective action comprising the Dirac-Born-Infeld (DBI) and Chern-Simons (CS) terms, S = S DBI + S CS . (2.10) The DBI action is where is the D8-brane tension, g mn is the nine-dimensional induced metric on the branes, F mn is the U(N f ) gauge field strength, and Str denotes the symmetrized trace over the flavor indices. In the limit of fixed brane geometry and weak gauge fields, the DBI action reduces to Yang-Mills, (2.12) Furthermore, the low-energy field modes we are interested in are constant along the S 4 sphere and the vector fields' directions are ⊥ S 4 . Consequently, we are going to dimensionally reduce the flavor gauge theory down to 5 dimensions: the 4 Minkowski dimensions x 0,1,2,3 , plus one coordinate z along the U-shaped line on the cigar. We find it convenient to choose a particular z coordinate that makes the 5D metric conformal although its relations to the u and x 4 coordinates of the cigar follow from the differential equations . (2.14) In the (x 0 , x 1 , x 2 , x 3 , z) coordinates, the 5D YM action for the flavor gauge fields becomes Near the bottom of the U-shaped flavor branes -which is the only region we are going to care about in this article - (2.17) The Chern-Simons term arises from the couplings of the gauge fields on the D8 brane to the bulk Ramond-Ramond field. In 9 dimensions After integrating over the S 4 and dimensionally reducing to 5D, the Chern-Simons term becomes In a particularly interesting case of 2 flavors, it is convenient to separate the U(2) gauge fields A M into their SU(2) components A M and the U(1) componentsÂ M . In terms of these components, the Chern-Simons action becomes (2.20) We shall see in a moment that the baryons and the multi-baryon systems have strong selfdual SU(2) magnetic fields F µν . 3 Thanks to the first term in this Chern-Simons action, the instanton number density acts as electric charge density for the abelian fieldÂ 0 ; the net electric charge of an instanton is Q el = N c /2.
Besides the U(N f ) gauge fields, the effective low-energy 5D theory also contains the scalar fields Φ a (x, z) describing the small fluctuations of the D8 branes in the transverse direction. For N f branes, the scalars form the adjoint multiplet of the U(N f ) gauge symmetry. The action for the scalar fields follow from the DBI action for the embedded metric g mn of the fluctuating branes. For the Φ(x, z) fields normalized to have similar kinetic energies to the vector fields, the scalar action looks like The details of the scalar potential V (Φ) = m 2 (z) × Φ 2 + a(z) × Φ 4 + · · · need not concern us here, what's important is the second term describing the backreaction of the gauge fields on the brane geometry. In the antipodal Sakai-Sugimoto model C(x) ≡ 0 and there is no backreaction because of a geometric symmetry, but in the non-antipodal models C(z) = 0 and the scalar fields Φ induced by the vector fields of the baryons lead to attractive nuclear forces [2]. The ratio of these attractive forces to the repulsive forces mediated by the abelian electric fields depends of the C(z) profile of the interaction term (2.22). For the Sakai-Sugimoto models so the net force is repulsive.
To see how that works, let us focus on the baryons in the Sakai-Sugimoto and other models of the holographic QCD.

Baryons as flavor instantons
In a compact space like S 4 , the net electric charge must vanish. and act as N c quarks.
We may put the D4 brane anywhere in space and anywhere on the cigar. However, the S 4 volume increases with the u coordinate, so the lowest-energy location of the D4 is the cigar's tip u = u Λ . At other location, the brane feels a gravity-like force pulling it down to the tip. However, the strings connected to the baryon vertex pull it towards the flavor D8 branes; in the non-antipodal models the D8 branes do not reach the cigar's tip, so the strings pull the baryonic vertex up from the tip towards the lowest point u 0 of the flavor branes. The 'tag of war' between the upward and downward forces on the baryon vertex determines its ultimate location. In some models, the forces reach equilibrium for the baryon vertex hanging on strings below the flavor branes [35], while in many other models -including the non-antipodal Sakai-Sugimoto -the string forces win and pull the baryon vertex all the way to the lowest point u 0 of the flavor branes [30].
In all such models, the baryonic vertex is a D p brane completely embedded in a stack of D p+4 flavor branes, so it is equivalent to zero-radius Yang-Mills instanton of the U(N f ) gauge symmetry on the flavor branes, and for N f > 1 it may be smoothly inflated to a finite-radius instanton [31]. In p + 5 dimensions of the flavor branes, this instanton is a fat D p -brane wrapping some compact cycle, but once we dimensionally reduce to 5 dimensions, the instanton becomes a finite-size particle. Thus, in the low-energy effective 5D theory of the holographic QCD, a baryon is realized as a finite-size instanton of the U(N f ) gauge theory.
Our research program of studying the instanton crystals as holographic duals of the cold nuclear matter stems from this realization of baryons as instantons in 5D, but it does not depend on any specific details of the Sakai-Sugimoto model. We use that model as a prototypical example, but any other model where the baryons are realized as instantons would be just as good for our purposes, for example the model of [26] for some 7-brane geometries [35], or the AdS 5 × S 1 model of [32] (the baryons of that model are studied in §6 of [30]).
In the baryon-vertex picture, each of the N c strings connecting the vertex to the flavor branes has electric charge 1/N f under the abelian U(1) subgroup of the U(N f ), so the whole baryon has abelian charge N c /N f . In the instanton picture, the same electric charge obtains from the Chern-Simons coupling between the abelian electric field and the non-abelian magnetic fields of the instanton.
where I(x) is the instanton number density of the magnetic fields. For N f ≥ 3, the Chern-Simons couplings also endow instantons with non-abelian electric charges. Altogether, where the net electric charge Q a el has both abelian and non-abelian components; in (2.28) Abelian or nonabelian, the like-sign electric charges repel each other; it is this Coulomb repulsion between different parts of the same instantons that prevents it from collapsing to a point-like D-brane.
However, the instantons do not grow large because the 5D gauge coupling decreases away from the z = 0 hyperplane: a large instanton would spread into regions of space where the coupling is weaker, and that would increase the instanton's energy. Instead, the equilibrium radius of the instanton scales like For a holographic model of a baryon, this radius is unrealistically small. Indeed, using the ρ meson's mass as a unit, the real-life baryon radius KK . Moreover, it raises the question of whether we may adequately describe so small an instanton using the DBI + CS action, or perhaps we need to include the higherorder stringy corrections. 4 On the other hand, assuming the DBI+CS description is OK, the radius as in eq. (2.29) allows for consistent expansion of the instanton's parameters in powers of λ −1 [5]. In particular, the leading contribution to the instanton's mass is O(λ × N c M KK ) while the corrections due to z-dependent gauge coupling and due to Coulomb self-repulsion To see how that works, consider a static instanton -a time-independent configuration of SU(N f ) magnetic fields, plus the electric fields induced by the CS couplings, and the scalar fields induced by the tr(ΦF F ) coupling to the magnetic fields. Since the canonically normalized couplings in 5D are O(λ −1/2 ), the leading contribution to the instanton's energy comes from the magnetic field. In the DBI approximation 5 Both g YM (z) and K(z) depend on the z coordinate on the distance scale O(1/M KK ), so for instantons of much smaller size we may start with the approximation of constant K and constant 5D gauge coupling. And in this approximation, the DBI energy is minimized by the magnetic fields that are exactly self-dual (in the 4D space of (x 1 , x 2 , x 3 , z)); moreover, the DBI energy of an instanton is equal to its Yang-Mills energy regardless of its radius and of the K parameter of the DBI action.
To determine the equilibrium radius of an instanton we need to calculate its energy to the next order in 1/λ expansion. To this order, we assume the magnetic fields to be exactly self-dual -which allows us to use the YM action instead of DBI even for small instantons -but the gauge coupling is z-dependent, and we also account for the electric and the scalar fields. For small instantons, we may approximate the 5D gauge coupling as terms would have on a small instanton. In a supersymmetric background, the instanton is BPS and its net mass is protected against stringy corrections, so the DBI action -or even the Yang-Mills action -gives the exact value. But what happens to small instanton in non-supersymmetric backgrounds is an open question. 5 In string theory, the complete action for scalar and vector fields living on a D-brane starts with the DBI and CS terms but also includes an infinite number of terms involving D n F mn and the higher derivatives. The DBI + CS action is merely a long-distance approximation that allows for strong gauge field tensions F mn as long as their gradients are small. It is not clear how well this approximation works for small instantons, but at present we do not have a better approximation.
for some numerical constants B and D; for the Sakai-Sugimoto model B = ζ/27π and D = (8ζ 3 − 5)/9ζ 2 , while other models may have different values. Consequently, the YM energy of the nonabelian magnetic fields evaluates to where a is the instanton's radius and Z is is the z coordinate of its center.
The electric potentials A a 0 couple to the F µνF µν products of the magnetic fields while the scalar potentials Φ a couple to the F µν F µν . For the self-dual magnetic fields, both potentials couple to the same source I(x, z) × Q a el , and the only difference is the coupling strength ratio C(z) (cf. eqs. (2.22) and (2.23)). For small instantons we may neglect the z-dependence of this ratio and let C(z) ≈ C(0) ≡ C, and as long as we do not go very far from the instanton (for distances r ≪ R KK ), we may also neglect the z-dependence of the gauge coupling. Consequently, both electric and scalar potentials become the 4D Coulomb potentials (2.34) The electric fields lead to repulsive forces between the charges while the scalar forces lead to the attractive forces. Altogether, the net Coulomb energy amounts to and for a single instanton of radius a it evaluates to (2.36) Note: for C < 1 the electric fields are stronger than the scalar fields and the net Coulomb energy of the instanton is positive -which makes for a net self-repulsive force that prevents the instanton from shrinking to zero radius. In models with C > 1 (assuming they exist), the scalar fields would be stronger, the net Coulomb energy would be negative, which means a self-attractive force rather than self-repulsive. In such a model, the instanton would shrink to zero radius and our approximations would not be valid.
Altogether, the net energy of an instanton amounts to KK ) both radius-dependent terms here are O(1/λ) corrections to the leading term. Minimizing the net energy, we find the equilibrium value of the instanton the instanton center is in equilibrium at Z = 0 (the bottom of the U-shaped flavor branes), and the instanton's mass is For the antipodal Sakai-Sugimoto model We may absorb 2 out of 3 model-dependent parameters B, C, D into a redefinition of the λ and M KK parameters of the effective 5D theory, for example For the static instantons or multi-instanton systems we may also get rid of the third modeldependent parameter C. Indeed, for the static systems Φ(x, z) = C × A 0 (x, z), so the only effect of the scalar fields is to reduce the net Coulomb force by a constant factor (1 − C 2 ). We may simulate this effect without any scalar fields by using different gauge couplings for the electric and magnetic fields, g 2 el = (1 − C 2 ) × g 2 mag , or equivalently by rescaling the time Consequently, the static instanton's energy becomes Besides the radius a and Z, the instanton has other moduli -the X 1,2,3 coordinates of the center (which corresponds to the baryon's coordinates in 3D) and 4N f − 5 orientation moduli in the SU(N f ) gauge algebra; and the net energy is degenerate in these moduli to all orders in 1/λ. For finite N c -even if it's very large -one should quantize the motion of the instanton in those moduli directions. Consequently, a holographic baryon acquires definite spin and isospin quantum numbers; for N f = 2 the baryons have I = J = 0, 1, 2, . . . , Nc 2 for even N c or I = J = 1 2 , 3 2 , 5 2 , . . . , Nc 2 for odd N c [5,30], and there are similar spin ↔ flavor multiplet relations for N f > 2.
However, in multi-baryon systems interactions between the baryons break the rotational and flavor symmetries of individual baryons, and only the overall SO(3) and SU(N f ) symmetries remain unbroken. In holography, the multi-instanton systems suffer from the same problem: the magnetic fields of multiple instantons interfere with each other, which spoils the degeneracy of the net energy with respect to orientations moduli of the individual instantons. In the large N c limit, this effect becomes more important than the quantum motion in the moduli space.
Consequently, in our research program into holographic nuclear matter and related issues, we stick to classical static instantons with definite classical orientations in space and in SU(N f ). From the quantum point of view, such instantons are superpositions of states with different spins and isospins (or rather SO(4) and SU(N f ) quantum numbers). We do not even care is the instantons are bosons or fermions -at our level of analysis it simply does not matter. Minimizing the classical energy of a classical multi-instanton system with respect to classical positions and orientations of all the instantons is already a very hard problemsee the next section for the highlights, -and we do not wish to delve into the quantization of moduli until we have the classical issues under control.

Multi-Baryon Systems
In the large N c limit, nuclear forces between the baryons are dominated by the static potentials. Holographically, a static system of A baryons corresponds to a time-independent configuration of the non-abelian magnetic flavor fields F a µν (x, z) (µ, ν = 1, 2, 3, z) of net instanton number A, accompanied by the Coulomb electric A a 0 (x, z) and scalar Φ a (x, z) potentials induced by their Chern-Simons and ΦF F couplings to the magnetic fields [2]. The whole configuration should minimize the net DBI + CS energy of the system subject to the constraint (2.46).
In the λ → ∞ limit, the DBI energy of the magnetic fields is O(λ) while the net effect of the electric and scalar fields is only O(1). Moreover, the magnetic fields are concentrated within O(λ −1/2 ) distance from the X 4 = 0 hyperplane, so to the leading order we may approximate the 4 + 1 dimensional spacetime as flat. Thus, where we neglect the z dependence of the g 2 YM and K = g 11 (z)/(2πα ′ ). Similar to the Yang-Mills energy of an A instanton system, this leading-order DBI energy is minimized by the self-dual configurations of the magnetic fields F a µν (x, z). In fact, all such self-dual configurations (of the same instanton number A) have the same leading-order energies The self-dual configurations form a continuous family parametrized by A × 4N f moduli, which correspond to the locations, radii, and SU(N f ) orientations of the A instantons. But the leading-order energy (2.48) does not depend on any of these moduli.
Fortunately, the sub-leading corrections to the net energy lift the degeneracy of the leading order, which provides for the O(λ 0 ) interactions between the baryons. To work out such interactions we need the degenerate perturbation theory for the magnetic field configurations and their energies. At first order of this perturbation theory, we (a) limit the F a µν (x, z) configurations to the degenerate minima of the leading-order energy function, i. e., to the self-dual magnetic fields; (b) calculate the O(λ 0 ) corrections ∆E for the energies of these configurations; (c) minimize the ∆E among the self-dual configurations. At the next order, we would calculate the O(λ −1 ) corrections ∆F a µν (x, z) to the magnetic fields -which would no longer be self-dual -and then use such corrections to calculate the net energy to the order O(λ −1 ). But in this article, we limit our calculations to the O(λ 0 ) interaction energies and zeroth-order magnetic fields, so we shall stop after step (c).
In more detail, to get O(λ 0 ) interactions between A baryons, we proceed as follows: 1. General self-dual magnetic field configurations obtain via ADHM construction [11] in terms of A × A and A × N f matrices obeying certain quadratic constraints. Our first task is to solve these constraints and write down the ADHM matrices in terms of the instantons' locations, radii, and orientations.
2. Given the ADHM matrices, we work out the instanton number density profile (2.49) and for N f > 2 also the non-abelian adjoint density 3. Next, we calculate the O(λ 0 ×N c M) corrections to the net energy of the system. There are two sources of such corrections: the z-dependence of the 5D gauge coupling, and the Coulomb electric and scalar potentials induced by the Chern-Simons and the ΦF F couplings to the magnetic fields, thus The z-dependent 5D gauge coupling changes the DBI energy of the magnetic fields by while the Coulomb energy depends on the N f . For N f = 2, the U(2) Chern-Simons and ΦF F terms couple the SU(2) magnetic fields to the U(1) Coulomb fields only. Consequently, the A a 0 and the Φ a fields are abelian and couple to the instanton density I(x, z), and their net energy is simply 4 + 1 dimensional Coulomb energy Chern-Simons and ΦF F terms couple the SU(N F ) magnetic fields to both abelian and non-abelian electric and scalar fields; the abelian Coulomb fields are sourced by the instanton density I(x, z) while the non-abelian fields are sourced by the adjoint density I a (x, z). Altogether, the net energy of these Coulomb fields is (2.54) 4. Steps 1, 2, and 3, give us ∆E as a function of baryons' locations, radii, and SU(N f ) orientations. Now we need to minimize the ∆E with respect to all these moduli.
This 4-step procedure is fairly straightforward for a few baryons -cf. calculations of the 2-body nuclear forces by Kim, Sin and Zahed [33] and by Hashimoto et al [34], -but it becomes prohibitively difficult for large numbers of baryons and outright impossible for infinite baryon crystals. At best, we can survey a small subspace of the A-instanton moduli space and try to minimize the ∆E over that subspace. For example, we may assume that all the instantons have the same radius, that their centers form a periodic lattice of some particular symmetry, and that the orientations of the instantons also form some kind of a periodic pattern; this gives us an ansatz for all the A × 4N f moduli in terms of just a few overall parameters, and we can try to calculate and minimize the ∆E as a function of these parameters. However, any such ansatz is likely to miss the true lowest-energy configuration of the system. Indeed, in condensed matter guessing the crystalline symmetry of some substance from the properties of the individual atoms is a game of chance with poor odds, and there is no reason why the instanton crystals should be any simpler. Moreover, even if we could somehow guess all the symmetries of the instanton crystal, actually working through steps 1-4 is impossible without additional approximations (besides λ ≫ 1). For example, in our previous paper [1] we needed the (lattice spacing) ≫ (instanton radius) approximation just to solve the ADHM constraints for the zigzag-shaped chain of instantons, and even for a straight chain we needed the spacing ≫ radius approximation to calculate the Coulomb energy ∆E C .
In this article we explore a shortcut around steps 1, 2, and 3. In section 3 we shall see that when the distances between the instantons are much larger than their radii, the interactions between the baryons are dominated by the two-body forces, for a manageably simple function F 2 of the two instanton's positions and orientations. Consequently, we can minimize the net interaction energy over the entire moduli space of the multi-instanton system using a simple numerical simulation: Starting from a random set of instanton positions and orientations, use the steepest descent algorithm to find the nearest local minimum of the net energy; repeat this procedure for different random starting points to find other local minima; eventually, we find all the local minima, compare their energies, and identify the global minimum. In our next paper [41] we shall use this algorithm to construct the 2D and 3D instanton lattices starting from scratch, i. e., from the F 2 function.
In the present article, we focus on the 1D lattice ( §4) and on the first step in the transition form 1D to 2D -the zigzag ( §5). In both cases, we assume the lattice geometry for the instanton positions, but we use the numeric simulation to find the lowest-energy pattern of their orientations. Once we know all the patterns that show up, we use them as ansatz's (with a few parameters) for which we calculate the net energy as analytic functions of the parameters. In this way, we map the boundaries between different orientation patterns much more accurately than we could do it in the numerical simulation.
To make the baryons form a 1D lattice instead of spreading out in three dimensions, we curve the x 2 and x 3 dimensions of the flavor brane similar to the curvature of the z coordinate. In terms of the effective 4 + 1 dimensional theory, this corresponds to the 5D flavor gauge coupling depending on the x 2 and the x 3 as well as the x 4 ≡ z, This gauge coupling acts as a harmonic potential for the instantons which pulls them towards the x 1 axis, so at low densities the instantons form a 1D lattice along the x 1 . At higher densities, the instantons push each other away from the x 1 axis and form more complicated 2D or 3D lattices, starting with the zigzag To make sure the transition from the straight chain to the zigzag happens for lattice spacings much larger than the instanton radius (which is required by the two-body force approximation), we assume This setting is similar to what we have used in our previous paper [1], except for a change of notations 6 and allowing for M 3 = M 4 . In §5 we shall see that M 3 = M 4 makes a big difference for the instanton orientation patterns in a zigzag: the lowest-energy patterns are abelian for M 3 ≈ M 4 but non-abelian for anisotropic M 3 < 0.94M 4 . Even the order of the phase transition from the straight chain to the zigzag depends on the M 3 /M 4 ratio: second-order for ratios > 0.725 but first-order for ratios < 0.725.

Two-Baryon versus Multi-Baryon Interactions
In real-life nuclear physics, besides the two-body nuclear forces due to meson exchanges, there are significant three-body forces, and presumably also four-body forces, etc., where n stands for the quantum numbers of the n th nucleon. Likewise, in the holographic nuclear physics interactions between multiple baryons include two-body forces and also threebody, four-body, etc., forces. Even in the classical infinite-mass limit where the holographic 6 Back in [1], we had so the zigzag was in the (x 3 , x 4 ) plane, with the long axis along the x 4 . Translating between the old notations and those of this paper, we have ( baryons become static instantons of the SU(N f ) gauge fields in 4+1 dimensions, the potential energy (due to g 5D = const and due to Chern-Simons interactions) of an A-instanton system has form with significant three-body, etc., terms.
What about the relative magnitudes of the 2-body, 3-body, 4-body, etc., interaction terms? When the baryons are packed cheek-by-jowl till their instanton cores overlap and merge, we expect all the n-body forces to have comparable strengths. But in the opposite low-density regime of baryons separated by distances much larger then their radii, the twobody forces dominate the interactions, while the multi-body forces are smaller by powers of (radius/distance) 2 . The purpose of this section is to prove this fact.
Before we delve into the technical issues, let us briefly summarize our argument. First, consider the simplest case of N f = 2 and intermediate-range distances between the baryons, which allows us to treat the 5D holographic space as approximately flat. For N f = 2 the holographic baryons are instantons of the SU(2) magnetic fields, which source the U(1) electric and scalar fields via CS and tr(ΦF µν F µν ) couplings. When those instantons are small and separated from each other by large distances |X m − X n | ∼ D ≫ a, their interactions come from two sources: 1. Direct Coulomb repulsion (electric − scalar) between nearly-point-like abelian charges in 4 + 1 dimensions, which is a manifestly two-body interaction.
2. The interference between the instantons changes the distribution of the instanton number density in space, which in turn changes the self-interaction energy (Coulomb and non-abelian) of each instanton by an amount comparable to (3.4). Within an instanton -i. e., at O(a) distance from some instanton's center X µ n -the interference from the other instantons should be relatively weak, so there should be some kind of a perturbation theory for it. At the first order of such perturbation theory, the ∆I arises from interference between the un-perturbed standalone-like instantons, so we expect it to be a sum of pair-wise interferences from the other instantons, and depends only on the instantons #n and #m -i. e., on their positions, radii, and orientations -but not on any other instantons. At the second-order, we expect to include the interference between the first-order ∆I and the additional instantons, so at this order we obtain 3-body effects, Likewise, the higher orders may involve more and more instantons, but the magnitudes of such high-order interference effects are suppressed by the higher powers of (a 2 /D 2 ). Now consider the Coulomb self-interaction of the instanton #n, where ρ n is the instanton's effective charge radius, A standalone instanton has 1 ρ 2 n = 4/5 a 2 n (3.14) but interference from the other instantons should change this radius by small amount of similar relative magnitude to ∆I/I standalone n , thus which changes the instanton's Coulomb self-interaction energy by (3.16) Note that this effect has a similar magnitude to the direct Coulomb repulsion (3.4) between the instances.
Moreover, the charge radius correction ∆ n is linear in the ∆I interference at x near the X n , hence in light of eq. (3.8), the leading-order contribution to the ∆ n is a sum of pair-wise interferences from the other instantons, thus where each ∆ n,m depends only on the instantons #n and #m. Consequently, the leading effect of the interference on the net Coulomb self-energy of all the instantons has form where the leading terms act as two-body interactions between the instantons.
Besides the Coulomb energy, the non-abelian energy is also affected by the interference between the instantons, (3.19) but this time the integral should be taken over the whole 4D space, including both the instantons and the inter-instanton space. Indeed, we shall see later in this section that at O(a) distance from an X n , ∆I ∼ 1 20) and both kinds of places make O(a 4 /D 2 ) contributions to the integral (3.19). Moreover, in both kinds of places, the leading terms in the ∆I(x) is a sum of independent two-instanton interference terms, Although our heuristic argument for such decomposition near instanton centers does not work in the inter-instanton space, we shall prove later in this section that the decomposition (3.21) works everywhere in the 4D space. Therefore, the non-abelian interactions between the instantons due to interference are also dominated by the two-body terms Thus far we have assumed N f = 2. For larger numbers of flavors, interference between magnetic fields of the SU(N f ) instantons is harder to describe mathematically, but its qualitative features for a ≪ D remain exactly as in eqs. (3.5) through (3.21). In particular, the effect of interference between small instantons separated by large distances on the density profile I(x) is dominated by the two-instanton interference terms as in eqs. (3.21) and (3.8-11). Consequently, the instantons' non-abelian and Coulomb interaction energies suffer small corrections from the interference that give rise to predominantly-two-body interactions between the instantons as in eqs. (3.18) and (3.22).
The qualitatively new aspect of N f ≥ 3 comes not from the SU(N f ) self-dual magnetic fields of the multiple instantons but from the electric and scalar Coulomb fields induced by the magnetic fields via Chern-Simons and tr(ΦF |muν F µν ) interactions. For N f = 2 such electric and scalar fields are purely abelian, A 0 , Φ ∈ U(1) ⊂ U(2), but for N f ≥ 3 they have both U(1) and SU(N f ) components. For a single instanton, all fields -magnetic, electric, and scalar -belong to the same Cartesian -so the Coulomb self-interaction of a standalone instanton works exactly as for N f = 2.
On the other hand, Coulomb interactions between multiple instantons whose magnetic fields span different SU(2) subgroups of the SU(N f ) become much more complicated.
Moreover, in a general multi-instanton background, the Coulomb fields A 0 (x) and Φ(x) do not behave as simply 1/r 2 . Instead, they satisfy gauge-covariant equations for the matrixvalued fields, and ditto for the scalar Φ(x). Consequently, the Coulomb force between two sources (such as instantons) depends on the magnetic fields in the space between them. In other words, for N f ≥ 3 the direct Coulomb interactions are no longer two-body; instead, the force between two instantons depends on the other instantons in the system.
Fortunately, for small instantons separated by large distances, the magnetic fields are concentrated in small volumes of 4-space, so their effect on the Coulomb field propagation becomes perturbatively weak: where b and c are adjoint indices of the U(N f ). The leading term here is the ordinary Coulomb propagator in 4D space, thus in the a ≪ D limit, the direct Coulomb interactions between small instantons are predominantly two-body, where overlap(n, m) ≥ 0 is the overlap between the SU (2) Fortunately, we do not need a hard calculation to see that at large distances from an instanton its magnetic fields are very weak. Indeed, even in flat 5D space the fields weaken with distance as A a µ ∼ a 2 /r 3 (in the IR-safe singular gauge), so at r ≫ a they are so weak that the field equations become effectively linear. In the curved space, we may decompose the weak 5D gauge fields into 4D mesonic fields, hence at large distances r 1/M from an instanton, its fields become Yukawa-like Consequently, at the location X µ n of any particular instanton, the background fields from the other instantons are very weak, and their effect on the instanton #n itself can be adequately accounted by the first-order perturbation theory. In other words, the effects of other instantons #m = n on the instanton #n are weak and add up linearly! For the self-interaction energy of the n th instanton, this means where the second term gives rise to two-body interaction energies We expect the two-body terms to be rather small -in addition to the usual 5D 1/|X m −X n | 2 factors they should carry Yukawa-like exponentials e −mr (or rather sums of such exponentials), but these are the leading interactions due to interference. To three-body or multibody interactions follow from higher-order perturbation by very weak fields (3.27), so they are much smaller that the two-body interactions.
As to the direct Coulomb interactions between the holographic baryons via electric or scalar fields, at large distances |X m − X n | ∼ 1/M they also decompose into sums of Yukawa forces. Moreover, since the scalar mesons generally have different masses from the vector mesons, the attractive potential due to scalars may have a different r dependence from the repulsive potential due to vectors. Thus, for a right model, the net two-body force between two holographic baryons may become attractive at large enough distances between them. But regardless of the model, the direct Coulomb interactions are always manifestly two-body for N f = 2, while for N f ≥ 3 the multi-body terms exist but become very small at large distances between the baryons. Thus, where the precise form of the potential V (2) (r) is model-dependent, but the two-body form of the direct interactions is quite universal.

⋆ ⋆ ⋆
At this point we have finished our heuristic arguments and about to start a more technical analysis of the interactions between the holographic baryons. To keep our presentation relatively simple, we assume just two quark flavors and intermediate-range distances between the holographic baryons -D ≫ a but D ≪ 1/M.
For N f = 2 the non-abelian gauge symmetry is SU(2), so the ADHM data for an Ainstanton configuration consists of a quaternionic array Y n (n = 1, . . . , A) and a symmetric quaternionic A × A matrix Γ mn . Equivalently, we may use A real instanton radii a n = |Y n |, A SU(2) matrices y n (equivalent to unimodular quaternions Y n /a n ) parametrizing the instanton's orientations, and 4 real symmetric matrices Γ µ mn whose diagonal elements Γ µ nn = X µ n are the locations of the instanton's centers. Finally, the off-diagonal matrix elements Γ µ m =n = α µ mn follow from the other parameters by solving the ADHM equations or equivalently η a µν Γ µ , Γ ν mn + a m a n × tr y † m y n (−iτ a ) = 0 (3.32) where η a µν is the 't Hooft's tensor mapping between the SU(2) gauge and the SU(2) L inside Spin (4) The ADHM data are somewhat redundant -an O(A) symmetry acting on all the Γ µ mn and Y n = a n y n does not change any physical properties of the multi-instanton data. This symmetry includes Z A 2 which flips the instanton orientations y n → −y n (independently for each n). It also includes small SO(A) rotations that change the off-diagonal elements by δα µ mn = ǫ mn (X m − X n ) µ + O(ǫ 2 ). To eliminate these rotations, the ADHM equations (3.32) for the off-diagonal elements should be combined with additional constrains (one for each m = n), for example ∀m = n : (X m − X n ) µ α µ mn = 0. (3.34) For large distances |X m − X n | ∼ D between the instantons, D ≫ a n , we may solve the ADHM equations (and the constraints (3.34) as a power series in a 2 /D 2 : Note that for each off-diagonal matrix element, the leading term α (1) µmn in this expansion depends only on the instantons #m and #n (i. e. , on the positions, radii, and orientations of only these two instantons) while the subleading terms α (2) µmn , α (3) µmn , . . . involve additional instantons.
Given the ADHM data Γ µ mn , a n , and y n of a multi-instanton system, the 4D instanton number density (3.5) can be obtained as Thanks to the double D'Alambertian in eq. (3.36), several moments of the instanton density may be obtained via integrating by parts: where T mn ≡ 1 2 a m a n tr y † m y n , (3.38d) The quadratic moment (3.38c) for µ = ν = 4 is particularly important since it gives us an exact formula for the non-abelian energy of the system directly in terms of the ADHM data, Obviously, the first sum on the last line here is the sum of individual instantons' potential energies due to their radii and locations (relative to the x 4 = 0 hyperplane) while the second sum comprises the interactions between the instantons, This formula is general and applies to any multi-instanton configuration, the problem is calculating the off-diagonal matrix elements α µ mn = Γ µ m =n . Fortunately, for the small instantons separated by large distances O(D) ≫ a, those off-diagonal matrix elements are given by the expansion (3.35.a-d) in powers of (a 2 /D 2 ). Moreover, for each m = n the leading term α while the multi-body interactions are O(N c λM 3 a 6 /D 4 ) or weaker. Specifically, where N mn is the 3-vector part of the unit 4-vector Now consider the Coulomb energy of the multi-instanton system. Alas, we cannot obtain this Coulomb energy directly from the ADHM data via integration by parts or some other clever trick but have to calculate it the hard way: Derive the instanton density profile from eq. (3.36), decompose I(x) into individual instantons plus interference terms, and finally perform the integral (3.44) for some approximation to the integrand.
Our starting point is the L mn (x) matrix (3.37) which governs the instanton density profile. For a system of small instantons separated by large distances, the diagonal elements of this matrix are much larger than the off-diagonal elements. Indeed, the diagonal matrix elements are given by This hierarchy between the diagonal and the off-diagonal matrix elements of L allows us to evaluate log det(L) as a power series in in the ratios of the off-diagonal to diagonal elements: log det(L) = log det L diag + log det 1 + L −1 diag L offdiag = tr log L diag + tr log 1 (3.47) Consequently, the net instanton density profile of the system may written as a sum of one-instanton, two-instanton, three-instanton, etc., terms:  which is precisely the density profile of a standalone instanton of radiusâ n . Note that this radius is slightly larger than the original instanton radius a n -this is a hidden effect of the interference from the other instantons m = n -while the visible effects of the interference come via the two-instanton terms I mn (x), the three-instanton terms I ℓmn (x), etc. Deriving explicit formulae for all these terms is a straightforward (albeit rather tedious) exercise, and the results are too cumbersome to present here. Instead, let us simply estimate the magnitudes of the multi-instanton terms.
In the inter-instanton space -i. e. at O(D) distances from all the instanton centers - L ℓm L mn L nℓ L ℓℓ L mm L nn ∼ a 6 D 6 , etc.
L ℓm L mn L nℓ L ℓℓ L mm L nn ∼ a 6 D 10 , . . . , (3.55) and therefore In the inter-instanton space ℓmn (x) ∼ a 6 D 10 , . . . . (3.56) As promised in eq. (3.21), the interference between small well-separated instantons is dominated by the two-instanton terms. In fact, in the inter-instanton space the two-instanton terms have similar magnitudes to the one-instanton terms, However, the entire I(x) in the inter-instanton space makes a negligible contribution to the net Coulomb energy of the system, Therefore, the net Coulomb energy of the multi-instanton system may be summarized as To work out the effects of interference on the charge radius ρ n of the n th instanton we need to estimate (and then evaluate) the interference terms I ℓm , I kℓm , etc., for x µ near the n th instanton's center X µ n . In this neighborhood, the estimates (3.56) do not work -or rather they work only for the interference terms between the other instantons #k, ℓ, m, . . . = n.

But the terms I
(2) mn , I ℓmn , etc., that do involve the n th instanton itself turn out to be much larger than (3.56) because of the 1/L nn (x) factor that happens to be O(1/a 2 ) instead of O(1/D 2 ). Moreover, this factor depends on x much more rapidly than all the other factors -on the scale of x ∼ a instead of x ∼ D -so taking the space derivatives of this factor produces much larger results than taking the same derivatives of the other factors. Thus and therefore etc., etc.
The first factor in all these interference terms evaluates to etc., etc.
Therefore, the two-instanton interference terms I 2 mn affect the Coulomb energy by the amount comparable to the direct repulsion (between point-like instantons) but the effect of interference terms involving three or more instantons at once is suppressed by extra powers of a 2 /D 2 and becomes negligible in the limit of small instantons separated by large distances. Since we have already verified similar behavior of the non-abelian energy of the multi-instanton system, this completes the proof of our theorem. Quod erat demonstrandum.
To conclude this section, we would like to derive explicit formulae for the two-instanton interference terms (near one of the two instantons) and hence the Coulomb interaction energy due to interference. Starting from eq. (3.61) we have the first factor on the RHS evaluated in eq. (3.63), so let us evaluate the second factor. For x µ = X µ m , eq. (3.46) gives us L mn (x = X I ) = L nm (x = X I ) = (X m − X n ) µ α µ mn + ℓ =m,n α µ ℓm α µ ℓn + a m a n × 1 2 tr y † m y n (3.68) where the first term on the RHS happens to vanish -cf. eqs. (3.34) -while the second term is O(a 4 /D 2 ), which is much smaller than the O(a 2 ) third term. Thus, L mn (x = X n ) = L nm (x = X n ) ≈ a m a n × 1 2 tr y † m y n (3.69) while and therefore It is easy to check that such two-instanton interference terms do not affect the net charge of the n th instantonnear Xn d 4 x I (2) mn (x) = 0 (3.73) -but they move the charge distribution close to the center X n since I (2) mn is positive for r n <â n but negative for r n >â n . Consequently, they decrease the effective charge radius ρ n , which corresponds to ∆ n > 0 and hence positive extra Coulomb energy. Specifically, (3.74) Note however that this ∆ n is a correction of the instanton's charge radius or (rather 1/ρ 2 n ) starting from the standalone radiusâ n instead of the original a n . Consequently, there is an additional correction The α µ mn are spelled out in eqs. (3.35); to the leading order and therefore (3.80) Finally, let us combine the non-abelian and the Coulomb energies of the system and re-organize the net energy into one-body, two-body, etc, terms: (3.82) (3.84) When comparing magnitudes of the non-abelian and the Coulomb terms here we have used . In fact, we may be more precise by noticing that for small instantons distant from each other, the one-body potential energies for the instanton radii a i are much larger than the two-body, etc., interactions between different instantons. Consequently, in the minimal-energy or near-minimal-energy configuration of the multi-instanton system, the instanton radii will be close to the equilibrium radius of a stand-alone instanton, where minimizes the one-body potential energy (3.82). Plugging this equilibrium radius into eq. (3.83), we obtain a simpler formula for the two-instanton interaction energy: (3.87) Note that the expression inside '[· · · ]' is always positive, so the two-body forces between the instantons are always repulsive, regardless of the instantons' SU(2) orientations. However, the orientations do affect the strength of the repulsion: two instantons with similar orientations repel each other 9 times stronger then the instantons at the same distance from each other but whose orientations differ by a 180 • rotation (in SO(3) terms) around a suitable axis. This fact will be at the core of our analysis of instanton crystals in subsequent sections.

Linear Chains of Instantons
Consider an infinite set of instantons arranged in a crystalline lattice. In light of eq. (3.87), the orientations of instantons at different lattice sites affects the system's energy just as much as the lattice geometry. In this section we shall see how this works for a simple onedimensional low-density lattice -i. e. an infinite chain of small instantons located along a straight line at X µ n = (nD, 0, 0, 0), n ∈ Z (4.1) where the lattice spacing D is much larger than the instanton radius a.
Before we go any further, we need a reason why the instantons would form a 1D chain rather then spread out into all 3 flat space dimensions like holographic baryons (or real-life baryons) in the compressed nuclear matter. In principle, we may sidestep this question by simply freezing the X µ n moduli of the instantons and allow only their orientation moduli y n to vary, but then it would be hard to place such a system in any kind of a physical context. Alternatively, we can make 5D gauge coupling depend on the x 2 and x 3 coordinates as well as the holographic dimension x 4 which creates an effective potential for the instanton centers X n , 3) At low instanton densities, this potential makes the instantons line up along the x 1 axis, while at higher densities, the repulsion between the instantons becomes stronger than this potential and they spread out into the other dimensions. We shall address such spreading out in the following section §5, but for the moment let us focus on the low-density 1D chains (4.1).
In a bigger holographic picture, a g 5 that depends on three coordinates x 2,3,4 indicates that the flavor branes are curved in all 3 of these directions, which is quite different from the usual holographic setup where the energy-scale coordinate U depends only on the x 4 . But in this article we are not going to explore the geometric aspects of such curvature. Instead, we shall simply use the gauge coupling (4.2) in an approximately flat space (for distances ≪ 1/M) as a tool to put instantons into a 1D lattice. In our next article [41], we shall make g 2 depend only on x 3 and x 4 to put instantons into a 2D lattice in the (x 1 , x 2 ) plane, and eventually go back to g 5 depending only on the x 4 to make 3D instanton lattices.
For simplicity, let us assume that the g 5 is much less sensitive to the x 2 and x 3 coordinates than the x 4 , thus M 2 2 , M 3 3 ≪ M 2 4 . In this limit, the M 2 2 and M 3 3 parameters give rise to the x 2 2 and x 3 3 terms in the one-body instanton potential (4.3), but their effect on the instanton radius a or the two-body forces between the instantons may be neglected (compared to the effect of the M 2 4 x 2 4 term). Therefore, the net two-body forces between the instantons remain approximately as in eq. (3.87).
Moreover, for the 1D lattice geometry (4.1) of the instanton centers, we have |X m −X n | 2 = D 2 × (m − n) 2 while N mn = (±1, 0, 0). Consequently, the net energy of the instanton chain as a function of the instantons' orientations y n may be summarized as To minimize this energy, each pair of instantons m and n wants to have y † m y n to be a linear combination of iτ 2 and iτ 3 ; in SO(3) terms, this corresponds to a relative rotation (between the 2 instantons) through a 180 • angle around some axis ⊥ x 1 . Alas, this cannot be achieved for all instanton pairs at once; indeed, if we minimize the energies of the (n, m) and the (n, ℓ) pairs, then the energy of the (m, ℓ) pair would be maximal instead of minimal: If y † n y m = i n 1 · τ & y † n y ℓ = i n 2 · τ for n 1 , n 2 ⊥ x 1 Then y † m y ℓ = ( n 1 · τ )( n 2 · τ ) = ( n 1 · n 2 ) + i( n 1 × n 2 ) · τ for n 1 × n 2 x 1 =⇒ tr 2 y † m y ℓ + tr 2 y † m y ℓ (−iτ 1 ) = 4( n 1 · n 2 ) 2 + 4( n 1 × n 2 ) 2 1 = 4 (maximum). for some angles φ n . Indeed, for any set of the angles φ n we have Eqs. (4.8) for various angles φ n define a big family of instanton configurations. Surprisingly, all these configurations have exactly the same net energies! Indeed, for all sets of φ n , all instanton pairs (m, n) with odd m − n have minimal energies while pairs with even m − n have maximal energies: (4.10) Consequently, regardless of the angles φ m , the net energy per instanton of the 1D lattice is for odd ℓ, 9 2 for even ℓ, (4.11) Now let's consider the instantons' orientations y n from the group-theoretical point of view. Each y n is an SU(2) matrix, but its overall sign is irrelevant, so we only care about the y n (modulo sign), which belongs to the SU(2)/Z 2 ∼ = SO(3). In a generic lowest-energy configuration (4.8) of the 1D chain, the instanton orientations span an SO(2) × Z 2 subgroup of the SO(3) corresponding to the rotational symmetries of a cylinder -rotations through arbitrary angles around the x 1 axis, and 180 • rotations around axes ⊥ x 1 . The y n alternate between the two types of rotations, but apart from that they generically do not follow any regular patterns.
However, the family (4.8) also contains some regular patterns in which the orientations y n (modulo sign) follow a repeating cycle of finite length p; moreover, the values of y n span a discrete subgroup of the cylindrical symmetry SO(2) × Z 2 . Here are some examples: • The anti-ferromagnetic chain, with 2 alternating instanton orientations: y even n = ±1, y odd n = ±iτ 3 . (4.12) In this configuration -which obtains for φ n ≡ 0 -the y n (modulo sign) span a Z 2 subgroup of the SO(2) × Z 2 .

⋆ ⋆ ⋆
All in all, the 1D lattice of instantons has a huge degenerate family (4.8) of lowest-energy configurations. However, this degeneracy obtains only in the M 2 , M 3 ≪ M 4 limit where the 5D gauge coupling (4.2) depends mostly on the holographic coordinate x 4 and only a little bit on the x 2 and x 3 . For the larger M 2 , M 3 ∼ M 4 , the degeneracy is lifted and there is a unique lowest-energy configuration of the 1D lattice, namely the link-periodic array (4.16) whose ϕ and θ parameters depends on the M 2 /M 3 ratio.
To see how this works we note that for M 2 ∼ M 3 ∼ M 4 , all 3 coordinates x 2,3,4 transverse to the 1D lattice play similar roles, and the 3 terms M 2 4 x 2 4 + M 2 3 x 2 3 + M 2 2 x 4 2 inside '(· · · )' in eq. (4.2) give rise to similar contributions to the non-abelian energy of a multi-instanton system: (where M 2 4 is the same as M 2 ). Consequently, the two-body interaction energy between the instantons becomes more complicated then in eq. and (4.21) Note that the coefficient of each a 2 n in the one-body potential (4.20) depends on the M 2 and the M 2 2 as well as the M 2 4 ≡ M 2 . Consequently, when we combine this potential with the Coulomb one-body potential for the instanton radius a, the equilibrium radius ends up smaller than in eq. (3.86), namely .

(4.22)
Plugging this radius into eq. (4.21) and combining with the Coulomb two-body interaction term from eq. (3.80), we obtain the net two-body force in the background with For instantons located along a straight line, this formula can be simplified a bit using N ν mn ≡ (±1, 0, 0, 0) and hence η a µν N ν mn tr y † m y n (−iτ a ) = ± tr y † m y n (−iτ 1 ) for µ = 4, ± tr y † m y n (−iτ 2 ) for µ = 3, ± tr y † m y n (−iτ 3 ) for µ = 2. (4.25) Consequently, the net interaction energy of the 1D lattice is where Q(m, n) def = 1 2 + tr 2 y † m y n + C 4 tr 2 y † m y n (−iτ 1 ) (4.27) Without loss of generality we assume M 4 ≥ M 3 ≥ M 2 and hence C 4 ≥ C 3 ≥ C 2 -otherwise we simply re-label the coordinates x 2,3,4 . Thus, the lowest-energy relative orientation of an instanton pair (m, n) is y † m y m = ±iτ 3 . If this orientation is unachievable, the next best bet is a linear combination of iτ 3 and iτ 2 . If that is also unachievable, the third-best choice would be a linear combination of all three iτ 1,2,3 ; in SO(3) terms, this corresponds to the 180 • rotation around a generic axis. Finally, rotations through other angles would be the least attractive option which the instantons #m and #n would rather avoid -unless they are forced to it by interactions with the other instantons between m and n.
In this setting, it is not intuitively obvious how to balance the energy cost of nearestneighbor interactions versus next-to-nearest neighbors and more distant instanton pairs. Instead of intuition, we have performed a computer experiment using a lattice of 200 SU(2) matrices y n . In each run, we start with the y n being independent random elements of the SU(2) group, and then let them evolve towards a minimum of the energy function (4.26) via the relaxation method. That is, we let the y n evolve with time according to where K is a constant mobility factor and the derivative with respect to an SU(2) matrix y n is defined as Each run ends when all the y n seem to converge to an equilibrium configuration and their derivatives (4.28) become very small. We made many runs for different C 3 and C 4 parameters, and here is what we have found: • In all equilibrium configurations for backgrounds with C 3 < C 4 , all twists y † n y n+1 between nearest neighbors have form y † n y n+1 = cos ψ n × iτ 3 + sin ψ n × iτ 2 + (tiny) × iτ 1 + (tiny) × 1, (4.30) where the tiny coefficients of 1 and iτ 1 are artefacts of imperfect convergence to equilibrium (they get smaller when we allow more time for convergence).
• The backgrounds with M 3 = M 4 =⇒ C 3 = C 4 are trickier due to symmetry between iτ 1 and iτ 2 components of the SU(2) matrices. (Or all three iτ 1,2,3 matrices for M 2 = M 3 = M 4 =⇒ C 2 = C 3 = C 4 .) In equilibrium configurations for such symmetric backgrounds, the nearest-neighbor twists y † n y n+1 generally do have iτ 1 components. However, all such components can be eliminated by a suitable global symmetry of the system, y n → y n × G, same G ∈ SU(2) for all y n , and then all the nearest-neighbor twists y † n y n+1 take the form (4.30).
• For most nearest-neighbor pairs, the angles ψ n in eqs. (4.30) take values ±ϕ (mod π) where ϕ depends on the C 2,3,4 parameters. Moreover, the ± sign here tends to alternate between odd and even n. However, some pairs do not follow these rules.
In general, our runs end up with patterns of (−1) n ψ n that look like this: Physically, this corresponds to the 1D lattice having two degenerate ground states, one with ψ n = +(−1) n ϕ and the other with ψ n = −(−1) n ϕ. In our numerical runs, some parts of the lattice converge to one of these states while other parts converge to the other state; altogether, we end up with several domains separated by walls. In a perfect simulation, the walls would eventually move towards each other and annihilate, so the whole lattice would end up in the same ground state. But this process is so slow that our numeric simulation stops before it barely begins, thus we always end up with multiple domains instead of a single ground state for the whole lattice.
⋆ The bottom line of our numerical simulations is that the ground state of the 1D instanton lattice always has a link-periodic instanton orientations (4.16) (up to a global symmetry, if any) for a periodicity angle ϕ that depends on ratios of parameters Once we know the ground state of the instanton orientations in a 1D lattice is linkperiodic, we may calculate the periodicity angle ϕ analytically. Indeed, applying eqs. (4.10) to the link-periodic array (4.16) of instanton orientations, we have: For even n − m, y † m y n = cos((n − m)ϕ) × i n−m + sin((n − m)ϕ) × i n+m+1 τ 1 =⇒ Q(m, n) = 1 2 + 4 cos 2 ((n − m)ϕ) + 4C 4 sin 2 ((n − m)ϕ) (4.32) where the ± sign is (−1) n . To calculate average interaction energy per instanton in the lattice, we should average Q(n, m) between odd and even n for fixed ℓ = n − m; indeed Thus, for even ℓ, Q(ℓ) = 5 2 + 2C 4 + 2(C 3 + C 2 ) cos(2ℓϕ), for odd ℓ, Q(ℓ) = 1 2 + 2(C 3 + C 2 ) − 2(C 3 − C 2 ) × cos(2θ) × cos(2ℓϕ), (4.35) and therefore where the sums are evaluated assuming |ϕ| < π Minimizing this expression WRT the ϕ and θ produces 4 degenerate minima, namely (4.37) however the last two minima are physically equivalent to the first two. Thus, in agreement with our computer 'experiments', the 1D instanton lattice has 2 degenerate ground states related by the ϕ → −ϕ symmetry. The value of |ϕ| also agrees with our 'experimental data': For all configurations in this family, the instanton orientations y n (modulo signs) repeat with period 4 while spanning the Klein groups Z 2 × Z 2 of 180 • rotations around 3 mutually ⊥ axes. For θ = − π 4 the y n are spelled out in eq. (4.13); for other values of θ we have similar cycles related by the U(1) ⊂ SU(2) symmetry.
Finally, for the very asymmetric background with M 2 2 ≪ M 2 3 , the two ground states (4.38) become indistinguishable as ϕ → ±0. In this limit, the instanton orientations form the anti-ferromagnetic order (4.12).

Instanton Zigzags
In this section we consider the most likely first step in the transition between 1D and 2D instanton lattices -a lattice with one infinite dimension x 1 while the other dimension x 2 has just two layers. Since the instantons repel each other, the two layers should be staggered in the x 1 direction, so the whole lattice looks like a single zigzag-shaped chain: In such a zigzag, the instanton #n has its center at X 1 n = n × D, X 2 n = (−1) n × ǫ, X 3 n = X 4 n = 0, (5.2) for some parameters D and ǫ; we shall refer to D as the lattice spacing and to ǫ as the zigzag amplitude.
To make the instantons form a zigzag lattice, we use the same background we have used in the last section to make a straight chain but increase the instanton density ρ = 1/D until the repulsion between the instantons pushes them away from their neighbors in some transverse direction. For M 2 < M 3 ≤ M 4 , the initial breakout from the straight line is going to be in the ±x 2 directions, with opposite signs for the nearest neighbors (to increase the distance between them), hence the zigzag geometry (5.2). For higher densities -and hence stronger repulsion between the instantons -we will get more complicated lattices: multiple layers in the x 2 direction, and eventually breaking out into the x 3 and x 4 dimensions. However, the first geometry right after the straight 1D lattice ought to be the zigzag.
Clearly, the above heuristic argument is not a proof. The proof -or disproof -will come from a computer simulation of multi-instanton systems where the instantons are allowed to move independently in all directions as well as to change their orientations. We shall describe such a simulation -and its results -in our next paper [41] in this project. Meanwhile, in the present paper we shall simply assume that the instanton centers from a zigzag-shaped chain (5.1) and focus on the instantons' orientations.
Since our formulae for the instanton interactions presume distances |X m − X n | much larger than the instanton radius a, we need the transition from a straight 1D chain to a zigzag to happen at a critical lattice spacing D c ≫ a, which calls for Indeed, consider an instanton that have moved a little bit away from the x 1 axis in the x 2 direction. The restoring force on this instanton due to x 2 dependence of the gauge coupling is while the repulsion from the neighboring instantons pushes it further away with the net force At the critical lattice spacing, the two forces balance each other, hence (5.6) At the same time, Consequently, eq. (4.23) can be rearranged as where Q z (m, n) = 1 2 + tr 2 y † m y n + C 3 a=1,2 Note that in general, the force between two instantons is not central -it depend not only on the distance between the instantons and their relative orientation y † m y n of their isospins, but also on the direction N mn of their separation X n − X n in space. However, in the background with M 3 = M 4 -which is the background we have used in our previous paper [1] (albeit in different notations) -C 3 = 1 2 , the last term in eq. (5.9b) goes away, and the force becomes central (but isospin-dependent).
The forces between instantons will be important in our next article [41] where we shall allow them to move in the x 1,2 plane seeking the lowest-energy lattice. For the moment, we simply assume that the instantons somehow form a zigzag of some particular lattice spacing D and amplitude ǫ and focus on their orientation moduli y n . Our immediate goal is to find the lowest-energy configuration of the orientations y n as a function of the ǫ/D and M 3 /M 4 ratios.
A priori, we do not know if the orientations -or rather the relative orientations y † n y n+1 of the neighbor instantons -follow any particular pattern or patterns. To find such pattern(s) we used numerical simulations: starting with completely random y n , we let them evolve towards a minimum of the energy function (5.9) according to eqs. (4.28); when the evolution seems to stop, we look at the relative orientations of the nearest neighbors y † n y n+1 to see if they form any recognizable pattern. Repeating this procedure for different combination of the ǫ/D and M 3 /M 4 parameters, we saw five distinct patterns of orientations (or rather four distinct patterns and one confused mess). Figure  Here is the color key to this figure: • The red dots on figure  • The yellow dots denote another abelian pattern (AB) in which all nearest neighbors differ by the same U(1) ⊂ SU(2) rotation, but now the rotation angle is < 180 • , same y † n y n+1 = exp i 2 φτ 3 for all n, 0 < φ < π. (5.11) • The blue dots denote a non-abelian link-periodic pattern (NA1) in which the relative rotation between nearest neighbors is always through a 180 • angle, but the direction of rotation alternates between two different axes in the (12) plane, one axis for the odd-numbered instantons and the other for the even-numbered. In SU(2) terms, for some A, B = 0 (A 2 + B 2 = 1).
• The green dots denote another non-abelian link-periodic pattern (NA2). Again, the relative rotation between nearest neighbors is always through a 180 • angle, but the direction of rotation alternates between two different axes. However, this time the two axes no longer lie within the (12) plane, thus where A, B, C all = 0 (A 2 + B 2 + C 2 = 1).
• Finally, the purple dots denote a confused non-abelian pattern NAX in which the nearest neighbors seem to differ by rotations through angles φ < 180 • around 2 alternating axes. In SU(2) terms, we see something like where A, B, C, D all = 0 (A 2 + B 2 + C 2 + D 2 = 1), but this picture is rather noisy: the coefficients A, B, C, D fluctuate from n to n and from run to run much stronger than their analogues in the other 4 patterns (5.10)-(5.13), and even the relative signs between A, B, C, D flip every few lattice sites -the long-range correlations are quite poor.
Physically, the disorder in the purple part of the phase diagram 5 indicates a nearby firstorder phase transition between incompatible orientation patterns. Near such a transition, our relaxation method fails to converge to a homogeneous phase that happens to have lowest energy for the (ǫ/D, M 3 /M 4 ) parameters in question. Instead, it leads to a mishmash of short domains of two nearly-degenerate phases; in fact, the domains are only a few lattice sites long, so they are significantly distorted by the boundary effects. The resulting mess is not even a local minimum of the zigzag's energy function, but once our simulation reaches this point, further evolution of the instanton orientations y n becomes very slow. Given very long time (and a better numeric algorithm) we would eventually see the domains of each phase growing longer and eventually merging into a single phase -but in our calculations we have simply run out of time and patience before this can happen.
The bottom line of our numeric analysis is that the instanton zigzag has 4 or 5 distinct orientation phases -we can not tell if NAX is a real phase or an artefact of mixed domains of other phases that we could not resolve. To answer that question -and also to get a better map of phase boundaries -we need to go back to analytic work.
Our starting point is the qualitative result of our numerical analysis: all 5 patterns (5.10)-(5.14) we have seen -or might have seen -can be fit into a single anzatz: for some angles φ, α, β. Indeed, the AF phase obtains for φ = π, α = π 2 ; the AB phase obtains for 0 < φ < π, α = π 2 , β = 0; the NA1 phase obtains for 0 < φ < π, α = 0, β = π 2 ; the NA2 phase obtains for 0 < φ < π, 0 < α < π 2 , β = π 2 ; and the NAX phase -if it exists at all -obtains for generic values of all three angles, 0 < φ < π, 0 < α < π 2 , 0 < β < π 2 . Given the ansatz (5.15), we may calculate the zigzag's net energy as an analytic function of the angles φ, α, β and parameters (ǫ/D, M 3 /M 4 ). The calculation is presented in Appendix A, and the result is for F = 3 2 + −1 + C 3 (cos 2 α cos 2 β + sin 2 α sin 2 β) × Σ 0 (φ) where In Appendix A we also minimize F -and hence the net energy per instanton -with respect to the orientation parameters α, β, and φ for fixed C 3 and ǫ/D. In particular, we show that the minimum never lies at generic values of the angles φ, α, β -which would correspond to the NAX phase. This means that the NAX phase does not really exist -its appearance in the purple region on the figure 5 is an artefact of poor convergence near a first-order phase transition (see below). are second-order. Likewise, the transition between the non-abelian NA2 and the abelian anti-ferromagnetic phase AF is also second-order.
On the other hand, the transitions between the other abelian AB phase and the nonabelian phases NA1 and NA2 are first-order. In fact, these are precisely the first-order transitions that confused our numeric analysis in the purple region of the figure 5 into appearance of the fifth orientation phase NAX.
Finally, the triple point at (M 3 /M 4 ) = 1, (ǫ/D) ≈ 0.38 between the AF, NA2, and AB phases is critical. At that point, the transitions between all 3 phases are second-order.
Indeed, in our previous paper [1] -where the analysis was limited to zigzags with M 3 = M 4 (in present notations) -we saw a second-order transition between the two abelian phases AF and AB.

⋆ ⋆ ⋆
Thus far, we have focused on the instanton orientations y n for a given zigzag geometry. That is, we have not only assumed that the instanton centers form a zigzag, but we have also treated the zigzag amplitude ǫ as an independent input parameter, just like the lattice spacing D or the M 2 , M 3 , M 4 parameters of the 5D gauge theory. Indeed, figure 6 shows the transitions between different orientation phases as function of independent, freely-adjustable parameters ǫ/D and M 3 /M 4 . But physically, the zigzag amplitude ǫ is a dynamical modulus whose value follows from minimizing the net energy of the multi-instanton system.
In fact, all instanton center coordinates X µ n are dynamical moduli, which raises two problems. First, in some situations the lowest-energy configuration of the instanton lattice may be more complicated than a zigzag or a straight chain; we shall address this issue in our next paper [41]. Second, for some lattice spacings D, a uniform lattice of any kind -a straight chain, a zigzag, or anything else -may be unstable against breaking into domains of different phases with different densities.
By way of analogy, consider a fluid governed by an equation of state such as Van-der-Waals. Formally, this equation allows a uniform fluid to have any density (up to some maximum). But in reality, at sub-critical temperatures one may have a low-density gas or a high-density liquid, but there are no uniform fluids with intermediate densities. If we constrain the overall volume V of some amount of fluid such that its average density would fall into the intermediate range, we would not get a uniform fluid; instead, part of the volume would be filled by the higher-density liquid while the other part by the lower-density gas. In the same way, if we fix the overall length L of the x 1 axis occupied by some large number N of instantons, their lowest-energy configuration is not necessarily a uniform lattice of some kind. Instead, for some average densities ρ = N/L we we would have L split into domains of two different lattices of different densities, one denser than N/L and the other less dense.
To keep the fluid uniform, one should control its pressure P rather than the volume V ; consequently, the preferred phase follows from minimizing the free enthalpy G = E−ST +P V rather than the free energy F = E −ST . Likewise, for the 1D lattice of instantons, we should control the net compression force F along the x 1 axis rather than the net length L or the lattice spacing D. Also, we should minimize the free enthalpy of the lattice G = E−ST +LF, but since we work at zero temperature all we need is the ordinary enthalpy H = E + LF. Equivalently, we may minimize the non-relativistic chemical potential (We focus on the non-relativistic chemical potentialμ = µ − M baryon because the relevant scale ofμ would be much smaller than the baryon mass. Indeed, according to eq. Thus, through the remainder of this section, we are going to impose a compression force F on the multi-instanton system, assume that the instantons form a uniform zigzag of some lattice spacing D and amplitude ǫ -or a uniform straight chain for ǫ = 0 -and vary D, ǫ and the orientation moduli (φ, α, β) of the zigzag to seek the minimum of the NR chemical potential (5.20) for any given combination of F and M 3 /M 4 .
We begin by changing variables from ǫ to ξ = ǫ/D, combining eqs. (5.20) and (5.17), and rewriting the result aŝ for F from eq. (5.18). Note that although we know F as an analytic function of all its arguments, and we know how to minimize it analytically with respect to the α and β angles (cf. Appendix A), the minimization with respect to the φ angle has to be done numerically, so we do not have an explicit formula for the F m (ξ; M 3 /M 4 ). Nevertheless, it is fairly easy to calculate numerically both the F m function itself and its derivative ∂F m /∂ξ.
Before we seek the global minimum of the chemical potential (5.21), let's find all the local extrema. Spelling out the extremality conditions we obtain and The last equation here has two branches of solutions: the zigzag branch ∂F m ∂ξ , (5.26) and the straight-line branch ξ = 0, any D > 0 (5.27) which follows from ∂F m /∂ξ = 0 for ξ = 0. Plugging these two solution branches into the compression force equation (5.24), we arrive at for the straight-line branch, for the zigzag branch. for the straight-line branch, (5.29a) and (F(ξ),μ(ξ)) for the zigzag branch (the red and green curves, according to the zigzag's orientation phase), thus The left plot here shows all the solutions while the right plot shows only the global minimum of the chemical potentialμ for any compression force F. We see that at F = F c the global minimum switches from one solution branch to another; physically, this corresponds to a first-order phase transition from a straight line to a zigzag.
For other values of the M 3 /M 4 ratio, we get more complicated (F,μ) plots than (5.30) For example, for M 3 /M 4 = 0.73 the zigzag branch goes back and forth several timesalthough this is very hard to see graphically without zooming on very small portions of the curve -so that for some values of the force F we get up to six different solutions of eqs. (5.28). Consequently, when we vary F there is a sequence of three phase transitions: First, a second-order transition from a straight line to a zigzag with the antiferromagnetic order of orientations. Second, a first order transition in which the zigzag amplitude drastically increases, the lattice spacing shrinks about 20%, and the orientation pattern changes from AF to the non-abelian NA2. Third, another first-order phase transition, this time with a smaller jump of the lattice geometry while the orientation pattern changes from the NA2 to another non-abelian phase NA1.   The stable straight-chain phase is colored pink, while the stable zigzag phases are colored red, yellow, blue, or green according to the instanton orientation pattern. Finally, the gray color denotes densities at which a uniform zigzag or straight chain would be mechanically or thermodynamically unstable.
Here are a few particularly noteworthy features of these diagrams: • Since we assume M 2 ≪ M 3 , the straight-chain phase always has the anti-ferromagnetic order of the instantons' orientations.
• The very first transition from the straight chain to a zigzag could be either first-order or second-order, depending on the M 3 /M 4 ratio: for (M 3 /M 4 ) < 0.725 the transition is first-order while for (M 3 /M 4 ) > 0.725 it's second-order. This difference is due to different orientation phases of the zigzag immediately after the transition: for (M 3 /M 4 ) > 0.725 the zigzag has the same antiferromagnetic order as the straight chain, which allows a second-order transition; but for (M 3 /M 4 ) < 0.725 the zigzag has a different orientation pattern -the non-abelian NA1 or NA2 -so the transition is first-order.
• The non-abelian phases NA1 and NA2 of the zigzag cover much larger areas of the phase diagrams 7(a,b) and 8 than the abelian phases AF and AB. In particular, at larger compression forces F -and hence larger chemical potentialsμ, larger densities, and larger zigzag amplitudes -the instanton orientations usually prefer the non-abelian patterns. Only the backgrounds with M 3 ≈ M 4 -such as the model we have analyzed in [1] -favor the abelian orientations.
• Figure 8 has gray areas at which an instanton zigzag with a uniform lattice spacing D = 1/ρ and a uniform amplitude ǫ (or a uniform straight chain for ǫ = 0) would be unstable against instantons' motion along the x 1 axis (the long direction of the zigzag). If we put N ≫ 1 instantons into a box of fixed length L = N/ρ and let them seek the lowest-energy configuration, they would organize themselves into domains of two different lattices with different lattice spacings and different amplitudes.
• The NA1 phase of the zigzag occupies two separate regions of the phase diagram separated by the region of the other non-abelian phase NA2. The phase transition between the lower-left NA1 region and the NA2 region is second-order, while the transition between the upper-right NA1 region and the NA2 is weakly first-order: the lattice spacing and the zigzag amplitude are discontinuous across the transition, but the discontinuity is very small and hard to see graphically on figure 8 or (5.30)-like plots.
• Likewise, the transition between the upper-right region of the NA1 phase and the abelian AB phases of the zigzag is weakly first-order. On the other hand, the transitions between the antiferromagnetic phases of the straight chain or zigzag and all the other zigzag phases is strongly first-order, with largish discontinuities of the lattice spacing and even larger discontinuity of the zigzag amplitude. But the zigzag is only the first step in the transition sequence from a 1D to a 2D instanton lattice. The subsequent steps -from a zigzag to a 3-layer lattice, to 4-layer lattice, etc.will be explored in our next paper [41].

Summary and outlook
In [1] we have investigated the crystalline structure of the holographic nuclear matter. The analysis done in that paper revealed the "popcorn transition" of a one-dimensional infinite chain of nucleons into a zigzag structure. We have also found out that the phases between adjacent nucleons are abelian. In this paper we continued the exploration of 1D instanton systems and their popcorn transition to the zigzag, but by generalizing the background to allow for M 3 = M 4 , we found a much richer phase structure, including two distinct nonabelian phases.
In §3 we proved that for low-density multi-instanton systems -in which the distances between the instantons are much larger that their radii -the interactions between the baryons are dominated by the two-body forces, while the 3-body forces, etc., are suppressed by powers of (radius/distance) 2 . The two-body-force approximation makes for much easier calculations, so in §4-5 we could put the instanton chains and zigzags into richer backgrounds than we had in our previous paper. In particular, we could use a less symmetric external potential for the baryons, by making the 5D gauge coupling vary as with three independent curvature parameters M 2 , M 3 , M 4 .
For the almost-holographic background M 2 , M 3 ≪ M 4 , we found that the 1D lattice of instantons has a highly degenerate family of ground states with many different patterns of instantons' orientations. Among them are periodic patterns in which instanton orientations span finite subgroups of the SU(2)/Z 2 : the anti-ferromagnetic chains spanning the Z 2 subgroup (two alternating orientations related by a 180 • twist), period = 4 chains spanning the Klein group Z 2 × Z 2 , and period = 2k chains spanning the prismatic groups Z k × Z 2 and the dihedral groups D 2k . There are also many link-periodic patterns in which the relative orientations y † n y n+1 of neighboring instantons are periodic but the orientations y n themselves are not periodic. But the vast majority of the degenerate ground states are not periodic at all. Raising the M 2 and M 3 parameters lifts the degeneracy. For For the zigzag-shaped chains of instantons we need M 2 ≪ M 3 ∼ M 4 to make sure the transition from a straight chain to a zigzag happens for lattice spacing ≫ instanton radius, but we may vary the M 3 /M 4 ratio. Consequently, we found 5 distinct phases shown in figures 7 and 8: The anti-ferromagnetic straight chain, the anti-ferromagnetic zigzag, the abelian link-periodic zigzag (period = 1), and two different non-abelian link-periodic phases (period = 2). Figures 7 and 8 also show a rather complicated phase structure: As we increase the compression force of the 1D lattice -and hence increase the instanton density and the chemical potential, -the zigzag amplitude and the orientation pattern go through a sequence of 1, 2, or 3 phase transitions, but the number of transition and their thermodynamic orders depend on the M 3 /M 4 ratio. In particular, the very first transition from a straight chain to a zigzag is first-order for M 3 < 0.725M 4 but second order for M 3 > 0.725M 4 . The subsequent transition(s) may also be first-order or second-order, depending on the M 3 /M 4 . This paper is part of a program focused on exploring the holographic nuclear matter, especially its low-temperature high-pressure phases. Thus far we have worked out the 1D instanton lattices and their zigzag deformations as a toy model for the transitions between 3D and 4D lattices. But many interesting questions remain open, and we would like to address them in our future papers; here is a short list: • The very next step in our program is to follow up the "popcorn transitions" from a 1D instanton chain into the second dimension beyond the zigzag phases through multiple lines of instantons towards a thick 2D lattice. We should also check that that the very first popcorn transition indeed goes from the straight chain to the zigzag rather than to some other configuration of the instantons.
• After that, we should set M 2 = 0 and study the infinite 2D lattices and their popcorn transitions into into the third dimension for M 3 ≪ M 4 . The techniques developed in this paper -using the two-body-forces approximation for a ≪ D -should work all the way up to infinite 3D lattices for M 2 = M 3 = 0. Alas, the ultimate popcorn transition from the 3D instanton crystal into the fourth space dimension probably happens for D ∼ a, so the two-body-forces approximation would not be accurate. But we ought to to try it anyhow, just to get a qualitative picture of the transition.
• Our program (starting with [1]) is based on the conjecture that the popcorn transition from a 3D instanton crystal to a 4D crystal is holographically dual to the transition of the large-N c nuclear matter from the baryonic phase to the "quarkyonic phase" [14] (Fermi liquid of quarks, with baryon-like excitations near the Fermi surface). Thus far, this conjecture was justified by qualitative arguments only, and we would like to make a quantitative argument by working out the details and the implication of the transition.
• In the skyrmion model of large-N c nuclear matter, Kugler and Shtrikman [19] showed that at low pressures, the lowest-energy configuration of is the FCC lattice of skyrmions with 4 different isospin orientations, while the high-pressure configuration is the simple cubic lattice of half-skyrmions. Later, Manton and Sutcliffe [13] described these lattices in terms of instantons on a T 4 torus. We would like to re-interpret their results in terms of instanton lattices to see where the Klebanov-Kugler-Shtrikman transition point (between whole-skyrmion and half-skyrmion lattices) corresponds to the popcorn transition from a 3D to a 4D lattice. Also, we would like to see how the chiral symmetry restoration in the high-density phase works in terms of the instanton lattices.
• Finally, we would like to see if any of our 1D, 2D, 3D lattices -or their generalizations -appear in condensed matter rather than in the holographic nuclear physics. We imagine a crystalline lattice of some high-spin atoms / ions / whatever whose orientations may form all kinds of complicated patterns -perhaps non-abelian patterns -depending on external parameters such as temperature, pressure, or doping. Hopefully, the methods developed in this article would be useful for studying such systems.
In this case, we would end up answering questions that have yet not been asked. and then we minimize that energy with respect to the orientation moduli α, β, and φ.

⋆ ⋆ ⋆
In the second part of this Appendix we minimize the zigzag's energy with respect to the orientation moduli φ, α, β. Specifically, we are going to minimize F (φ, α, β) in two stages: At first we hold φ fixed and minimize WRT α and β, and then we allow φ to vary and seek the overall minimum.
For the first stage it is convenient to change variables α → a = cos α and β → b = sin β, −1 ≤ a, b ≤ +1.
(A. 22) Then in terms of (φ, a, b) the F function becomes a polynomial in a and b with φ-dependent coefficients, where A(φ) = 3 2 − Σ 0 (φ) + 5 2 Σ 1 + 2Σ 2 (φ), Consequently, F becomes independent of b, but that's OK because the β angle does not affect the instantons' relative orientations y † m y n for φ = π, cf. eqs. (A.5) and (A.12). In other words, β or b becomes an irrelevant variable at φ = π -just like the geographic longitude becomes an irrelevant coordinate at latitudes ±90 • . As to the a variable, B > 0 means that the minimum is at a = 0, i. e., α = ± π 2 . Altogether, we have y † n y n+1 = ±iτ 3 (A. 26) which means the anti-ferromagnetic phase AF of the instanton orientations. The energy of this phase corresponds to Similarly, for φ = 0 we have Σ 0 = Σ 5 = 0 while Σ 2 = +Σ 1 and Σ 4 = +Σ 3 , hence This time, F becomes independent of a, but that's OK since the α angle becomes irrelevant to the instanton's relative orientations y † m y n for φ = 0, cf. eqs. (A.5) and (A.12). As to the b variable, C < 0 calls for maximal b 2 = sin 2 β, thus b = ±1 and β = ± π 2 . Altogether we obtain y † n y n+1 = ±iτ 1 (A. 29) and hence a new anti-ferromagnetic phase AF ′ . However, this new anti-ferromagnetic phase AF ′ has a higher energy than the old anti-ferromagnetic phase AF - To see this, consider any straight line through the center of the square. Along such a line the ratio a/b is fixed, so F becomes a quadratic polynomial of radius 2 = a 2 + b 2 with a non-positive coefficient of the (radius 2 ) 2 term. As a function of an un-bounded real variable such a polynomial does not have any minima at all. But since the radius 2 variable is bounded -between zero at the center and the maximum at the edge of the square -F may have a minimum (or minima) at the center or/and at the edge. Applying this argument to all lines through the center, we find that the minimum or minima of F may lie at the square's boundaries a = ±1 or b = ±1, or at the center a = b = 0, but not at any other interior point.
Let us consider such possible minima in more detail: • The minimum at a = b = 0 -i. e. α = ± π 2 , β = 0 -leads to y † n y n+1 = cos and hence the abelian AB phase of the instanton orientations. Local stability of this minimum requires C > 0 and BC > J 2 , and its energy corresponds to F (AB) = A(φ). (A.33) • Minimum at the boundary b = ±1 (but generic −1 < a < +1) -i. e. β = ± π 2 and generic α -leads to y † n y n+1 = (−1) n cos φ 2 × (±iτ 1 ) + sin φ 2 × cos α × (iτ 2 ) + sin α × (iτ 3 ) (A. 34) which corresponds to the non-abelian orientation phase NA2. The specific location of such a minimum is b = ±1, a = ± J B − K (A. 35) and its local stability at |a| < 1 requires The energy of this minimum is . (A.37) • In principle there could be a similar minimum at the other boundary a = ±1 (but generic b). Specifically, such a minimum would happen at a = ±1, b = ± J C − K (A. 38) and its local stability at |b| < 1 would require However, thanks to B > C these two requirements contradict each other, so this type of minimum does not happen.
The phase diagram on figure 6 (on page 67) was obtained by repeating this numerical