Crystals of superconducting Baryonic tubes in the low energy limit of QCD at finite density

The low energy limit of QCD admits (crystals of) superconducting Baryonic tubes at finite density. We begin with the Maxwell-gauged Skyrme model in (3+1)-dimensions (which is the low energy limit of QCD in the leading order of the large N expansion). We construct an ansatz able to reduce the seven coupled field equations in a sector of high Baryonic charge to just one linear Schrodinger-like equation with an effective potential (which can be computed explicitly) periodic in the two spatial directions orthogonal to the axis of the tubes. The solutions represent ordered arrays of Baryonic superconducting tubes as (most of) the Baryonic charge and total energy is concentrated in the tube-shaped regions. They carry a persistent current (which vanishes outside the tubes) even in the limit of vanishing U(1) gauge field: such a current cannot be deformed continuously to zero as it is tied to the topological charge. Then, we discuss the subleading corrections in the 't Hooft expansion to the Skyrme model (called usually \mathcal{L}_{6}$, $\mathcal{L}_{8}$ and so on). Remarkably, the very same ansatz allows to construct analytically these crystals of superconducting Baryonic tubes at any order in the 't Hooft expansion. Thus, no matter how many subleading terms are included, these ordered arrays of gauged solitons are described by the same ansatz and keep their main properties manifesting a universal character. On the other hand, the subleading terms can affect the stability properties of the configurations setting lower bounds on the allowed Baryon density.


Introduction
Exact analytic results on the phase diagram of the low energy limit of QCD at finite density and low temperatures are extremely rare (it is often implicitly assumed that they are out of reach of the available techniques). This fact, together with the non-perturbative nature of low energy QCD, is one of the main reasons why it is far from easy to have access to the very complex and interesting structure of the phase diagram (see [1]- [4], and references therein) with analytic techniques.
One of the most intriguing phenomena that arises in the QCD phase diagram at very low temperatures and finite Baryon density, is the appearance of ordered structures like crystals of solitons (as it happens, for instance, in condensed matter theory with the Larkin-Ovchinnikov-Fulde-Ferrell phase [5]). From the numerical and phenomenological point of view, ordered structures are expected to appear in the low energy limit of QCD (see, for instance, [6]- [10], and references therein). The available analytic results have been found in (1 + 1)-dimensional toy models and all of them suggest the appearance of ordered structures 1 of solitons.
Even less is known when the electromagnetic interactions arising within these ordered structures are turned on. Analytic examples of crystals of gauged solitons with high topological charge in (3 + 1)dimensions in the low energy limit of QCD would reveal important physical aspects of these ordered phases. The only available analytical examples 2 are derived either in lower dimensions and/or when some extra symmetries (such as SUSY) are included (see [12], and references therein).
We search for analytic solutions (despite the fact that these questions can be addressed numerically, as the previous references show) because a systematic tool to construct analytic crystals of gauged solitons can greatly enlarge our understanding of the low energy limit of QCD: the analytic tools developed here below disclose novel and unexpected phenomena. The gauged Skyrme model [13], [14], which describes the low energy limit of QCD minimally coupled with Maxwell theory at the leading order in the 't Hooft expansion [15]- [18] (for two detailed reviews see [19] and [20]), will be our starting point. Using the methods introduced in [21], [22], [23] and [24] we will construct analytic gauged multi-soliton solutions at finite Baryon density with crystalline structure and high topological charge. These crystals describe ordered arrays of superconducting tubes in which (most of) the topological charge and total energy are concentrated within tube-shaped regions 3 . They carry a persistent current (which vanishes outside the tubes) which cannot be deformed continuously to zero as it is tied to the topological charge.
These regular superconducting tubes can be considered as explicit analytic examples of the superconducting strings introduced in [27]. The spectacular observable effects that such objects could have (see [28], and references therein) together with the fact that these objects can be constructed using natural ingredients generate a huge interest both theoretically and phenomenologically. However, until now, there are very few explicit analytic (3 + 1)-dimensional examples built using only ingredients arising from the standard model. In fact, the present superconducting tubes appear in the low energy limit of QCD minimally coupled with Maxwell theory.
Then we move to the subleading correction to the Skyrme model in the 't Hooft expansion (see [29]- [32] and references therein). Although one could believe that such complicated corrections could destroy the nice analytic properties of the crystals of superconductive Baryonic tubes, we will show that no matter how many subleading terms are included, these ordered arrays of gauged solitons are described by the very same ansatz and keep unchanged their main properties manifesting a clear universal character. On the other hand, the stability properties of these crystals of superconducting Baryonic tubes can be affected in a non-trivial way by the subleading terms in the 't Hooft expansion.
The paper is organized as follows. In the second section the general field equations will be derived and the definition of topological charge will be introduced. In the third section, the ansatz which allows to solve analytically the field equations (no matter how many subleading terms are included) in the ungauged case in a sector with high topological charge will be discussed. Then, it will be explained how this ansatz can be generalized to the gauged case with the inclusion of the minimal coupling with the Maxwell field. The physical properties of these gauged crystals and their universal character will be analyzed. We will also study how the subleading terms can affect the stability properties of the configurations setting lower bounds on the allowed Baryon density. In the last sections, some conclusions and perspectives will be presented.

The gauged generalized Skyrme model
The starting point is the action of the U(1) gauged Skyrme model in four dimensions, which corresponds to the low energy limit of QCD at leading order in the 't Hooft expansion: where K and λ are the Skyrme couplings, d 4 v is the four-dimensional volume element, g is the metric determinant, m is the Pions mass, A µ is the gauge potential, ∇ µ is the partial derivative and σ i are the Pauli matrices. In Eq. (1), L corr represents the possible subleading corrections to the Skyrme model which can be computed, in principle, using either Chiral Perturbation Theory (see [33] and references therein) or the large N expansion [34], [35]. The expected corrections have the following generic form and so on [29], where the c p (p ≥ 6) are subleading with respect to K and λ.
A natural question arises here: in which sense these correction terms are subleading with respect to the original Skyrme model? First of all, we would like to remark that this question does not affect directly the present construction since the analytic method presented here allows to construct exact solutions no matter how many "subleading terms" are included (as it will be shown in the following sections). However, from the physical point of view, the above question is very interesting.
In principle, as remarked in [30], [31] and [32], one should expect generically higher-derivative terms of the chiral field U in the low-energy limit of QCD. Due to the fact that each term is larger in canonical dimension, dimensional constants to the same power minus four must go with each of them. These constants are expected to be proportional to the mass scale of the degrees of freedom integrated out of the underlying theory. Therefore, as long as the energy scales being probed are much smaller than the lowest mass scale of a state that was integrated out, the higher-derivative expansion may make sense and thus converge. Another intuitive argument is due to 't Hooft and Witten (in the classic references [15], [34] and [35]) which argued that in the large-N limit, QCD becomes equivalent to an effective field theory of mesons and which the higher order terms with respect to the non-linear sigma model (NLSM henceforth) action are accompanied by inverse power of N (where here N is the number of colors).
The analysis here below will clarify that no matter how many further subleading terms are included, the pattern will never change. In particular, using precisely the same ansatz discussed in the next section for the SU(2) -valued Skyrmionic field, the field equations for the generalized Skyrme model with all the corrections included always reduce to a first order integrable 4 ordinary differential equation.
Moreover, if the minimal coupling with Maxwell theory is included, the gauged version of the field equations for the SU (2) Without such a reduction, even the numerical analysis of the electromagnetic properties of these (3+1)dimensional crystals of superconducting tubes would be a really hard task (which up to now, has not been completed to the best of our knowledge). While, with the present approach, the numerical task to analyze the electromagnetic properties of these crystal of superconducting tubes is reduced to a linear Schrödinger equation with an explicitly known potential.

Field equations
The field equations of the model are obtained varying the action in Eq. (1) w.r.t. the U field and the Maxwell potential A µ . To perform this derivation it is useful to keep in mind the following relations where δ U denotes variation w.r.t the U field, and Here we have used and we have defined From the above, the field equations of the gauged generalized Skyrme model turns out to be together with where the current J µ is given by

Energy-momentum tensor and topological charge
Using the standard definition we can compute the energy-momentum tensor of the theory under consideration with T U(1) µν the energy-momentum tensor of the Maxwell field According to Eq. (8), a direct computation reveals that The topological charge of the gauged Skyrme model is given by [14], [36]: Note that the second term in Eq. (10), the Callan-Witten term, guarantees both the conservation and the gauge invariance of the topological charge. When Σ is space-like, B is the Baryon charge of the configuration.

Crystals of superconducting Baryonic tubes
In this section we will show that the gauged generalized Skyrme model admits analytical solutions describing crystals of superconducting Baryonic tubes at finite density.

The ansatz
Finite density effects can be accounted for using the flat metric defined below: where 4π 3 L 3 is the volume of the box in which the gauged solitons are living. The adimensional coordinates have the ranges 0 ≤ r ≤ 2π , 0 ≤ θ ≤ π , 0 ≤ φ ≤ 2π.
From Eqs. (13) and (14) the topological charge density reads and therefore, as we want to consider only topologically non-trivial configurations, we must demand Now, the problem is to find a good ansatz which respect the above condition and simplify as much as possible the field equations. A close look at Eq. (5) (see Appendix II for its explicit form in terms of α, Θ and Φ) reveals that a good set of conditions is A suitable choice that satisfies Eqs. (15) and (16) is the following [23]: Additionally some other useful relations are satisfied by the above ansatz, namely

Solving the system analytically
The identities in Eqs. (16) and (18) satisfied by the ansatz in Eq. (17) greatly simplify the field equations keeping alive the topological charge. This can be seen as follows.
Firstly, a direct inspection of the field equations reveals that all the terms which involve sin 2 Θ are always multiplied by ∇ µ Φ∇ µ Φ so that all such terms disappear.
Secondly, since Θ is a linear function, in all the terms in the field equations ∇ µ Θ∇ µ Θ becomes just a constant.
Thirdly, since the gradients of α, Θ and Φ are mutually orthogonal (and, moreover, ∇ µ Φ is a light-like vector), all the terms in the field equations which involve Fourthly, the above three properties together with Eq. (18) ensures that two of the three field equations of the generalized Skyrme model are identically satisfied (see Appendix I and Appendix II).
It is also worth to emphasize that the four properties listed here above are true no matter how many subleading terms (L 10 , L 12 and so on) are included in the generalized Skyrme action. For the above reasons, the three non-linear coupled field equations of the generalized Skyrme model in Eq.
(5) with the ansatz in Eq. (17) are reduced to the following single ODE for the profile 5 α: This is already a quite interesting fact in itself. Moreover, the above analysis clearly shows that it will remain true even including further subleading term. What is really remarkable is that Eq. (19) can always be reduced to a first order ODE: which is explicitly solvable in terms of generalized Elliptic Integrals [37]. Here E 0 is a positive integration constant and X = dX dr . Therefore Eq. (20) implies that, with the ansatz defined in Eq. (17), the field equations are integrable and reducible to the following quadrature 6 : It is also worth to emphasize that the integration constant E 0 can be chosen in such a way that, first of all, α never changes sign (which is a necessary condition for stability) and, secondly, the topological charge is B = np (as we will show in the following subsection).
Quite surprisingly, these very intriguing properties of the ansatz are not destroyed by the inclusion of the minimal coupling with Maxwell field. The coupling of the generalized Skyrme model with the Maxwell theory is introduced replacing the partial derivatives acting on the SU(2)-valued scalar field U with the following covariant derivative 5 It is interesting to note that the terms in the field equations arising from L6 in the generalized Skyrme model vanish identically due to the properties of the ansatz in Eqs. (16), (17) and (18). On the other hand, such a term can affect the stability properties of the solutions, as we will see below. 6 The identities in Eqs. (16) and (18) satisfied by the ansatz in Eq. (17) ensures that (even if the subleading corrections L10, L12 and so on are included) the ansatz in Eq. (17) will always reduce the three coupled non-linear field equations to a single first order ODE for the profile α.
A straightforward computation shows that the above replacement in Eq. (23) is completely equivalent to the replacement here below (in terms of α, Θ and Φ) It is worth to emphasize that D µ Φ determines the "direction" of the electromagnetic current (as it will be discussed below).
Obviously, when the derivative is replaced with the Maxwell covariant derivative (as defined in Eq. (23) or, equivalently, in Eq. (24)), in the field equations of the gauged generalized Skyrme model many new terms will appear which couple the SU(2) degrees of freedom with the U(1) gauge potential A µ . Thus, one may ask: Which is the best choice of the ansatz for the gauge potential A µ which keeps as much as possible the very nice properties of the ansatz of the SU (2)-valued scalar field in Eqs. (15), (16) and (18) which allowed the complete analytic solutions in the previous case?
In order to achieve this goal, it is enough to demand The above conditions determine that the Maxwell potential A µ must be of the form [23]: From the expressions of L µ (see Appendix I) one can see that, despite the explicit presence of A µ in the U(1)-covariant derivative, the three field equations of the gauged generalized Skyrme model still reduce to the Eq. (19). The reason is that all the potential terms which, in principle, could couple the SU(2)-valued scalar field U with A µ in the field equations actually vanish due to the identities in Eqs. (16), (18), (25) and (26) satisfied by the choice of our ansatz (that is why we have chosen the ansatz in that way). One can verify easily that the four Maxwell equations are reduced to the following single PDE: where Ω(α) is given by Note also that Eq. (28) can be written as a periodic two-dimensional Schrödinger equation Therefore, with the ansatz defined in Eqs. (12), (17) and (27)  A µ = 0, the current does not vanish. Such a residual current cannot be deformed continuously to zero, and the reason is that the only way to "kill" it would also kill the topological charge but, as it is well known, there is no continuous transformation which can change the topological charge. We will return to this very important issue when we address the superconducting nature of the Baryonic tubes.

Boundary conditions and Baryonic charge
From Eq. (10) we can compute the energy density of the configurations presented above, which turns out to be The above expression can be written conveniently as and one can check that the topological charge becomes Assuming the following boundary condition for u and α and taking into account that q is an odd integer, the topological charge becomes It is worth to stress here that (unlike what happens in the case of a single Skyrmion in a flat space without boundaries, when the boundary conditions are just dictated by the condition to have finite energy) when finite density/volume effects are taken into account the choice of the boundary conditions is not unique anymore. A very important requirement that any reasonable choice of boundary conditions must satisfy is that the integral of the topological density (which, of course, by definition is the topological charge itself) over the volume occupied by the solutions must be an integer. If this condition is not satisfied, the configurations would not be well defined. Hence, the boundary conditions should be fixed once and for all within the class satisfying the requirement described here above: our choice is the simplest one satisfying it. Now one can note that, according to Eqs. (21) and (22), the integration constant E 0 is fixed in terms of n through the relation It is easy to see that the above equation for E 0 always has a real solution as the integrand interpolates from very small absolute values (when E 0 is very large in absolute value) to very large (when E 0 is such that the denominator can have zeros). Hence, one can always find values of E 0 able to satisfy Eq. (33).

Baryonic crystals at finite density and its superconducting nature
From Eq. (9) one can compute the energy density E of the configurations defined in Eqs. (12), (17) and (27), and this turns out to be wherẽ It is interesting to note that, despite the fact that the term L 6 does not contribute to the field equations (as it has been already emphasized), it does contribute to the energy-momentum tensor. In order to have a positive definite energy density a necessary condition is c 6 ≥ 0.
On the other hand, the U(1) current in Eq. (7), in the ansatz defined by Eqs. (12), (17) and (27), with Ω(α) defined in Eq. (29). From the expression of the current in Eq. (35) (see Appendix I for the explicit form of the components of the current) the following observations are important.
1) The current does not vanish even when the electromagnetic potential vanishes (A µ = 0).
2) Such a "left over" is maximal where the energy density is maximal and vanishes rapidly far from the peaks as the plots show (see Fig. 1 and Fig. 2).
3) J (0)µ cannot be turned off continuously. One can try to eliminate J (0)µ either deforming α and/or θ to integer multiples of π (but this is impossible as such a deformation would kill the topological charge as well) or deforming Φ to a constant (but also this deformation cannot be achieved for the same reason). Note also that, as it happens in [27], Φ is defined modulo 2π (as the SU (2) valued field U depends on cos Φ and sin Φ rather than on Φ itself). This implies that J (0)µ defined in Eq. (35) is a superconducting current supported by the present gauged tubes. Moreover, these properties are not spoiled by any of the higher order corrections, parameterized by c p .
The plots of the energy density and the current clarify the physical interpretation of the present multi-solitonic configurations. In Fig. 1 and Fig. 2 we have chosen K = 2, λ = 1, c 6 = c 8 = 1 5 , m = 0 and q = p = 1. The components of the electric and magnetic fields can be also computed and are Figure 1: From left to right we can see the energy density E, the magnetic field, the electric field and the current J µ for a configuration with Baryon charge B = 1. The electric and magnetic fields vanish at the peaks of the energy density while the current takes its maximum value.
given by The current is concentrated in the tube-shaped regions defined where (most of) the E is contained, and vanishes outside the tubes. The maximum values of E and the current coincide in the lattice, which is periodic in r, θ and perpendicular to the φ direction, along which the strings exist.

About the existence of exact crystals and the universality of the ansatz
In the previous sections we have shown that the low energy limit of QCD supports the existence of crystals of superconducting Baryonic tubes. Of course, this result is very technical in nature and, a priori, it is not clear whether or not one could have expected the appearance (in the low energy limit of QCD) of topological defects supporting superconducting currents. Here we will give an intuitive argument which justifies why one should have expected the existence of such defects.
The first necessary (but, in general, not sufficient) condition that must be satisfied in order to support the existence of superconductive currents in a relativistic context is the existence of a massless excitation which can be coupled consistently to a U(1) gauge field (see the pioneering paper [27]).
According to Eqs. (13) and (14) the SU(2) valued Skyrme field U describes the dynamical evolution of three scalar degrees of freedom α, Θ and Φ which are coupled through the non-linear kinetic terms typical of Skyrme like models (see Eq. (50) which is the explicit expression of the Skyrme action in terms of α, Θ and Φ: of course such action is equivalent to the usual one written implicitly in terms of the SU(2) valued field U ). This fact hides a little bit which is the "best candidate" to carry a superconductive current since our intuition is built on models where the interactions appear in potential terms (like in the Higgs model or in the Ginzburg-Landau model) and not in "generalized kinetic terms" as in the present case. So, the question is: how can we decide a priori which whether or not there is an excitation able to carry a persistent current? In other words, which one of the three degrees of freedom α, Θ and Φ associated to the SU(2) valued scalar field U can be a carrier of a superconductive current? In what follows we will detail the intuitive arguments that lead us to consider Φ as the most natural choice.
To illustrate our argument, let us first consider the simpler and well known case of two scalar fields Ψ i , with i = 1, 2, interacting with a quartic potential in a SO(2) invariant way, In order to disclose which degree of freedom is a natural candidate to carry a chiral current, we can write

then, Eq. (36) becomes
From Eq. (37) it is clear that R can not be chiral field because of the presence of a non-trivial potential term that only depends on R and generates a natural mass scale in the dynamics of R (there should be no characteristic mass scale in a superconducting current). On the other hand, χ (which represents a phase and so is defined only modulo 2π) describes excitations along the valley of the potential, and, consequently, is a more suitable candidate to carry a superconductive current. Of course, all of this is well known in the analysis of the Higgs and Goldstone modes, but this short review helps to identify the correct chiral field in our case in which there is no potential to look at (as the interactions happen in non-linear kinetic-like terms). Moreover, the above Lagrangian can also be naturally coupled to a U(1) gauge field as follows: It is easy to see that the U(1) current J µ arising from the above action is proportional to J µ ∼ (∂ µ χ − eA µ ). These are part of the main ingredients of [27] to build topological defects supporting superconducting currents. Hence, the fingerprints to identify the degree of freedom (call it for convenience Φ * ) suitable to carry the superconducting current are, firstly, that such a degree of freedom Φ * should only appear with a kinetic term in the full action of the theory (as χ in Eq. (37) here above) and without any other explicit non-linear term involving Φ * itself 7 . Secondly, the coupling of the theory with a U(1) gauge theory should only affect the kinetic term of the field Φ * . Clearly, the above requirements allow to identify χ as the field Φ * candidate to be a carrier of a superconducting current in the above example.
What happens in the present case? Obviously, in the Skyrme case there is no interaction potential which is responsible for the interactions term: in all the Skyrme-like models the non-linear behavior is related with generalized kinetic terms. Nevertheless, as can be seen in Eq. (50), the α and Θ fields have explicit non-linear interaction terms, all of them proportional to sin 2 α and/or sin 2 Θ.
Hence (although in the present case there is no potential which clearly allows to identify the "proper valleys and Goldstone modes"), it is clear that neither α nor Θ can be the analogue of χ in the previous example. The reason is that if one sets α to a generic constant value Θ will still have non-linear interaction terms and viceversa if one sets Θ to a generic constant value α will still have non-linear interaction terms 8 . Consequently, the fields α and Θ are rather similar to the field R than to the field χ in the previous example. The field Φ on the other hand, has precisely the same characteristics as the field χ in the previous example: it is only defined modulo 2π (since U depends on Φ only through sin Φ and cos Φ) and moreover it appears in the action only with the corresponding kinetic term. Thus, if you set the other fields α and Θ to generic constant values, then the field 7 Thus, when one set to constant values all the other degrees of freedom of the full action (but Φ * itself) then Φ * behaves as a massless field. Indeed, this is the case for the field χ in the above example. 8 Here, "generic constant values" mean different from nπ as otherwise there would be no kinetic term at all in the action (as it happens when one sets R = 0 in the action in Eq. (37)). Thus such values are not relevant for the present analysis.
Φ can behave as a chiral massless field. Thus, a priori, one should have expected that also in the (generalized) Skyrme model superconducting Baryonic tubes should be present. Moreover, the minimal coupling of the (generalized) Skyrme model(s) with the Maxwell theory is defined by the following . From the viewpoint of the α, Θ and Φ the minimal coupling rule is completely equivalent to change in the action only the derivatives of Φ as follows Hence, also from the viewpoint of the interaction with the Maxwell theory, the field Φ is the analogue of χ in the previous example and this is exactly what we need, according to Witten [14], to have a superconducting current.
As a last remark, at a first glance, one could also argue that the presence of a mass term for the Pions should destroy superconductivity. In fact, this is not the case since, in terms of α, Θ and Φ, the mass term for the Pions is m 2 π (1 − cos α), and it only affects α (in the same way as a mass term in the previous example would only affect R but would not set a mass scale for χ).
These are the intuitive arguments which strongly suggest a priori that it certainly pay off to look for superconducting solitons in the generalized Skyrme model(s) and that Φ should be the superconducting carrier.
Furthermore, by requiring that ∇ µ Φ∇ µ Φ vanishes (as it is expected for chiral fields), the field equations are enormously simplified (see Eqs. (51), (52) and (53) on Appendix II). This simplification occurs not only on the Skyrme model case, but even if higher derivative order terms are considered, as we have already discussed.

Stability analysis
One of the most intriguing results of the present framework is that the physical properties of these superconducting Baryonic tubes remain the same no matter how many subleading terms are included in the generalized Skyrme model 9 . In other words, these topologically non-trivial configurations are almost "theory independent". As it has been already emphasized, this happens since the ansatz defined in Eq. (17) works in exactly the same way without any change at all no matter how many higher order terms are included in the generalized Skyrme action. In particular, the field equations will always be reduced to a single integrable ODE for α and the corresponding configurations will describe superconducting tubes. Hence, the present topological gauged solitons are likely to be a universal feature of QCD as they stay the same at any order in the large N expansion.
To give a flavor of why such a property is so surprising, let us consider the 't Hooft-Polyakov monopole [39], [40]. The stability of these configurations in the Georgi-Glashow model is of course very well understood. However, if one deforms even slightly the theory, the properties of the 't Hooft-Polyakov monopole are going to change as well (see, for instance, [41] and references therein). To give just an example: in [41] the authors considered a very natural correction to the Georgi-Glashow model which leads to a non-spherical deformation of the 't Hooft-Polyakov monopole (so that, in particular, the ansatz for the 't Hooft-Polyakov monopole must be changed accordingly). Consequently, the shape of non-Abelian monopoles is also going to change when these types of deformations of the Yang-Mills theory are included. On the other hand, the superconducting Baryonic tubes constructed here keep their properties at any order in the large N expansion. To the best of authors knowledge, these are the first examples of "universal" gauged solitons in the low energy limit of QCD described by an ansatz able to survive to all the subleading large N corrections. Indeed, the subleading corrections to the generalized Skyrme model will only change slightly the plot of α(r) keeping unchanged the plots and the properties of the superconducting currents and of the energy density (see Figure 3). Here below we write the field equation for α(r) with corrections up to order L 12 together with the plots of the energy density of the superconducting tubes in the sector with Baryonic charge B = 1 in Figure 4.
For this we have chosen K = 2, λ = 1, c 6 = c 8 = c 12 = 1 5 , m = 0 and q = p = 1. The field equations are given by that can be written again as a first order equation where in this case Note also that Eq. (40) can be seen as a cubic polynomial in the variable z = α 2 which allows, once again, to reduce the complete field equations to a simple quadrature of the form: where χ(α, E 0 ) is the positive real root of the cubic polynomial in z = α 2 defined in Eq. (40). The integration constantẼ 0 always allows such polynomial to have positive real roots.

Perturbations on the profile
A remark on the stability of the above superconductive tubes is in order. In many situations, when the hedgehog property holds (so that the field equations reduce to a single equation for the profile) the most dangerous perturbations 10 are perturbations of the profile which keep the structure of the  hedgehog ansatz (see [42], [43] and references therein). In the present case these are perturbations of the following type: which do not change the Isospin degrees of freedom associated with the functions Θ and Φ. A direct computation reveals that the linearized version of Eq. (19) around a background solution α 0 (r) of Baryonic charge B = np always has the following zero-mode: ξ (r) = ∇ r α 0 (r). Due to the fact that the integration constant E 0 (defined in Eqs. (20), (21) and (22)) can always be chosen in such a way that ∇ r α 0 (r) never vanishes, the zero mode ξ (r) has no node so that it must be the perturbation with lowest energy. Thus, the present solutions are stable under the above potentially dangerous perturbations. Although this is not a complete proof of stability, it is a very non-trivial test.

Electromagnetic perturbations
A very useful approach to study the stability of the superconducting Baryonic tubes is to perform electromagnetic perturbations on the effective medium defined by the topological solitons. This is a good approach in the 't Hooft limit, since in the semiclassical interaction Photon-Baryon, the Baryon is essentially unaffected due to the Photon has zero mass (see [44]).
The complete stability analysis requires to study the most general perturbations of the solutions defined in Eqs. (17), (21) and (22). This is a very hard task even numerically as it involves a coupled system of linear PDEs, therefore in practical terms, consider only electromagnetic perturbations greatly simplifies the stability analysis and allows to reveal very relevant features of the superconducting tubes, as we will see immediately.
Here we will analyze the simplest non-trivial case which is related to the role of the subleading corrections in the 't Hooft expansion to the Skyrme model of the sixth order. As it has been already emphasized in the previous sections such sixth order term does not even appear in the equation for the profile (while it does enter in the corresponding Maxwell equations). This is very interesting since it shows that, despite the universal character of the present crystals of gauged solitons (which are almost unaffected by the subleading terms), their stability properties may depend explicitly on the subleading terms themselves. Also for simplicity reasons, we will set m to zero.
At first order in the parameter ε the Maxwell equations become where V and Ω(α) are defined in Eqs. (29) and (31)  and ξ 2 must depend on the temporal coordinate. But, according to the previous equations if ξ 1 and ξ 2 depend on time these functions must also depend on the coordinate φ, that is We can assume that By consider the Fourier transformation in the coordinate φ we get an equation for ξ i in the form is the Fourier transform of c i , and the non-vanishing eigenvalue k = l/(2π) is the wave-number along the φ-direction, with l a non-vanishing integer.
According to Duhamel's principle (see, for instance, [45] and references therein), an inhomogeneous equation for a function W = W (x, t) of the form with M a non-negative operator and initial conditions W (·, 0) = ψ 1 , ∂ t W (·, 0) = ψ 2 , has the following general solution In our case, to ensure that the perturbed Maxwell equation can be solved we need to demand that Since V depends on α and the square of its derivative α = α (E 0 ), defined via Eq. (20), one can find the following upper bound to the potential V : Then, a necessary condition 11 for a positive defined effective potential V eff in Eq. (42) is x = L 2 , (the most restrictive case being the one with l 2 = 1 as it is easier to satisfy the above inequality when l 2 is large). The above inequality set an upper bound on the allowed values of L (which is the same as a lower bound on the allowed values of Baryon densities): x The above inequality is equivalent to together with the obvious condition that Thus, in the range of parameters in which L 2 Max is positive (which always exists) the conditions on L 2 is Thus, at a first glance, from Eq. (19) one could think that the presence of the c 6 , which for energetic considerations must be positive (see Eq. (34)), do not play any role in the perturbation of the system. However, it is quite interesting to see that this term can in fact affect the stability of the system determining the maximum allowed value of the size of the box in which the superconducting strings are confined. We hope to come back on the physical properties of these gauged crystals in the low energy limit of QCD in a future publication.

Conclusions and perspectives
The Maxwell-gauged Skyrme model in (3 + 1)-dimensions together with all the subleading corrections in the 't Hooft expansion admit configurations describing ordered arrays of Baryonic superconducting tubes where (most of) the Baryonic charge and total energy is concentrated in tube-shaped regions.
The corresponding current cannot be deformed continuously to zero as it is tied to the topological charge. Quite remarkably, no matter how many subleading terms are included, these ordered arrays of gauged solitons are described by the very same ansatz and keep their main properties manifesting a sort of universal character. The similarity with the plots obtained numerically in the analysis of nuclear spaghetti phase is quite amusing [9]. These results open the unprecedented possibility to analyze these complex structures with analytic tools which are able to disclose novel features which are difficult to analyze with many body simulations. On the other hand, the subleading terms in the

Appendix II: Reducing the Skyrme equations
In this appendix we will show how and why the Skyrme equations (using the ansatz defined by Eqs. (12), (13), (14), and (17)) are reduced to just one integrable ODE for the soliton profile α = α(r). To see this fact it is possible to take two paths, as we detail below.
In order to make this reduction clear, in this appendix we will consider the action in Eq. (1) without the higher order terms and without the coupling with Maxwell theory, i.e. we will deal only with the usual Skyrme action The reason for doing this is that the mechanism which makes the present strategy successful with the usual Skyrme model works in exactly the same way when the higher order terms are included.
The most direct way to see that the Skyrme equations are reduced to just one equation with the ansatz defined by Eqs. (12), (13), (14), and (17) corresponds to write the action in Eq. (49) explicitly in terms of the functions α = α(x µ ), Θ = Θ(x µ ) and Φ = Φ(x µ ), according to Eqs. (13), (14). In this parameterization Eq. (49) becomes Now, varying the action in Eq. (50) w.r.t the functions α, Θ, Φ, in a long but direct calculation, we get to the following set of equations: The equations system written above are completely equivalent to the system in Eq. (5) when the parameterization in Eqs. (13) and (14) is considered. For instance, one can check that with this parametrization the original spherical hedgehog ansatz of Skyrme himself [13] reads where x is the radial coordinate of flat space-time metric in spherical coordinates. If one plugs the ansatz in Eqs. (54) and (55) into the field equations (51), (52) and (53) one can see that, first of all, the field equations for Θ and Φ are identically satisfied and, secondly, that the remaining field equation for α reduces to the well known equation for the profile of the spherical Skyrmion. This is the defining characteristic of the hedgehog ansatz and is equivalent to the statement that the field equations reduce consistently to just one ODE for the profile α.
On the other hand, one could be dissatisfied with the above method to establish the hedgehog property since in most textbook such a property is defined by looking at the SU(2) valued field without using the explicit parametrization in terms of the fields α, Θ and Φ. Here we offer another method to derive the hedgehog property which perhaps is more familiar to many of the readers. This method uses the properties of the normalized Isospin vector n i of the generalized hedgehog ansatz in Eq. (51), (52) and (53), is using.
Let us remind that the most general parametrization for the Skyrme field is n 1 = cos Θ sin Φ , n 2 = sin Θ sin Φ , n 3 = cos Θ , and the Maurer-Cartan form reads where the generators of the SU (2) group, t i , satisfy One can check that, for the ansatz in Eq. (17), the n i vectors satisfy the following eigenvalue equation It is worth to note that the original ansatz for the spherical Skyrmion in Eqs. (54) and (55) satisfies a similar property (but with a different Σ which depends explicitly on the radial coordinate: thus, we can say that in this sense the present generalized hedgehog ansatz is simpler than the usual spherical hedgehog ansatz). Eq. (57) is a very important result since it allows to reduced all the Skyrme system to just only one equation for the soliton profile thanks to a very nice factorization property of the complete field equations (such a factorization is the matrix version of-and completely equivalent to-the property discussed above which is responsible for the fact that the three Skyrme field equations (51), (52) and (53) for α, Θ and Φ reduce to just one equation for α).
One can directly check that the components of the Maurer-Cartan tensor L µ defined above are = n k ∇ µ α + 1 2 sin(2α)∇ µ n k + ijk sin 2 (α)n i ∇ µ n j .
Hence, one can see that the divergence ∇ µ L k µ of the L k µ tensor in Eq. (60) (that corresponds to the NLSM field equations) is factorized into the Isospin vector n k itself (which obviously never vanishes) time a factor which depends on α. Consequently, in the NLSM case, such a factor is nothing but the equation for the profile α. Hence, the choice of α and of the Isospin vector n k in Eq. (51), (52) and (53) reduces the three field equations of the NLSM ∇ µ L k µ = 0 to just α + 1 2 Σ sin(2α) = 0 .
The factorization of the divergence ∇ µ L k µ of the L k µ tensor in Eq. (60) is the matrix form of the property stated here above that the field equations (51), (52) and (53) reduce to just one equation for the soliton profile α: however, this "matrix form" of the hedgehog property can be more familiar to most of the readers so that, for pedagogical reasons we have included it here in the present discussion.
Once again, we see that the hedgehog property is not related at all with the spherical symmetry and nice ansatz can be constructed even at finite density and without spherical symmetry. Even more, the present non-spherical hedgehog, useful to describe multi-solitonic solutions at finite Baryon density is actually simpler than the spherical hedgehog (which describes one Skyrmion since the function Σ in Eq. (60) is constant, as one can see from Eq. (57)).
One may wonder whether this nice factorization property survives when the Skyrme term (and, in fact, also the higher order corrections terms mentioned in the main text) is included. In order to see that this is indeed the case, one can proceed as follows. In fact, the commutator in Eq. (59) can be written as where we have defined S ν µ = L i µ L ν i = ∇ µ α∇ ν α + sin 2 (α)∇ µ n b ∇ ν n b , so that S = L i µ L µ i = ∇ µ α∇ µ α − Σ sin 2 (α) .
This analysis clearly shows that the ansatz defined in Eqs. (12), (13), (14) and (17) reduces the Skyrme equations to a single equation for the soliton profile thanks to the factorization property mentioned here above. It is straightforward to show that the same derivation is still valid when the higher order terms of the generalized Skyrme model are included. To the best of authors knowledge, this is the first complete discussion of the equivalence of these two different viewpoints on the hedgehog ansatz.