Non-Abelian vortex lattices

We perform a numerical study of the phase diagram of the model proposed in \cite{Shifman:2012vv}, which is a simple model containing non-Abelian vortices. As per the case of Abrikosov vortices, we map out a region of parameter space in which the system prefers the formation of vortices in ordered lattice structures. These are generalizations of Abrikosov vortex lattices with extra orientational moduli in the vortex cores. At sufficiently large lattice spacing the low energy theory is described by a sum of $CP(1)$ theories, each located on a vortex site. As the lattice spacing becomes smaller, when the self-interaction of the orientational field becomes relevant, only an overall rotation in internal space survives.


Introduction
Non-Abelian vortices, first constructed in [2] and [3] [4], are Abrikosov-Nielsen-Olesen (ANO) vortices which support additional orientational (non-Abelian) moduli on their world-sheets. They have been widely studied in the literature (see, for example, [5] [6] [7] [8] [9] [10] and references therein) as candidates for the vortices responsible for the dual confinement mechanism [11] [12] [13], as originally proposed in [14] [15]. Although the original models in which they appeared contained a relevant degree of complexity (for example most models involved supersymmetry) recently a particularly simple extension of ANO's original model (based on Witten's superconducting string model [16]) was proposed which was shown to also contain them [1]. In fact, the general idea behind this model was successfully applied to many solitonic solutions with the same outcome: the condensation of a scalar field in the core of the solutions leading to orientational degrees of freedom [17] [18] [19] [20] [21] [22] [23]. More rencently, the model was also used in a holographic setup to find non-Abelian vortices in a dual 2+1 dimensional superconductor [24]. The appearance of the non-Abelian degrees of freedom has been attributed to distinct sources. In [25], mainly focusing on cosmic string applications, these were thought to originate from a dark matter sector. In [20] [22] however, a more condensed matter approach was taken and the superconducting system with such additional directional degrees of freedom was shown to have a close analogue in liquid crystals. In particular it was argued that this particular kind of superconductor can be thought of as a superconducting liquid crystal state, with the non-Abelian degrees of freedom playing the part of the director. It remains unclear whether this kind of superconductor can be realised in nature, and if it has any relation to the confining phase of QCD. Therefore, this research can be thought of as a toy model for many physical setups. In a cosmic setup it predicts the existence of periodic arrays of cosmic strings, stabilised by dark matter condensates. From a condensed matter point of view it describes a new phase of matter in which a superconducting liquid crystal form periodic vortex solutions. Finally, from a particle physics perspective it describes the lowest energy configurations of non-Abelian vortices, candidates to be responsible for confinement of quarks.
The action enjoys a local U (1) gauge symmetry in the ψ sector and a global O(3) symmetry in the χ i sector. We will take χ i to be a real triplet field. In recent related work, the additional χ sector was interpreted as dark matter [25]. The standard quartic potential in the ψ sector breaks the U (1) symmetry spontaneously. In general there is an interesting range of vacua in this model, the vacuum equations ∂ ψ V = 0 and ∂ χ V = 0 lead to several branches of solutions with the second and the third branch coalescing at the special point v = µ. We will be interested in the second branch of vacuum solutions which is only a global minimum if Vacua in which the χ field condenses are interesting and have been important in specific studies of cholesteric non-Abelian vortices [20] [22] but will not be treated here. The condensation of the scalar field "Higgses" the photon and gives it a mass defining the penetration depth of the superconductor as d ≈ (1/m A ). The scalar field is also massive with mass This mass defines the coherence length ξ ≈ (1/m ψ ) of the superconductor. The mass of the χ field is The coupling potential between the two scalar fields, with the assumption that v > µ, leads to χ i not condensing in the vacuum. However, whenever ψ vanishes (in general whenever its value is less than µ), the potential destabilises the χ field which therefore condenses.
In particular this happens in the core of an ANO vortex formed by ψ, which is a well known solution of the χ i = 0 theory. This mechanism leads to the existence of non-Abelian vortices, as shown in [1]. These vortices are lower in mass than the ANO vortices, so that the condensation of the additional field lowers their energy. The low energy theory describing the orientational gapless excitations is given by a CP (1) non-linear sigma model on the vortex world-sheet. This is most easily seen by the pattern of global symmetry breaking of the χ i sector, SU (2)/U (1) → CP (1), the remaining U (1) symmetry group corresponding to rotations in the plane defined by the direction of the χ i field in internal space (this will be apparent later when we discuss the ansatz). It is well known however that, for the case of superconductors of type II for which m ψ > m A , the above system with χ i = 0 has an energy minimizing periodic solution describing a lattice of ANO vortices, each carrying a unit of magnetic flux, the so-called "Abrikosov lattice". This point of transition between superconducting types, and in general whether the condensation of the additional core field, leading to orientational degrees of freedom, supports periodic structures, has not yet been investigated in this model and is essential in order to understand its general behavior over the whole region of parameter space. It is the purpose of this paper to numerically map out the precise nature of the superconducting transition in this model, classifying the superconductor type. The result of the numerical study is Figure  13, which shows that in the presence of the χ field the condition for type II superconductors changes and provides a whole view of the phase transition line plotted against the relevant parameters of the system. In the region marked type II, the system allows periodic vortex structures supporting non-Abelian moduli in their cores, which we also find and present numerically.
It is important for the rest of the paper to give a short review of Abrikosov's solution for the lattice near the superconducting critical point and of its extension to the full range of magnetic fields (see [26] for an application to color magnetism). We do this below.
The set of coupled equations of motion derived from (1) are (we apply no gauge-fixing condition yet) and the energy-momentum tensor functional reads We can generally parametrise the fields in order to gauge away the phase of the scalar field. Therefore we pick a parametrization of the form ψ = f (x, y)e iφ(x,y) , where f (x, y) and φ(x, y) are real, and A µ = Q µ + 1 e ∂ µ φ . From here on we fix the gauge so that ∂ µ A µ = 0 and consider only static solutions with A 0 = Q 0 = 0. We also define the dimensionless Note that in terms of these parameters global vacuum stability requires for type II superconductors we require a < 1.
We pick an ansatz for the χ i field which points in one direction only in the internal space, Rescaling the dimensionless fields to bef = f /v ,χ = 2β µ 2 χ andQ = Q/v, the dimensionless equations of motion become (we drop all tildes from the fields or derivative operators) For the time being let us ignore the χ sector (this is equivalent to setting b = 0). Then, following Abrikosov, we can find periodic solutions to the above equations in the vicinity of f ≈ 0, which is the point at which superconductivity is destroyed close to the critical magnetic field. In order to do so begin by expanding the fields as the series where is a small parameter denoting the deviation of the applied magnetic field B to the critical one B c where superconductivity is completely destroyed, = (B c − B)/B c . At first order the Q b will give us the applied magnetic field, hence we take this as Q b = (−By, 0, 0), and B is a constant. Abrikosov showed that, in this background and at first order, a periodic solution for the scalar field equation exists for which where θ 1 is the first elliptic theta function and x 1 , x 2 and y 2 are parameters which determine the lattice structure. The surface area of the lattice cell is simply S = x 1 y 2 . For a square lattice y 2 = x 1 and x 2 = 0, for a triangular lattice instead y 2 = √ 3 2 x 1 and x 2 = 1 2 x 1 . This result can alternatively be written as a Fourier series as where The common procedure to extend the lattice solution to the full range of magnetic fields is to use this solution as a seed to determine the full non-linear structure of the vortex lattice, for an arbitrary magnetic flux. It can be understood as a way to fix the total flux (or equivalently vortex number density) in the action minimization problem, as explained below. Note that for a fixed flux B there are infinitely many solutions corresponding to lattices given by parameters x 1 and y 2 (nature then presumably chooses these values so as to minimise the energy).
To find the lattice solution for any value of the magnetic flux (and not just those values close to the critical field) the strategy is the following: one must first fix the total fluxB passing through the integration domain and the lattice structure (given by x 1 and y 2 ) and then solves in this background for the local gauge field given by Q µ and the scalar field f . To do this consider the modified action where G µν is the field strength of an external U (1) gauge field G µ whose magnetic flux is where S here denotes the integration domain (which may comprise more than a single unit cell of the lattice). Therefore G denotes an external constant field, while F µν (A) denotes the internal electromagnetic response of the superoconductor. This leads to the equations of motion It is convenient to use the field variable Q → Q b i + G i , that is to split the gauge field into the external gauge field plus an internal one 1 , so that the equations become Given a fixed fluxB we must now look for periodic lattice solutions of the above equations with the requirement that Q b i has to satisfy a zero flux requirement. A particularly clever way to achieve this is to pick the external gauge field of the form [27] where ω A (B) is given in terms of the original Abrikosov lattice solution for the scalar field (28) but with B replaced byB (this is what we meant when saying that this solution was used as a seed in the full non-linear problem). In this definition κ = d/ξ is the ratio of the magnetic penetration depth to the coherence length. For type II superconductors κ ≥ 1/ √ 2. In terms of the units used in this paper we have that so that the critical a is a = 1 corresponding to the point where the scalar and gauge field masses are equal. In terms of these units the critical magnetic fields are B c2 = κ and corresponding to the upper critical magnetic field (field above which superconductivity is completely destroyed) and lower critical fields (approximate field value below which a vortex doesn't want to form) respectively. Note that in these units the critical fields, and as shown later also the flux quantum, depend on κ and hence on a. This is very important and it amounts to a choice of convention. In order to discuss results at fixed flux quanta we change conventions in a later section.
This definition for G allows one to pick the lattice structure by effectively choosing x 1 and y 2 in the definition of ω A . The external gauge field G i is then singular at the vortex positions. This definition satisfies where Φ 0 =B/n v (n v being the number of vortices in the cell) is the single vortex flux and δ (2) denotes the two dimensional delta function. The external flux requirement (32) then trivially follows. In the above r i denote the vortex positions in the lattice with coordinates (x i , y i ).
Note that this definition does not restrictB to be in the vicinity of the critical flux B c . It rather defines a vortex lattice structure with the appropriate magnetic field singularities located at the vortex cores, supporting a total fluxB. The delta function originates from the fact that close to the vortex cores ω A ≈r 2 , wherer = r − r 0 and r 0 is the coordinate position of the vortex core, so that G ≈ ẑ× r r 2 and we have the important result that Therefore, for a particularly chosen lattice at fixedB, once x 1 is specified and its relation to x 2 and y 2 is fixed (and therefore ω A (B) is known), this can be used to determine the local gauge field Q b i appearing in equations (37) from a numerical minimization problem. In terms of the parameters of our model, the external total flux over the integration domain is given byB = 2πn/κ, therefore the magnetic field carried by each unit cell is B = 2π x 1 y 2 κ where n counts the number of lattice cells in our numerical integration domain.
Some model solutions are shown in figure 1. In these figures B is the magnetic field of Q i whose total flux isB. The solutions represent Abrikosov Lattices for a generic fluxB per unit cell of both square and triangular geometries. The energy of these solutions is discussed in a later section.

Numerical procedure
The equations are solved numerically using a relaxation procedure. The derivatives are discretized using a second order central finite difference method. The accuracy of the procedure is O(10 −6 ). The results are further checked using a Newton-Rhapson method, which reproduced the same results. We impose periodic boundary conditions at the borders of the lattice and demand that the scalar field f vanish at the vortex centers.

Non-Abelian Vortex Lattices
It is the purpose of this paper to investigate the extension of the above setup to include the condensation of the additional χ field. At large enough lattice spacing the condensation of the χ field should resemble that of isolated non-Abelian vortices. However, as the lattice spacing is reduced, these self interactions will become relevant and the structure of the lattice must change.

Non-Abelian vortex -vortex forces
Let us discuss briefly the relevant aspects of the non-Abelian vortex interactions which are important for the rest of the paper.
The vortex-vortex interaction can most easily be calculated by an analysis of the fardistance behaviour of the fields and an assumption of large vortex separation (we will follow the discussion in [30]). We restrict our considerations to the case of largely separated non-Abelian vortex orientations (a general case can be found in [31]). If we consider the case of an isolated non-Abelian vortex, we must go back to equations (15) - (17) and, switching to cylindrical coordinates with radial variable ρ use the ansatz all of which are understood in a-dimensional units. Here n e is the quantum of flux carried by each vortex. With these ansatz the field equations reduce to Note that we have rescaled the field Q → Q/ √ a compared to (37), this rescaling is important in discussing the phase diagram in a later section 2 . This change of convention is important as it removes the dependence on a in the flux quanta. Therefore it allows us to change the parameter a appearing in the equations of motion without changing the flux of the vortices. Throughout the paper we will use the first convention when displaying solutions to the lattice system (this is done in order to remain in contact with [27]) and the latter convention when discussing single vortex results (such convention is adopted in [1], for example). Either convention is especially simple in discussing both aspects separately, and hence we choose to switch between both when needed.
The magnetic field is B = 1 r Q , with denoting differentiation w.r.t r. Linearising the fields around their leading behaviours far from the vortex cores, the resulting leading behaviours are given by the usual modified Bessel functions where c i are integration constants, which have leading order expansions at large ρ of the form wherec i are in general different constants from c i (these are obtained from numerical integration of the vortex profiles). As expected, the scalar field χ has a general behaviour which is very similar to the scalar field f in terms of its exponential decay. The magnetic field is B = 1 r Q , with denoting differentiation w.r.t r. The corresponding expression for the energy of the isolated vortex is where and The vortex-vortex interaction energy E int at large vortex separations, is generally calculated by subtracting from the total energy of a vortex-vortex configuration the energy of the two isolated vortices. To do so one must first take a vortex-vortex solution ansatz. The good ansatz is given by where the subscript on the fields indicates they are the field profiles of vortex 1 and 2 respectively. Therefore, to conserve the right topological properties one must take the product of the scalar field profile f . Then, at large vortex separations, it is sufficient to consider only the far field behaviours of the fields in the energy, where we understood the field profiles to be those given by equations (53) in the above. Then, inserting these expressions into (54) and subtracting the energies of each isolated vortex we find, at leading order, Using the field expansions (53), assuming similar profiles for all the χ i , and the techniques outlined in [32] this expression can be integrated to give where,ĉ is the constant vector of constants appearing from the integration of χ i and s is the vortex separation assumed to be large. If χ = 0, in the standard Abrikosov vortex case, when a = 1, at the so called critical point,c 1 =c 2 and therefore E int = 0. In the presence of the χ field, which is after all adding a scalar sector, we have an additional force channel between the vortices. When the vortices have parallel internal orientations, thenĉ·ĉ > 0 and the force is attractive. This attractive mechanism was also observed separately in a numerical study of Skyrmions [19]. When the vortices are antiparallel,ĉ ·ĉ < 0 and the force is repulsive. When they are completely transverse instead, the χ field interaction vanishes. In particular, it is no longer true that at a = 1 critical vortices have zero interaction energy (we expect that at this order howeverc 1 =c 2 should remain true). This means that the condition on type II vortices and the point of criticality is generally expected to be different than before. In field terms, there is no BPS condition at a = 1 even though it doesn't exclude the possibility that there is one since E int might still vanish in general for specifically chosen values of the parameters. This should now be determined numerically from the full field profiles. We will not perform a numerical investigation on the values of the parametersc i in this paper as we generally solve the full field equations. However, this would be an interesting avenue of research per se.
Note also that the change in nature of the interaction force between non-Abelian vortices with respect to their internal orientation is indicative of the possibility of alternative (as in besides the parallel orientation case) lattice structures. We will devote a section of this paper to the exploration of such structures.

Solutions
Now we wish to find the full 2D solutions which describe lattices of non-Abelian vortices. We therefore switch back to the conventions of (37). In the presence of the χ field the equations we must solve are We use an identical numerical procedure as outline above in order to find the solutions. Figures 2 and 3 show solutions at large lattice spacing (of square and triangular geometries respectively). These represent a lattice of isolated non-Abelian vortices, with the χ field condensing only in the core of the Abrikosov vortices and quickly decaying to zero outside. As we bring the vortices closer together by decreasing the lattice spacing, self-interactions of the χ field become important. The field gradually and smoothly lifts and becomes nonvanishing over the whole lattice (see figure 4 and figure 5). These solutions do not represent a lattice of non-Abelian vortices, in this case the χ field is non zero over the whole lattice, which implies the internal orientational moduli are delocalised from the vortex cores.

Anti-parallel configurations on the square lattice
In addition to the solutions we have found in the previous subsection, we have considered more general orientations of the χ field on the square lattice. Since we deal with a numerical      relaxation procedure, we must make sure we start with an initial configuration which is close to an actual solution. To this end, it is difficult to imagine a case where the parallel configuration from lattice site to lattice site would not be the most stable state. However, it may be possible that such alternative configurations could achieve some meta-stability for certain cases of the lattice spacing and constants. For example, we are unable to rule out the possibility that the anti-parallel configuration of χ fields (analogous to an anti-ferromagnet) is metastable. In figures (7) and (8) we show the solutions for the f , χ, and B fields for two values of the lattice spacing. These are generated using a similar ansatz for the f and B fields as the cases above, and an anti-parallel configuration for χ.
In order to test the stability of these anti-parallel configurations we may consider small perturbations of the χ field, and observe the response. For the case of solution shown in (7) we find that small perturbations leave χ disordered from site to site, with directors pointing in random directions. This is as expected since the large spacing between lattice sites prevents interaction of the χ field localized at each site. We conclude that in this particular case, no configuration of χ directors is preferred. This is simply a confirmation that, for large lattice spacing (or small vortex interaction) the system has independent orientational moduli on each vortex site (see the following section).
On the other hand, for closer lattice spacings, such as the case presented in (8), we find that a small perturbation of the initial anti-parallel configuration, leads to an instability. In this case, the χ directors reorient in the parallel configuration. In figure (9) we illustrate the field configuration at various time steps in the minimization procedure. We generate the initial perturbation of χ by introducing a small wave form in χ with Fourier coefficients generated randomly. We conclude that the anti-parallel configuration is not stable in this case either.
Of course we cannot argue on general grounds that stable anti-parallel configurations for χ may exist for particular ranges of the parameter space and lattice spacings. For all cases we have considered we find that the either the parallel or the disordered configuration for the χ field are the only stable solutions, depending on lattice spacing.

Searching for other meta-stable configurations
In an effort to be more complete with the solution space for the χ field configuration, we attempted to find additional configurations for the χ lattice. A priori it is hard to imagine any relaxation seed besides the parallel or anti-parallel configurations, therefore we resorted to a more statistical approach based on random initial orientations per lattice site, repeating this many times for different configurations. This was done by Fourier decomposing the χ field with randomly generated coefficients. The minimization procedure was then carried out with this initial ansatz for χ, and was continued until convergence was achieved. This procedure was repeated several times for several different values of the parameter space. In figure (10) we show the relaxation procedure for a particular initial χ configuration. In all cases considered, the χ configuration relaxed to either a disordered state, or a parallel configuration depending on the spacing between lattice sites. As previously mentioned, this is no surprise when interaction strength between lattice sites are considered.
We mention that this analysis is far from complete, and the possibility of parameter ranges where other metastable configurations appearing cannot be ruled out. In order to find them however, it is necessary to start with an intelligent seed and hence a clear idea of what these might look like.

Energy
In this section we present the numerical results regarding the energetics of our solutions. We compare the energies of square and triangular lattice configurations and, more importantly, those with and without the condensation of the χ field in or around the vortex cores. In the dimensionless units used, the energy functional is where   Here and In the absence of the χ sector it is a well-known result that the triangular lattice has lower energy per unit cell compared to the square one. We have indeed confirmed this numerically. For example, for the solutions shown in figure (1), comparing the energies per unit cell we obtained where the subscript on the energy indicates the geometry of the lattice. The energy of the triangular geometry is lower, as is well known. This result holds in general for all the ranges of parameters considered in this paper, we present some representative data in the table shown in figure 6.
A more interesting and new comparison is that between Abrikosov lattices and non-Abelian vortex lattices. The energy for the solution shown in Figure 2, for example, is E 2 = 7.991. The corresponding energy without χ is E χ=0 2 = 7.999. Similarly, the solution of figure 3 has energy E 3 = 7.988 and that without χ (that of the corresponding pure Abrikosov lattice) E χ=0 2 = 7.997. Therefore, once again, even in the presence of the additional field, the triangular lattice has least energy per unit cell. Note that it is energetically favourable for the lattice to nucleate the χ field in the vortex cores. This is not surprising, these solutions represent an ideal non-Abelian vortex lattice with well separated non-Abelian vortices. These vortices have lower energy than the corresponding ANO vortices [23], a lattice of them should therefore also have lower energy if these don't interact. Interestingly the energy is lower but not significantly so, this was also found in [23] which allowed an approximate kink solution to be constructed.
In this system the non-Abelian vortex lattices of triangular geometries are therefore the lowest energy solutions, in the limit of large separation of vortices. This result holds for all numerical values of the parameters we investigated, and seems to be a general result of this setup. We present some numerical data of representative solutions in the table shown in figure 11.
When the non-Abelian vortices are more tightly packed, for smaller lattice spacing, we must verify whether the solutions with a delocalised χ field over the whole lattice are actually energetically preferred over solutions in which χ vanishes. We numerically checked that this was the case, some representative values are shown in the table in figure 12 (these lattice spacing values correspond to the solutions shown in figures 4 and 5). Therefore, for all lattice spacings, even when the χ field does not vanish outside of the vortex cores, the system always prefers to have the orientational field than not.

Phase diagram of ideal non-Abelian vortex lattices
As shown in section (2.1) the condition of type II superconducting vortices changes in the presence of the χ field. In particular, we must check for what values of the parameters it is energetically convenient for a non-Abelian vortex of flux Φ to form a lattice of n vortices of flux Φ/n. Only then will the lattice formation be preferred for this type of superconductor. This section is devoted to carrying out such a numerical analysis and therefore mapping out the phase space of this system. This will determine the kind of superconductor that exists (type I or type II) as a function of parameters.
Our strategy is the following: we solve the isolated vortex ODEs for non -Abelian vortices as a function of the parameters a and b for a single vortex of 2 flux quanta and calculate its energy, we then compare this energy to that of 2 vortices of a single flux quantum, at the same parameter values assuming no interaction (large separation). If the energy of the latter is smaller, then the system will want to be in a vortex lattice state and is a superconductor of type II. As discussed previously, the transition to this kind of superconductivity for the standard Abrikosov system (corresponding to χ = 0 or equivalently b = 0), happens at a = 1, when the gauge boson and scalar field masses are equal. This is commonly known as the BPS point, where the vortex interaction energy vanishes. For this comparison, we want to keep the flux constant and vary only the parameters of the equations. Therefore we will switch back to the rescaled convention adopted in section (2.1). In this convention, the quantum flux is Φ = 2π √ 2n, and importantly is independent of a (this in turns makes the upper critical field also independent of a).
To study isolated non-Abelian vortices we must go back to equations (48) -(50) and solve them with the boundary conditions Our comparison is therefore between energies of solutions with n = 2 and n = 1 as we vary a and b. The result of this analysis is shown in figure 13. We can distinguish two main sections of this plot. The section marked II in this plot corresponds to a region in which the energy of two separated vortices of a single flux is less than that of a single vortex of flux 2. This region corresponds to type II superconductivity. We see immediately that for b = 0, which corresponds to the usual Abrikosov vortex with χ = 0, the critical point of this phase is at a = 1, as discussed in section (37). As we vary b we enter the non-Abelian vortex solutions. Here we find that there is a whole dome of parameter space of type II superconductive behaviour. We can label this as non-Abelian type II superconductivity. Beyond this region, in the region on the plot called I, the high flux single vortex solution is energetically preferred and the system is no longer of type II. It appears that as we increase b, and make the system increasingly more non-Abelian (in the sense of increasing the core value of χ), the critical point of transition in the a parameter is lowered, i.e. the extra scalar field means criticality is achieved at higher Higgs scalar masses, and no longer when this is equal to the gauge boson mass. This was the result anticipated in the discussion of section (37), but it is now mapped numerically. The additional scalar channel provides an attractive force between the vortices, hence the χ field pushes the system towards type I. Note that the upper value of b ≈ 2 here is actually physical in the sense that above this line we violate the vacuum condition (20) and would enter a phase in which the χ field wants to condense in the vacuum. We repeated this analysis for several values of c and β and find similar results for the shape of the critical line. The main difference being that the physical limit line on the b axis changes as one changes c or β.
This diagram illustrates neatly where the solutions discussed in section 2 (which adopted the other convention) actually exist. They are representatives of the energy minimizing solutions, at fixed flux, that exist in the region of type II non-Abelian superconductivity.   Finally, as an aside, we re-instanted the full 2D convention and also used the isolated vortex to compare our numerical lattice solutions with the 1D radial plots. First we compare the lattice solutions with single vortices carrying the total flux of the unit cell, this is shown in figure 14. In this figure the single vortex (solid line) carries 4 flux quanta, while the isolated vortices of the lattice each carry just one. In figure 15 we show the isolated vortex (solid line) carrying the same flux as a vortex of the lattice (not the whole flux of the lattice). As seen by the figure there is a remarkable agreement between the isolated and lattice solutions demonstrating that our numerical procedure for the 2D solver reproduces the radial ODE solutions well.

Low Energy theory
The solutions represent lattices of non-Abelian vortices. In the core of the vortices, where ψ = 0, the χ field condenses and breaks the O(3) global symmetry down to U (1), the rotations in the (1,2)-plane of internal space. If the lattice spacing is sufficiently large for the vortices to be well-separated, as per the solutions shown in figures 2 or 3, the low energy theory will be a lattice of CP (1) non-linear sigma models localised on each vortex site. In order to see this let us consider a small neighbourhood around a vortex core (defined approximately as a region of area in which χ is non-vanishing) in the (x, y)-plane and take the ansatz where S i (t, z) are orientational fields depending on the world sheet coordinates and satisfying S i S i = 1. Inserting this ansatz into the action gives the low-energy action for the orientational moduli as with k = (t, z), and α = µ 2 2β v dxdyχ 2 , where the integral runs over the neighbourhood of the vortex site. The overall action is then the sum of all the vortex site contributions so that where the index i runs on over each vortex site and α i refers to the integral α performed over connecting neighbourhoods of each vortex site. It is important to stress that for this kind of lattice each core moduli are independent to rotate on each vortex site since the χ fields are well localised.
When the lattice spacing is lowered, such that the χ field delocalises from the vortex cores (see figure 4), only part of these orientational degrees of freedom survive as gapless excitations. These are the global rotations of the internal orientation throughout the whole lattice, namely locked rotations in internal space of all vortex sites. In this case each vortex site is not independent to rotate freely in internal space, the only gapless excitation is for each vortex to lock onto each other and rotate equally throughout the lattice. Then, the low energy theory is described by the CP (1) model where the integral now runs over the whole lattice. This is similar to what happened for Skyrmions in [19].

Conclusions
This paper extends the previously known features of the Abrikosov string supporting non-Abelian moduli proposed in [1]. In particular, it concludes the investigation on the nature of the superconducting regimes in it, with a complete map of superconducting type in parameter space. This model is a toy model of several physical systems, ranging from dark matter considerations [25] to confinement in QCD, wherever solitonic solutions of this type play an important role [13] [33] [34]. Previous results in this model showed that such isolated vortices exist, and some studies on the properties of these vortices were made [1] [23]. This paper extends the investigation to the full analysis of what type of superconductivity is involved, and what are the real low lying energy solutions the system adopts, in the case of parallel orientation. In particular, we first provide analytical evidence that the additional scalar field sector, providing the orientational degrees of freedom, acts as an additional attractive channel between isolated vortices. In turn, this modifies the usual a = 1 BPS point of no vortex interaction energy. We then demonstrate numerically the existence of periodic arrays of parallely oriented non-Abelian vortices, both in the limit of large separations and when they are tightly packed, in square and triangular geometries. We map out the phase diagram of this superconductor and deduce the critical line of phase transition between type 2 and type 1 superconductivity in parameter space. When the system is in the type 2 region, the lowest energy solutions at fixed external flux and lattice spacing correspond to triangular non-Abelian vortex lattices. When the lattice spacing is large, the low energy theory of this system is described by an array of CP (1) theories described by the local independent rotations of the orientational degrees of freedom at each lattice site. When the spacing is reduced, the self interactions of such degrees of freedom becomes important and the only remaining gapless excitation is a locked rotation of all lattice sites. This is similar to a ferromagnetic phase transition, in which independent spins couple into a single orientation under an applied magnetic field. With this in mind it would be interesting to allow all possible "spin" orientations of the vortices, rather than restricting to parallel ones.
We have considered alternate possible configurations of the spin lattice. For the "antiferromagnetic" case (anti-aligned spins per vortex site) we determined that the configuration is unstable in the case of small lattice spacing. We also tried starting with random initial orientations and repeating the relaxation procedure with the hope that the system would spontaneously find alternative spin configurations. No such configurations were found. Based on this analysis we can conclude that the parallel orientation lattice is the most stable solution at small lattice spacings. However, we admit that this analysis is far from concluded (alternative spin configurations might arise, for example, for lattice geometries besides the triangular or square case). A more complete investigation of other potential stable configurations will have to wait for future projects.