Crystals of gauged solitons, force free plasma and resurgence

We show that the (3+1)-dimensional gauged non-linear sigma model minimally coupled to a U(1) gauge field possesses analytic solutions representing gauged solitons at finite Baryon density whose electromagnetic field is a Force Free Plasma. These gauged solitons manifest a crystalline structure and generate in a very natural way persistent currents able to support Force Free Plasma electromagnetic fields. The trajectories of charged test particles moving within these configurations can be characterized. Quite surprisingly, despite the non-integrable nature of the theory, some of the perturbations of these gauged solitons allow to identify a proper resurgent parameter. In particular, the perturbations of the solitons profile are related to the Lam\'e operator. On the other hand, the electromagnetic perturbations on the configurations satisfy a two-dimensional effective Schrodinger equation, where the soliton background interacts with the electromagnetic perturbations through an effective two-dimensional periodic potential. We studied numerically the band energy spectrum for different values of the free parameters of the theory and we found that bands-gaps are modulated by the potential strength. Finally we compare our crystal solutions with those of the (1+1)- dimesional Gross-Neveu model.


Introduction
One of the most important types of plasma in plasma physics is the so-called Force Free Plasma (FFP henceforth), and several reasons justify this statement. First of all, FFP are extremely relevant in many astrophysical situations. For instance, objects with extremely intense magnetic fields like pulsars are typical sources of FFP. Secondly, FFP can be characterized in a very elegant way with a non-linear set of PDEs for the electromagnetic field. This is the so-called force free condition, The above equation is usually the starting point of the theoretical analysis on FFP (see for instance [1], [2] and references therein). From one point of view, the above system of equations is very convenient since many features of FFP can be analyzed without taking into account the electromagnetic sources (which are absent in Eq. (1)). On the other hand, this implies that often the very important physical question about which are the actual concrete sources able to generate such FFP is somehow neglected.
In the present paper, as it will be explained in the following sections, we identify a very natural and concrete source of FFP in the low energy limit of Quantum Chromodynamics (QCD). FFP play also a key role in the so-called Blandford-Znajek process [3], mechanism that allows to extract rotational energy from a spinning black hole. This happens when the electric potential becomes larger than the threshold of electron-positron pair creation, so that the black hole will be surrounded by a FFP in such a way that the rotation energy can be radiated away. The pioneering theoretical analysis in [4], [5], [6], [7], [8], [9], [10], [11], [12] and [13] shed considerable light on this process.
In recent years relevant analytic examples of FFP have been constructed (see, for instance, [14], [15], [16], [17], [18], [19], [20], [21], [22], [23] and references therein), however, many of the available examples (which analyze in detail solutions of the system in Eq. (1)) leave open the issue to identify the actual sources of the force-free electromagnetic field. Due to the interesting astrophysical applications of FFP it is of great importance to find also analytic examples of FFP in which the sources of the electromagnetic field can be traced back directly to observed particles of the standard model. One of the goals of the present paper is to construct explicitly such solutions. In the present approach, the SU (2) non-linear sigma model (NLSM henceforth), which is the low limit of QCD describing Pions dynamic, provides with explicit and topologically non-trivial sources of FFP. This is in agreement with the fact that FFP must occur for objects made of Hadronic matter (such as pulsars).
The NLSM is one of the most interesting effective field theories. Besides the well-known relation with the low energy limit of QCD, such a model is also very useful in statistical mechanics, in the theory of the quantum hall effect, in the analysis of superfluids and so on (see [24], [25], [26] and [27]). It is worth to emphasize that the NLSM on flat space-time does not possess static solutions with non-trivial topological charge due to Derrick's scaling argument [28]. A well-known method to avoid this no-go theorem was found by Skyrme (even before the actual paper of Derrick): the so-called Skyrme term [29] together with the NLSM allows the existence of static solitonic solutions with finite energy and topological charge, which represents the Baryonic charge, called Skyrmions (see [30], [31], [32], [33] and references therein). However, the original arguments in [30] to prove that the topological charge must be identified with the Baryonic charge do not require explicitly the presence of the Skyrme term. Moreover, the Skyrme term is not the only way to avoid the Derrick's argument. The Derrick's scaling argument can be avoided by constructing time-dependent ansatz for the matter fields with the characteristic that the corresponding energy-density is stationary. Hence, as we will show in the next sections, the minimal ingredients to build FFP together with the corresponding sources are already present in the NLSM minimally coupled with Maxwell field.
The present strategy is based on the techniques developed in [34], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], that have allowed to construct exact NLSM solutions with non-trivial topological charge at finite density in [50], [51], [52], [53], [54] and [55]. The gauged solitons constructed in [52] and [53] live at finite Baryon density and most of the energy and the topological charge are contained within tube-shaped regions which are regularly spaced 1 . These structures are expected in the description of cold and dense nuclear matter as a function of topological charge, commonly called nuclear pasta, and therefore are quite relevant in the phase diagram of the low energy limit of QCD (see [58] and [59]). Although until very recently the nuclear pasta phase was considered to be a very hard nut to crack from the analytic viewpoint, there are very strong observational evidences supporting it (see [60], [61], [62], [63] and references therein). The comparison of the energy density and Baryon density contour plots in [51], [52] and [55] with the ones in, for instance, [62] is very encouraging. Not only the analytic plots of energy density and Baryon density in [51], [52] and [55] are very close to the phenomenological ones in [60], [61], [62], [63], the present analytic framework also allows the explicit computation of relevant quantities (such as the computation of the shear modulus of the nuclear lasagna in [54] which is close to the value of the shear modulus of nuclear lasagna found in [64]). Here we will show that the electromagnetic field generated by the solitons crystals found in [51], [52] (which describe quite well many features of the nuclear spaghetti phase) is of force free type.
This connection between FFP and the nuclear pasta phase is a quite surprising result which (to the best of our knowledge) has not been discussed previously in the literature. As it will be discussed in the next sections, this allows among other things to discuss the trajectories of charged particles close to these gauged solitons.
Another novel result of the present analysis is that some of the main physical properties of these gauged solitons, which are natural sources of FFP, manifest a "resurgent character ". Resurgence is one of the few theoretical techniques which can help to make sense of the ubiquitous factorially divergent series appearing in perturbation theory [65]. Such strategy works as follows: First of all, one must use Borel summation in the complex-g plane (g being the coupling constant of the theory) to handle the given factorially divergent series so that the divergent series becomes a finite expression. At a first glance, it seems that we have just traded one problem for another (maybe worse) problem since the Borel sum possesses ambiguities which manifest themselves along suitable lines in the complex g-plane (for nice reviews see [66] and [67]). However, when the models under investigation possess suitable topological charges labelling non-trivial non-perturbative sectors, the perturbative expansions in these topologically non-trivial sectors (which, usually, are also factorially divergent) allow a remarkable rescue. In these relevant theories, the ambiguities in the Borel summation of the perturbative expansions in the topologically non-trivial sectors cancel the ambiguities of the perturbative sector giving rise to a well-defined result [68], [69]. Starting from [70] in recent years there has been a great revival of these ideas (see, for instance, [66], [67], [71], [72], [73] and references therein) with applications in Quantum Mechanics, topological strings, and integrable Quantum Field Theory in low dimensions (see [71], [74], [75], [76], [77], [78], [79] and references therein), however there are relatively few applications of resurgence in non-integrable models. It is therefore surprising that in the analysis of FFP interacting with Hadronic matter (which is a sort of prototype of non-integrable systems), operators which can be analyzed with such a tool (such as the Lamé operator) appear. In fact, one dimensional periodic quantum systems like Mathieu and Lamé equations have well known resurgence structures, and in the usual cases which one is the suitable "resurgent parameter" is pretty clear from the beginning [80], [81], [82]. On the other hand, in the present situation in which gauged solitons at finite Baryon density generate a FFP, there are many potential expansion parameters and, a priori, it is not obvious which one is the "most suitable" from the resurgent viewpoint. We will see that our analysis shed light on this issue as well.
This paper is organized as follows: In Section 2 the gauged NLSM is introduced. In Section 3 we show that the gauged NLSM can be solved analytically to obtain topological multi-solitons solutions.
In Section 4 we show that the electromagnetic field generated by the multi-solitons is a FFP, and we study its main physical features, including the trajectory of charged particles in this background. In Section 5, in order to test the stability of the solutions, we perform perturbations on the solitons as well as on the electromagnetic field. In Section 6 we draw a relation between the perturbations on the solutions and the theory of resurgence. Section 7 is devoted to the comparison of our solutions with the (1+1)-dimensional crystals of the Gross-Neveu model. In the final section some conclusions will be drawn.

The gauged non-linear sigma model
The action of the gauged NLSM minimally coupled with the Maxwell theory is where g is the determinant of the metric, m is the Pions mass and K is the coupling constant of the NLSM. Here A µ is the gauge potential, ∇ µ denotes the partial derivative and σ i are the Pauli matrices. The covariant derivative is defined by The field equations, obtained varying the action in Eq. (2) w.r.t the fields U and A µ , are where the current, J µ , in Eq. (6) is given by The energy-momentum tensor of the model is withT µν the electromagnetic energy-momentum tensor, On the other hand the topological charge is where the topological charge density ρ B is defined by The above expression was constructed in [91] (see also the pedagogical analysis in [92]), where the second term in Eq. (10) is added in order to guarantee both the conservation and the gauge invariance of the topological charge. When Σ is space-like, B turns out to be the Baryonic charge of the configuration.

Gauged solitons at finite Baryon density
In this section we show that the gauged SU (2)-NLSM admits analytic solutions describing crystals of topological solitons.

The ansatz
As it has been already emphasized, the analytic description of cold and dense nuclear matter as a function of the topological charge is a very interesting but quite challenging theoretical problem (see [58], [59] and references therein). Since one of the main goals of the present analysis is precisely to understand how gauged solitons react to the presence of a finite volume containing them, we need to analyze the model in the following flat metric where the adimensional coordinates {r, θ, φ} have the ranges The above metric represents a box with different lengths on each side, so that V = 4π 3 L r L θ L, is the volume of the box in which the solitons are confined.
For the U field we adopt the standard parametrization of an element of SU (2), that is where 1 2 is the 2 × 2 identity matrix and According to Eqs. (9), (10), (13) and (14) it follows that the topological charge density is Here, from Eq. (15) one can see that to have topologically non-trivial configurations we must demand that dα ∧ dΘ ∧ dΦ = 0 . On the other hand, as we want to construct analytical solutions, it is necessary to have a good ansatz that significantly reduces the field equations in Eqs. (5) and (6). The approach developed in [51], [52] and [55] lead to the following ansatz for the gauged solitons as well as for the electromagnetic potential It is a direct computation to verify that the ansatz defined in Eqs. (16) and (17) possesses the "hedgehog property". Such a property, which in the past was usually only associated to spherically symmetric solitons, corresponds to the following desirable feature: one would like to reduce the coupled system of non-linear field equations defining the solitons to only one non-linear equation for the profile keeping alive the topological charge. In the present case, quite remarkably, this property holds despite the lack of spherical symmetry and despite the minimal coupling with the Maxwell gauge field, as we will see below.

Analytic crystal of topological solitons
Using the ansatz in Eqs. (11), (13), (16) and (17) one can check directly that, first of all, the three coupled field equations of the gauged NLSM reduce to a single ODE for the profile α, namely which can be reduced to a first order equation of the form with E 0 an integration constant. The above equation can be solved explicitly in terms of Inverse Elliptic Functions and is reduced to the following quadrature A necessary condition for stability is that α does not change sign, therefore we must require, according to Eq. (20), that On the other hand, with the ansatz defined in Eqs. (11), (13), (16) and (17) the Maxwell equations are reduced to just one linear Schrödinger-like equation with an effective two-dimensional potential, that is where we have defined Note that the Maxwell equation in Eq. (22) can be solved once the soliton profile is determined, whose solution is implicit in Eq. (20).

Topological charge and energy density
From Eqs. (9), (13), (15), (16) and (17) one can see that the natural boundary conditions in the present case (see [52] for a detailed analysis) are the following The above boundary conditions in Eq. (24) determine that the topological charge turns out to be Now, from Eqs. (12), (20) and (25) we see that the integration constant E 0 is fixed in terms of n through the relation The above equation always has a solution according to Eq. (21), which can be found numerically.
On the other hand, the energy density of the configurations presented above is given by It is also convenient to define the "reduced energy density" ε r in the r-direction obtained integrating Eq. (27) along the coordinates θ and φ, that is This reduced energy density represents the effective energy density (namely, the energy per unit of length) in the r-direction, and can be compared directly with 1-dimensional models possessing solitons crystals. In the particular case in which both the Maxwell gauge potential and the Pions mass vanish one gets 4 On the force free nature of the electromagnetic field In this section we will show that the electromagnetic field generated by the gauged solitons discussed in the previous sections satisfies the FFP condition. To the best of the authors knowledge, this is the first analytic example in which a FFP emerge in a natural way.

Emergent force free plasma
The force free condition is realized in many relevant examples of plasmas such as in the solar corona [1] or close to rotating neutron stars and black holes [93], [94], [95]. As it was already mentioned the FFP plays an important role in the Blandford-Znajek process [3] but, even in such process, the FFP emerge as a suitable choice that "provide a reasonable approximation to the time-averaged structure of the magnetosphere", and then, the black hole parameters are estimate values in order to satisfy the force-free condition.
It can be readily seen that the Maxwell field surrounding our gauged solitons satisfies the force-free condition in Eq. (30). First, note that the U (1) gauge potential in Eq. (17) can be written as so that the field strength is found to be It follows that the divergence of the field strength, identified with the current, is which is manifestly orthogonal to the field strength F µν . Therefore we see that the force-free condition in Eq. (30) is automatically satisfied.
The electromagnetic field lines (for suitable choice of the parameters) are presented in Fig. 2, where we have used the following boundary conditions for the Maxwell potential u(r, 0) = u(r, π) = 0 , u(0, r) = u(2π, 0) = 0 .
At this point it is important to emphasize that, although several analytic solutions to the force-free  Figure 1: Electric field, magnetic field and the current of two gauged solitons, with n = 2, m = 0, and p = q = 1. The electric and magnetic fields vanish in the center of the tubes while the current is completely contained inside these. Maxwell equations could be found (for instance, using the approaches in [14], [16], [17], [18], [20], [21] and [22]) most of these analytic solutions to the force-free electrodynamics do not discuss the corresponding sources. Instead, the source of the FFP electromagnetic field found here can be clearly identified with the Hadronic degrees of freedom described by the NLSM.

Trajectory of charged particles on the FFP
Here we will draw some qualitative plots of the trajectories of charged test particles (like electrons) on the FFP generated by the gauged solitons.
First of all, it follows from Eq. (17) that the electric and magnetic fields generated by the gauged solitons are where the components read For a test particle of charge q e and velocity v, the Lorentz force generated by the FFP acting on the particle is with Eq. (34) is a set of three coupled differential equations, namely The above equations can be numerically integrated to know the trajectory of test particles. We show the trajectory of a single particle for times 0 < t < 1, for 1 < t < 50 and for 490 < t < 500 in Fig. 3.
In these figures one can notice that the charged test particle oscillates in the r − θ plane while moving along the axis of the tube. It is also worth to emphasize that this motion is quite different from the usual trajectory of a charged test particle in the magnetic field of a thin straight wire. Indeed, in the case of a thin solenoid there is only the θ−component of the magnetic field (and no electric field) while in the present case there is also a radial component and both the electric and the magnetic fields are non-vanishing.

Perturbative analysis
A complete perturbative analysis of these solutions and the corresponding FFP is beyond the scope of the present paper (since such an analysis would involve seven coupled PDE in the non-trivial background provided by the soliton crystal itself). This task is very hard even from a numerical viewpoint. Hence, in this section, we will analyze two special types of perturbations of the Hadronic field and of the FFP which have both interesting physical meaning and allows some analytic control.

SU (2) perturbations
A relevant perturbation of these gauged solitons is the one which keeps the hedgehog structure intact, namely, which keep alive the hedgehog property defined in the previous section. In the present context such perturbations are defined by the following small deformations of the original ansatz in Eqs. (13) and (16): with 0 < ε 1. Such perturbations (and, in particular, the little disturbance F 1 of the profile α(r)) are very likely to be the "smallest energy perturbations" 2 of the Hadronic degrees of freedom. For this reason the linear operator which determines the spectrum of these perturbations has a very important role.
The field equations in Eq. (5) at first order in ε for the perturbations defined in Eq. (38) read where we have defined First of all, one can notice that F 2 (r) can be readily integrated giving rise to On the other hand the situation for F 1 (r) is different. We will write down such linear operator in the case m = 0 since in this way the analysis is cleaner, however a non-vanishing m will not change the results and will only make the analytic formulas a bit more involved. In the m = 0 case the soliton profile (see Eqs. (18) and (19)) is given by 2 The reason is that it takes "a little effort" to perform a radial deformation of α (see, for instance, the discussion of the "hedgehog ansatz" in [96] and [97]) as compared, for instance, to deformations of the Isospin degrees of freedom Θ and Φ. For instance, the classic reference [33] showed that the low energy "Isospin perturbations" have a gap. Hence, there are examples in the literatures in which the unstable negative energy-modes are precisely in this sector while the perturbations of the Isospin degrees of freedom have positive energies.
where am(x; k) is the the Jacobi Amplitude [98] and k is the elliptic parameter (related to the elliptic modulus m in [98] as k ≡ √ m). From here on we are going to use the parameter k instead m. The periodic structure appears due to the well-known properties of the function am(r; k), namely sin (am(r; k)) = sn(r; k) , cos (am(r; k)) = cn(r; k) , where sn(r; k), cn(r; k) and dn(r; k) are the Jacobi Elliptic Functions [98].The function sn 2 (r; k) is periodic, with period 2K(k), where K(k) = π/2 0 dθ/ 1 − k sin 2 θ is the elliptic-quarter period.
Going back to the linearized equations, the equation for F 1 (r) in Eq. (39) in the case m = 0 takes the form Quite interestingly, using the properties of Jacobi Amplitudes in Eqs. (43), (44) and (45), the operator ℘ (0) (which determines the stability of the present gauged solitons system under the perturbations in Eq. (38)) can be cast in the form of a Lamé Operator: Consequently, the spectrum of the family of operators ℘ (0) (or ℘ (m) ) defined above is very important not only in relation with the stability analysis of the present gauged solitons but also because such spectrum encodes relevant physical informations about the band spectrum of these gauged solitons.
As far as the stability of this gauged solitons-FFP system one need to study the following eigenvalues problem so that the stability condition under this family of perturbations is It can be seen that the condition in Eq. (49) can be satisfied. Indeed, it is easy to find a zero mode just taking where α 0 is a background solution of Eqs. (41) and (42)  Moreover we can also deduce that Ψ 0 has no node since ∂ r α 0 does not vanish. Consequently, standard theorems in Quantum Mechanics ensure that Ω 2 = 0 is the minimal eigenvalue and all the other are positive. In fact, Eq. (48) encodes many more information since it has to do with the spectrum of the lowest energy perturbations of the present gauged solitons. In other words, the spectrum of ℘ (0) (or ℘ (m) ) will tell us the "phonons" of the system. In the next section we will analyze such spectrum in the case m = 0 in which the results are cleaner and directly related with the theory of resurgence. The results with m = 0 are very similar but more difficult to analyze with resurgence since such spectrum involves two relevant parameters (this case belongs to the so-called parametric resurgence).
Last but not least, this family of operators is especially well-suited for the resurgence analysis as it has been shown in [77], [80] and [81]. One of the main results of this section is that the proper resurgent parameter g (in the m = 0 case) is Although this result is very simple (one just needs to compare Eq. (47) with those in [77], [80] and [81]) it is very non-trivial. It is relevant to note that the effective resurgent parameter 3 g does not involve the coupling constant of the theory, namely K. Instead, the direct computation discussed in this section shows that the suitable parameter is in Eq. (50), which depends on the odd integer q as well as on the "asymmetry ratio" It is worth to emphasize the similarity of the spectrum of the "phonons" defined by the operator in Eq. (48) with the ones of the crystal of kinks in [99] as well as in the analysis of the small fluctuations of solitons crystals in integrable field theory models in (1+1) dimensions. In Section 7 we will elaborate more on this comparison.

Electromagnetic perturbations
An effective technique to analyze the stability of solutions presented here is to consider the gauged solitons as an effective background medium on which the electromagnetic perturbations propagate. This is an excellent approximation especially in the 't Hooft limit because in the semiclassical Photon-Baryon scattering, the Baryon is basically not affected since the Photon has zero mass (for a detailed discussion see [100]). Consequently, the Photons perceive the gauged soliton as an effective medium while the solitons are not affected by the small fluctuations of the electromagnetic field. Correspondingly, in this section we will consider the SU (2) degrees of freedom as fixed to be the background solution.
As it has been already emphasized, the complete stability analysis not only is completely out of reach from analytical methods but also numerically is a very hard task. Nevertheless, the perturbations we are considering here encode relevant features of the gauged solitons and of the corresponding FFP, as we will see below.
Let us consider the following perturbations on the Maxwell potential with 0 < ε 1. At first order in the parameter ε and for the perturbation defined in Eq. (51) the Maxwell equations in Eq. (6) become Since we want to test linear stability we need to check the (absence of) growing modes in time. This implies that ξ 1 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 For simplicity we will assume that so that the above equations system is reduced to with V defined in Eq. (23). Then, performing the Fourier transformation in the coordinate φ and time t, we obtain an eigenvalue equation forξ i , with The non-vanishing eigenvalue in the effective potential in Eq. (55) is the wave-number along the φ-direction, with l a non-vanishing integer. A sufficient condition ensuring linear stability under the perturbation defined in Eqs. (52) and (53) is the requirement The obvious reason is that the above condition implies that the eigenvalues of the operator O S defined as are positive and, consequently, the parameter ω in the time-Fourier transform of the perturbations is real (so that there are no growing modes in time).
However, this condition is not necessary in the sense that it could be possible for V eff to be "slightly negative" keeping, at the same time, the eigenvalue ω 2 in Eq. (54) real and positive. On the other hand, the mathematical task to find a sharp characterization of the stability in this sector can be very complicated from the viewpoint of functional analysis in the case of Schrodinger-like potentials which depend in a non-trivial way on two (or more) spatial coordinates. Hence, here we will only consider the criterion in Eq. (57) which is more than enough to provide a qualitative picture. We hope to come back on this interesting mathematical issue in a future publication.
The requirement in Eq. (57) imposes an upper bound for the length of the box. This can be seen as follows. The "less favorable case" (from the stability viewpoint) corresponds to the least possible value for (k 3 ) 2 in Eq. (56) when the positive part of the effective potential is the smallest. Correspondingly, In terms of the Baryon density, the above bound implies that the present configurations are viable when the Baryon density is of the order of 1 Baryon per f m 3 or higher (this range is well within the range of validity of the NLSM). In fact, we expect that a more refined mathematical analysis would give a better bound showing that these solutions are viable even at lower densities.
6 Some comments on the fluctuations As it has been described in the previous sections, the perturbations for the SU (2) field are governed by operators with a well-established resurgence character such as the Lamé operator. The singularities analysis (using, for instance, the Borel-Pade approximation) of these operators in the Borel plane of the effective coupling constant is crucial in order to understand the resurgent character of these type of perturbations. This will be done in this section.
On the other hand, although resurgence techniques in one-dimensional quantum mechanical systems are very well tested, in the case of two-dimensional (and, indeed, higher dimensional) potentials the situation is far less clear. Since the electromagnetic perturbations satisfy an effective twodimensional Schrödinger equation, the spectrum of the electromagnetic perturbations must be determined numerically. In this section we will obtain numerical results relating the electromagnetic perturbations with the following two-dimensional Schrödinger equation 6.1 WKB Analysis for the SU (2) perturbations According to the previous section about the perturbations on the U field, we need to study a Schrödinger equation of the form where the potential is given by We show the energy spectrum of our configurations as a function of the coupling constant g in Fig. 4 (see Eq. (50)). For this purpose we employ the WKB method to obtain an expansion in g 2 → 0, This is an standard procedure in which one solve the resulting Ricatti equation using a power series ansatz for S(x) as well as for the energy, that is Here we have used the BenderWu package [101] to compute the WKB expansion obtaining a perturbative asymptotic expansion of the ground state energy (we will not consider higher level states in this paper).

Numerical approach for the perturbed Maxwell equation
Here we present numerical solutions for the two-dimensional eigenvalue problem in Eq. (58) using the Finite Difference Method, considering a two-dimensional grid in the (r, θ)-plane of N r × N θ grid of 70 × 70 on a domain delimited by the ranges 0 < r < 2π and 0 < θ < π. We have used the Python Library Linear algebra (scipy.linalg) for eigenvalue problems. For the numerical implementation we will consider for simplicity L r = L θ = L. In this case the gauged NLSM field equations and the perturbed Maxwell equations read    The boundary conditions for the profile α(r) according to Eqs. (24), (25) are α are α(2π) − α(0) = nπ, so that the parameter n in the boundary conditions is equal to number of "peaks" of the potential (or energy density) in the r-direction of the lattice, as it is shown in Figure 6). The potential is fully determined by fixing the parameters q, K and m. In Figure 7 we plot the energy eigenvalues as function of L. Making the comparison with the one dimensional Mathieu system [77] (where the band-gaps are modulated by ) one can note that the parameter L plays the role of an effective .
This theory exhibits a phase diagram with a crystalline structure similar to what is expected in QCD, which is invisible to standard perturbation theory being only accessible to the non-perturbative 1/N expansion [88], [89], [90]. Using an auxiliary field, the GN field equations can be written as a Hartree-Fock-Dirac equation (see, for instance, [83] and [84]), that is The self-consistent analysis of the previous references allows to assume that S(x) is a Lamé potential where E = a 2 2 ω 2 + κ 2 .
The eigenvalue term E depends on the constants a and , which are fixed by the mean density. The elliptic parameter κ is fixed by minimizing the ground state energy density It is relevant to note the close similarity of the results in the references [84], [85], [86], [87], [88], [89], [90] and [99] in the case of the solitons crystals in (1+1)-dimensions with the (3+1)-dimensional gauged solitons discussed here. Not only the energy density of the solitons crystal of the GN model and of the Sine-Gordon model in [99] are very similar to the energy density of the present gauged solitons. This similarity is especially clear comparing the form of the reduced energy density defined in Eqs. (28) and (29) with, for instance, the corresponding expression in [99]. Also, the spectrum of the fluctuations of these (1+1)-dimensional solitons crystals in [84], [85], [86], [87], [88], [89], [90] and [99] are determined by a one-dimensional Schrödinger operator of the same Lamé family as the one in Eq. (47). This supports very strongly the existence of non-homogeneous condensates in the low energy limit of QCD in (3+1)-dimensions.

Conclusions and perspectives
In this paper we have shown analytically that topologically non-trivial gauged solitons of the (3+1)dimensional gauged non-linear sigma model at finite Baryon density are natural sources of Force Free Plasma. For these multi-solitons most of the total energy and the topological charge are gathered within tube-shaped regions while the magnetic field lines go around the tubes. Our explicit analytical solutions allow to discuss the trajectories of charged test particles moving close to these gauged solitons and also identify the proper resurgent parameters through the analysis of the perturbations.
In particular, the perturbations of the solitons profile are related to the Lamé operator with a suitable resurgent parameter. On the other hand, the electromagnetic perturbations of the above system satisfy a two-dimensional effective Schrödinger equation, where the soliton's background interacts with the electromagnetic perturbations through an effective periodic potential in two spatial dimensions. Also we have studied numerically the band energy spectrum for different values of the free parameters of the theory and we found that bands-gaps are modulated by the potential strength. Finally, we have shown that the crystal solutions constructed here are qualitatively very similar to those of the Gross-Neveu model in (1+1) dimensions, which strongly support the existence of non-homogeneous condensates in the low energy limit of QCD in (3+1)-dimensions.