Many Phases of Generalized 3D Instanton Crystals

Nuclear matter at large number of colors is necessarily in a solid phase. In particular holographic nuclear matter takes the form of a crystal of instantons of the flavor group. In this article we initiate the analysis of the three-dimensional crystal structures and the orientation patterns for the two-body potential that follows from holographic duality. The outcome of the analysis includes several unexpected results. We perform simulations of ensembles of O(10000) instantons whereby we identify the lattice structure and orientations for the different values of the weight factors of the non-Abelian orientation terms in the two-body potential. The resulting phase diagram is surprisingly complex, including a variety of ferromagnetic and antiferromagnetic crystals with various global orientation patterns, and various"non-Abelian"crystals where orientations have no preferred direction. The latter include variants of face-centered-cubic, hexagonal, and simple cubic crystals which may have remarkably large or small aspect ratios. The simulation results are augmented by analytic analysis of the long-distance divergences, and numerical computation of the (divergence free) energy differences between the non-Abelian crystals, which allows us to precisely determine the structure of the phase diagram.


Introduction and Summary
It is believed that real life N c = 3 QCD cold nuclear matter forms a quantum liquid, but for large N c it becomes a crystalline solid [1]. This structure of nuclear matter has been intensively investigated in the context of Skyrme theory [2][3][4]. A more modern approach to nuclear matter at large N c is based on the use of holography. In holography mesons are open string configurations with endpoints on flavor branes. Baryons have the structure of a baryonic vertex connected with N c strings to flavor branes. It was shown in [5] that the baryonic vertex can be realized in terms of a wrapped D p brane over a p-cycle with a RR flux of N c . Such structures have been analyzed in a confining holographic background in [6]. A prototype model that admits this structure is the so called Witten-Sakai-Sugimoto model (WSS) [7,8]. The holographic dual of the pure glue theory is based on near horizon limit of D 4 branes compactified along one space dimension with anti-periodic boundary conditions for the fermions. The quark sector is mapped into stacks of D 8 andD 8 probe flavor branes that are asymptotically separated and merge in the IR region. This model was then generalized to include non-antipodal flavor branes [9]. It was shown in [10], using energy considerations, that in that model the baryonic vertex is immersed in the flavor branes. In such a case the baryon can be viewed as a U (N f ) instanton [11] in a curved Euclidean 4D built from the three ordinary space dimensions and the holographic dimension. At leading order in the strong coupling limit these instantons can be approximated by the BPST instantons. Thus holographic nuclear matter takes the form of a lattice of BPST instantons. A key ingredient for the structure of the latter is the interaction between the instantons. It was shown in [1] that the interaction in the WSS model is repulsive at any separation distance. For long distances between the baryons the repulsion is due to the fact that the lightest isoscalar vector, whose exchange yields repulsion, is lighter than the lightest scalar that yields attraction. In the near and intermediate zones the interaction between two instantons of the WSS model is purely repulsive. By using the generalized WSS model instead of the original model, the severeness of the problem can be reduced. As was shown in [1], only in the generalized model (and not in the original WSS) there is, in addition to the repulsive force, also an attractive one due to an interaction of the instantons with a scalar field that associates with the fluctuation of the embedding, though the ratio of the attractive to repulsive potential can never exceed 1/9. In [12] it was shown that there is another holographic model [13] with a dominance of the attraction at long distances, together with a small binding energy, as one encounters in nature.
Thus, the crystalline of holographic nucleons has to be a multi-instanton solution of the ADHM equations. Motivated by the interesting behavior of skyrmions at high density and being equipped with the new methodology of holography, first steps in determining the structure of the cold nuclear matter were made in a context of a 1D toy model in [14] and [15]. For the straight 1D lattices, it was found that the orientations are 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. In addition we discovered a formation of zigzagshaped lattices, where instantons are popping up to the holographic direction. For these configurations we detected 4 distinct orientation phases -the anti-ferromagnet, another Abelian phase, and two non-Abelian phases. We further discussed the phase diagram of 2D instantons [16].
Related studies appeared in . In particular, [18] and [19] have analyzed the phase diagram of configurations with finite baryon density, treating the instantons as point-like sources. A completely different approach was taken in [19,[40][41][42], where the instanton structure was replaced by a homogeneous non-Abelian field, motivated by dense instanton configurations smeared over spatial directions.
In analogue to the 1D toy model, realistic holographic models of nuclear matter are expected to involve 3D and 4D crystals and a transition between them with increasing density. In this article, we will only discuss the 3D crystals. But solving even the 3D ADHM equations is a tremendous task. Instead in this paper we use an approximation of a two body instanton interaction which is adequate for crystals of low density. In [15] and in [24] the two body potential between two instantons was derived. We write down a more general potential depending on two parameters α and β, E 2 body (m, n; α, β) = 2N c (1 + 2α + 2β)λM 1 |X n − X m | 2 × 1 2 + α tr 2 y † m y n + β tr 2 y † m y n (−i N mn · τ ) , (1) such that the potential discussed in the earlier reference is obtained for α = β = 1. Here in this research work we generalize the potential by letting α and β vary. The notation in this potential will be explained in detail below. At the moment it is enough to know that the SU (2) matrices y m and y n carry the information on the orientations of the instantons n and m. In this way we consider different weights for the orientation terms in the potential.
In this paper we determine the crystal structure and the instanton orientation patterns that minimize the total energy of the system. In order to understand our results, it is important to notice that the two-body interaction (1) gives rise to a strong long-distance (IR) divergence in three dimensions. The potential due to a configuration of (large) size R behaves as drr 2 E 2 body (r) ∼ dr r 0 ∼ R so it is linearly divergent. Therefore long distance effects are strongly enhanced, leading to drastically different results in comparison to the results for lower dimensions [15,16]. In particular, much of the phase structure is determined through the comparison of the coefficients of the divergent term for various orientation structures. In several of these phases, minimization of the energy leads to nontrivial long distance correlations between the instanton orientations: for example, we obtain orientation structures which are spherical at long distances.
Because of the divergence we need a long-distance cutoff in all our simulations and computations, so we need to be particularly careful that our results are independent of the cutoff. Moreover, in case of simulations, the long distance interactions lead to clustering of the instantons at the surface of the simulation volume if we only set a hard wall cutoff which forces the instantons to stay within the volume. The removal of this undesired effect necessitates the use of a smooth external force which pushes the instantons towards the center of the simulation. Because of the long-distance divergence, however, there is no obvious "correct" choice for the external force for all configurations that we consider. A simple choice which is applicable to all configurations that we encounter is to take a force which sets the (locally averaged) instanton density to be constant. For most of the phase diagram this force matches with the (regularized) force due to a homogeneous density of instantons outside the simulation volume, i.e., the force due to the instantons left out of the simulation. Notice also that all our simulation results turn out to be insensitive to the precise choice of the force.
We are now ready to summarize our results. The results are obtained by using three different methods of analyzing the lattices of instantons: (i) We perform simulations of ensembles of instantons of the order of O(10000) instantons subjected to the two body interactions between any two of them. In addition, in order to reduce boundary effects, we also added an external force which, as discussed above, sets the average density of instantons to be constant. Using these simulations were were able to determine the geometry and orientation structures of the instanton collection with the lowest total energy.
(ii) We determine the orientation patterns based on the behavior of the two body potential for far apart instantons, i.e., the long-distance divergence discussed above. In this range of distances between the instantons we take the continuum limit and ignore the lattice geometry. We compute the total energy of the system as a function of α and β and in particular determine the associated lowest energy configuration.
(iii) We compute the energy difference between various pairs of geometrical and orientation structures by which we eventually determine the phase diagram. The computations of energy differences yield finite results since the divergences are cancelled out in the differences.
The main outcome results of the simulations are the following for the basic cases: • In the non-orientation case α = β = 0 the crystal structure is face-centered-cubic (fcc).
• The basic oriented case α = β = 1 is a face-centered tetragonal lattice with a large aspect ratio, 1 i.e., fcc with one direction rescaled, breaking the cubic symmetry. The aspect ratio c is large: we find c ≈ 2.467. That is, the instantons form clearly separated layers with two-dimensional square lattice structure. A sample of this structure is shown in fig. 2. The two dimensional layers have antiferromagnetic structure. The orientations between the layers repeat in cycles of two, as seen from fig. 2. Overall, there are therefore four distinct and linearly independent orientations so that the set of orientations does not single out any direction. We call this class of orientation structures "non-Abelian".
We have also revealed the structure of the phase diagram of the ensembles of instantons as a function of α and β in (1). The results are mostly based on simulations, with the boundaries between the phases refined by using the methods (ii) and (iii). We identify the following phases: • Oriented non-Abelian phase. This phase contains crystals in the same class with the case α = 1 = β, i.e., the orientations span the whole four dimensional orientation space. We find this phase in the region given by the conditions α > β, β > −1/8, and α > γ 1 β with γ 1 ≈ 0.0575. The phase is further divided into subphases having different lattice and orientation structures, which can be classified in the following classes: (a) Tetragonal/cubic (fcc or fcc-related) lattices with "standard" orientation pattern (b) Tetragonal/cubic (fcc or fcc-related) lattices with "alternative" orientation pattern (c) Hexagonal lattices (d) Simple cubic/tetragonal lattices.
In the classes (a), (b), and (d) the crystals have patterns of four distinct orientations, whereas the hexagonal lattices have six orientations (details will be given in Sec. 5). We find crystals in all of theses classes also with nontrivial aspect ratios, i.e., with the spatial structure rescaled in one coordinate directions. We also identify regions in the phase diagram where the crystal is exactly fcc (aspect ratio is equal to one) within the class (a), and regions where the crystal is exactly simple cubic, class (d). Moreover we find a line where the crystal is exactly body-centered-cubic (bcc) for the class (b) (aspect ratio is 1/ √ 2 in fcc units). The two basic cases (non-oriented α = 0 = β and oriented α = 1 = β) are included in the non-Abelian phase as special limiting cases.
We find that the orientation structure is locally antiferromagnetic: in domains of small size, only two different orientations appear. Globally the orientations are arranged spherically. The detailed description of the global structure will be given below. The spatial lattice has a strong tendency of forming spherical shells. The two-dimensional lattice structure depends on α and β. For α = O (1) and β = O (1), the spherical shells have the structure of two-dimensional square lattice (with antiferromagnetic orientations). As α and β decrease, the lattice changes to rhombic. For small values, α = O (0.1) and β = O (0.1), the layers are closer to triangular and the spatial structure is locally close to bcc.
• Spherical ferromagnetic phase. In the sector β > Γα > 0, we find a pattern which is locally fcc and ferromagnetic: all nearby orientations are (almost) the same. However globally the orientations form a spherical structure. The simulations indicate that the lattice structure is independent of α and β within the phase.
• Ferromagnetic egg phase. In the remaining corner of the phase diagram, i.e., when γ 2 β < α < γ 1 β and β > −1/8, we find a special ferromagnetic crystal. It is again locally fcc, and the interior of the sample orients as a global ferromagnet, but there is crust of finite width where the orientations obey a global spherical alignment.
The methods that we have developed in this paper for the analysis of holographic nuclear matter, namely, the long distance analysis of the potential, determining the stable crystalline and orientation structure by performing simulations and the computations of energy differences, could be applied to other non-Abelian lattices.
The paper is organized as follows: In the next section we spell out the basic setup of the lattices of instantons. In particular we write down the two body interactions between the instantons, their generalization and the symmetries of the system. In section §3 we describe the outcome of the simulations. We start with the unoriented case. Next we present the results of the standard orientation of the unmofified model of α = β. Oriented non-Abelian crystals with α ≥ β come next. For the range of 0 < α < β an oriented (anti) ferromagnetic spherical phase was found and finally global ferromagnetic phases at α < 0. Section §4 is devoted to the computation of the potential and total energy at long distances for the different types of orientations. The phase diagram of the non-Abelian crystals is determined in section §5. We conclude and write down several open questions in §6. Two appendices were added. In the first we determine the potential for a ball of homogeneous matter and in appendix B we present numerical details of the setup of the simulations and on the computations of the energy differences.

The basic three dimensional setup
We consider a three-dimensional system of point-like SU (2) instantons so that all instantons are located at the same value of the holographic coordinate. Therefore each instanton is characterized in terms of the position X = (X 1 , X 2 , X 3 ) and the unimodular quaternion (y 1 , y 2 , y 3 , y 4 ) carrying the information on the orientations of the instantons. The basis vectors in the quaternion space are denoted as usual: 1 = (1, 0, 0, 0), i = (0, 1, 0, 0), j = (0, 0, 1, 0), and k = (0, 0, 0, 1). We also use the representation of the quaternions as a 2×2 matrices, where τ j are the Pauli matrices. Then unimodularity is equivalent to the condition tr(y † y) = 2 4 k=1 (y k ) 2 = 2. The ensemble is described in terms of the pairs ( X n , y n ) with n = 1, 2, 3, . . . N .

Two body forces between instantons
In holographic nuclear physics, like in real life, besides the two-body nuclear forces due to meson exchanges, there are significant three-body and probably also higher number-body forces, where n stands for the quantum numbers of the n th nucleon.
In the low-density regime of baryons separated by distances much larger then their radii, the two-body forces dominate the interactions, while the multi-body forces are smaller by powers of (radius/distance) 2 . The dominance of the two body forces for low density was proven in [15]. The two body interaction was also determined there by combining the non-Abelian and the Coulomb energies of the system and re-organizing the net energy into one-body, two-body, etc, terms.
We use a two-body potential motivated by the Witten-Sakai-Sugimoto model [8]. Assuming that the separation |X m − X n | between the instantons n and m satisfies where λ is the coupling constant and M is the Kaluza-Klein scale, the two-body potential in this model takes a simple form [15,24]: where we defined the unit vector The lower limit in (4) arises from the instanton size and the upper limit from neglecting curvature corrections in the AdS geometry. Notice that the range of validity is parametrically large in the strong coupling limit λ → ∞, where the dual description in terms of classical gravity is reliable. The normalization in (5) was chosen such that taking the average over the orientations of either of the instantons (with the obviously chosen measure) gives 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. The total energy of the system (including only the dominant two-body terms) is the sum over all pairs:

Generalizations of the defining instanton interactions
It is natural to consider a generalization of the two-body potential which includes both the oriented interaction (5) and the orientation independent potential of (7). We therefore define where the normalization was chosen such that after averaging over the orientations, (7) holds. Notice that α = 0 = β gives the unoriented potential and α = 1 = β gives back (5). Let us then comment on how the parameters α and β affect the interaction in physical terms. The α coefficient multiplies a term where the dependence on orientation is independent of the spatial interactions. For positive (negative) α perpendicular (parallel) spins of nearby neighbors are preferred, and the effect increases with increasing |α|. The interaction involving β is a bit more complicated since there is nontrivial coupling between the orientations and directions in coordinate space. Picking an instanton with unit orientation (y = 1), positive (negative) β means that the orientation y of a neighboring instanton which is perpendicular (parallel) to the spatial link between the two instantons is preferred. (This with the understanding that the orientations are mapped to the directions in the coordinate space in the standard fashion: for example, the orientation y = i is parallel to the x axis in coordinate space.) We will comment on how these observations are reflected in the structure of the final phase diagram in Sec. 5.

Symmetries of the interaction
The two body potential (5) and its generalizations (9) for various values of α and β are invariant under several transformations.
• Since the potential dependence on the space coordinates is via X n − X m and N ij it is obvious that it is invariant under space translation of the locations of the instantons; for any constant three dimensional vector A.
• Whereas the non-orientation part and also for the case of β = 0 the potential is invariant under global rotations. Since N ij is not invariant under rotations the full potential is invariant only under rotation which is accompanied with a rotation in the y space described below.
• Changing the sign of either y m or y n does not have an effect on the energy. That is, the orientations matter only up to a sign.
• The y n s reside on an S 3 and hence are invariant under SU (2) L × SU (2) R symmetry transformations related to the orientation structure, one of which is (in general) coupled to the spatial rotations. First, notice that the interactions only depend on the combination y † m y n which we call the twist. Therefore, obviously a left multiplication of all quaternions y n by a unimodular quaternion u leaves the potential invariant: • The other SU R (2) rotation is coupled to spatial rotations as follows. Under rotation by an angle θ around the axis n, the orientation dependent factor in (9) transforms as We therefore identify the action of the rotations to the quaternions as the following right multiplication: which, when applied together with the spatial rotation, leaves the energy invariant.
• In the special case β = 0, spatial rotations are decoupled from the orientations; the vectorial transformation (13) and spatial rotations are good symmetries separately.

Simulations of the various instanton lattices
We have studied the three dimensional crystal structures arising from the above interactions by using numerical simulations. Starting from a random initial condition with N = O(10000) of instantons in a spherical cavity, we use an algorithm to relax the system towards the state of lowest energy. In more detail, we use a setup which roughly corresponds to adding masses and drag forces for the instantons, and simulating their dynamics as the system converges to a final state with low energy. We also add a compensating "external" force, which (as we explained in the introduction) prevents the instantons from clustering at the edge of our simulation volume. To be precise, we choose a force which sets the density of instantons, averaged locally over small distances, to be constant (see Appendix A). This means that we compute the force field generated by the sample in the continuum limit, i.e., by approximating the instanton configuration by homogeneous matter of constant density, and apply the force of the same magnitude but opposite direction on the instantons during the simulations. We note that for many of the crystals which we find (including the global ferromagnetic and non-Abelian crystals defined below) the force that sets the instanton density to constant is apparently the unique natural choice for the compensating force. This is true for the crystals where the orientation structure is simple in the continuum limit; see Appendix A.2 for a detailed analysis. We have also tried modifying the external force, and the crystals found in the simulations are to large extent insensitive to such modifications. Moreover, convergence to crystals of good enough quality for our purposes typically requires adding a pseudo-temperature in the form of random kicks to the instantons that would slow down the relaxation process. The setup for the simulations is discussed in more detail in Appendix B.
Notice that the simulation results contain finite size effects which may not be highly suppressed even for 10000 instantons. In particular, compared to simulations in lower dimensions [15,16], the simulation volume is effectively smaller compared to the boundary: For example, the distance of a random instanton to the boundary is ∼ N 1/d where d is the dimension, and is reduced as d increases. This enhances the dependence on boundary conditions and other boundary effects. This is seen in the simulation results, e.g., as the emergence of unwanted layer structures which are aligned with the surface of the volume of the simulation, and only appear near the surface. Therefore it is important to check that the results are independent of the number of instantons included in the simulation. Moreover the choice of initial conditions (even if they are random) as well as the pseudo-temperature may cause bias towards some of the crystals we find in the results. It is also important to compare the results with those obtained by using other methods in sections 4 and 5. We typically find good agreement, as we will discuss below.
Originally, we studied the standard unoriented case, given by the potential (7), and the standard instanton interaction, given in (5). Then it is natural to also consider the interaction that interpolates between these two, i.e., the choices with α = β in the more general interaction of (9). However, as it turns out, the line α = β is a critical line, that is, a phase boundary on the (α, β) phase diagram determined by the general interactions. Therefore we decided that it is best to extend the studies to the whole (α, β) plane.
There are however some quite obvious constraints on the values of α and β. Namely, for a pair (n, m) of instantons there are orientation choices where the factors tr 2 (· · · ) in (9) take independently any possible value, i.e., a value within the range from zero to four. That is, if α < −1/8 or β < −1/8 the factor in the square brackets can become negative, so that the interaction turns attractive (at all distances). Such an interaction will lead to a collapse of the system, and therefore must be excluded from the phase diagram. That is, we restrict our studies to the region α > −1/8 and β > −1/8. In addition, we have decided to exclude from our analysis the potentials where α or β is much larger than one, so that the orientation dependent terms would dominate the interaction.
We will then discuss our simulation results. We start by the unoriented and oriented standard cases which are given by the points α = 0 = β and α = 1 = β on the phase diagram, respectively. Then we will describe the results in the rest of the plane, taking into account the constraints discussed above.

Unoriented
Let us start from the unoriented potential, given in (7). The convergence of the simulations is far from optimal for the unoriented potential: we obtain fcc structures, but domains are very small. However clear domains of any other regular structures are rare. In some instances, chunks of crystal closer to bcc than to fcc were found. This suggests that fcc and bcc are very close in energy, which may explain in part why it is difficult to identify the lattice with lowest energy. We will confirm this in Sec. 5, and show that the fcc structure indeed has the lowest energy. A sample of a simulation result is shown in fig. 1. This sample is a carefully chosen region of a simulation with 10000 instantons. It represents roughly the cleanest fcc structure that we could obtain. Notice that this figure, as well as the other figures in this section, may be viewed as an interactive 3D version if the pdf is opened in Adobe Reader or Adobe Acrobat (version 9 or higher).

Oriented "standard" instanton interactions (α = β = 1)
Let us then discuss our results for the original interaction (5). In this case, the crystal with lowest energy could be identified with high certainty, and the result is somewhat surprising. Namely, the spatial structure is a face-centered tetragonal lattice with a fairly large aspect ratio, i.e., fcc with one direction rescaled, breaking the cubic symmetry. The lattice can be equally described as body-centered tetragonal, which is related to the face centered version by a rescaling of the aspect ratio by √ 2 (and a rotation by 45 degrees of the unit cell). The aspect ratio c is indeed large: we find 2 c ≈ 2.35 (using the face-centered conventions). That is, the instantons form clearly separated layers with two-dimensional square lattice structure. A sample of this structure is shown in fig. 2. This sample was taken from a well-converged simulation with 5000 instantons.
The orientation structure could also be extracted reliably. The two dimensional layers have anti-ferromagnetic structure. The orientations between the layers repeat in cycles of two, as seen from fig. 2. If we choose the unit cell of the lattice to be aligned with the axes such that z is the asymmetric direction, the orientations can be chosen to be 1 and k in one set of the layers, whereas the other set has the directions i and j (where we used the standard (1, i, j, k) basis for the quaternions). Recall that the orientations are determined only up to signs and left multiplication by a constant unimodular quaternion. These four orientations are shown as different colors in fig. 2.
Finally, we remark that the standard instanton interactions lie on a special line of the phase diagram of general interactions, which depends in the parameters α and β. That is, Figure 2: A sample of a simulation result for the oriented two-body interaction of 5d instanton dynamics (i.e. with α = β = 1). The lattice is tetragonal with non-Abelian orientation structure. The colors show the four different orientations of the instantons. The figure was created by using the Jmol software [45]. An interactive 3D version of the figure can be viewed with Adobe Reader.
there is a phase transition on the line α = β on the (α, β)-plane, and above this line the structure is drastically different (as we will detail below) from the result of fig. 2 while below the line they are similar. The fact that the results for α = β are continuously connected to the lower phase may be dependent on finite size effects and other details of the simulation (e.g., the pseudo-temperature). Actually in simulations we find that the critical line is at small positive β − α. Despite the point being on the critical line, the convergence on this line is much faster than for the unoriented case and the resulting crystal have typically good quality.

Oriented non-Abelian crystals (α ≥ β)
Let us then consider the simulation results for the generalized interaction (9), varying values of the parameters α and β. We will here only sketch the phase diagram based on the simulation results, and the final result for the diagram will be presented in Secs. 4 and (5) by using the other methods. The orientation of fig. 2 structure falls in a wider class which we call "non-Abelian" crystals. In this class, there is no two-dimensional subspace which would contain all the instanton orientations. The set as such does not single out any preferred direction. We start by exploring the region where non-Abelian crystals are found. Based on the simulation results, this means that α β, α 0, and β > −1/8.
The phase of the non-Abelian crystals may be further divided into smaller subphases. We observe that the simulation results show most variability near "critical" lines given by α = β, β = 0, and β = −1/8. Outside these critical regions, we find plain fcc when β > 0 and plain simple cubic when β < 0. Examples of simulation results for both of these lattices are shown in fig. 3. The orientation pattern for the fcc lattice is the natural one with four orientations: it is actually the same as in fig. 2, but just the aspect ratio is much smaller, equal to one in this case. For the cubic lattice the orientation pattern is also the natural pattern. That is, for both these lattices the way the orientations are distributed on the instanton lattice does not single out any direction, so that the symmetry of the lattice protects the value of the aspect ratios (which equal one for both lattices). In other words, one can choose a set of three orthogonal axes such that the lattice is invariant under 90 degree rotations around all the axes.
Near the critical lines, however, other kind of lattices are found. We first discuss the line with β = α which has the richest phase diagram. As we remarked above, this line is actually the boundary of the non-Abelian phase, but due to finite size (or other simulation) effects our simulations converge to non-Abelian crystals also on this line.
For 0.5 α < 1 the structure is the tetragonal lattice described in Sec. 3.2 and shown in fig. 2, except that the aspect ratio slightly decreases with decreasing α. At β = α ≈ 0.5 we encounter a first order phase transition to another kind of lattice, which appears to be dominant for 0.2 α 0.5. This is a hexagonal lattice, which has the hexagonal close packed (hcp) structure but smaller aspect ratio (in the direction perpendicular to the 2D layers with triangular structure). There are in total six different orientations, which alternate between the 2D triangular layers. Let us take the layers to be aligned with the xy plane and separated in the z direction. Then one set of the triangular layers has orientations in the plane spanned by 1 and k, which are related to each other by rotations by 120 degrees. The layers with different orientation structure have orientations in the plane spanned by i and j, similarly related through rotations by 120 degrees in the plane. We give the precise description of the structure in Sec. 5. A sample of this lattice is shown in fig. 4 (top left), which is taken from a simulation with 3000 instantons. Notice that for this lattice, the simulation result is particularly clean for the spatial structure (while there are some minor defects in the orientations). We used a coloring scheme which is different from the other plots to show all the six orientations.
At β = α ≈ 0.2 the simulations suggest that there is another phase transition, interestingly, to another fcc related structure with different value for the aspect ratio, which is now c ∼ 0.7. That is, the lattice is close to bcc. A sample of this lattice is shown in fig. 4 (top right), and it is related to those in fig. 2 and fig. 3 (left) by rescaling of one of the coordinates.
We have also checked what happens for β = α > 1. Simulations indicate that there is at least one phase transition: instead of layers with anti-ferromagnetic 2D square lattices, one obtains layers of anti-ferromagnetic hexagonal lattices at high values of α, see fig. 4 (bottom left). It is difficult to pinpoint the value of the phase transition precisely based on the simulation results, which often have domains of both lattices, but it appears to be near α = β = 2. Let us then discuss the simulation results near β = 0. We encounter another structure, which is shown in fig. 4 (bottom right). The spatial structure is close to bcc for β = 0 exactly, but the orientations are arranged in a different way than in fig. 4 (top right). When β = 0 but small, we observe that the aspect ratio varies relatively fast (increasing with β).
The remaining critical line is β = −1/8. Our simulations in this region converge slowly, but suggest the existence of simple tetragonal phases, i.e., lattices of fig. 3 (right), but with nontrivial aspect ratio. We will return to this in Sec. 5. Figure 5: anti-ferromagnetic spherical lattices found for 0 < α < β. Left: A sample with two-dimensional square lattice layers from a simulation with 15000 instantons at α = 2/3 and β = 1. Right: A sample with two-dimensional rhombic layers from a simulation with 15000 instantons at α = 0.1 and β = 0.15.

Oriented (anti) ferromagnetic spherical phase (0 < α < β)
Next we discuss the results in the region 0 < α < β. In this region the results are drastically different from those at α > β: there is a very strong tendency for the orientations to align with the position of the instanton, which leads to a global spherical orientation structure. We note that simulations in these "spherical" phases converge to the final state over ten times faster than for the non-Abelian phases found for β ≤ α which suggests that the spherically arranged ground states are energetically strongly favored. The simulation results agree with analytic analysis of the long distance behavior we will carry out in the next section: in all simulations for 0 < α < β, we find that the orientations are arranged spherically. The spherical phase further divides into two regions, where the orientation structure is either locally antiferromagnetic (constraining to small regions, two orientations appear, with closest distance links being between different orientations) or ferromagnetic (in small regions, all orientations are aligned). We verify that there is a phase transition from (local) anti-ferromagnetic to (local) ferromagnetic order at β = Γα with the value of Γ between 2.2 and 2.3. This is close to the analytic prediction Γ ≈ 2.14 from Sec. 4, and the small deviation is probably a finite size effect and/or a bias due to the simulation setup.
In the anti-ferromagnetic phase we observe a very strong tendency of the system to converge to almost exactly spherical shells. While there is some tendency to form spherical shells also in the other phases, at least near the boundary of the spherical cavity in our simulations, in the anti-ferromagnetic phase this effect is clearly stronger than in the other phases. That is, it is plausible that the spherical shells observed in the other phases are essentially boundary defects, but in the anti-ferromagnetic phase, even for the largest simulations with 40000 instantons which we have run, the whole system converges rapidly to easily recognizable spherical structures.
Within the anti-ferromagnetic phase, the simulation result also strongly depends on α and β. For large values, α 2 + β 2 0.5 2 , the spherical layers are clearly separated and form two-dimensional anti-ferromagnetic square lattices, see fig. 5 (left) for an example. Due to finite size effects it is difficult to pin down the preferred local crystal structure, but we conjecture that it is simple tetragonal so that the different layers prefer to be precisely aligned with exactly same orientation pattern. The separation between the layers slightly grows with α and β. Closer to the origin on the (α, β)-plane, the layers are rhombic instead of exactly rectangular, see fig. 5 (right). At very small values of the parameters, i.e. for α 2 +β 2 0.1 2 , the two-dimensional layers are closer to triangular, and the distance between the layers is comparable to the distances between nearest-neighbor instantons within the layers. The crystal structure is locally close to bcc, while the orientations remain antiferromagnetic. It is difficult to pin down the precise preferred orientation structure in this region as the orientation dependent terms in the interaction potential are weak.
We notice that the simulation results in the anti-ferromagnetic phase therefore have some resemblance to those in the non-Abelian phase at α = β: At α 2 + β 2 ∼ 1 we find clearly separated square lattices in both phases, and at small α 2 + β 2 we find bcc-like structures in both phases. Hexagonal lattices, which are found in the non-Abelian phase, are however absent in the anti-ferromagnetic phase.
Finally, we have also carried out simulations in the ferromagnetic phase (β > Γα > 0). As we pointed out above, the results are in agreement with the analytic long distance analysis: the orientations are spherically aligned. There is some tendency of forming spherical structures but this is considerably weaker than in the anti-ferromagnetic phase. It is again difficult to determine the local crystal structure with high certainty but it is plausible that it is fcc: locally the simulation results are similar to fig. 1. The simulation results are essentially independent of α and β within this phase.
Notice also that, as we remarked above, for the spherical crystals it is not obvious what one should choose as the compensating external force (see detailed discussion in Appendix A.2). In order to make sure that our results are consistent, we have also run simulations for different choices of potentials. For example, instead of the potential which enforces constant instanton density, we have used potentials which are due to an IR regulated potential obtained by integrating a continuum of spherically arranged instantons outside the cavity of simulation. Changes in the external force did not lead to any changes in the local crystal structure, or remove the global spherical arrangement of the orientations. Furthermore, since the spherical structures are surprising, we have checked that the results are insensitive to the boundary conditions by varying the size of the ensemble and the shape of the cavity of simulation. We were able to carry out simulations with up to about 40000 instantons without any essential change in the results. Naturally, we cannot exclude that some changes occur for even larger ensembles.

Global ferromagnetic phases (α < 0)
The line of α = 0 marks another phase transition. Namely, for negative α (and with α > −1/8), we find a simple global ferromagnetic structure, where all instantons have the same orientation to a good precision. The two-body potential is then proportional to that of the same lattice is then fcc, and simulation results are similar than that shown in fig. 1.
There is however a small corner in the region where both α and β which shows highly unexpected behavior: the orientation structure is not uniform. The region where this happens is given by 0.5β α 0 and β > −1/8. Most of the simulation result obeys a global ferromagnetic structure, but there is a crust that is locally ferromagnetic but globally spherical. The crust is relatively thin (typically having thickness of about 10% of the radius of the simulation) but it scales with the system size and therefore is present even in the continuum limit. We discuss this structure in more detail in the next section.

The instanton energy at long distances
In this section we calculate the net self-interaction energy in a large spherical piece (radius R) of the instanton matter of density ρ. Due to linear IR divergence of 1/distance 2 twobody potential in 3D, the net IR energy -which is independent on the instanton lattice geometry -scales like E IR ∼ ρ 2 R 4 , much larger than the geometry-dependent effect whose net contribution scales like ∆E ∼ ρ 2 R 3 a where a is the lattice spacing. In this section we focus on the leading IR effects only, so we take the continuum limit a → 0 and ignore the lattice geometry -it can be cubic, hexagonal, tetragonal, whatever, it does not matter.
However, unlike the lattice geometry, the orientation pattern of the instantons does affect the IR energy of the lattice. Or rather, what matters is the spectrum of instanton orientations in each small but macroscopic volume d 3 x of the instanton matter; but the IR limit does not care for the specific locations of those orientations within the instanton lattice. The spectrum of orientations matters because the two-body force between the instantons depends on their orientations. For convenience we set the prefactor of eq. (1) to one.
so we have where y 1 and y 2 are the two instantons' orientations in SU (2) notations and n 12 is the unit vector in the direction of x 2 − x 1 . For our purposes, we treat the two instantons' locations x 1 and x 2 as continuous variables, and average the Q 12 over the orientations of instantons found within small but macroscopic volumes d 3 x 1 and d 3 x 2 . Given such an average Q 12 , the net IR energy of the instanton matter of macroscopic density ρ(x) is Since our numeric simulations used instanton matter of approximately uniform density confined to a spherical cavity of large radius R, we are going to use similar setting in this section: Uniform ρ(x) = const while all space integrals are limited to a ball of radius R. Consequently, the net IR energy of the instanton matter becomes For future reference, let us explicitly calculate the IR energy of the baseline case of constant Q 12 = 1: This baseline energy serves as a convenient unit of IR energies for various orientation phases of the instanton matter. In particular, for the phases having uniform Q 12 = const, we immediately have Through the rest of this section, we shall calculate the Q 12 as a function of x 1 and x 2 and hence the net IR energy for all the orientation patterns we have observed in our simulations, namely the non-Abelian patterns, the spherical ferromagnetic and antiferromagnetic patterns, the global ferromagnetic pattern, and the ferromagnetic 'egg' (which obtains for some negative α and β). For completeness sake, we also include the global antiferromagnetic pattern we have not observed but suspected might show up in some obscure corner of parameter space.

Non-Abelian patterns
Let's start with the non-Abelian patterns involving 4, 6 or perhaps more different instanton orientations in any lattice cell. In terms of the S 3 group manifold of the SU (2), the net 4D quadrupole momentum of such orientations should have no traceless part. In the 4-orientation case y = 1, iτ 1 , iτ 2 , iτ 3 (so that the corresponding quaternions are 1, i, j, k) vanishing of the traceless quadrupole moment is obvious; we have checked that it also vanishes for the 6-orientation pattern of a hexagonal lattice and we believe this is a general case for any non-Abelian orientation pattern. In terms of the instanton orientations y i themselves, this means that averaging over a complete periodic cell of the lattice -and hence over any macroscopic volume -one gets average over y i of tr 2 (y i Y † ) is 1 for any given SU (2) matrix Y.
In the context of (9), Y can be either y of a distant instanton or (in · τ ) y , hence Since this averaged coefficient of the 1/r 2 force is constant, it follows that the net IR energy is

Global ferromagnetic pattern
Now consider the global ferromagnetic pattern in which all the instantons have the same orientation y = const. Let's assume no macroscopic domains, so that y stays the same over the whole big ball of the instanton matter. In this case y 1 y † 2 = 1 for any two instantons, hence Q 12 = 1 2 + 4α.
Again, we have a constant Q (although a different constant from the NA patterns), hence the net IR energy is

Global anti-ferromagnetic pattern
Next, the global anti-ferromagnetic pattern in which instantons have 2 alternating orientations y a and y b related by (the SU (2) representation of) a 180 • rotation around some axis N, that is y b = ±i(N · τ ) y a . By the global AF pattern we mean there are no domains so the two orientations y a and y b are the same everywhere in the big ball of the instanton matter, thus Averaging between these two cases, we get where n 12 is the direction of the line between the two instantons. For the purpose of calculating the net IR energy, we may further average the Q over the directions of the x 2 and of the x 1 and hence over the directions of the n 12 , thus and therefore

Spherical anti-ferromagnetic pattern
Now consider the anti-ferromagnetic spherical shells that we found in the simulations for β > α. In the anti-ferromagnetic case, the microscopic orientation pattern is anti-ferromagnetic, but macroscopically the alternate orientations y a and y b vary from place to place such that near x the flip axis N points in the direction of x, thus Note: taking 2 lattice steps in some direction, an AF orientation y a should turn into the y b and then back into y a , but since the n x unit vectors for the two steps are slightly different, we find that the y a orientations at the double-step separation are not exactly the same.

Spherical ferromagnetic pattern
Finally, consider the ferromagnetic spherical shell pattern in which the orientation is locally ferromagnetic -all the instantons in a neighborhood have similar orientation y, -but globally y(x) slowly varies from place to place. In particular, in the pattern seen in the simulations, up to a global SU (2) rotation. For this pattern, we need no averaging but Q 12 = Q 12 depends on the angle θ 12 between the directions of the instantons' locations x 1 and x 2 . Specifically, tr(y(x 1 )y † (x 2 )) = ±2 cos θ 12 while tr y(x 1 )y † (x 2 )(in 12 τ ) = ±2(n 1 × n 2 ) · n 12 = 0, hence The net IR energy is rather simple:

Summary of homogeneous orientation patterns
Thus far, we have calculated the IR energy of all the homogeneous orientation patters we have seen in our simulations: any non-abelian pattern: global ferromagnetic pattern: global anti-ferromagnetic pattern: spherical anti-ferromagnetic shells: for ν = 4 log(2) − 1 3 ≈ 0.59, spherical ferromagnetic shells: Actually, we have never seen the global anti-ferromagnetic pattern in our simulations, but we have included it just in case. Now let's compare these IR energies for various ranges of α and β and find which pattern has the lowest IR energy. Note however that the UV stability of the instanton matter requires the two-body 1/r 2 to be repulsive for all orientations and hence α ≥ − 1 8 and β ≥ 1 8 . Also, in this section we shall focus on the parts of the parameter space where α ≥ 0, or β ≥ 0, or both. The remaining part of the parameter space where both α and β are negative will be explored in the next subsection 4.7 once we calculate the energy of the in-homogeneous ferromagnetic egg pattern.
Here is the summary of our results: • The non-Abelian orientation patterns win -i. e. have lower IR energy that any other pattern -for α > β > 0 and also for α > 0 while β < 0 (but β ≥ − 1 8 ).
-Note that there are several different non-Abelian patterns, but they all have degenerate IR energies. In the next section 5 we shall compare their energies beyond the IR limit and find the specific winning non-Abelian patterns for each range of α and β.
× Finally, the global anti-ferromagnetic pattern never wins: for any allowed α and β some other pattern has lower IR energy than the global AF patter.
To conclude this summary, let me emphasize the phase boundaries between the different IR patterns: The boundary between the non-Abelian and the spherical antiferromagnetic patterns lies along the positive β = α line. In our simulations, this boundary seems to lie at small but positive β − α, but this is probably due to finite-sample-size or finitetemperature effects.
The boundary between the spherical anti-ferromagnetic and ferromagnetic shells lies at positive β = Γα line. In simulations, this boundary seems to lie at slightly higher β/α ratios between 2.2 and 2.3, but again this is probably due to finite-sample-size or finite-temperature effects.
Finally, the boundary between the global and the spherical-shell ferromagnetic patterns lie along the α = 0 line, and that's exactly what we saw in our simulations.

Ferromagnetic egg pattern
Thus far, we have calculated the IR energies of all the homogeneous orientation patters we have seen in out simulations. However, for some values of negative α and negative β the simulations produced the in-homogeneous egg-like pattern where the central part of the spherical cavity -the yolk of the egg -is a ferromagnetic lattice of constant orientation y, while the outer part of the cavity -the white of the egg -is made of ferromagnetic shells where the instanton orientation y slowly changes with the spherical coordinates. Moreover, all the y(θ, φ) in the 'egg's white' are perpendicular to the uniform y of the 'yolk'. Specifically, the orientations are for 0 < r < ξR, y = 1, for ξR < r < R, y = in · τ (43) where R is the cavity's radius and ξR is the boundary between the 'yolk' and the 'egg white' orientation patterns. Typically 0.73 < ξ < 0.87, so both the 'yolk' and the 'white' patterns take a significant fraction of the overall cavity volume. For such ferromagnetic egg pattern, the Q 12 for for two instantons at x 1 and x 2 does not need macroscopic averaging, but its value depends not only on the directions n 1 and n 2 of the two instanton positions but also on their radii r 1 and r 2 : • For r 1 , r 2 < ξR -i. e. when both instantons are in the egg's yolk, (44) similar to the global ferromagnetic pattern.
• For r 1 , r 2 > ξR -i. e. when both instantons are in the egg's white, similar to the ferromagnetic spherical shell pattern.
• But when one of the instanton is in the yolk while the other is in the egg's whitesay, for r 1 < ξR < r 2 , - For simplicity, we assume that both the yolk and the white parts of the ferromagnetic egg have equal and uniform instanton densities, ρ(x) = const throughout the cavity. Nevertheless, for the above formulae for the Q 12 , the energy integral becomes rather complicated and takes Mathematica to evaluate. Even after several rounds of simplification, we get a mess of dilogarithmic functions, where f (ξ) = π 2 3 + 48 − π 2 12 ξ 4 + 2(1 − ξ 2 ) 2 ar tanh(ξ) and g(ξ) = − π 2 + 16 2 Note: although these formulae for the f (ξ) and g(ξ) are rather formidable, numerically plotting them as functions of ξ shows that for all ξ between 0 and 1, and also f (ξ) + g(ξ) > 2 and f (ξ) + Γg(ξ) > π 2 4 .
In light of these bounds, the ferromagnetic egg pattern never wins the lowest energy competitions against the homogeneous patterns for α ≥ 0 or β ≥ 0 or both α, β ≥ 0. In these regimes, the winner is always a homogeneous patters discussed in the previous section $4.6, namely, non-Abelian, AF spherical shell, FM spherical shell, or global FM. But the regime of both α, β < 0 requires a more careful consideration. To find the winning pattern in the α, β < 0 regime, we minimize the ferromagnetic egg's energy (48) as a function of the relative yolk size ξ and then compare the minimum to the other two competing patterns in this regime, namely the non-Abelian pattern and the global ferromagnetic pattern. For simplicity, we first take a difference and then plot h as a function of ξ for different values α/β. If the we see the whole h(ξ) curve lying above the h = 0 axis, then the lowest-energy pattern is non-Abelian rather than the the ferromagnetic egg. On the other hand, if the curve dips below h = 0, then the lowest energy is some kind of a ferromagnetic pattern: Global FM if the minimum obtains for ξ = 1, spherical FM shells if the minimum obtains for ξ = 0, and FM egg if the minimum lies for ξ strictly between 0 and 1.  Figure 6 shows the actual plots h(ξ) for several values of the α/β ratio. We see that for (α/β) < 0.0575 the whole curve lies above h = 0, so the lowest energy pattern is non-Abelian. For (α/β) = 0.0575 the curve's minimum touches zero -see the left half of figure 7 for the zoomed up plot -which indicates a first-order phase transition, -and for the higher (α/β) ratios the minimum goes negative, so the lowest energy pattern is ferromagnetic. Moreover, as long as (α/β) 0.5 there is a unique minimum at ξ < 1 so the winning pattern is the FM egg rather than the global FM. As we increase the α/β ratio, the minimum moves to larger ξ, and for (α/β) ≈ 0.5 the curve develops two local minima, one at ξ ≈ 0.85 (FM egg) and the other at ξ = 1 (global FM). At (α/β) = 0.543 -see the right half of figure 7 for the zoom up, -the two minima become degenerate, which The bottom line is: in the regime of negative α and negative β, specifically − 1 8 < α, β < 0, there are three IR-different orientation patterns, depending on the α/β ratio: • The globally ferromagnetic pattern for − 1 8 < α < 0.543β. • The ferromagnetic egg pattern with relative yolk radius 0.73 < ξ < 0.87 for 0.543β < α < 0.0575β.

Summary
To summarize this whole section, the diagram of all IR-different orientation patterns for all allowed values of α, β > − 1 8 is shown on fig. 8.

The phase diagram of the non-Abelian crystals
We then continue to explore the phase diagram of the instanton crystals with the generalized interaction of eq. (9) depending on the parameters α and β. We have a more detailed look of the phase diagram in the non-Abelian phase, i.e., the magenta area of fig. 8. This phase is interesting because (unlike in the other phases) simulations suggest that various different crystal structures appear in this phase, giving rise to a highly nontrivial phase diagram. Moreover, these crystals are "regular" so that they leave some of the (discrete) translation symmetries intact, unlike the structures found in the spherical (anti-)ferromagnetic and mixed ferromagnetic phases. Therefore they admit a precise mathematical description, which makes further analysis possible. We proceed as follows. Based on the simulation results, we identify the precise structure (including spatial lattice and orientations) of all crystals appearing in the non-Abelian phase. This allows us to generate numerically very large perfect crystal domains, having O 10 7 instantons. Then we can compute the energies of all these crystals to a high precision, and compare them in order to draw the phase diagram.

Non-Abelian crystal structures
By carrying out simulations in Sec. 3, we have identified several types of crystals as the ground states of the system. Up to changes in aspect ratios, they can be categorized into four different classes (excluding the hexagonal layers of fig. 4 (bottom left) which is found at large α and β; we restrict here to α and β of at most O (1)). The unit cells for these classes are shown in fig. 9, and the precise analytic definitions are the following: (a) Tetragonal/cubic (fcc) lattices with "standard" orientation pattern. By this we mean fcc with the natural pattern made of four orientations (1, i, j, k), and crystals which are obtained from this by rescaling the coordinate system in one spatial direction, so that the configuration depends on one aspect ratio. The crystals in fig. 2, fig. 3 (left), and fig. 4 (top right) fall into this class.
The explicit definition is as follows. We take a unit cell defined by the following three vectors where e j are the unit vectors for the spatial coordinate system and c is the aspect ratio.  We then place the following four instantons The whole lattice is then obtained by carrying out shifts with linear combinations of a j with integer coefficients. In fig. 9 (a) we show the defining instantons and all other instantons which touch the unit cell. In our conventions, fcc is obtained for c = 1 (which is the case shown in the figure) and body-centered-cubic (bcc) for c = 1/ √ 2. In simulations, we have observed phases both with a wide range of values of c, including c < 1 (often close to the bcc value), c = 1 exactly, and c > 1. These crystals are denoted with blue color in the plots below.
(b) Tetragonal/cubic lattices with "alternative" orientation pattern. This system is otherwise the same as (a), i.e. the spatial structure is obtained from fcc by scaling and the orientations are (1, i, j, k), but there is an additional shift in the orientations in alternating layers of the unit cells in z direction. The crystal in fig. 4 (bottom right) falls into this class.
Precise definition of the lattice structure is as follows. The units cell now corresponds to two unit cells of the fcc lattice, The cell is then filled with the following eight instantons (see fig. 9 (b)) 3 The whole lattice is again obtained by shifting with linear combinations of a j with integer coefficients. We have found this lattice in simulations for β ≈ 0 and with aspect ratio c < 1, being typically close to the bcc value of 1/ √ 2. These crystals are denoted with red color in the plots below.
(c) Hexagonal lattices. We have obtained as end points of the simulation lattices that are related to hexagonal-close-packed (hcp) lattice by a coordinate rescaling perpendicular to the triangular planes. The orientation pattern includes six different orientations which form two equilateral triangles. The crystal in fig. 4 (top left) falls in to this class.
The precise definition can be taken to be the following. The unit cell is defined in terms of Notice that for hexagonal systems the unit cell is not rectangular. Moreover our definitions do not follow the standard in the literature; this is because due to the orientation structure it is convenient to define a unit cell which is larger in the e 1 and e 2 directions than what is typically used for hexagonal lattices. The hcp lattice has c = 2 √ 2/3 ≈ 0.9428 in our conventions.
The orientations and locations of the defining six instantons are given by (see fig. 9 (c)) (64) In simulations we have seen a variant of the lattice with the aspect ratio around c ≈ 0.5. This kind of hexagonal lattice is found in nature in the tungsten carbide compound (which has c ≈ 0.57). The comparison of energies done in this section reveals another variant, having c close to the hcp value, is dominant in a narrow region on the (α, β)-plane. Hexagonal crystals are denoted with green color in the plots below.
(d) Simple cubic/tetragonal lattices. Finally we have also seen non-Abelian simple cubic lattices in the region β < 0. The crystal in fig. 3 (right) falls in to this class. The precise definition for these lattices can be taken to be the following. The unit cell is defined in terms of a 1 = 2 e 1 , a 2 = 2 e 2 , a 3 = 2c e 3 .
Here we chose the units such that the shortest distances between the instantons in the directions e 1 and e 2 will be equal to one, but a factor of two was added due to the orientation structure, which will repeat in cycles of two for each spatial directions. The orientations and locations of the defining eight instantons are given by (see fig. 9 (d)) In simulations, we have only found lattices with c = 1 exactly. However, based on the analysis of the minimum energies of the various lattices, we expect that cubic lattices with c = 1 (both lattices with c > 1 and with c < 1) have the lowest energy near the critical value β = −1/8. All the simple cubic or tetragonal lattices are marked with yellow in the phase diagrams.

The phase diagrams
In order to carry out the detailed analysis of the non-Abelian phase, we computed numerically the energies of large lattices (with O 10 7 instantons) for the various crystal structures. See Appendix B for details of the computation. We included all the crystals of fig. 9 and scanned over all reasonable values of aspect ratios for each of them. First, before going to the non-Abelian structure, we verify that the lowest energy configuration for the unoriented interaction of (7) is indeed fcc. The fcc, hcp, bcc, and simple cubic lattices are all local minima, with the energies satisfying E(fcc) < E(hcp) < E(bcc) E(simple cubic). We also checked that the dominant crystal for the oriented, standard instanton interaction of (5) is the tetragonal (fcc-related) crystal having a large aspect ratio, which is shown in fig. 2. For the aspect ratio we find c = 2.467, a number slightly larger than that estimated from the simulations, c ≈ 2.35.
We then draw the phase diagram of the non-Abelian lattices on the (α, β)-plane. The result for the phase diagram is shown in fig. 10. In the figure the blue, red, green, and yellow colors indicate crystals (a), (b), (c), and (d) of fig. 9, respectively. The brightness of the color is given by the aspect ratio, with brighter (darker) shades corresponding to lower (higher) values. The largest solid blue domain is fcc, i.e., aspect ratio is c = 1, and the large yellow domain is simple cubic. In all the other domains, the aspect ratio varies with α and β.
We continue by analyzing the phase diagram close to the critical lines already observed in the simulations in Sec. 3, starting from the results near the diagonal, 0 < α = β < 1.
We show the phase structure in this region in fig. 11, where we switched to angular coordinates and show the "radial" coordinate in log scale. In this area only the crystals (a) (fcc) and (c) (hexagonal) of fig. 9 appear, but the aspect ratios vary. In the light-blue domain the aspect ratio is close to the bcc value except for very close to the transition line to fcc and the tail at largest values of α within this domain. In the hexagonal green domain, the aspect ratio is close to the value of tungsten carbide (denoted by WC in the plot) at large α and β. As α decreases there is however a crossover to region where the aspect Figure 11: The details of the phase diagram near α = β. We use angular coordinates and show the horizontal axis in logarithmic scale. Notation as in fig. 10.
ratio is close to the hcp value of 2 √ 2/3 ≈ 0.94. The crossover takes place at α ≈ 0.15. In the remaining dark blue tetragonal domain, the aspect ratio is large and increases with β, reaching the value of c ≈ 2.467 at the physical point α = 1 = β. Therefore in this region the structure is most naturally viewed as aligned planes of anti-ferromagnetic square lattices with alternating orientations, as we also saw in simulations in Sec. 3. There are additional details in the area where the (near) bcc, fcc and hexagonal phases meet, which are poorly visible in fig. 11. To resolve the details, we show a zoom into this region in fig. 12. The blue regions in this plot are all versions of the same crystal, (a) in fig. 9, but having slightly different aspect ratios close to the fcc value c = 1. Since the aspect ratios are close, they are shown with similar colors. To make the various phases better visible we added stripes in the plot. All phase boundaries are first order transitions, but the transitions between the near fcc phases are weak. In particular the triple point is critical: transitions along generic lines intersecting the point are of second order.
Results for the energy differences and aspects ratios along the line α = β are shown in fig. 13. The plots on the top row show the energy difference ∆E between various phases  Figure 13: The non-Abelian crystals for 0 < α = β < 1. Top row: the energy difference of different phases with respect to the fcc configuration as a function of α. Bottom row: the aspect ratios as a function of α. Blue colors indicate cubic/tetragonal lattices, fig. 9 (a), and green colors are hexagonal, fig. 9 (c). In more detail, the solid blue, dashed light blue, and dotted dark blue curves are the results for the fcc, (near) bcc, and tetragonal large-c lattices, respectively, whereas the dot-dashed green curves show the result for the hexagonal lattices. Dominant phases are shown as thick curves and subdominant as thin curves. and the fcc configuration (with aspect ratio c = 1) as a function of α = β, in linear scale in the left plot and in log-log scale in the right plot. The "units" of ∆E depend on the system size (see Appendix B). The bottom row plot show the aspect ratio as a function of α in the various phases. Along the line α = β there are two first order phase transitions, as we can also see from fig. 11. In the unoriented limit of very small α, the fcc phase (with aspect ratio c = 1 exactly) is the ground state. At α ≈ 0.0019 we find a first order phase transition to a hexagonal phase (dot-dashed green curves). This phase transition is so close to the origin that it is not visible in the top left plot but it is identified as the point where the dot-dashed green curve hits zero on the top right log-log plot. At α ≈ 0.41 there is another transition, to the tetragonal lattice with large aspect ratio (dark blue dotted curves). This is the ground state up to α ∼ 2, including the physically interesting point α = 1. It is surprise that we find hexagonal ground state at α ∼ 0.1, since the simulations of Sec. 3 would converge to a bcc structure in this region. Notice however that the energy difference between the hexagonal structure and the bcc (light blue dashed curves) is extremely slim. Apparently our simulations are not sensitive enough to discern this difference.
Notice also that the crossover in the hexagonal phase from approximate hcp lattice to tungsten carbide lattice at α ≈ 0.15 is clearly seen in the bottom plot: the aspect ratio in the hexagonal phase varies rapidly, starting from c ∼ 0.9 at low α and ending at c ∼ 0.5 around α = 0.4. At the same time the energy in the top left plot shows a smooth kink around α ≈ 0.15 for the hexagonal phase. The value of the aspect ratio for the tungsten carbide compound appearing in nature, c ≈ 0.57, is crossed at α ≈ 0.21.
Recall that in the fcc configuration (solid blue curve) c = 1 exactly due to the enhanced symmetry, but e.g. for the bcc configuration (dashed blue curves in fig. 13 and the light blue domain in fig. 11) the aspect ratio equals only approximately the bcc value 1/ √ 2, and depends mildly on α and β. We have verified numerically, however, that c → 1/ √ 2 for the bcc configuration as α → 0 and β → 0, which is natural as the effect of the orientations is suppressed in this limit. Similarly in the hexagonal phase c → 2 √ 2/3 as α → 0 and β → 0, which is the value for hcp lattice in our conventions. Notice that neither the bcc configuration nor the hcp configuration is the ground state at small values of α and β but they still exists as subdominant local minima, as seen from fig. 13. Figure 14: The phase diagram for the non-Abelian crystals near β = 0 in angular coordinates. Notice that the horizontal axis is in logarithmic scale. Notation as in fig. 10 We go on discussing the phase diagram in the vicinity of the line β = 0. As already seen from fig. 10, in this region the tetragonal lattice with "alternative" orientation pattern marked with red color (lattice (b) in fig. 9) dominates. We have zoomed into this region in fig. 14 by using angular coordinates, in analogy to fig. 11. We see that the tetragonal phase (red domain) is divided, at small α, into two phases separated with a first order transition line ending at a critical point. The lower (upper) phase has c ≈ 1/ √ 2 (c ≈ 1) so the lattice is close to bcc (fcc), respectively. We have verified numerically that (apart from the section in the near fcc phase at very small α) c = 1/ √ 2 on the β = 0 line (the horizontal axis in the plot) so the lattice is exactly bcc on this line.
Finally, we have a closer look at the region of negative β in fig. 15. First, we notice that the transition line between the tetragonal phase (red domain) with alternative orientation pattern, and the simple cubic lattice (yellow domain) is almost precisely horizontal; we find that the transition takes place at β ≈ −0.011. Second, there are two additional phases near the critical line β = −1/8. These phases have simple tetragonal structure, i.e., they have the orientation pattern of fig. 9 (d) but nontrivial aspect ratios. The phases were not found in simulations in Sec. 3, but the analysis of energies shows that they are favored with respect to the simple cubic in the vicinity of the critical line. The left phase (light yellow, smaller α) has c < 1 and the right phase (dark yellow, larger α) has c > 1. All the three yellow phases are separated by first order transitions, which become weak in the vicinity of We also notice that, as expected, the orientation dependence of the final phase diagram reflects the nature of the interaction terms (which was discussed in Sec. 2.2). First of all, as we already see from Fig. 8, negative α leads to ferromagnetic order, whereas for large positive α we encounter non-Abelian order. This is expected as positive (negative) α in the potential (9) favors different (the same) orientations. The somewhat more complicated effect of the β term can also be discerned. For example the typical non-Abelian phases for β > 0 in Fig. 10 are the fcc or tetragonal phases (a). From Fig. 9 we see that for nearest neighbor pairs, the twist y † m y n is perpendicular to the link N mn between the instantons (assuming the standard mapping from orientations to spatial directions). For the phases found at negative β, i.e., simple cubic or simple tetragonal phases (d), the twist of any pair of nearest neighbors is parallel to the link between them.

Conclusions and open questions
The analysis of non-Abelian lattices and in particular those associated with large N c nuclear matter is in an infant stage and there are many questions to further investigate. Here we enlist some of them.
• The basic two body instanton interaction is with α = β = 1. In this paper we have discovered a very rich structure of the phase diagram that followed from generalizing the interaction by varying the values of α and β. On top of the simulations described in this paper we have also performed several simulations in the region of α > 1 or/and β > 1 and found interesting structure. It will be interesting to further study this region.
• The results of this paper are based on the generalized two instanton potential given in (9). A natural question is to what extent is this potential indeed general. In particular what is the form of the potential in the Skyrme and other nuclear matter models. A simple modification of the potential is achieved by using 1 r d instead of 1 r 2 . • Another obvious generalization is transforming the flavor symmetry U (2) → U (N f ) and in particular the eight-fold symmetry N f = 3. In the opposite direction one would like to explore also the breaking of the isospin symmetry.
• The systems discussed in this paper are at zero temperature. It will be interesting to consider instead lattices of instantons (Torons) at finite temperature.
• We assumed here that the instantons were essentially point-like. Taking into account effects due to their finite size is a demanding task. Some nontrivial configurations are known however for the Skyrme models (see, e.g., [46][47][48]).
• The simulations performed in the paper have been made for instantons that were placed in a spherical cavity and a compensating "external" force was exerted on them to ensures that the sample has a constant density of instantons at large scales. We have made certain checks about the dependence of the final results on these two features and found very mild dependence mainly in the outer shells of instantons. This dependence on the structure of the cavity and the form of the external forces deserves further more elaborated study.
• In recent years there has been an effort to apply holography to the study of neutron starts and their mergers. It will be interesting to understand if and how taking crystal structures rather than the liquids change the picture of neutron stars. Moreover, our analysis included the spherical fermionic and anti-fermionic structure. These phases may be relevant to the nuclear structure of neutron stars that obviously are spherically symmetric.
• Dressing the instantons with electric charges and differentiating between neutral and charged nucleons is a very important factor in real nuclei and nuclear matter. It is interesting to redo the analysis of the crystalline structure for chunks of holographic neutrons and protons.
• This research work dealt with uniform nuclear matter. It will be interesting to explore in what physical systems such uniformity can be realized. Needless to say that the study of non-uniform ensembles of instantons should also be of great interest.

A Potentials for a ball of homogeneous matter
In this Appendix we compute the potential of a spherical homogeneous configurations (i.e., ball of radius R) of constant density, corresponding to the continuum limits of the various phases encountered in the simulations. The setup is therefore similar to that analyzed in Sec. 4: The lattice geometry does not matter. The geometry-independent IR potential of a probe instanton scales like V IR ∼ ρR, while the geometry-dependent contributions scale like ∆V ∼ ρa, where a is the lattice spacing. Therefore the geometry-dependent terms vanish in the continuum limit a → 0. Moreover, we are taking an average over the orientations locally. More precisely, as explained in Sec. 4, it is enough to consider the average of the orientation dependent expression Q 12 of eq. (15). These potentials are important because they are also used as an input for our numerical simulations: an external force is required for the (locally averaged) simulation results to have constant density. This amounts to the interactions between the instantons in the simulations and those outside the simulated ball (either exactly or up to some additional factors, depending on the phase). We discuss the interpretation of these potentials in more detail below.

A.1 Potentials for the various orientation patterns
In general, the IR potential for a probe instanton located at point x 1 (or rather average potential for the instantons in a small neighborhood of x 1 ) inside a large ball of instanton matter of uniform density ρ is Further integrating over the location x 1 of the probe within the ball of instanton matter would give us the net IR energy of all the instantons that we have analyzed in Sec. 4, specifically But in this Appendix, we are interested in the potential (72) and the consequent bulk force which we would need to cancel by fiat in order to get a uniform instanton density in our numeric simulations. Before we focus on specific orientation patterns, let's consider the baseline case of Q 12 = 1 and hence This baseline potential is spherically symmetric, and it is easy to check that it is positive for all r 1 < R and monotonically decreases with the radius r 1 = |x 1 |, see fig. 16 (left). Consequently, the bulk force (74) always points away from the center, The force is shown in 16 (right). Notice that it diverges logarithmically as r → R. Given this baseline case, let's focus on the orientation patterns we have studied in section 4: • Non-Abelian patterns. In this case we found that Q 12 = 1 2 + α + β in Sec. 4. Since this averaged coefficient of the 1/r 2 force is constant, it follows that the IR potential is • Global ferromagnetic pattern. Again, we have a constant Q 12 = 1 2 + 4α (although a different constant from the NA patterns), hence the IR potential is • Global anti-ferromagnetic pattern. For the global anti-ferromagnetic pattern, This time, the averaged Q is not constant but direction dependent, so the IR potential for the probe instanton at x 1 is not going to be spherically symmetric. We therefore start by averaging the expression (77) over the direction of the x 2 but not over the x 1 . This gives us where θ 12 is the angle between the x 1 and x 2 directions while θ 1 is the angle between the x 1 and the direction N of the AF flip. The θ 1 dependence of eq. (78) is a combination of a constant and a quadrupole, so the IR potential also has a spherically symmetric and a quadrupole terms, The radial dependence of this quadrupole term is plotted on figure 17. As we can see, V AF quadrupole is nearly constant within a ±5% range, so it does not contribute much to the radial force. Instead, it cases a 1/r force in the longitudinal direction n θ , • Spherical anti-ferromagnetic pattern. In this case we found that Q 12 = 1 2 +α+β + (α−β)×cos θ 12 − 2β× (r 1 − r 2 cos θ 12 )(r 2 − r 1 cos θ 12 ) Note that this Q 12 depends on the angle θ 12 between the x 1 and x 2 directions but is invariant under simultaneous rotation of both x 1 and x 2 . Consequently, after integrating over the x 2 we get a spherically symmetric IR potential for the instantons at x 1 . Specifically, where which is drawn in fig. 18 is positive for all r 1 , while the ∆ 2 V happens to vanish, Note that for α = β the potential for the spherical AF setup is equal to that of the NA setup.
• Spherical ferromagnetic pattern. For this pattern, Q 12 = 1 2 + 4α cos 2 θ 12 . Note that Q depends only on the relative angle between the directions of x 1 and x 2 , so integrating over the x 2 gives us a spherically symmetric IR potential for an instanton at x 1 . Unfortunately, the explicit formula for this potential is rather messy, The potential is drawn in fig. 19 and the corresponding force is given by which is drawn in fig. 20. A.2 Long-distance behavior of the potentials Finally, let us comment on the long-distance cutoff dependence of these potentials. This is also linked to the fact that we are using the potentials to determine the external forces applied to the instantons in the simulations that we carry out. Actually, the natural choice for the external potential would be the potential of a homogeneous distribution (having a constant density which equals the average density of instantons in the simulation) outside the volume of simulation. However, due to the long-distance divergence of the two-body interaction, such a potential is not a priori well-defined, but depends on the long-distance cutoff.
For some of the potentials discussed above there is only one reasonable definition of the external potential, at least up to an irrelevant constant. These are the patters for which Q 12 is constant, i.e., the non-Abelian and global ferromagnetic crystals. In this case the long-distance divergence for the homogeneous distribution is that of the plain 1/r 2 interaction. Then the result is independent of the precise definition of the cut-off, with the (reasonable) assumption that it does not introduce a "dipole" term (for example the cutoff is invariant as x → −x). With such a cutoff and at constant density, the force due to homogeneous matter inside the ball of radius R is exactly the opposite of the force due to matter outside the ball. This fact is reflected in the large R expansion of the base potential, i.e., interpreting the value of R as the long-distance cutoff. The only cutoff dependent piece is an irrelevant constant. For the crystals having nontrivial long-range correlations between orientations, such as the spherical antiferromagnetic pattern, the above argument does not hold. Indeed, by expanding the potential at large R we find that V IR SAF (r) = 4πρ so that even with a simple spherical cutoff there is a nontrivial potential term depending on R. We also note that the force due to this potential is singular at the origin: it diverges as ∼ log r. Obviously such a divergent force field cannot be generated by any distribution of instantons outside the sphere, let alone a homogeneous distribution. Similar issues appear in the spherical ferromagnetic and ferromagnetic egg phases. Therefore the potentials which we use in the simulation, which sets the density of the instantons to (roughly) constant, are lacking a clear physical interpretation in these phases.

B Numerical details
The numerical analysis of the configuration consists of two steps: First, we carry out simulations for various choices of the two-body potential, starting from a random initial condition and using an algorithm to seek for a minimum of total energy. We use a setup which roughly corresponds to adding masses and drag forces for the instantons, and simulating their dynamics as the system converges to a final state with low energy. We will explain this in more detail below. The likely structure of the ground state can then be extracted from the final simulation results. The various crystals and the phase structure arising from this analysis are presented in Sec. 3. Second, for final states that are regular lattices, we can generate large configurations and compute their energies numerically directly without doing simulations. The energies can then be used to draw the phase diagram of the system as a function of the parameters of the interaction potential. The results from this analysis are used in Sec. 5 where we draw the final phase diagram of the non-Abelian crystals.

B.1 Setup for simulations
We took N instantons inside a three dimensional ball with radius R, starting from random initial positions X k and random initial orientations y k , and used an algorithm to minimize the total energy of the system. In order to reduce boundary effects, we also added an external force which sets the density of the instantons (averaged in a small neighborhood) to constant (see Appendix A). For the algorithm searching for the minimum of the energy, we implemented both a simple direct steepest descent method (first order formalism) and a second order formalism, which amounts to considering a dynamics for the instantons with a mass and a drag force and keeping track on their momentum. The second order formalism was seen to be clearly more efficient than the first order formalism.
We briefly sketch of the second order formalism which we use. For simplicity, we discuss mostly the relaxation of the positions (rather than orientations) of the instantons. Relaxing the system is an iterative procedure, which aims at minimizing the total energy where the two-body potential is given in (9) and the external potential is chosen such that the (locally averaged) density of the instantons is constant. The external potentials are given in Appendix A for all crystals we encountered. For the non-Abelian phase, where we did most of our simulations, it equals the interaction energy between the instanton n and the unoriented continuum outside the ball. It is given (up to an irrelevant constant and normalization) by the base potential in (75), where the normalization was fixed by the requirement that density of the continuum matches with the average density of instantons inside the ball. Let us then go to the details of the iterative algorithm which we used to minimize the energy. For each step in the procedure we compute a pseudo force as where the estimate in the denominator essentially averages over orientation dependence, i.e., we plug in ∇ 2 (k) E tot est = n =k 6 |X n − X k | 2 E 2 body (n, k) + Here ∇ (k) acts on X k , and the first term is estimated by using the Laplacian for the unoriented potential (α = 0 = β). We divide by this estimate because then the force (93) roughly measures the distance of the instanton k from its equilibrium position when we are close to a local minimum. With this normalization it is easier to optimize the algorithm for the search of the minimum than when taking the force to be exactly the gradient of the energy. A simple (first order) gradient descend would now simply iterate the positions simultaneously for all the instantons by using where "old" and "new" refer to the positions of the instanton before and after the iterative update, respectively. A good guess for the coefficient κ = O(1) after normalizing the pseudo force as in (93). But as promised above, we implement a second order method. For this, we define a momentum P k for each instanton which is initially set to zero. Then the iterative update is defined as where κ P and κ F are parameters which were tuned by trial and error. They correspond to adding masses and dissipation to the system. Notice that for κ P = 0 one gets back the first order formalism, and one needs to have 0 ≤ κ P < 1; the limit of κ P → 1 from below is the limit of no dissipation. It turned out that good values for κ P are relatively close to one (meaning small dissipation); typically κ P = 0.99 would lead to good convergence, i.e., a crystal with low total energy with a relatively low number of iteration steps. For κ F , values close to the largest possible (where the system would still converge) were found to be optimal. In order to determine the optimal parameters (which slightly depend on both the system size and ground state crystal) we found it useful to monitor the total energy (91) while running the simulation.
The above algorithm will typically produce final state which is suggestive of crystal structure but has too many defects to definitely determine the crystal with lowest energy. In order to improve the quality, we add one more ingredient: a pseudo temperature. That is, we replace (93) by where T is the pseudo temperature and each component of the random force F rand was taken from a uniform distribution between −1 and 1. A strategy which would produce crystal domains of good quality is the following. First let the system relax with zero temperature T = 0, repeating the update of (96) and (97), until the convergence stops so that the positions or orientations of the instantons no longer change. Then set T to a value, found by trial and error, which just below the critical value which would completely melt the lattice structure. Then repeat the iteration procedure of (96) and (97) a large number of times, monitoring the total energy E tot , until it no longer decreases. As the last step, set the temperature again to zero and let the system to relax to a final state.
Finally let us comment on the relaxation in the orientation space. To implement this, we study the total energy in the vicinity of y k by defining E tot ( , k) = E tot | y k →y k e i · τ (99) and then define a pseudo torque by where for the estimate we used simply The rest of the implementation is fully analogous to the position space approach. In particular we also implemented a second order algorithm for the orientation space even though this turned out to be much less important for convergence than the implementation in the position space. The complete algorithm then updates the orientations and (isospin) angular momenta of the instantons simultaneously with the positions and momenta of Eqs. (96) and (97).
We also remark that the external potential term in (91) depends on the final ground state crystal, which is not known before running the simulation. Therefore we first assumed the "standard" formula (92), and if the simulation would converge to a spherical ferromagnetic and antiferromagnetic crystal, we repeated the simulations taking the correct form for the external force and checked that the system converges to the same crystal.

B.2 Computation of the energy differences
Here we discuss how we compute the precise phase diagram for the non-Abelian crystals. Once the precise structure of the crystals has been identified, as we have done in Sec. 5.1, it is in principle simple to compute numerically the energies of large samples of different crystals and compare them to draw the phase diagram. However, the size of the system is limited by computing resources, and because the two-body interaction leads to IR divergences for setups that are homogeneous at large scales, numerical noise due to the long distance cutoff in practice prevents a direct numerical comparison of the energies.
To overcome this issue we did the following. First, for all non-Abelian crystals which we have encountered, it is enough to compare the interaction energy of a single instanton with the rest of the sample. We choose this instanton to have unit orientation and to lie at the origin. We can then create large samples of various crystals including instantons up a radius R with fixed density of instantons (which also include the instanton at the origin), and compute the energy E(R) = n>1 E 2 body (n, 1) where the index 1 refers to the instanton at the origin. However even at the largest values of R which we can handle, the function E(R) contains essentially random noise coming from how exactly the cutoff r < R happens to select the instantons in the crystal structure. An efficient way to suppress this noise is the following. Instead of keeping R fixed, we choose its value from a probability distribution (e.g. Gaussian) around some large value and with the deviation δ much larger than the spacing between the instantons. We then take an average over R defined by this probability distribution. Choosing samples of O 10 7 instantons then turns out to be enough to remove essentially all noise from the results.
This approach boils down to computing the energy between the chosen instantons and the rest of the sample with a smooth large distance cutoff. That is, instead of (9) we use E 2 body (m, n; α, β; ∆, δ) = 2N c (1 + 2α + 2β)λM 1 |X n − X m | 2 × 1 2 + α tr 2 y † m y n + β tr 2 y † m y n (−i N mn · τ ) in (102) so that the interactions are cut off at distance ∆ with deviation δ. At N tot = 10 7 instantons, generated in a ball with radius R, good choices turn out to be ∆ = 0.7R and δ = 0.04R. One can show that the IR divergent term in the total energy is independent of the precise orientation structure and therefore exactly the same for all non-Abelian lattices at constant instanton density. The orientation dependent term is well convergent in the IR. We compute the energy differences between various non-Abelian lattices so that the divergent terms cancel, possibly up to terms depending on the details of the IR regulator. It is therefore enough to check that the results are independent of the parameters of the IR regulator to verify that the regulator-dependent terms are suppressed and the results are reliable. While it is presumably possible to do this analytically by using methods similar to those of Appendix A, we have only carried out a direct numerical check. Indeed we have verified numerically that the energy differences are insensitive to changes in ∆ and δ to a very good precision when the parameters are varied in a reasonable range. Moreover we checked that applying an ellipsoidal cutoff instead of spherical, or changing the number of instantons does not change the results.