Stability of Classical Chromodynamic Fields

A system of gluon fields generated at the earliest phase of relativistic heavy-ion collisions can be described in terms of classical fields. Numerical simulations show that the system is unstable but a character of the instability is not well understood. With the intention to systematically study the problem, we analyze a stability of classical chromomagnetic and chromoelectric fields which are constant and uniform. We consider the Abelian configurations discussed in the past where the fields are due to the single-color potentials linearly depending on coordinates. However, we mostly focus on the nonAbelian configurations where the fields are generated by the multi-color non-commuting constant uniform potentials. We derive a complete spectrum of small fluctuations around the background fields which obey the linearized Yang-Mills equations. The spectra of Abelian and nonAbelian configurations are similar but different and they both include unstable modes. We briefly discuss the relevance of our results for fields which are uniform only in a limited spatial domain.


I. INTRODUCTION
Soon after the discovery of asymptotic freedom [1,2], when quantum chromodynamics was recognized as an underlying theory of strong interactions, the stability of various configurations of classical chromodynamic fields was investigated [3][4][5][6][7][8]. These studies, which revealed that numerous configurations are actually unstable, were not performed with a specific application in mind, rather it was about better understanding the newborn theory.
Today, classical chromodynamics is often used as an approximation of quantum theory. We are interested in the early phase of relativistic heavy-ion collisions experimentally studied at RHIC and the LHC. Within the Color Glass Condensate (CGC) approach, see e.g. the review articles [9,10], color charges of valence quarks of the colliding nuclei act as sources of long wavelength chromodynamic fields which can be treated as classical because of their large occupation numbers. The system of non-equilibrium gluon fields created in the nuclear collision is called glasma. At the earliest moment of the collision, the glasma is dominated by the chormoelectric and chromomagnetic fields parallel to the beam axis and later on transverse fields show up.
It has been found numerically [11,12], see also [13], that there is an unstable exponentially growing mode in the course of glasma's evolution. The mode was identified as the Weibel instability [14] which is well known in physics of electromagnetic plasma. A presence of the chromodynamic Weibel instability in relativistic heavy-ion collisions was first argued in [15] and further on studied in detail, see the review [16].
The Weibel instability occurs when charged particles with anisotropic momentum distribution interact with the magnetic field generated by the particles. As explained in detail in [16], there is an energy transfer from the particles to the field which causes its exponential growth. In glasma there are no particles but high-frequency modes of classical fields are often treated as particles [11,12].
The problem of unstable glasma was studied in the series of papers [17][18][19][20] where a particular attention was paid to strong longitudinal chormoelectric and chromomagnetic fields generated at the earliest phase of nuclear collisions. It was suggested [17][18][19][20] that the unstable mode found in the numerical simulation [11,12] is not the Weibel but rather Nielsen-Olesen instability [21] which occurs when spin 1 charged particles circulate in a uniform magnetic field. There was considered [19] a possible role of the vacuum instability due to strong electric field which according to the Schwinger mechanism [22] causes a spontaneous generation of particle-antiparticle pairs from vacuum. A stability of oscillatory chromomagnetic fields was also studied [23] in the context of glasma and a coexistence of the Nielsen-Olesen instability with the phenomenon of parametric resonance was found.
We intend to clarify what are the unstable modes of evolving glasma. Since the simulation [11,12] was preformed in terms of classical fields we study a stability of classical field configurations. We start with the simplest case of constant and uniform chromomagnetic and chromoelectric fields. The fields which are truly constant and uniform are obviously an idealization but our results are relevant for fields which are approximately constant and uniform in a limited space-time domain. Here we focus on a specific aspect of nonAbelian theory which has not been explored yet. The constant and homogeneous chromoelectric and chromomagnetic fields can occur due to the potentials which are of single color and as in electrodynamics linearly depend on coordinates. We call such configurations Abelian. However, the fields can be also generated by the multi-color non-commuting potentials and then we have genuinely nonAbelian configurations. We note that the Abelian and nonAbelian configurations are physically inequivalent as they cannot be related to each other by a gauge transformation. It was also proved [24] that there are only these two gauge-inequivalent configurations which produce the space-time uniform chromodynamic fields.
Stability of the Abelian configurations of constant and uniform chromoelectric and chromomagnetic fields was studied in [5,6] and later on repeatedly analyzed, see e.g. [18][19][20]23]. However, the nonAbelian configurations seem to be more relevant for glasma. The point is that the chromoelectric and chromomagnetic fields E a , B a from the earliest phase of the collisions are generated along the beam axis z in a nonAbelian manner [25,26]. Specifically, the E a , B a fields occur due to the transverse pure gauge potentials of initial nuclei where f abc are the structure constants of the SU(N c ) group and zij is the totally asymmetric tensor. We are aware of only one paper [7] where the stability of nonAbelian uniform configuration of chromomagnetic field was briefly discussed. A presence of an unstable mode was indicated but a complete spectrum of modes was not derived.
We perform a comparative study of linear stability of Abelian and nonAbelian configurations of constant and homogeneous chromomagnetic and chromoelectric fields. We are mostly interested in the nonAbelian configurations but for a completeness of our study we repeat with minor refinements the stability analyses of Abelian configurations presented in [5,6]. Throughout our whole study we use the background gauge while the axial gauges (different for the chromomagnetic and chromoelectric configurations) were applied in [5,6]. Using one gauge facilitates comparisons of various cases.
We note that comparative analyses of the Abelian and nonAbelian configurations of uniform chromoelectric and chromomagnetic fields can be found in [24] and [27]. A motion of a classical particle was shown to be significantly different in the two cases [24]. Quantum matter fields of spin 0 and 1/2 also behave differently in the background of Abelian and nonAbelian chromodynamic fields [24,27]. However, the self-interaction of nonAbelian fields, which is of our main interest, was not studied in [24] and [27].
Our paper is organized as follows. In Sec. II we present the linearized Yang-Mills equations in the background gauge which are subsequently used in stability analyses. In Secs. III and IV we discuss, respectively, Abelian and nonAbelian configurations of the constant homogeneous chromomagnetic field. Secs. V and VI are devoted analogously to the constant homogeneous chromoelectric field. Our study is closed in Sec. VII. After summarizing our considerations, we briefly discuss the relevance of our results for fields which are uniform only in a limited spatial domain. Finally, we outline a perspective for further research.
Throughout the paper the indices i, j = x, y, z and µ, ν = 0, 1, 2, 3 label, respectively, the Cartesian spatial coordinates and those of Minkowski space. The signature of the metric tensor is (+, −, −, −). The indices a, b = 1, 2, . . . N 2 c −1 numerate color components in the adjoint representation of SU(N c ) gauge group. We neglect henceforth the prefix 'chromo' when referring to chromoelectric or chromomagnetic fields. Since we study chromodynamics only, this should not be confusing.

II. LINEARIZED CLASSICAL CHROMODYNAMICS
The Yang-Mills equations written in the adjoint representation of the SU(N c ) gauge group are where D ab µ ≡ ∂ µ δ ab − gf abc A c µ , j ν a is the color current and the strength tensor is The electric and magnetic fields are given as where ijk is the Levi-Civita fully antisymmetric tensor. We assume that the potentialĀ µ a solves the Yang-Mills equation (1) and we consider small fluctuations a µ a around A µ a . So, we define the potential A µ a (t, r) ≡Ā µ a (t, r) + a µ a (t, r), such that |Ā(t, r)| |a(t, r)|. Assuming that the background potentialĀ µ a satisfies the Lorentz gauge condition ∂ µĀ µ a = 0 while the fluctuation potential a µ a that of the background gaugeD ab µ a µ b = 0, whereD ab µ ≡ ∂ µ δ ab − gf abcĀc µ , the Yang-Mills equation linearized in a µ a can be written as The background gauge appears particularly convenient for our purposes because different color and space-time components of a a µ are mixed only through the tensorF µν b which enters Eq. (6). In case of other gauges, e.g. the Lorentz gauge ∂ µ a µ a = 0, the mixing is more severe. Throughout our analysis, which includes Abelian and nonAbelian configurations of magnetic and electric fields, we use the background gauge which facilitates comparisons of various cases. However, our further considerations are limited to the SU(2) gauge group when f abc = abc with a, b = 1, 2, 3.

III. ABELIAN CONFIGURATION OF MAGNETIC FIELD
The constant homogeneous magnetic field along the axis x occurs for a potentialĀ i a which is known from the Abelian theory. Specifically,Ā a (t, r) = δ a1 (0, 0, yB), where B is a constant and r = (x, y, z). Then, using Eqs. (3), one finds E a (t, r) = 0 and B a (t, r) = δ a1 (B, 0, 0). We also note that the only non-vanishing components of the strength tensor areF zy The Abelian configuration solves the Yang-Mills equations (1) with vanishing current j µ a . The nonAbelian terms disappear because there is only one color component. We also note that the chosen potential satisfies the Lorentz gauge condition.
WhenĀ µ a (t, r) = δ a1 (0, 0, 0, yB), the equation (6) One sees that the color component a µ 1 decouples from the remaining two color components and it satisfies the free equation of motion. So, the functions a µ 1 represent free waves which we do not consider anymore. Defining the functions the equation (7) provides The equations of T ± and X ± have the diagonal form. To diagonalize the equations of Y ± and Z ± one defines the functions which allow one to change the equations (11) -(12) into − 2igBy∂ z ∓ 2gB + g 2 B 2 y 2 W ± = 0.
We assume that the functions a µ a depend on t, x, z as e −i(ωt−kxx−kzz) . Since the functions should be real, only their real parts are of physical meaning. Now, the equations (9), (10) and (14), (15) − We note that one obtains exactly the same equations (17), (18) and (19) using the temporal axial gauge a 0 a = 0 which was applied in Refs. [5,6].
Since the eigenenergy Schrödinger equation of harmonic oscillator can be written as where m is the oscillator mass, E its energy andω is the frequency of the corresponding classical oscillator, one observes that Eqs. (16) -(19) coincide with Eq. (20) under the following replacements where d = 0 for Eqs. (16), (17) and d = ∓2gB for Eqs. (18), (19).
Since E =ω(n + 1/2) with n = 0, 1, 2, . . . , the frequency squared ω 2 is for Eqs. (16), (17) and for Eqs. (18), (19). It should be stressed that although we refer to the Schrödinger equation to find the frequencies (22), (23) the solutions are purely classical -the Planck constant does not show up in the final formulas. The frequencies squared (22), (23) are 'quantized' because the fluctuation field a µ a is assumed to be limited everywhere. This is analogous to the requirement that a wave function, which solves the Schrödinger equation, is normalizable.
One sees that ω 2 0 ≥ 0 and ω 2 + ≥ 0 for any n but ω 2 − = −gB + k 2 x for n = 0 and consequently, it is negative for k 2 x < gB. Then, there are unstable modes of U − and W + which grow as e γt with γ ≡ gB − k 2 x . This is the wellknown Nielsen-Olesen instability [21]. The unstable modes are paired with the overdamped modes which decrease in time as e −γt .
Let us now discuss a character of the solutions of Eqs. (16) - (19). The potentials a µ a are assumed to depend on time as e iωt but to see how a given combination of a µ a evolves in time, one should consider only the real parts of a µ a .
Modes T ± and X ± The modes T ± and X ± are stable. Assuming that T + = 0 while T − = X ± = U ± = W ± = 0, one finds that T + represents the wave which rotates in the two-dimensional color space spanned by the colors 2 and 3. The mode T − is similar but it rotates in the opposite direction than T + . The modes X ± behave as T ± .
Modes U ± and W ± The modes U + and W − are always stable. When U + = 0 and U − = X ± = W ± = 0, one finds that U + represents the wave which rotates in both two-dimensional 2-3 color and y-z coordinate spaces. There is analogous situation with the stable modes W − . The modes U − and W + can be stable or unstable. If the modes are stable, they are similar to U + and W − . In case of unstable and overdamped modes U − and W + , which depend on time as e γt and e −γt , the small field wave does not rotate neither in color nor in coordinate space.

IV. NONABELIAN CONFIGURATION OF MAGNETIC FIELD
A nonAbelian configuration ofĀ i a which produces a constant homogeneous magnetic field B a = δ a1 (B, 0, 0) can be chosen asĀ where the Lorentz index µ numerates the columns and the color index a numerates the rows. The potential (24), which obviously satisfies the Lorentz gauge condition, does not solve the Yang-Mills equation (1) with j µ a = 0. Instead one gets  Following [7], we assume that the current, which enters the Yang-Mills equation, equals the left-hand side of Eq. (25). Then, the potential (24) solves the Yang-Mills equations (1). The equation of motion of the small field a µ a (6) is found to be where A ≡ B/g. Assuming that a µ a (t, x, y, z) = e −i(ωt−k·r) a µ a , where k = (k x , k y , k z ) and r = (x, y, z), Eqs. (26) are changed into the following set of algebraic equationsM t B a t = 0, and Since the equations (27), (28) and (29) are all homogeneous, they have solutions if which are the dispersion equations.
where k ≡ |k| and k T ≡ k 2 y + k 2 z . The solutions are One observes that ω 2 ± (k) and Computing the determinant of the matrixM yz B as one finds that the dispersion equation detM yz B = 0 is cubic in x ≡ ω 2 and it reads with We note that because of the square in the determinant (36) each solution of the cubic equation is doubled.
As well known, see e.g. [28], all three roots of a cubic equation can be found algebraically. Since the coefficients a 0 , a 1 , a 2 are real, the character of the roots depends on a value of the discriminant ∆ = 18 a 0 a 1 a 2 − 4 a 3 2 a 0 + a 2 1 a 2 2 − 4 a 3 1 − 27 a 2 0 .
One distinguishes three cases: • if ∆ > 0, the roots are real and distinct; • if ∆ = 0, the roots are real and at least two of them coincide; • if ∆ < 0, one root is real and the remaining two are complex.
The discriminant (41) with the coefficients (38), (39), (40) is computed as As seen, ∆ > 0 and there are three distinct real solutions of the equation detM yz B = 0. The real solutions of the cubic equation (37) can be written down in the Viète's trigonometric form [28] x n = 2 −p 3 cos 1 3 arccos 3q 2p where n = 1, 2, 3 and These formulas assume that p < 0 and that the argument of the arccosine belongs to [−1, 1]. These conditions are guaranteed as long as ∆ > 0 which is the case under consideration. We show ω 2 n (k) with n = 1, 2, 3 as a function of k 2 for Θ = π/2 in Fig. 1 and ω 2 n (k) as a function of Θ for k 2 = gB in Fig. 2. The angle Θ defines the orientation of the wave vector k with respect to the magnetic field along the axis x. Therefore, k T = k sin Θ. One observes that ω 2 1 (k) and ω 2 2 (k) are everywhere positive and the corresponding modes ±ω 1 (k) and ±ω 2 (k) are stable. There is a domain shown in the left panel of Fig. 3 where ω 2 3 (k) is negative. For comparison we show in the right panel of Fig. 3 the domain of instability of the Abelian configuration discussed in Sec. III. The Abelian mode depends only on k x = k cos Θ. One observes that the domain of instability of the Abelian configuration extends to infinity for a wave vector which is perpendicular to the magnetic field. In case of nonAbelian configuration, the domain of instability is limited for any orientation of the wave vector. We note that with the pure imaginary unstable modes which exponentially grow in time there are paired overdamped modes which exponentially decay in time.
The waves represented by the solutions of the equation detM yz B = 0 rotate in the y-z plane, as the analogous solutions of the Abelian configuration, but the rotation in the color space is not in the two-dimensional subspace but in the three-dimensional space.

V. ABELIAN CONFIGURATION OF ELECTRIC FIELD
The constant homogeneous electric field along the axis x occurs for a potentialĀ i a which is known from the Abelian theory. Specifically,Ā µ a (t, r) = δ a1 (−xE, 0, 0, 0), where E is a constant and r = (x, y, z). Then, using Eqs. (3), one finds E a (t, r) = δ a1 (E, 0, 0) and B a (t, r) = 0. We also note that the only non-vanishing elements ofF µν a arē F x0 1 = −F 0x 1 = E. The chosen potential solves the Yang-Mills equations (1) with vanishing current. The nonAbelian terms disappear because there is only one color component. We also note that the chosen potential satisfies the Lorentz gauge condition.
WhenĀ µ a (t, r) = δ a1 (−xE, 0, 0, 0), the equation (6) of a µ a becomes One sees that the color component a µ 1 decouples from the remaining two color components and it satisfies the free equation of motion. So, the functions a µ 1 represent free waves which we do not consider any more. Defining the functions Eqs. (45) provides the equations The equations of Y ± and Z ± have a diagonal form while the equations of T ± and X ± are diagonalized using Then, Eqs. (47) and (48) provide Assuming that the functions G ± , H ± , Y ± , Z ± depend on t, y, z as e −i(ωt−kyy−kzz) we find We note that one obtains exactly the same equations (54), (55) and (56) using the axial gauge a z a = 0 which was applied in Ref. [5].
Since the eigenenergy Schrödinger equation of inverted harmonic oscillator can be written as one sees that Eqs. (54) -(57) coincide with Eq. (58) under the following replacements As discussed in detail in [29], there are no normalizable solutions of the Schrödinger equation of inverted harmonic oscillator which reflects the fact that the solutions run away either to plus or minus infinite. In this sense the configuration of constant electric field is genuinely unstable.

VI. NONABELIAN CONFIGURATION OF ELECTRIC FIELD
A nonAbelian configuration ofĀ i a , which produces a constant homogeneous electric field E a = δ a1 (E, 0, 0), can be chosen asĀ where the Lorentz index µ numerates the columns and the color index a numerates the rows. The potential (60), which obviously satisfies the Lorentz gauge condition, does not solve the Yang-Mills equation (1) with j µ a = 0. Instead one gets Following [7], we assume that the current, which enters the Yang-Mills equation, equals the left-hand side of Eq. (61). Then, the potential (60) solves the Yang-Mills equations (1). The equation of motion of the small field a µ a (6) is found to be where A ≡ E/g. Assuming that a µ a (t, x, y, z) = e −i(ωt−k·r) a µ a , where k = (k x , k y , k z ) and r = (x, y, z), Eqs. (62) are changed into the following set of algebraic equationsM tx E a tx = 0, M y E a y = 0, whereM and Since the equations (63), (64) and (65) are all homogeneous, they have solutions if which are the dispersion equations.
A. Equations detM y E = 0 and detM z E = 0 Let us start with the equation detM y E = 0. Computing the determinant of the matrix (67) as the dispersion equation detM y E = 0 is again the cubic equation (37) but the coefficients are The discriminant equals 2 (k) and ω 2 3 (k). One could think that the curves of ω 2 2 (k) and ω 2 3 (k) computed for Θ = 0 and shown in Fig. 4 cross each other. However, when the curves are computed for Θ = 0.05 and shown in Fig. 5 one sees that the curves instead only approach each other. For bigger values of Θ the curves ω 2 2 (k) and ω 2 3 (k) are well separated. The phenomenon of mode coupling is discussed in detail and explained in §64 of the textbook [30].
One observes in Figs. 4 -7 that ω 2 3 (k) can be negative. Then, there is a pair of pure imaginary modes, one is unstable and one is overdamped. We show ω 2 3 (k) as a function of k 2 and Θ in the left panel of Fig. 8. In the right panel of Fig. 8  Let us now discuss the equation detM tx E = 0. We start with the simple special k = 0 when the determinant of the matrix (66) is computed as The solutions of the equation detM tx E = 0 are: the double solution ω 2 = 0 and double solutions  which give the mode frequencies where φ = arctg One observes that the modes ω (+,+) and ω (−,−) are unstable as their imaginary parts are positive. The general case of k = 0 is much more complicated than the k = 0 case, but there is an important simplification. The determinant of the matrix (66), which is the polynomial of ω 2 of order 6, appears to be the square of the polynomial of order 3 that is So, we have again the cubic equation (37) in ω 2 with the coefficients We note that because of the square in the determinant (79) each solution of the cubic equation is doubled. The discriminant ∆ with the coefficients (80), (81) and (82) equals It can be either positive or negative as shown in Figs. 9 and 10. When ∆ > 0 the solutions of the cubic equation are real and can be expressed in the Viète's trigonometric form (43). When ∆ < 0 there are one real and two complex solutions which can be found using the Cardano formula [28]. The solutions are The solutions in both regions ∆ < 0 and ∆ > 0 are shown in Figs. 11-14. Figs. 11-13 present the real solutions from the domain of ∆ > 0 combined with the real parts of the complex solutions from ∆ < 0. The solutions are shown as functions of k 2 for three values of Θ = 0.05, π/4, π/2. In Fig. 14 we show the imaginary parts of the complex solutions from the domain ∆ < 0. The whole spectrum of the solutions of the equation detM tx E = 0 is far not simple. We first note that the real parts of ω 2 1 (k) and ω 2 2 (k) are equal to each other (Reω 2 1 (k) = Reω 2 2 (k)) in the region ∆ < 0 while the imaginary parts shown in Fig. 14 are of the opposite sign (Imω 2 1 (k) = −Imω 2 2 (k)). One sees in Figs. 11-13 that the real part of low momentum solution ω 2 1 (k) or ω 2 2 (k) from the region ∆ < 0 bifurcates into two real solutions ω 2 1 (k) and ω 2 2 (k) from the region ∆ > 0. The bifurcation occurs at ∆ = 0. The solution ω 2 3 (k) is real in both regions ∆ < 0 and ∆ > 0 and it is smooth at ∆ = 0, as seen in Figs. 11-13. Fig. 11 shows that, as in case of the equation detM y E = 0 discussed in Sec. VI A, the solution ω 2 2 (k) is coupled to ω 2 3 (k). There are two unstable modes which occur as square roots of ω 2 1 (k) and ω 2 2 (k) from the region ∆ < 0. These are the modes with positive imaginary parts which at k = 0 are equal to ω (+,+) and ω (−,−) from Eq. (77). Fig. 14 shows that the imaginary parts are maximal at k = 0 and equal to Imω (+,+) = Imω (−,−) . So, the fastest modes grow as e γt with γ = Imω (+,+) ≈ 0.57 √ gE.

VII. SUMMARY, DISCUSSION AND OUTLOOK
Linear stability of classical magnetic and electric fields, which are constant and uniform, have been studied. We have considered the Abelian configurations where the fields are generated by a single-color potential linearly depending on coordinates and the nonAbelian configurations where the fields occur due to constant non-commuting potentials of different colors. While the potentials of Abelian configurations solve the sourceless Yang-Mills equations the nonAbelian configurations require non-vanishing currents to satisfy the equations. All four cases have been analyzed using the same background gauge condition. The Abelian and nonAbelian configurations are physically inequivalent and indeed the spectra of eigenmodes of small fluctuations around the background fields are different though similar. The spectra include unstable modes in all cases.
One asks how our analysis changes when the fields are uniform not in the infinite coordinate space but only in a limited domain of a size L? Assuming that the domain is a cube centered at r = 0 and demanding that the real potentials a µ a vanish at the edge of the cube, the wave vectors k x , k y , k z should be replaced as (k x , k y , k z ) −→ (2l x + 1, 2l y + 1, 2l z + 1) π L , FIG. 13. ω 2 n (k) as a function of k 2 for Θ = π/2. FIG. 14. Imaginary parts of ω 2 1 (k) and ω 2 2 (k) as functions of k 2 for Θ = 0, π/4, π/2.
where l x , l y , l z are integer numbers. Consequently, a spectrum of eigenmodes becomes discrete and some unstable modes can disappear. As an illustrative example we consider the unstable mode of the Abelian configuration of magnetic field which is ω 2 = −gB + k 2 x . After the replacement (88) it becomes ω 2 = −gB + π 2 (2l x + 1) 2 /L 2 . One sees that ω 2 < 0 at least for l x = 0 if gBL 2 > π 2 . So, the instability exists if the field is sufficiently strong and uniform over a sufficiently big domain. Since all unstable modes we found occur at small wave vectors there are analogous conditions for existence of instabilities of the fields which are uniform only in a limited spatial domain.
We intend to extend our analysis to the case of parallel electric and magnetic fields which are both present at the same time. Such a situation is expected at the earliest phase of relativistic heavy-ion collisions described in the framework of Color Glass Condensate. To make our considerations more relevant for relativistic heavy-ion collisions we plan to analyze stability not of the fields which are constant and uniform but rather the fields which are invariant under Lorentz boosts along the collision axis.