# On instabilities of perturbations in some homogeneous color-electric and
-magnetic backgrounds in $SU(2)$ gauge theory.

Divyarani C. Geetha  Faculty of Science and Technology, University of
Stavanger, 4036 Stavanger, Norway  Anders Tranberg  Faculty of Science and
Technology, University of Stavanger, 4036 Stavanger, Norway

(March 2024)

###### Abstract

We consider the instabilities of field perturbations around a homogeneous
background color-electric and/or -magnetic field in SU(2) pure gauge theory.
We investigate a number of distinct cases of background magnetic and electric
fields, and compute the dispersion relations in the linearised theory,
identifying stable and unstable momentum modes. In the case of a net
homogenous non-abelian $B$-field, we compute the non-linear (quadratic and
cubic) corrections to the equation of motion, and quantify to what extent the
instabilities are tempered by these non-linearities.

##  I Introduction

The thermalisation process of non-abelian gauge fields has been the subject of
intense study, in particular in the context of heavy ion collisions (HIC) (see
Mrow3  for a comprehensive review). In the initial stages of the collision,
anisotropic particle distributions along and perpendicular to the beam line
may create large anisotropic classical fields lappi1 ; lappi2 ; lappi3 ;
epelbaum1 ; epelbaum2 . In this background, which is not necessarily
homogeneous, perturbations may under certain conditions grow
(semi-)exponentially. These (plasma) instabilities may result in long-
wavelength fluctuations with very large occupation numbers, that may impact
the equilibration process as well Baier ; Mrow1 ; Mrow2 ; Mrow3 . It turns out
that fast hydrodynamization may likely be understood from kinetic theory with
input from perturbation theory kurk1 ; kurk2 ; kin1 ; kin2 ; kin3 ; kin4 ;
schlichting1 ; schlichting2 ; schlichting3 , but a thorough understanding of
also the non-perturbative dynamics is still of great interest.

Instabilities in anisotropic backgrounds are well-known from both abelian and
non-abelian gauge theory Randrup . While a U(1) theory produces linear field
equations of motion for the gauge field instabilities, which could then in
principle grow very large, for QCD the equations of motion are non-linear and
the instabilities could turn out to be short-lived and irrelevant. The
abelianization phenomenon Arnold  implies that there will may be directions in
field space, where non-linear contributions vanish, and where the instability
could continue unhindered for some time. Ultimately, this will depend on the
relative magnitude of the non-linearities, the range of modes that become
unstable and the self-coupling. Substantial work has been performed
analytically and numerically on plasma instabilities in QCD and QCD-like
theories (see for instance sim1 ; sim2 ; sim3 ; sim4 ; sim5 ), both in a HIC
context and for simplified models.

In the present work we will explore the instabilities of gauge field
perturbations around homogeneous and constant E and B fields in pure classical
SU(2) gauge theory. For simplicity, we will ignore any fermions coupled to the
gauge fields, which may feel and enhance the anisotropy. Our approach is
inspired by bazak1 ; bazak2 , where the authors investigated the dispersion
relations of perturbations in the linearized theory, in a variety of gauge
field backgrounds. Part of the present work is an extension of their analysis,
while later parts attempts to go beyond the linear approximation.

###  I.1 The equations of motion order by order

The classical equation of motion for pure SU(2) gauge theory reads

$\displaystyle D^{ab}_{\mu}F^{\mu\nu,b}=j^{\nu,a},$ (1)

where the covariant derivative is

$\displaystyle
D_{\mu}^{ab}=\partial_{\mu}\delta^{ab}-g\epsilon^{abc}A_{\mu}^{c},$ (2)

and the field strength tensor is

$\displaystyle
F^{\mu\nu,a}=\partial^{\mu}A^{\nu,a}-\partial^{\nu}A^{\mu,a}+g\epsilon^{abc}A^{\mu,b}A^{\nu,c}.$
(3)

The gauge field $A_{\mu,a}$ is a 12-component vector, corresponding to colour
indices $a=1,2,3$ and Lorentz indices $\mu=0,1,2,3$. Because of gauge
invariance, there is some redundancy in these degrees of freedom, and we will
make use of this point shortly. The metric signature is taken to be $(+---)$,
so that $x_{0}=x^{0}=t$ and $-x_{i}=x^{i}=(x,y,z)$, and
$\partial^{0}=\partial_{0}=\partial/\partial t$,
$-\partial^{i}=\partial_{i}=\partial/\partial x^{i}$. Then $E^{i,a}=F^{i0,a}$,
while $B^{i,a}=\frac{1}{2}\epsilon^{ijk}F^{kj,a}$.

We will assume that the gauge field may be decomposed into a background field
and a fluctuation, with the notation

$\displaystyle A^{\mu}_{a}(x,t)=\bar{A}^{\mu}_{a}(x,t)+h^{\mu}_{a}(x,t),$ (4)

and with an assumption that the fluctuations are in some sense ”small”,
$h_{a}^{\mu}\ll\bar{A}^{\mu}_{a}$. It then makes sense to expand order by
order in $h^{\mu}_{a}$.

#### Background field

The equation of motion for the background field (expanding to zero’th order in
the fluctuation $h^{\mu}_{a}$) is simply

$\displaystyle\bar{D}^{ab}_{\mu}\bar{F}^{\mu\nu,b}=j^{\nu,a},$ (5)

where the bar on $\bar{D}^{ab}_{\mu}$ and $\bar{F}^{\mu\nu,b}$ indicate that
they are expressed only in terms of the background field $\bar{A}$. We have
assumed that the external current enters in this background field equation in
its entirety. In practice, we will in the following state the background field
$\bar{A}^{\mu}_{a}$ and simply infer what the current $j^{\nu}_{a}$ is
required to be to generate such a field. In principle, this current may be
taken to arise from some distribution of colour-charged particles, but we will
not discuss this further.

#### Linear in fluctuations

To linear order in fluctuations, one finds

$\displaystyle\partial_{\mu}\partial^{\mu}h^{\nu}_{a}+g\epsilon_{abc}\left(2\partial^{\nu}\bar{A}_{\mu}^{b}h^{\mu}_{c}+2\partial_{\mu}\bar{A}^{\nu}_{c}h^{\mu}_{b}+2\bar{A}^{\mu}_{b}\partial_{\mu}h^{\nu}_{c}\right)$
$\displaystyle-g^{2}\left(2\bar{A}_{\mu}^{c}\bar{A}^{\nu}_{a}h^{\mu}_{c}-2\bar{A}^{\nu}_{c}\bar{A}_{\mu}^{a}h^{\mu}_{c}+\bar{A}_{\mu}^{c}\bar{A}^{\mu}_{c}h^{\nu}_{a}-\bar{A}_{\mu}^{c}\bar{A}^{\mu}_{a}h^{\nu}_{c}\right)=0.$

We adopt the Lorenz gauge $\partial_{\mu}\bar{A}^{\mu}_{a}=0$ for the
background field, and have imposed the background gauge fixing condition on
the fluctuations, $\bar{D}^{ab}_{\mu}h^{\mu,b}=0$, where $\bar{D}_{\mu}^{ab}$
again denotes the covariant derivative in terms of the background field
$\bar{A}_{\mu}^{b}$ only. In a more compact notation, this may be rewritten as

$\displaystyle\left[g^{\mu\nu}(\bar{D}_{\rho}\bar{D}^{\rho})_{ac}+2g\epsilon_{abc}\bar{F}_{b}^{\mu\nu}\right]h_{\nu,c}=\mathcal{M}^{\mu\nu}_{ac}[12\times
12]h_{\nu,c}=0.$

As we will see, in momentum space this turns into a $[12\times 12]$ matrix
equation for $h^{\mu}_{a}$, which has unstable modes for certain choices of
$\bar{A}^{\mu}_{a}$. Some of these correspond to non-zero homogeneous colour-
electric and/or -magnetic fields, and we will in section II consider a fairly
broad subset. Given $\bar{A}^{\mu}_{a}$, one may compute the dispersion
relations for these modes.

#### Quadratic in fluctuations

At early times, when $h^{\mu}_{a}$ is truly $\ll\bar{A}^{\mu}_{a}$, the linear
approximation is valid and unstable modes grow exponentially. However, non-
linearities will eventually impact the evolution. At second order in
fluctuations, we have

$\displaystyle
g\epsilon_{abc}\left(2\partial_{\mu}h^{\nu}_{c}h^{\mu}_{b}+\partial^{\nu}h_{\mu}^{b}h^{\mu}_{c}\right)$
$\displaystyle-g^{2}\left(2\bar{A}^{\mu}_{c}h_{\mu}^{c}h^{\nu}_{a}-2\bar{A}^{\mu}_{a}h_{\mu}^{c}h^{\nu}_{c}+\bar{A}^{\nu}_{a}h_{\mu}^{c}h^{\mu}_{c}-\bar{A}^{\nu}_{c}h_{\mu}^{c}h^{\mu}_{a}\right).$

which upon taking the expectation values for the fluctuations, and truncating
at quadratic order, naturally adds to the background equation of motion (5).
Defining

$\displaystyle C^{\mu\nu}_{ab}=\langle
h^{\mu}_{a}h^{\nu}_{b}\rangle,\quad\bar{j}^{\nu}_{a}=g\epsilon_{abc}\langle
2\partial_{\mu}h^{\nu}_{c}h^{\mu}_{b}+\partial^{\nu}h_{\mu}^{b}h^{\mu}_{c}\rangle,$
(9)

we have schematically

$\displaystyle\bar{D}^{ab}_{\mu}\bar{F}^{\mu\nu,b}-g^{2}\bar{M}^{\nu\mu}_{ab}[C]\bar{A}_{\mu,b}=j^{\nu}_{a}-\bar{j}^{\nu}_{a}.$
(10)

where $\bar{M}$ is some matrix involving the $C$ correlators in (9). This in
effect changes both the background equation of motion and the effective
current, and if $j^{\nu}_{a}$ is constant will eventually ruin the assumption
of a constant background electric/magnetic fields. We will investigate the
effect of this non-linear back reaction below.

#### Cubic in fluctuations

Finally, the equation of motion has cubic contributions in the fluctuations,
which read

$\displaystyle-g^{2}\left[h^{\nu}_{a}h_{\mu}^{c}h^{\mu}_{c}-h^{\mu}_{a}h_{\mu}^{c}h^{\nu}_{c}\right],$
(11)

and which again upon taking expectation values naturally add to the linear
equation for the modes as

$\displaystyle\left[g^{\mu\nu}(\bar{D}_{\rho}\bar{D}^{\rho})_{ac}+2g\epsilon_{abc}\bar{F}_{b}^{\mu\nu}-g^{2}\hat{M}^{\mu\nu}_{ac}(C)\right]h_{\nu,c}=0.$

The matrix $\hat{M}$ is expressed in terms if the $C$ correlators, providing a
corrected linear equation for the fluctuations $h^{\mu}_{a}$. This gives rise
to new time-dependent dispersion relations, altering the pattern of
instabilities.

###  I.2 This work

In section II, we will consider a fairly general set of homogeneous,
anisotropic field backgrounds. For three of these we will in section A compute
the dispersion relations of the fluctuation modes explicitly, and establish
for which momentum regions the fluctuations are unstable. The remaining cases
are more involved, and we will simply sketch how one would go about finding
dispersion relations for them in Appendix B, relying heavily on bazak1 ;
bazak2 . In section III we choose one of the three cases and compute the
eigenmodes explicitly, in order to compute the $C$ correlators. This will in
turn allow us to compute the cubic correction matrix $\hat{M}$ and the
quadratic corrections $\bar{M}$ and $\bar{J}$, from which we will be able to
study the evolution of the instabilities. We provide some discussion,
conclusions and suggestions for future work in section IV.

##  II Anisotropic backgrounds with constant, homogeneous electric/magnetic
fields

###  II.1 Electric field only

We will first consider a constant, homogeneous colour-electric field only, and
set up our procedure and notation. Firstly, without loss of generality, we can
choose the electric field in the $i=1$, $a=1$ direction. We hence require that

$\displaystyle\bar{F}^{i0}_{a}=E^{i}_{a}=\delta^{i1}_{a1}E.$ (13)

Since

$\displaystyle\bar{F}^{i0}_{a}=\partial^{i}\bar{A}^{0}_{a}-\partial^{0}\bar{A}^{i}_{a}+g\epsilon_{abc}\bar{A}^{i}_{b}\bar{A}^{0}_{c},$
(14)

such an electric field may correspond to a multitude of choices for
$\bar{A}^{\mu}_{a}$ (some of which may be gauge-equivalent). We choose to
restrict ourselves to two cases: One where only the linear-in-$A$ (”Abelian”)
derivative part of $\bar{F}$ is non-zero. And one where only the quadratic-
in-$A$ (”non-Abelian”) part is non-zero. A homogeneous, time-independent
Abelian electric field may then be written (in a compact notation we will use
repeatedly for the 12-component objects such as $\bar{A}^{\mu}_{a}$,
$j^{\mu}_{a}$), as

$\displaystyle\bar{A}^{\mu}_{a}$
$\displaystyle=\left[\left(\begin{array}[]{c}\bar{A}^{0}_{1}\\\
\bar{A}^{0}_{2}\\\
\bar{A}^{0}_{3}\end{array}\right)\left(\begin{array}[]{c}\bar{A}^{1}_{1}\\\
\bar{A}^{1}_{2}\\\
\bar{A}^{1}_{3}\end{array}\right)\left(\begin{array}[]{c}\bar{A}^{2}_{1}\\\
\bar{A}^{2}_{2}\\\
\bar{A}^{2}_{3}\end{array}\right)\left(\begin{array}[]{c}\bar{A}^{3}_{1}\\\
\bar{A}^{3}_{2}\\\ \bar{A}^{3}_{3}\end{array}\right)\right]$ (40)
$\displaystyle=\left[\left(\begin{array}[]{c}-x\,C_{1}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}-t\,C_{2}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\right],$

so that

$\displaystyle E^{1}_{1}\equiv E=C_{1}+C_{2}.$ (42)

This corresponds to a vanishing current $j^{\nu,a}=0$. Similarly a homogeneous
time-independent non-Abelian electric field may be written as

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ D_{2}\\\
D_{1}\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{4}\\\
D_{3}\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\right],$ (55)

to find

$\displaystyle E^{1}_{1}\equiv E=g(D_{4}D_{1}-D_{3}D_{2}).$ (56)

This corresponds to a current with four non-vanishing components:

$\displaystyle j^{\nu}_{a}=-gE\left[\left(\begin{array}[]{c}0\\\ D_{3}\\\
-D_{4}\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{1}\\\
-D_{2}\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\right].$ (69)

###  II.2 Magnetic field only

For a homogeneous magnetic field we can again choose the $\mu=1$, $a=1$
direction, and write

$\displaystyle\frac{1}{2}\epsilon^{ijk}F^{kj}_{a}=B_{a}^{i}=\delta^{i1}_{a1}B.$
(71)

An Abelian background field can be parametrized as

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}-z\,C_{4}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}y\,C_{3}\\\ 0\\\
0\end{array}\right)\right],$ (84)

so that

$\displaystyle B^{1}_{1}=B=C_{3}+C_{4}.$ (85)

This again corresponds to a vanishing current $j^{\nu,a}=0$. The non-Abelian
realisation is

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{6}\\\
D_{5}\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{8}\\\
D_{7}\end{array}\right)\right],$ (98)

with

$\displaystyle B^{1}_{1}=B=g(D_{5}D_{8}-D_{6}D_{7}).$ (99)

This again corresponds to a current with four non-vanishing components:

$\displaystyle j^{\nu}_{a}=-gB\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{7}\\\
-D_{8}\end{array}\right)\left(\begin{array}[]{c}0\\\ -D_{5}\\\
D_{6}\end{array}\right)\right].$ (112)

These expressions are minor generalisations to bazak1 , where the authors
consider Abelian $E$ with $C_{1}=E$, $C_{2}=0$, Abelian $B$ with $C_{3}=B$,
$C_{4}=0$, non-Abelian $E$ with $D_{2}=-D_{3}=\sqrt{E/g}$, $D_{1}=D_{4}=0$ and
non-Abelian $B$ with $D_{5}=D_{8}=\sqrt{B/g}$, $D_{6}=D_{7}=0$. In bazak2  the
same authors in addition considered arbitrary $D_{2},D_{3}$, and
$D_{5},D_{8}$, respectively. We see, that there are a host of other
possibilities for choosing the $C_{i}$ or the $D_{i}$ for a given $E$ or $B$.
It will not be possible to keep full generality in the following, but we will
further investigate some of the combinations in terms of the resulting
dispersion relations.

###  II.3 Electric and magnetic fields, combined

The natural generalisation is to consider the combination of an electric field
$E^{1}_{1}$ together with a magnetic field that is/is not aligned in space
and/or colour. Allowing for all combinations of Abelian/non-Abelian
$\bar{A}^{\mu}_{a}$ one might expect 16 distinct cases. Some combinations are
however not possible since additional components of magnetic and electric
field are sourced. The task is therefore to identify vector potentials with
the property that they generate only $E^{1}_{1}$ and then some combination of
$B^{i}_{a}$, enforcing that they are either of Abelian or non-Abelian type,
and that $E$ and $B$ are constant in space and time. One then finds a limited
set of options, which we will briefly review in the following.

#### Abelian $E$ and Abelian $B$

Choosing a vector potential of the form

$\displaystyle\bar{A}^{\mu}_{a}=\left[\begin{pmatrix}-x\,C_{1}\\\ \\\ \\\ 0\\\
0\\\ \end{pmatrix}\begin{pmatrix}-t\,C_{2}-y\,C_{16}\\\ +z\,C_{9}\\\ \\\ 0\\\
0\\\ \end{pmatrix}\begin{pmatrix}-z\,C_{4}\\\ +x\,C_{15}\\\ \\\ 0\\\ 0\\\
\end{pmatrix}\begin{pmatrix}-x\,C_{10}\\\ +y\,C_{3}\\\ \\\ 0\\\ 0\\\
\end{pmatrix}\right],$

gives Abelian $E$\- and $B$-fields aligned in colour, with

$\displaystyle E^{1}_{1}=C_{1}+C_{2},$ $\displaystyle
B^{i}_{a}=g\left[\left(\begin{array}[]{c}C_{3}+C_{4}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}C_{9}+C_{10}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}C_{15}+C_{16}\\\ 0\\\
0\end{array}\right)\right].$ (123)

Explicitly setting $\bar{A}^{0}_{1}=0$ ($C_{1}=0$), this opens up a few more
options:

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}-t\,C_{2}-y\,C_{16}+z\,C_{9}\\\
-y\,C_{18}+z\,C_{11}\\\
-y\,C_{20}+z\,C_{13}\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\right],$ (137)

which gives Abelian $E$\- and $B$-fields orthogonal in space:

$\displaystyle E^{1}_{1}=C_{2},$ $\displaystyle
B^{i}_{a}=g\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}C_{9}\\\ C_{11}\\\
C_{13}\end{array}\right)\left(\begin{array}[]{c}C_{16}\\\ C_{18}\\\
C_{20}\end{array}\right)\right].$ (148)

#### Abelian E and non-Abelian B

A vector potential of the form

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}-tC_{2}\\\ D_{4}\\\
D_{3}\end{array}\right)\left(\begin{array}[]{c}D_{11}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}D_{12}\\\ 0\\\
0\end{array}\right)\right],$ (162)

gives Abelian $E$\- and non-Abelian $B$-fields orthogonal in space and colour:

$\displaystyle E^{1}_{1}=C_{2},$ $\displaystyle
B^{i}_{a}=g\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{3}D_{12}\\\
-D_{4}D_{12}\end{array}\right)\left(\begin{array}[]{c}0\\\ -D_{3}D_{11}\\\
D_{4}D_{11}\end{array}\right)\right].$ (172)

#### Non-Abelian E and non-Abelian B

Finally, we may choose

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ 0\\\
D_{1}\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{4}\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
D_{5}\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
D_{7}\end{array}\right)\right],$ (186)

which gives non-Abelian $E$\- and $B$-fields orthogonal in colour, but aligned
in space.

$\displaystyle E^{1}_{1}=gD_{1}D_{4},$ $\displaystyle
B^{i}_{a}=g\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}D_{4}D_{7}\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}-D_{4}D_{5}\\\ 0\\\
0\end{array}\right)\right].$ (196)

The combination

$\displaystyle\bar{A}^{\mu}_{a}=\left[\left(\begin{array}[]{c}0\\\ D_{2}\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
D_{3}\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{6}\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{8}\\\
0\end{array}\right)\right],$ (210)

is equivalent.

In summary, we have seen that we can find:

  * •

Abelian $E$ or $B$, or non-Abelian $E$ or $B$ with a somewhat generalised
parametrisation of what was reported in bazak1 .

  * •

Space-aligned, colour-aligned $E$ and $B$: Possible as an Abelian/Abelian
combination. $C_{1}\neq 0$, $C_{3}\neq 0$ is the special case reported in
bazak1 .

  * •

Space-aligned, colour-orthogonal $E$ and $B$: Can not be realised in our
setup.

  * •

Space-orthogonal, colour-aligned $E$ and $B$: Possible as Abelian/Abelian and
non-Abelian/non-Abelian combinations.

  * •

Space-orthogonal, colour-orthogonal $E$ and $B$: Possible as Abelian/Abelian
and Abelian/non-Abelian combinations.

We also see that there are several computationally distinct cases, depending
on whether the vector potential is itself space- and/or time-dependent. We
will press on with the three purely non-Abelian cases (E only, B only, E and
B), where $\bar{A}^{\mu}_{a}$ is homogeneous and constant. We refer to
Appendix B for a brief discussion of the issues involved for the Abelian
vector potentials.

###  II.4 Non-Abelian B, E=0

Let us return to the case of a non-Abelian B-field and zero E-field. In
4-dimensional momentum space, the linear equation for the modes of the
fluctuations reduce to a set of 12 coupled linear algebraic equations for
$h^{\mu}_{a}$, which however split into 3 independent sections.

$\displaystyle\mathcal{M}^{\mu\nu}_{ac}\left[12\times
12\right]h_{\nu}^{c}=\left(\begin{array}[]{cccc}M_{00}&0&0&0\\\
0&M_{11}&0&0\\\ 0&0&M_{22}&M_{23}\\\ 0&0&M_{32}&M_{33}\\\
\end{array}\right)\left(\begin{array}[]{c}h^{0}_{1,2,3}\\\ h^{1}_{1,2,3}\\\
h^{2}_{1,2,3}\\\ h^{3}_{1,2,3}\\\ \end{array}\right),$ (219)

where in addition $M_{00}=M_{11}=M_{22}=M_{33}$ are all hermitian, and
$M_{23}=M_{32}^{\dagger}$, so that the whole matrix $M$ is hermitian. Hence
the equations for $h^{0}_{1,2,3}$ ($M_{00}$) and $h^{1}_{1,2,3}$ ($M_{11}$)
are identical

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}+D_{5}^{2}+D_{6}^{2}+D_{7}^{2}+D_{8}^{2}&2i(D_{5}k_{y}+D_{7}k_{z})&-2i(D_{6}k_{y}+D_{8}k_{z})\\\
-2i(D_{5}k_{y}+D_{7}k_{z})&\Box_{k}+D_{5}^{2}+D_{7}^{2}&-(D_{5}D_{6}+D_{7}D_{8})\\\
2i(D_{6}k_{y}+D_{8}k_{z})&-(D_{5}D_{6}+D_{7}D_{8})&\Box_{k}+D_{6}^{2}+D_{8}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{0,1}_{1}\\\
h^{0,1}_{2}\\\ h^{0,1}_{3}\\\ \end{array}\right)=0,$ (227)

while the $h^{2,3}_{1,2,3}$-sections mix ($M_{22,23,32,33}$)

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}+D_{5}^{2}+D_{6}^{2}+D_{7}^{2}+D_{8}^{2}&2i(D_{5}k_{y}+D_{7}k_{z})&-2i(D_{6}k_{y}+D_{8}k_{z})\\\
-2i(D_{5}k_{y}+D_{7}k_{z})&\Box_{k}+D_{5}^{2}+D_{7}^{2}&-(D_{5}D_{6}+D_{7}D_{8})\\\
2i(D_{6}k_{y}+D_{8}k_{z})&-(D_{5}D_{6}+D_{7}D_{8})&\Box_{k}+D_{6}^{2}+D_{8}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{2}_{1}\\\
h^{2}_{2}\\\ h^{2}_{3}\\\
\end{array}\right)+\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&-2B\\\
0&2B&0\end{array}\right)\left(\begin{array}[]{c}h^{3}_{1}\\\ h^{3}_{2}\\\
h^{3}_{3}\\\ \end{array}\right)=0,$ (241)

and

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}+D_{5}^{2}+D_{6}^{2}+D_{7}^{2}+D_{8}^{2}&2i(D_{5}k_{y}+D_{7}k_{z})&-2i(D_{6}k_{y}+D_{8}k_{z})\\\
-2i(D_{5}k_{y}+D_{7}k_{z})&\Box_{k}+D_{5}^{2}+D_{7}^{2}&-(D_{5}D_{6}+D_{7}D_{8})\\\
2i(D_{6}k_{y}+D_{8}k_{z})&-(D_{5}D_{6}+D_{7}D_{8})&\Box_{k}+D_{6}^{2}+D_{8}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{3}_{1}\\\
h^{3}_{2}\\\ h^{3}_{3}\\\
\end{array}\right)+\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&2B\\\
0&-2B&0\end{array}\right)\left(\begin{array}[]{c}h^{2}_{1}\\\ h^{2}_{2}\\\
h^{2}_{3}\\\ \end{array}\right)=0.$ (255)

To somewhat mitigate notational clutter we have absorbed $g$ into $\bar{A}$,
$gD_{i}\rightarrow D_{i}$, and defined the box operator in momentum space
$\Box_{k}=-\omega^{2}+k_{x}^{2}+k^{2}_{T}$, with
$k_{T}^{2}=k_{y}^{2}+k_{z}^{2}$. We note that the $D_{5,6,7,8}$ do not only
appear in the combination $B=D_{5}D_{8}-D_{6}D_{7}$, but depend on the $D_{i}$
individually.

Figure 1: Non-Abelian B: The three eigenvalues in the 0, 1 sectors for ${\bf
k}=0$, general $D_{i}$ (top) and for $\sqrt{B}=D_{6}/\lambda=-D_{7}\lambda$
(bottom).

#### Non-Abelian B: 0 and 1 sector

Figure 2: Non-Abelian B: The three eigenvalues in the 0, 1 sectors for
$\sqrt{B}=D_{6}=-D_{7}$, with $k_{T}=0$ (top) and $k_{x}=0$ (bottom).

Figure 3: Non-Abelian B: The three eigenvalues in the 0, 1 sectors for
$\sqrt{B}=D_{6}=-D_{7}$, with non-zero $D_{5}=\sqrt{nB}$, $n=1$. For non-zero
$k_{x}$ (top), $k_{y}$ (middle) and $k_{z}$ (bottom).

The condition on $\omega$ to satisfy the equation is that
$\textrm{Det}[M_{00}]=0=\textrm{Det}[M_{11}]$, which reduces to the cubic
equation

$\displaystyle\Box_{k}^{3}+2\alpha_{B}\Box_{k}^{2}+\beta_{B}\Box_{k}+\gamma_{B}=0,$
(257)

with

$\displaystyle\alpha_{B}$ $\displaystyle=$
$\displaystyle(D_{5}^{2}+D_{6}^{2}+D_{7}^{2}+D_{8}^{2})\geq 2B,$ (258)
$\displaystyle\beta_{B}$ $\displaystyle=$
$\displaystyle\alpha_{B}^{2}+B^{2}-4\left((D_{5}k_{y}+D_{7}k_{z})^{2}+(D_{6}k_{y}+D_{8}k_{z})^{2}\right),$
$\displaystyle\gamma_{B}$ $\displaystyle=$ $\displaystyle
B^{2}(\alpha_{B}-4k_{T}^{2}).$ (260)

Let us first consider the zero-momentum case $k_{T}=k_{x}=0$, so that
$\omega^{2}=-\Box_{k}$, and

$\displaystyle-\left(\omega^{2}\right)^{3}+2\alpha_{B}\left(\omega^{2}\right)^{2}-\left(\alpha_{B}^{2}+B^{2}\right)\omega^{2}+B^{2}\alpha_{B}=0.$

The solutions for $\omega^{2}$ only depend on the combination $\alpha_{B}/B$,

$\displaystyle\frac{\omega^{2}}{B}=\frac{\alpha_{B}}{B},\quad\frac{\omega^{2}}{B}=\frac{1}{2}\left(\frac{\alpha_{B}}{B}\pm\sqrt{\frac{\alpha_{B}^{2}}{B^{2}}-4}\right),$
(262)

and they are all real and positive (see Figure 1, top). For instance in the
case $D_{5}=D_{8}=0$, $D_{6}=-D_{7}=\sqrt{B}$, we find $\alpha_{B}/B=2$ so
that the eigenvalues reduce to $\omega^{2}=(2B,B,B)$ . For the more general
case $D_{6}=\sqrt{B}\lambda$, $-D_{7}=\sqrt{B}/\lambda$ for some $\lambda$
bazak1 , $\alpha_{B}/B=\lambda^{2}+\lambda^{-2}\geq 2$, and the solutions are
again real and positive (see Figure 1, bottom). As a sanity check, we note
that when $B=0$, the eigenvalues are $\omega^{2}=(0,\alpha_{B},\alpha_{B})$.
These are real and positive, and only zero when all the $D_{i}$ are zero.

Figure 4: Non-Abelian B: The lowest eigenvalues in the 0, 1 sectors for
$\sqrt{B}=D_{6}=-D_{7}$, with non-zero $D_{5}=\sqrt{nB}$, $n=0,1,2,3,4$. For
non-zero $k_{y}$ (top) and $k_{z}$ (bottom).

For non-zero momentum $k_{T},k_{x}$, the explicit realisation of $D_{5,..,8}$
becomes more important. Again for the case $D_{6}=-D_{7}=\sqrt{B}$ one finds
($\alpha_{B}/B=2$, $\beta_{B}/B^{2}=5-4k_{T}^{2}/B$,
$\gamma_{B}/B^{3}=2-4k_{T}^{2}/B$),

$\displaystyle\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)^{3}+4\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)^{2}$
$\displaystyle+\left(5-4\frac{k_{T}^{2}}{B}\right)\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)+2-4\frac{k_{T}^{2}}{B}=0,$

which has the solutions (see Figure 2)

$\displaystyle\frac{\omega^{2}}{B}=1+\frac{k_{x}^{2}}{B}+\frac{k_{T}^{2}}{B},\quad\frac{\omega^{2}}{B}=\frac{3}{2}+\frac{k_{x}^{2}}{B}+\frac{k_{T}^{2}}{B}\pm\frac{1}{2}\sqrt{1+16\frac{k_{T}^{2}}{B}}.$

These are real and positive for all $B,k_{x},k_{T}$, and hence the 0 and 1
components of the fluctuations are stable for this choice of $D_{5,..,8}$.

For fixed $B$, a 3-dimensional space of realisations of $D_{5,..,8}$ is
available, which we will not fully explore here. We will however briefly
consider again the case $D_{6}=-D_{7}=\sqrt{B}$, fixing $D_{8}=0$, allowing us
to choose any value for $D_{5}$ without changing $B$. We choose
$D_{5}=\sqrt{nB}$, and vary $n$. In Figure 3 we show the eigenvalues for
$n=1$, and their dependence on $k_{x,y,z}$ (respectively, keeping the other
two constant). We note how the $k_{T}=\sqrt{k_{y}^{2}+k_{z}^{2}}$-dependence
in Figure 2 is resolved into different dependence on $k_{y}$ and $k_{z}$. We
can further examine the $n$-dependence of the lowest eigenvalue, as function
of $k_{y}$ and $k_{z}$, shown in Figure 4. We see a non-trivial dependence,
but the eigenvalue remains real and positive.

#### Non-Abelian B: 2 and 3 sector

Figure 5: Non-Abelian B: The three double eigenvalues in the 2, 3 sectors for
${\bf k}=0$, general $D_{i}$ (top) and for
$\sqrt{B}=D_{6}/\lambda=-D_{7}\lambda$ (bottom).

Figure 6: Non-Abelian B: The three eigenvalues in the 2, 3 sectors for
$\sqrt{B}=D_{6}=-D_{7}$, with $k_{T}=0$ (top) and $k_{x}=0$ (bottom).

For the 2-3 component, the requirement $\textrm{Det}[M_{22,23,32,33}]=0$
amounts to (compare to (257))

$\displaystyle\left[\Box_{k}^{3}+2\alpha_{B}\Box_{k}^{2}+(\beta_{B}-4B^{2})\Box_{k}+\gamma_{B}-4B^{2}\alpha_{B}.\right]^{2}=0,$
(265)

We may focus on the expression inside the square bracket, and infer that any
solution is a double solution. For $k_{T},k_{x}=0$ this becomes

$\displaystyle-(\omega^{2})^{3}+2\alpha_{B}(\omega^{2})^{2}-(\alpha_{B}^{2}-3B^{2})\omega^{2}-3B^{2}\alpha_{B}=0,$

which has solutions

$\displaystyle\frac{\omega^{2}}{B}=\frac{\alpha_{B}}{B},\quad\omega^{2}=\frac{1}{2}\left(\frac{\alpha_{B}}{B}\pm\sqrt{\frac{\alpha_{B}^{2}}{B^{2}}+12}\right).$
(267)

One of these eigenvalues is negative for all $\alpha_{B}/B$ (see Figure 5,
top). In particular for $D_{5}=D_{8}=0$, $D_{6}=-D_{7}=\sqrt{B}$ we find
$\omega^{2}=(3B,2B,-B)$, while the parametrization $D_{6}=\sqrt{B}\lambda$,
$-D_{7}=\sqrt{B}/\lambda$ gives the negative eigenvalue for $\omega^{2}$ in
the interval $[-1,0]$ (see Figure 5, bottom). Again, when $B=0$,
$\omega^{2}=(0,\alpha_{B},\alpha_{B})$.

Figure 7: Non-Abelian B: The lowest (negative) eigenvalue in the 2, 3 sectors
for $\sqrt{B}=D_{6}=-D_{7}$, with non-zero $D_{5}=\sqrt{nB}$, $n=0,1,2,3,4$.
For non-zero $k_{x}$ (top) $k_{y}$ (middle) and $k_{z}$ (bottom).

For non-zero momentum, we consider again $D_{6}=-D_{7}=\sqrt{B}$, for which we
need to solve

$\displaystyle\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)^{3}+4\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)^{2}$
$\displaystyle+\left(1-4\frac{k_{T}^{2}}{B}\right)\left(-\frac{\omega^{2}}{B}+\frac{k_{T}^{2}}{B}+\frac{k_{x}^{2}}{B}\right)-6-4\frac{k_{T}^{2}}{B}=0,$

This is a general cubic equation, which may be solved straightforwardly, to
reveal (Figure 6) that the (double) unstable mode becomes stable as $k_{x}$
and/or $k_{T}$ increase. Explicitly, for $k_{T}=0$, the unstable region is
$|k_{x}|<1$, while for $k_{x}=0$ it is
$|k_{T}|<\sqrt{(3-\sqrt{8})^{1/3}+(3+\sqrt{8})^{1/3}}\simeq 1.53$.

We can again extend our analysis somewhat, by allowing $D_{5}=\sqrt{n}B$,
while keeping $D_{8}=0$. In Figure 7 we show the negative eigenvalue and its
momentum dependence (nonzero $k_{x}$, $k_{y}$ or $k_{z}$, respectively), for
different values of $n=0,1,2,3,4$. We see that larger $n$ reduces the region
of instability, and we find that for $n\rightarrow\infty$, $\omega^{2}({\bf
k=0})\rightarrow 0$.

Hence, in the presence of a constant, homogeneous non-Abelian magnetic field,
for this choice of $D_{5,..,8}$, of the 12 sets of momentum modes present, two
are potentially unstable, growing as

$\displaystyle\propto\exp{\left(\pm|\omega|t\right)},$ (269)

with an $\omega$ given by the two lowest, degenerate solutions to (265)
(exemplified by the $-$ solution in the special case of (267)). Furthermore,
the region of instabilibity is centered around the origin (${\bf k}=0$), and
that for large enough ${\bf k}$, these modes are again stable. Hence only a
finite IR region of momentum modes are unstable, and this region is not only
dependent on $B$, but also on the concrete choice of $D_{5,..,8}$, exemplified
here by varying $D_{5}$.

In section III we will further find the eigenmodes of the system for non-
Abelian B only, with $\sqrt{B}=D_{6}=-D_{7}$, $D_{5}=D_{8}=0$, compute the
quadratic and cubic corrections to the evolution equations and solve them
including these corrections.

###  II.5 Non-Abelian E, B=0

For a non-Abelian background E-field and zero B-field, the linear equation for
the Fourier modes again splits into 3 independent sections. The section for
$h^{2}_{1,2,3}$ ($M_{2}$) and $h^{3}_{1,2,3}$ ($M_{3}$) are now identical and
decoupled from each other, and read (absorbing again $g$ as $gD_{i}\rightarrow
D_{i}$)

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}-D_{1}^{2}-D_{2}^{2}+D_{3}^{2}+D_{4}^{2}&-2i(D_{1}\omega-
D_{3}k_{x})&2i(D_{2}\omega-D_{4}k_{x})\\\ +2i(D_{1}\omega-
D_{3}k_{x})&\Box_{k}-D_{1}^{2}+D_{3}^{2}&D_{1}D_{2}-D_{3}D_{4}\\\
-2i(D_{2}\omega-
D_{4}k_{x})&D_{1}D_{2}-D_{3}D_{4}&\Box_{k}-D_{2}^{2}+D_{4}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{2,3}_{1}\\\
h^{2,3}_{2}\\\ h^{2,3}_{3}\\\ \end{array}\right)=0.$ (276)

while the $h^{0,1}_{1,2,3}$-sections mix ($M_{00},M_{01},M_{10},M_{11}$)

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}-D_{1}^{2}-D_{2}^{2}+D_{3}^{2}+D_{4}^{2}&-2i(D_{1}\omega-
D_{3}k_{x})&2i(D_{2}\omega-D_{4}k_{x})\\\ 2i(D_{1}\omega-
D_{3}k_{x})&\Box_{k}-D_{1}^{2}+D_{3}^{2}&D_{1}D_{2}-D_{3}D_{4}\\\
-2i(D_{2}\omega-
D_{4}k_{x})&D_{1}D_{2}-D_{3}D_{4}&\Box_{k}-D_{2}^{2}+D_{4}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{0}_{1}\\\
h^{0}_{2}\\\ h^{0}_{3}\\\
\end{array}\right)+\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&-2E\\\
0&2E&0\end{array}\right)\left(\begin{array}[]{c}h^{1}_{1}\\\ h^{1}_{2}\\\
h^{1}_{3}\\\ \end{array}\right)=0,$ (290)

and

$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}-D_{1}^{2}-D_{2}^{2}+D_{3}^{2}+D_{4}^{2}&-2i(D_{1}\omega-
D_{3}k_{x})&2i(D_{2}\omega-D_{4}k_{x})\\\ 2i(D_{1}\omega-
D_{3}k_{x})&\Box_{k}-D_{1}^{2}+D_{3}^{2}&D_{1}D_{2}-D_{3}D_{4}\\\
-2i(D_{2}\omega-
D_{4}k_{x})&D_{1}D_{2}-D_{3}D_{4}&\Box_{k}-D_{2}^{2}+D_{4}^{2}\end{array}\right)\left(\begin{array}[]{c}h^{1}_{1}\\\
h^{1}_{2}\\\ h^{1}_{3}\\\
\end{array}\right)+\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&-2E\\\
0&2E&0\end{array}\right)\left(\begin{array}[]{c}h^{0}_{1}\\\ h^{0}_{2}\\\
h^{0}_{3}\\\ \end{array}\right)=0.$ (304)

#### Non-Abelian E: 2 and 3 sector

Figure 8: Non-Abelian E: The three eigenvalues in the 2, 3 sectors for
$\sqrt{E}=D_{1}/\lambda=D_{4}\lambda$. Figure 9: Non-Abelian E: Regions of
non-zero imaginary part of the eigenvalues in $\alpha_{+}/\alpha_{-}$ space.
Overlaid the curve $\sqrt{E}=D_{1}/\lambda=D_{4}\lambda$ (red) and curves when
further adding non-zero $D_{3}=\sqrt{nB}$, $n=0.5$, $1$ (green).

The condition on $\omega$ is now that
$\textrm{Det}[M_{22}]=0=\textrm{Det}[M_{33}]$, which is again a cubic equation
of the form

$\displaystyle\Box_{k}^{3}+2\alpha_{-}\Box_{k}^{2}+\beta_{E}\Box_{k}+\gamma_{E}=0,$
(306)

where still $\Box_{k}=-\omega^{2}+k_{x}^{2}+k_{T}^{2}$,
$k_{T}^{2}=k_{y}^{2}+k_{z}^{2}$, and

$\displaystyle\alpha_{-}$ $\displaystyle=$ $\displaystyle-
D_{1}^{2}-D_{2}^{2}+D_{3}^{2}+D_{4}^{2},$ (307) $\displaystyle\alpha_{+}$
$\displaystyle=$ $\displaystyle D_{1}^{2}+D_{2}^{2}+D_{3}^{2}+D_{4}^{2},$
(308) $\displaystyle\beta_{E}$ $\displaystyle=$
$\displaystyle\alpha_{-}^{2}-E^{2}-4\left((D_{1}\omega-
D_{3}k_{x})^{2}+(D_{2}\omega-D_{4}k_{x})^{2}\right),$
$\displaystyle\gamma_{E}$ $\displaystyle=$
$\displaystyle-E^{2}(\alpha_{-}-4(k_{x}^{2}-\omega^{2})),$ (310)

keeping in mind that $E=D_{1}D_{4}-D_{2}D_{3}$. For $k_{T}=k_{x}=0$, we
obtain111Note that because $\omega$ appears not only in the $\Box_{k}$, this
is different from the $B\neq 0$, $E=0$ case.

$\displaystyle-(\omega^{2})^{3}+2\alpha_{+}(\omega^{2})^{2}-(\alpha_{-}^{2}+3E^{2})\omega^{2}-E^{2}\alpha_{-}=0.$

There is now one real solution for $\omega^{2}$, while the other two solutions
are in general non-real. Hence already in the $2$ and $3$ sectors, there are
unstable modes in the presence of finite $E$. However, for the special case
$D_{1}=D_{4}=\sqrt{E}$, $D_{2}=D_{3}=0$, the eigenvalues are
$\omega^{2}=(0,E,3E)$, while for the generalisation $D_{1}=\sqrt{E}\lambda$,
$D_{4}=\sqrt{E}/\lambda$, one of the eigenvalues is negative in the interval
$\lambda\in[0,1]$ and the other two are real and positive (see Figure 8). It
appears that the curve traced out by $D_{1}D_{4}=E$ in the space of
$D_{1,..,4}$ is quite special, avoiding the regions where the eigenvalues are
non-real. In Figure 9, we show the $\alpha_{+}/\alpha_{-}$ plane, and the
contours of non-zero imaginary value of the eigenvalues. Overlaid in red, the
curve traced out at $\lambda$ is varied, which precisely avoids those regions.
Overlaid in green, the curves traced out by adding non-zero $D_{3}=\sqrt{nE}$,
$n=0.5,1$. As another sanity check, we note that for $E=0$,
$\omega^{2}=(0,\alpha_{+}\pm\sqrt{\alpha_{+}^{2}-\alpha_{-}^{2}})$, which is
real and positive, but non-zero whenever the $D_{i}$ are.

Figure 10: Non-Abelian E: The three eigenvalues in the 2, 3 sectors for
$\sqrt{E}=D_{1}=D_{4}$, as a function $k_{x}$ (top) and $k_{T}$ (middle), and
the region in the $k_{x}/k_{T}$-plane, where the lowest eigenvalue is negative
(bottom).

For non-zero momentum, we immediately specialize to $D_{1}=D_{4}=\sqrt{E}$,
$D_{2}=D_{3}=0$, so that we need to solve

$\displaystyle\left(-\frac{\omega^{2}}{E}+\frac{k_{T}^{2}}{E}+\frac{k_{x}^{2}}{E}\right)\times$
$\displaystyle\qquad\left[\left(-\frac{\omega^{2}}{E}+\frac{k_{T}^{2}}{E}+\frac{k_{x}^{2}}{E}\right)^{2}-4\left(\frac{\omega^{2}}{E}+\frac{k_{x}^{2}}{E}\right)-1\right]$
$\displaystyle\qquad\qquad\qquad-4\left(\frac{\omega^{2}}{E}-\frac{k_{x}^{2}}{E}\right)=0.$

The solutions are real, and one of them is negative for some regions of
$(k_{x},k_{T}^{2})$, as shown in Figure 10.

#### Non-Abelian E: 0 and 1 sector

Figure 11: Non-Abelian E: Regions of non-zero imaginary part of the
eigenvalues in the 0, 1 sector in $\alpha_{+}/\alpha_{-}$ space (shaded).
Overlaid the curve $\sqrt{E}=D_{1}/\lambda=D_{4}\lambda$ (red) and curves when
further adding non-zero $D_{3}$ (green).

We compute $\textrm{Det}[M_{00,01,10,11}]=0$ to find (compare (306)):

$\displaystyle\left[\Box_{k}^{3}+2\alpha_{-}\Box_{k}^{2}+\left(\beta_{E}+4E^{2}\right)\Box_{k}+\gamma_{e}+4E^{2}\alpha_{-}\right]^{2}=0,$
(313)

and so we may consider the expression inside the square, and treat each
solution as a double root. We again first consider the case $k_{x}=k_{T}=0$,
for which we have

$\displaystyle-(\omega^{2})^{3}+2\alpha_{+}(\omega^{2})^{2}-(\alpha_{-}^{2}+7E^{2})\omega^{2}+3E^{2}\alpha_{-}=0.$

Again, we have one real solution and two solutions that are non-real in part
of the parameter space, shown in Figure 11 as a function of $\alpha_{+}$ and
$\alpha_{-}$. Overlaid is the now-familiar curve parametrised by
$D_{1}=\sqrt{E}\lambda$, $D_{4}=\sqrt{E}/\lambda$, $D_{2}=D_{3}=0$ which lies
fully in the region of non-real eigenvalues. For $\lambda=1$, the eigenvalues
are $\omega^{2}=(0,(2\pm i\sqrt{3})E)$. The real part of two eigenvalues is
always positive, while the third has negative real part whenever $\alpha_{-}$
is negative.

Figure 12: Non-Abelian E: Real part and imaginary part for the three (double)
eigenvalues in the 0,1 sector, as a function of $k_{x}$ (top two) and $k_{T}$
(bottom two), for $k_{T,x}=0$, respectively.

For non-zero momentum and $D_{1}=D_{4}=\sqrt{E}$, we need to solve

$\displaystyle\left(-\frac{\omega^{2}}{E}+\frac{k_{T}^{2}}{E}+\frac{k_{x}^{2}}{E}\right)\times$
$\displaystyle\qquad\left[\left(-\frac{\omega^{2}}{E}+\frac{k_{T}^{2}}{E}+\frac{k_{x}^{2}}{E}\right)^{2}-4\left(\frac{\omega^{2}}{E}+\frac{k_{x}^{2}}{E}\right)+3\right]$
$\displaystyle\qquad\qquad\qquad-4\left(\frac{\omega^{2}}{E}-\frac{k_{x}^{2}}{E}\right)=0.$

which provides one real solution and two solutions that are in general non-
real. In Figure 12 we show the real and imaginary parts of these eigenvalues
as a function of $k_{x}$ (top) and $k_{T}$ (bottom), respectively. We see that
the eigenvalues are real and positive everywhere except a finite region near
the origin (an ellipse in the $k_{x}-k_{T}$-plane).

And so the non-Abelian electric field only, of the set of 12 momentum modes
present, several are potentially unstable, depending on the precise choice of
$D_{1,..,4}$. For the special case of $D_{1}=D_{4}=\sqrt{E}$, the 0,1 sector
has a finite momentum region with two unstable modes near the origin with non-
real eigenvalues, while in the 2, 3 sector there are two unstable modes in
certain momentum regions with real, negative eigenvalues (imaginary $\omega$).

###  II.6 Non-Abelian E and B

Figure 13: Non-Abelian E and B: Regions in the E-B plane, where one or more of
the first three (double) eigenvalues has negative real part (top) or non-zero
imaginary part (bottom).

Including both an electric and a magnetic field allows us to specify two
directions in space, and hence completely break isotropy. We recall the non-
Abelian/non-Abelian configuration

$\displaystyle\bar{A}^{0}_{3}=D_{1},\quad\bar{A}^{1}_{2}=D_{4},\quad\bar{A}^{3}_{3}=D_{7},$
(316)

which gives

$\displaystyle E^{1}_{1}=gD_{1}D_{4},\quad B^{2}_{1}=gD_{4}D_{7}.$ (317)

The linear equation in momentum space is again a 12-by-12 matrix, but now the
sectors are coupled. Introducing a shorthand, where we absorb $g$ as before
and define

$\displaystyle M_{\rm diag}=$
$\displaystyle\left(\begin{array}[]{ccc}\Box_{k}-D_{1}^{2}+D_{4}^{2}+D_{7}^{2}&2i\left(D_{7}k_{z}-D_{1}\omega\right)&-2iD_{4}k_{x}\\\
-2i\left(D_{7}k_{z}-D_{1}\omega\right)&\Box_{k}-D_{1}^{2}+D_{7}^{2}&0\\\
2iD_{4}k_{x}&0&\Box_{k}+D_{4}^{2}\end{array}\right),$ (321)

and

$\displaystyle M_{E}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&-2E\\\
0&2E&0\end{array}\right),\quad M_{B}=\left(\begin{array}[]{ccc}0&0&0\\\
0&0&2B\\\ 0&-2B&0\end{array}\right),$ (329)

We have schematically

$\displaystyle\left(\begin{array}[]{cccc}M_{\rm diag}&M_{E}&0&0\\\
-M_{E}^{T}&M_{\rm diag}&0&M_{B}\\\ 0&0&M_{\rm diag}&0\\\ 0&M_{B}^{T}&0&M_{\rm
diag}\\\ \end{array}\right)\left(\begin{array}[]{c}h^{0}_{1,2,3}\\\
h^{1}_{1,2,3}\\\ h^{2}_{1,2,3}\\\ h^{3}_{1,2,3}\\\ \end{array}\right)=0,$
(339)

This couples the 0,1 and 3 sectors, while the 2-sector is decoupled. Note that
this is not simply combining the non-Abelian E and the non-Abelian B that we
considered above.

We write $D_{1}=E/D_{4}$ and $D_{7}=B/D_{4}$, in which case the determinant of
the entire matrix may be conveniently written as

$\displaystyle\frac{1}{D_{4}^{16}}\left[D_{4}^{4}\Box_{k}^{3}+2D_{4}^{2}\alpha_{3}\Box_{k}^{2}+\beta_{3}\Box_{k}+\gamma_{3}\right]^{2}\times$
$\displaystyle\qquad\qquad\left[D_{4}^{4}\Box_{k}^{3}+2D_{4}^{2}\alpha_{3}\Box_{k}^{2}+(\beta_{3}-4D_{4}^{4}(\alpha_{3}-D_{4}^{4}))\Box_{k}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.+\gamma_{3}-4D_{4}^{2}\alpha_{3}(\alpha_{3}-D_{4}^{4})\right]^{2}=0,$

with

$\displaystyle\alpha_{3}$ $\displaystyle=$ $\displaystyle
B^{2}-E^{2}+D_{4}^{4},$ (342) $\displaystyle\beta_{3}$ $\displaystyle=$
$\displaystyle\alpha_{3}(\alpha_{3}+D_{4}^{4})-D_{4}^{8}-4D_{4}^{6}k_{x}^{2}-4D_{4}^{2}(Bk_{z}-E\omega)^{2},$
$\displaystyle\gamma_{3}$ $\displaystyle=$ $\displaystyle
D_{4}^{2}\left[\alpha_{3}(\alpha_{3}-D_{4}^{4})-4D_{4}^{2}(\alpha_{3}-D_{4}^{4})k_{x}^{2}\right.$
(344)
$\displaystyle\qquad\qquad\qquad\qquad\left.-4D_{4}^{2}(Bk_{z}-E\omega)^{2}\right],$

which reduces to (307) in the limit $B=0$ ($D_{7}=0$). We see that the 12
roots split up into 6 double roots, 3 for each term.

Figure 14: Non-Abelian E and B: Regions in the E-B plane, where one or more of
the second set of three (double) eigenvalues has negative real part (top) or
non-zero imaginary part (bottom).

We can again specialise to zero momentum, and we will rescale all quantities
with the appropriate power of $D_{4}$ ($\omega/D_{4}$, $B/D_{4}^{2}$,
$E/D_{4}^{2}$, $\alpha_{3}/D_{4}^{4}$ and so on, which amounts to setting
$D_{4}=1$) to find the requirements

$\displaystyle-(\omega^{2})^{3}+2(B^{2}+E^{2}+1)(\omega^{2})^{2}$
$\displaystyle\quad\quad-((B^{2}-E^{2})^{2}+3B^{2}+E^{2}+1)\omega^{2}$
$\displaystyle\quad\quad\quad\quad\quad+(B^{2}-E^{2})(B^{2}-E^{2}+1)=0,$ (345)

and

$\displaystyle-(\omega^{2})^{3}+2(B^{2}+E^{2}+1)(\omega^{2})^{2}$
$\displaystyle\quad\quad-((B^{2}-E^{2})^{2}-B^{2}+5E^{2}+1)\omega^{2}$
$\displaystyle\quad\quad\quad\quad\quad-3(B^{2}-E^{2})(B^{2}-E^{2}+1)=0,$
(346)

These provide 6 distinct (double) eigenvalues for $\omega^{2}/D_{4}^{2}$.
Figure 13 shows the regions in E-B space, where one or more has negative real
part (top) or non-zero imaginary part (bottom). Figure 14 show the same for
the second set of three eigenvalues (solutions of (II.6)). Ten of twelve
eigenmodes are subject to instability in some part of momentum space.

The solutions for non-zero momentum now involves $k_{x}$, $k_{y}$ and $k_{z}$
independently, in addition to $E$ and $B$, and it is perhaps not so helpful
attempting to display the eigenvalues in all generality. But for instance for
$k_{y}\neq 0$, we simply get the same relations but with
$\omega^{2}\rightarrow\omega^{2}-k_{y}^{2}$. Hence all solutions for
$\omega^{2}$ are shifted by $+k_{y}^{2}$, so that unstable, but real-$\omega$
solutions become stable for large $k_{y}$. For $k_{x}\neq 0$, a similar shift
appears, but there are additional terms involving $k_{x}^{2}$ as well in both
equations

$\displaystyle+4k_{x}^{2}(\omega^{2}-k_{x}^{2}-B^{2}+E^{2}).$ (347)

Something new happens for $k_{z}\neq 0$, where in addition to the shift
$\omega^{2}\rightarrow\omega^{2}-k_{z}^{2}$, we pick up terms linear in
$\omega$, $k_{z}$, and proportional to $BE$

$\displaystyle+4(B^{2}k_{z}^{2}-2BEk_{z}\omega)(\omega^{2}-k_{z}^{2}-1),$
(348)

so that the (relative) signs of $E$, $B$, $k_{z}$ matter for determining
$\omega$.

We will not pursue this further, but simply note that there is much fun to be
had in this large parameter space.

##  III Computing the quadratic and cubic contributions for non-Abelian B

###  III.1 Setting up the matrix structures: Mode equation

Having completed our survey of realisations of instabilities, we now return to
Eq. (I.1), which implies that the equation of motion for the fluctuations
$h^{\mu}_{a}$ including self-interaction at leading order takes the form

$\displaystyle\left[g^{\mu\nu}(\bar{D}_{\rho}\bar{D}^{\rho})_{ac}+2g\epsilon_{abc}\bar{F}_{b}^{\mu\nu}\right]h_{\nu,c}-g^{2}\hat{M}^{\mu\nu}_{ab}(C)h_{\nu,b}=0,$

where the matrix $\hat{M}^{\mu\nu}_{ab}$ is expressed in terms of the
expectation values of the fluctuation fields themselves

$\displaystyle C^{\mu,\nu}_{a,b}=\langle
h^{\mu}_{a}(x,t)h^{\nu}_{b}(x,t)\rangle=\int\frac{d^{3}k}{(2\pi)^{3}}\langle
h^{\mu}_{a}(k,t)(h^{\nu}_{b}(k,t))^{\dagger}\rangle.$

For the case of a non-Abelian B-field, the matrix structure is as follows (in
the 2-3 sector)

$\displaystyle\hat{M}^{i}_{j,ab}=\left[\begin{array}[]{cccccc}-\hat{M}_{11}&0&0&0&0&0\\\
0&-\hat{M}_{22}&-\hat{M}_{23}&0&0&-\hat{M}_{26}\\\
0&-\hat{M}_{32}&-\hat{M}_{33}&0&-\hat{M}_{35}&0\\\ 0&0&0&-\hat{M}_{44}&0&0\\\
0&0&-\hat{M}_{53}&0&-\hat{M}_{55}&-\hat{M}_{56}\\\
0&-\hat{M}_{62}&0&0&-\hat{M}_{65}&-\hat{M}_{66}\end{array}\right],$ (357)
(358)

where

$\displaystyle\hat{M}_{11}=C^{3,3}_{2,2}+C^{3,3}_{3,3},\quad\hat{M}_{22}=C^{3,3}_{3,3}+C^{3,3}_{1,1},$
(359)
$\displaystyle\hat{M}_{33}=C^{3,3}_{1,1}+C^{3,3}_{2,2},\quad\hat{M}_{44}=C^{2,2}_{2,2}+C^{2,2}_{3,3},$
(360)
$\displaystyle\hat{M}_{55}=C^{2,2}_{3,3}+C^{2,2}_{1,1},\quad\hat{M}_{66}=C^{2,2}_{1,1}+C^{2,2}_{2,2},$
(361)

and

$\displaystyle\hat{M}_{23}=\hat{M}_{32}$
$\displaystyle=-C^{3,3}_{3,2}=-C^{3,3}_{2,3},$ (362)
$\displaystyle\hat{M}_{56}=\hat{M}_{65}$
$\displaystyle=-C^{2,2}_{3,2}=-C^{2,2}_{2,3},$ (363)
$\displaystyle\hat{M}_{26}=\hat{M}_{62}$
$\displaystyle=2C^{2,3}_{2,3}-C^{2,3}_{3,2},$ (364)
$\displaystyle\hat{M}_{35}=\hat{M}_{53}$
$\displaystyle=2C^{2,3}_{3,2}-C^{2,3}_{2,3}.$ (365)

We have at this point already implemented that several correlators vanish
identically, and we have for simplicity by hand set to zero all correlators
involving the modes in the 0-1 sector, since they are stable. The 0-1 sector
of the matrix has the same structure as the 2-3 sector, but without the
sector-mixing components $\hat{M}_{26,62,35,53}$. In our particular setup, we
in addition find that $\hat{M}_{11}=\hat{M}_{44}$,
$\hat{M}_{22}=\hat{M}_{33}=\hat{M}_{55}=\hat{M}_{66}$,
$\hat{M}_{23}=\hat{M}_{32}=\hat{M}_{56}=\hat{M}_{65}$ and
$\hat{M}_{35}=\hat{M}_{53}=-\hat{M}_{26}=-\hat{M}_{62}$, as there are only
four distinct correlators to take into account,
$C^{2,2}_{1,1},C^{2,2}_{2,2},C^{2,3}_{2,3},C^{2,2}_{2,3}$.

Our task is then to compute these time-dependent correlators, and solving
again the eigenvalue equations ($\textrm{Det}[\mathcal{M}-g^{2}\hat{M}]=0$) to
find corrected, time-dependent values of $\omega^{2}$. We will do this below,
but already at this stage, we note that in the limit where
$\hat{M}\gg\mathcal{M}$, the corrected system now has three double eigenvalues
$\hat{M}_{11},\hat{M}_{22}\pm\sqrt{\hat{M}_{23}^{2}+\hat{M}_{26}^{2}}$, of
which one is negative whenever
$\hat{M}_{23}^{2}+\hat{M}_{26}^{2}>\hat{M}_{22}^{2}$. This means that (some
of) the instabilities remain even at late times, provided the background field
remains constant. In the 0-1 sector, all modes remain stable, since
$\hat{M}_{26}\rightarrow 0$ and as we will see, $\hat{M}_{22}>\hat{M}_{23}$.

###  III.2 Setting up the matrix structures: Background equation

Similarly, we can express the quadratic contribution (10) in terms of
correlators,

$\displaystyle\bar{D}^{ab}_{\mu}\bar{F}^{\mu\nu,b}=j^{\nu}_{a}-\bar{M}^{\nu\mu}_{ab}[C]\bar{A}_{\mu,b}.$
(366)

Again, in our specific setup the space-derivative terms in the correction to
the current $\bar{j}^{\nu}_{a}$ vanishes because the correlators are space-
independent. The time-derivative terms also vanish because they either act on
a stable correlator (in the 0-$\mu$ sector, which we neglect) or because of
antisymmetrization. We are therefore left with the matrix $\bar{M}(C)$ which
multiplies $\bar{A}$, and we find that this matrix has the same structure as
$\hat{M}$. Since there is no original current in the 0-1 sector, it is
consistent to set $\bar{A}^{0,1}_{a}=0$ throughout, and we will consider it no
further. In the 2-3 sector, however, we would have to solve the full Yang-
Mills equation with a time-dependent current.

We will not attempt this here, but simply estimate at which time the
assumption of constant $\bar{A}$ is likely to break down. We quantify this as
when the correction is larger than the original current

$\displaystyle\bar{M}^{\nu\mu}_{ab}[C]\bar{A}_{\mu,b}>j^{\nu}_{a}.$ (367)

Recall, that

$\displaystyle j^{\nu}_{a}=-gB\left[\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ D_{7}\\\
-D_{8}\end{array}\right)\left(\begin{array}[]{c}0\\\ -D_{5}\\\
D_{6}\end{array}\right)\right].$ (380)

with $D_{6}=-D_{7}=\sqrt{B/g}$, and $D_{5}=D_{8}=0$. The correction simply
becomes

$\displaystyle\left(\bar{M}\bar{A}\right)^{\nu}_{a}=-gB\left[\left(\begin{array}[]{c}0\\\
0\\\ 0\end{array}\right)\left(\begin{array}[]{c}0\\\ 0\\\
0\end{array}\right)\left(\begin{array}[]{c}0\\\ J_{7}\\\
-J_{8}\end{array}\right)\left(\begin{array}[]{c}0\\\ -J_{5}\\\
J_{6}\end{array}\right)\right].$ (394)

with

$\displaystyle
J_{6}=-J_{7}=\sqrt{\frac{B}{g}}\frac{g\left(C^{33}_{11}+C^{33}_{33}-3C^{23}_{23}\right)}{B},$
$\displaystyle J_{5}=-J_{8}=\sqrt{\frac{B}{g}}\frac{g\,C^{22}_{23}}{B},$ (396)

so that the form of the total current is unchanged, and will still generate a
$B$ field in the $\mu=1$, $a=1$ direction. This $B$-field is however now time
dependent, and the time-dependence of $\bar{A}^{\mu}_{a}$ could also
potentially source an Abelian E-field. This means that as soon as either of
the correlators is larger than $B/g$, our approximation breaks down. We will
now estimate when that happens.

###  III.3 Eigenvectors, quantization and correlators

Figure 15: Non-Abelian B: The correlators $C^{22}_{11}$, $C^{22}_{22}$,
$C^{22}_{23}$ and $C^{23}_{23}$ (top to bottom) as a function of $k_{y}$ for
time $gBt=1,2,3$ (red, blue, green). Renormalised (full lines) and
unrenormalised (dotted lines).

The 2-3 sector of the linear system has 6 eigenvalues for each set of momenta
$k_{x,y,z}$, $\omega_{i}^{2}(k)$, corresponding to 6 eigenvectors
$\phi_{i}(k,t)$ which are solutions to

$\displaystyle\partial_{t}^{2}\phi_{i}(k,t)=-\omega_{i}^{2}(k)\phi_{i}(k,t),\quad
i=1,...,6,$ (397) $\displaystyle\rightarrow$
$\displaystyle\phi_{i}(k,t)=\alpha_{k}^{i}e^{i\omega_{i}(k)t}+\beta_{k}^{i}e^{-i\omega_{i}(k)t}.$
(398)

In order to define the state in which to compute the expectation values, we
will assume that the system is prepared in a vanishing background
$\bar{A}^{\mu}_{a}=0$ for $t<0$, and that $\bar{A}^{\mu}_{a}$ is
instantaneously turned on at $t=0$, remaining constant for $t>0$. Then we may
write for $t<0$

$\displaystyle\phi_{i}(k,t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega^{0}_{k}}}\left(a_{k}^{i}e^{i\omega_{k}^{0}t}+(a_{-k}^{i})^{\dagger}e^{-i\omega_{k}^{0}t}\right),$
(399) $\displaystyle\partial_{t}\phi_{i}(k,t)$ $\displaystyle=$
$\displaystyle\frac{i\omega_{k}^{0}}{\sqrt{2\omega^{0}_{k}}}\left(a_{k}^{i}e^{i\omega_{k}^{0}t}-(a_{-k}^{i})^{\dagger}e^{-i\omega_{k}^{0}t}\right),$
(400)

where we note that for $\bar{A}^{\mu}_{a}=0$, all the eigenvalues coincide
$\omega_{i}(k)\rightarrow\omega^{0}_{k}=|{\bf k}|$.

For $t>0$, we may write in general

$\displaystyle\phi_{i}(k,t)$ $\displaystyle=$
$\displaystyle\alpha_{k}^{i}e^{i\omega_{k}^{i}t}+\beta_{k}^{i}e^{-i\omega_{k}^{i}t},$
(401) $\displaystyle\partial_{t}\phi_{i}(k,t)$ $\displaystyle=$ $\displaystyle
i\omega_{k}^{i}\left(\alpha_{k}^{i}e^{i\omega_{k}^{i}t}-\beta_{k}^{i}e^{-i\omega_{k}^{i}t}\right).$
(402)

Matching the expressions at $t=0$, we can express $\alpha_{k}^{i}$ and
$\beta_{k}^{i}$ in terms of $a_{k}^{i}$ and $a_{-k}^{i}$ as,

$\displaystyle\alpha_{k}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega^{0}_{k}}}\left[a_{k}^{i}\left(1+\frac{\omega_{k}^{0}}{\omega_{k}^{i}}\right)+(a_{-k}^{i})^{\dagger}\left(1-\frac{\omega_{k}^{0}}{\omega_{k}^{i}}\right)\right],$
$\displaystyle\beta_{k}^{i}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega^{0}_{k}}}\left[a_{k}^{i}\left(1-\frac{\omega_{k}^{0}}{\omega_{k}^{i}}\right)+(a_{-k}^{i})^{\dagger}\left(1+\frac{\omega_{k}^{0}}{\omega_{k}^{i}}\right)\right].$

Choosing our initial state as the vacuum at $t<0$, $a_{k}^{i}|0\rangle=0$ and
using the standard commutation relations
$[a_{k}^{i},(a^{j}_{k^{\prime}})^{\dagger}]=(2\pi)^{3}\delta^{ij}\delta^{3}(k-k^{\prime})$,
we can straightforwardly compute expectation values of any combination of
$\alpha_{k}^{i}$, $\beta_{k}^{i}$.

The fluctuation fields $h^{\mu}_{a}(k)$ in the 2-3 sector are linear
combinations of the eigenmodes $\phi_{i}(k)$, written in terms of a unitary
matrix $U^{ij}(k)$. With ($\mu=2,3$, $a=1,2,3$ $\rightarrow$ $1,..,6$) so that

$\displaystyle\langle h^{i}(k)(h^{j}(k))^{\dagger}\rangle$
$\displaystyle=\langle U_{ir}\phi^{r}(k)[U_{js}\phi^{s}(k)]^{\dagger}\rangle$
(405)
$\displaystyle=U^{ir}[U^{\dagger}]_{rj}\langle\phi^{r}(k)[\phi^{r}(k)]^{\dagger}\rangle,$

since $\langle\phi^{r}(k)[\phi^{s}(k)]^{\dagger}\rangle\propto\delta^{rs}$. We
wish to compute the correlators

$\displaystyle C^{\mu,\nu}_{a,b}$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}\langle
h^{\mu}_{a}(k)h^{\nu}_{b}(-k)\rangle$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}\frac{\langle
h^{i}(k)[h^{j}(k)]^{\dagger}+[h^{i}(k)]^{\dagger}h^{j}(k)\rangle}{2}.$

Which follow from combining (398, LABEL:eq:matching, 405), and using these we
can construct the non-linear correction to the linear mode equation (358).

One final complication is that the diagonal correlators are (IR and UV)
divergent, and so in order to gauge the effect of the instability, we will in
practice compute

$\displaystyle C^{\mu,\nu}_{a,b;ren}=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}\bigg{[}\langle
h^{\mu}_{a}(k,t)h^{\nu}_{b}(-k,t)\rangle-\langle
h^{\mu}_{a}(k,0)h^{\nu}_{b}(-k,0)\rangle\bigg{]}.$

In Figure 15 we show the four distinct (renormalised) correlators in momentum
space at times $t=1,2,3$ (in units where $gB=1$). We see that the
$C^{22}_{22}$ and $C^{23}_{23}$ are substantially larger than the other two.

###  III.4 Corrected dispersion relations

In Figure 16, we show $M_{11}$, $M_{22}\pm\sqrt{M_{23}^{2}+M_{26}^{2}}$, which
was our late-time estimate for the eigenvalues in section III.1. The lowest
eigenvalue is indeed negative throughout, and so for any momentum, for late
enough time, instabilities remain if they are not switched off for other
reasons.

Figure 16: Non-Abelian B: The three eigenvalues $M_{11}$,
$M_{22}\pm\sqrt{M_{2}3^{2}-M_{26}^{2}}$. Figure 17: Non-Abelian B: The
corrected eigenvalues $\omega_{i}$ as a function of time.

We can do better, and simply compute the eigenvalues at $k=0$ of the whole
matrix $\mathcal{M}-g^{2}\hat{M}$, as a function of time. Figure 17 shows the
6 eigenvalues in the 2-3 sector as a function of time. The initially unstable
mode is briefly stabilized around $gBt=3.5$, but by $4.7$ the late time
behaviour takes over and the eigenvalue becomes negative again.

###  III.5 Corrected current

Figure 18: Non-Abelian B: The quantities $J_{5,6,7,8}/\sqrt{B/g}$ in time.

In a similar way, we can insert the correlators into the effective current.
Figure 18 shows the quantities $J_{5,6,7,8}$ (396) in time. We see that
$J_{6}=-J_{7}$ grow fastest, and become order 1 (our criterion for $g=B=1$)
around $gBt=4-5$. At this time, we can then no longer assume a constant
background field. Whether this leads to stability is beyond our approach here,
but may be readily studied using numerical simulations.

##  IV Conclusions

To summarize, we have investigated the instabilities of field fluctuations in
pure SU(2) classical gauge theory with a constant homogeneous background E
and/or B field. We first in section II presented a fairly general set of
background vector potentials, generalising the setup of bazak1 ; bazak2 . We
then in section A focused on non-Abelian E and B field backgrounds, and
derived the dispersion relation that the fluctuations obey. Non-real instances
of frequencies $\omega(k)$ signal unstable modes in the linear theory, and we
confirmed that these arise for homogeneous E fields, B fields and when
combining the two. As an aside, we also found that the instability does not
simply depend on the B or E-field, but on the vector field components
$\bar{A}^{\mu}_{a}$ separately.

Of particular interest is to understand to what extent these instabilities
shut themselves off, due to the non-linearities of the theory. To that end, we
computed the non-linear (Gaussian) corrections to both the background field
equation and the linear fluctuation equation. We found that the corrections
certainly influence the instabilities of the fluctuations, but that rather
than turning them off, they can make them stronger. This however assumes that
the background field stays constant and homogeneous throughout, which is
spoiled by the non-linearities in the background field equation. We estimated
the time when the background field is no longer a constant non-Abelian B field
to be $gBt\simeq 4$, but further details of how this comes about, and the
explicit effect on the instabilities is likely better adressed numerically.

Although the context of this type of instabilities is the thermalisation of
the plasma in HICs, anisotropic gauge fields are not unique to QCD.
Anisotropic conditions also appear in cosmology, for instance in the context
of first order phase transitions. In that case, the walls of nucleated bubbles
sweep through the ambient plasma, so that immediately in front and behind the
advancing wall, the plasma is anisotropic (see for instance tranberg ; GW1 ;
GW2 ). Although it is not clear to what extent net E and B fields are
produced, it would be worth investigating whether an anisotropic out-of-
equilibrium background could potentially source similar instabilities.

The obvious continuation of this work is to make a full investigation the
other realisations of background E and B fields presented in Appendix B. This
could also readily be compared to direct numerical simulations of the pure
gauge SU(2) system, going beyond the Gaussian approximation to the non-
linearities. Further generalising to SU(3) gauge theory involves writing down
a $\mathcal{M}(32\times 32)$ matrix equation and solving for the dispersion
relations and then unstable modes there. This is ongoing work, but we note
that the results for SU(2) presented here naturally arise as a subsector in
the SU(3) theory.

## Acknowledgments

We thank Aleksi Kurkela and Germano Nardini for enlightening discussions and
Alexander Rothkopf for many useful comments. This work is funded by the grant
Project. No. 302640 funded by the Research Council of Norway.

##  Appendix A Linear equations of motion in momentum space

If the background field $\bar{A}^{\mu}_{a}$ has no space dependence, we can
simply Fourier transform in three dimensions and write

$\displaystyle
h^{\mu}_{a}(x,y,z,t)=\int\frac{dk_{x}dk_{y}dk_{z}}{(2\pi)^{3}}e^{-ik_{x}x-ik_{y}y-ik_{z}z}h^{\mu}_{a}(k_{x},k_{y},k_{z},t).$

Then the equation of motion in momentum space reads

$\displaystyle\mathcal{M}(\partial_{t}^{2},\partial_{t},k_{x},k_{y},k_{z},\bar{A})h^{\mu}_{a}=0,$
(409)

which can be solved as a set of coupled linear differential equations. One may
further Fourier transform in time, in the sense that

$\displaystyle
h^{\mu}_{a}(k_{x},k_{y},k_{z},t)=\int\frac{d\omega}{2\pi}h^{\mu}_{a}(k_{x},k_{y},k_{z},\omega)e^{i\omega
t},$ (410)

in which case one may substitute $\partial_{t}\rightarrow i\omega$ in (409).
The expectation is that for some of the modes (for certain values of
$k_{x},k_{y},k_{z},\bar{A}$) $\omega_{k}$ has a non-zero imaginary part, and
the mode grows (or decays) exponentially in time.

If the background field depends on a spatial coordinate, we can only Fourier
transform in the remaining coordinates (say $y,z$), to find

$\displaystyle
h^{\mu}_{a}(x,y,z,t)=\int\frac{dk_{y}dk_{z}}{(2\pi)^{2}}e^{-ik_{y}y-ik_{z}z}h^{\mu}_{a}(x,k_{y},k_{z},t),$

and we are led to a relation of the form

$\displaystyle\mathcal{M}(\partial_{t}^{2},\partial_{t},\partial_{x}^{2},\partial_{x},k_{y},k_{z},\bar{A})h^{\mu}_{a}=0.$
(412)

This implies that the field profile is not a product of plane waves in $x,y,z$
but that the $x$-dependence is some other function to be solved for
explicitly. Similary, if the background field depends on two spatial
coordinate, we can only Fourier transform in the remaining coordinates (say,
$z$), to find

$\displaystyle
h^{\mu}_{a}(x,y,z,t)=\int\frac{dk_{z}}{(2\pi)}e^{-ik_{z}z}h^{\mu}_{a}(x,y,k_{z},t).$
(413)

Neither the $x$\- and $y$-dependence can then be expressed as plane waves.
Finally, if $\bar{A}^{\mu}_{a}$ is a function of $t$, also the Fourier
transform in time fails, and the field has a non-trivial time evolution. We
again refer to bazak1  and Appendix B for a discussion on this.

##  Appendix B Space- or time-dependent vector potentials

For completeness, we now briefly present a few cases of Abelian vector
potentials, and how one might go about solving the linearised equations of
motion.

###  B.1 Abelian Configuration of B field ($C_{3}$ $=$ $B$ $\neq 0$)

Consider a potential $\bar{A}^{\mu}_{a}$ in the $\mu=3$ direction with a
constant homogeneous magnetic field in the $i=2$ direction and in the $a=1$
color direction,

$\bar{A}^{\mu}_{a}=\begin{bmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}By\\\ 0\\\ 0\end{pmatrix}\end{bmatrix}.$ (414)

The linearized Y-M equation becomes

$\displaystyle\Box
h^{\mu}_{a}-2gBy\epsilon^{ab1}\partial_{z}h^{\mu}_{b}-2gB\epsilon^{ab1}(\delta^{\mu
y}h^{z}_{b}-\delta^{\mu z}h^{y}_{b})$
$\displaystyle-g^{2}B^{2}y^{2}\epsilon^{ac1}\epsilon^{cb1}h^{\mu}_{b}=0.$

Let us define the functions

$\displaystyle T^{\pm}$ $\displaystyle=h^{0}_{2}\pm ih^{0}_{3},$
$\displaystyle X^{\pm}$ $\displaystyle=h^{x}_{2}\pm ih^{x}_{3},$
$\displaystyle Y^{\pm}$ $\displaystyle=h^{y}_{2}\pm ih^{y}_{3},$
$\displaystyle Z^{\pm}$ $\displaystyle=h^{z}_{2}\pm ih^{z}_{3}.$

Modifying Eq. (B.1) and then solving for each modes in Fourier space, we
obtain the dispersion relations such as

$\left(-\omega^{2}+k_{x}^{2}+(k_{z}\mp
gBy)^{2}-\frac{d^{2}}{dy^{2}}\right)T^{\pm}(y)=0,$ (416)
$\left(-\omega^{2}+k_{x}^{2}+(k_{z}\mp
gBy)^{2}-\frac{d^{2}}{dy^{2}}\right)X^{\pm}(y)=0,$ (417) $\left(-\omega^{2}\pm
2gB+k_{x}^{2}+(k_{z}-gBy)^{2}-\frac{d^{2}}{dy^{2}}\right)U^{\pm}(y)=0,$ (418)
$\left(-\omega^{2}\mp
2gB+k_{x}^{2}+(k_{z}+gBy)^{2}-\frac{d^{2}}{dy^{2}}\right)W^{\pm}(y)=0.$ (419)

where

$\displaystyle U^{\pm}$ $\displaystyle\equiv Y^{+}\pm iZ^{+},$ $\displaystyle
W^{\pm}$ $\displaystyle\equiv Y^{-}\pm iZ^{-}.$

The above equations resemble a Schrödinger equation for a harmonic oscillator,
and have a discrete spectrum of frequencies. We can solve for the frequencies
for $W^{+}$ and $U^{-}$, to find

$\omega_{0}^{2}=2gB\left(n+\frac{1}{2}\right)+k_{x}^{2}.$ (420)

and

$\omega_{\pm}^{2}=2gB\left(n+\frac{1}{2}\right)\pm 2gB+k_{x}^{2}.$ (421)

for discrete $n=0,1,2,\ldots$ and continuous $k_{x}$.

It is clear that $\omega_{0}^{2}$ and $\omega_{+}^{2}$ are always positive
$(\geq 0)$ for any values of $n$, but $\omega_{-}^{2}$ is negative whenever
$n<\frac{1}{2}-\frac{k_{x}^{2}}{2gB}$. That means there are unstable modes of
$U^{-}$ and $W^{+}$ which grow exponentially. The remaining modes $T^{\pm}$
and $X^{\pm}$, $U^{+}$ and $W^{-}$ are stable for all $k_{x}$. We see that the
dependence of $k_{z}$ is absorbed in a redefinition of the $y$-coordinate.

###  B.2 Abelian Configuration of E field ($C_{1}$ $=$ $E$ $\neq 0$)

Consider instead a potential $\bar{A}^{\mu}_{a}$ in the $\mu=0$ direction with
a constant homogeneous electric field in the $i=1$ direction and in the $a=1$
color direction

$\bar{A}^{\mu}_{a}=\begin{bmatrix}\begin{pmatrix}-Ex\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix}\end{bmatrix}.$ (422)

The linearized Y-M equation becomes

$\displaystyle
g^{\mu\nu}\left[(\partial_{\rho}\delta^{ad}-gf^{ade}\bar{A}^{e}_{\rho})(\partial^{\rho}\delta_{dc}-gf_{dcf}\bar{A}^{\rho}_{f})h^{c}_{\nu}\right]$
$\displaystyle+2gf^{abc}\bar{F}^{\mu\nu}_{b}h^{c}_{\nu}=0.$

Solving Eq. (B.2) for the modes in Fourier space give us dispersion relations
such as

$\left(k_{y}^{2}+k_{z}^{2}\pm
2igE-(\omega+gEx)^{2}-\frac{d^{2}}{dx^{2}}\right)G^{\pm}(x)=0,$ (424)
$\left(k_{y}^{2}+k_{z}^{2}\mp 2igE-(\omega-
gEx)^{2}-\frac{d^{2}}{dx^{2}}\right)H^{\pm}(x)=0,$ (425)
$\left(k_{y}^{2}+k_{z}^{2}-(\omega\pm
gEx)^{2}-\frac{d^{2}}{dx^{2}}\right)Y^{\pm}(x)=0,$ (426)
$\left(k_{y}^{2}+k_{z}^{2}-(\omega\mp
gEx)^{2}-\frac{d^{2}}{dx^{2}}\right)Z^{\pm}(x)=0.$ (427)

where the functions $G^{\pm}$ and $H^{\pm}$ are defined as

$\displaystyle G^{\pm}$ $\displaystyle\equiv T^{+}\pm X^{+},$ $\displaystyle
H^{\pm}$ $\displaystyle\equiv T^{-}\pm X^{-}.$

The above equations again resemble an inverted Schrödinger equation, an
”upside-down” harmonic oscillator. This does not have any normalizable
solutions.

###  B.3 Abelian E and Abelian B ($C_{1}$ $\neq 0$ $\neq$ $C_{3}$)

Consider a potential $\bar{A}^{\mu}_{a}$ corresponding to a constant
homogeneous electric field in the $x$ direction and a constant homogeneous
magnetic field in the $y$ direction,

$\bar{A}^{\mu}_{a}=\begin{bmatrix}\begin{pmatrix}-Ex\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\ 0\end{pmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}By\\\ 0\\\ 0\end{pmatrix}\end{bmatrix}.$ (428)

The linearized Y-M equation in this case becomes

$\displaystyle\left((k_{z}-gBy)^{2}-(\omega+gEx)^{2}\pm
2igE-\frac{d^{2}}{dx^{2}}-\frac{d^{2}}{dy^{2}}\right)$ $\displaystyle\times
G^{\pm}(x,y)=0,$ $\displaystyle\left((k_{z}+gBy)^{2}-(\omega-gEx)^{2}\mp
2igE-\frac{d^{2}}{dx^{2}}-\frac{d^{2}}{dy^{2}}\right)$ $\displaystyle\times
H^{\pm}(x,y)=0,$ $\displaystyle\left((k_{z}-gBy)^{2}-(\omega+gEx)^{2}\pm
2gB-\frac{d^{2}}{dx^{2}}-\frac{d^{2}}{dy^{2}}\right)$ $\displaystyle\times
U^{\pm}(x,y)=0,$ $\displaystyle\left((k_{z}+gBy)^{2}-(\omega-gEx)^{2}\mp
2gB-\frac{d^{2}}{dx^{2}}-\frac{d^{2}}{dy^{2}}\right)$ $\displaystyle\times
W^{\pm}(x,y)=0.$

Since there is explicit dependence on $x$ and $y$ we are left with a partial
differential equation in both coordinates. The equations (B.3) $-$ (B.3) can
be further solved variable separation. For example, considering Eq. (B.3), we
can write

$U^{\pm}(x,y)=U^{\pm}_{E}(x)U^{\pm}_{B}(y).$ (433)

We can split them into two equations such as

$\left[Q_{U}^{\pm}-\frac{d^{2}}{dx^{2}}-g^{2}E^{2}\left(\frac{\omega}{gE}+x\right)^{2}\right]U^{\pm}_{E}(x)=0,$
(434) $\left[-Q_{U}^{\pm}\pm
2gB-\frac{d^{2}}{dy^{2}}+g^{2}B^{2}\left(\frac{k_{z}}{gB}-y\right)^{2}\right]U^{\pm}_{B}(y)=0.$
(435)

where $Q_{U}^{\pm}$ is the separation constant. Equation (435) is another
Schrödinger equation with the following replacements, $Q_{U}^{\pm}\mp 2gB$
$\rightarrow$ $2m\varepsilon$, $gB$ $\rightarrow$ $m\bar{\omega}$,
$\frac{k_{z}}{gB}$ $\rightarrow$ $y_{0}$. Since the oscillator energy
$\varepsilon_{n}$ $=$ $(n+\frac{1}{2})\bar{\omega}$, $n$ $=$ $0,1,2,3,...$

$Q_{U}^{\pm}=gB(2n+1)\pm 2gB.$ (436)

Equation (434) is an inverted Schrödinger equation with the following
replacements, $Q_{U}^{\pm}$ $\rightarrow$ $-2m\varepsilon$, $gE$ $\rightarrow$
$m\bar{\omega}$, $\frac{\omega}{gE}$ $\rightarrow$ $x_{0}$. It has no
normalizable solutions.

Similary, separating variables one can write

$W^{\pm}(x,y)=W^{\pm}_{E}(x)W^{\pm}_{B}(y).$ (437)

It can also be split into two equations such as

$\left[Q_{W}^{\pm}-\frac{d^{2}}{dx^{2}}-g^{2}E^{2}\left(\frac{\omega}{gE}-x\right)^{2}\right]W^{\pm}_{E}(x)=0,$
(438) $\left[-Q_{W}^{\pm}\mp
2gB-\frac{d^{2}}{dy^{2}}+g^{2}B^{2}\left(\frac{k_{z}}{gB}+y\right)^{2}\right]W^{\pm}_{B}(y)=0.$
(439)

where $Q_{W}^{\pm}$ is the separation constant. Similarly one can solve the
remaining two equations.

###  B.4 Abelian E and non-Abelian B ($C_{2}$ $=$ $E$, $D_{3}$ $=$
$\lambda\sqrt{B/g}$ and $D_{12}$ $=$ $\frac{1}{\lambda}\sqrt{B/g}$)

Consider a potential $\bar{A}^{\mu}_{a}$ corresponding to a constant
homogeneous electric field in the $x$ direction and constant homogenous
magnetic field components in the $x$ and $z$ directions

$\bar{A}^{\mu}_{a}=\begin{bmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}-tC_{2}\\\ 0\\\
D_{3}\end{pmatrix}\begin{pmatrix}0\\\ 0\\\
0\end{pmatrix}\begin{pmatrix}D_{12}\\\ 0\\\ 0\end{pmatrix}\end{bmatrix}.$
(440)

In this case, the linearized Yang-Mills equation reduces to

$\displaystyle\Box
h^{\mu}_{a}+2g\epsilon^{a1b}(D_{12}\partial_{z}-tC_{2}\partial_{x})h^{\mu}_{b}$
$\displaystyle\quad\quad+2g\epsilon^{a3b}D_{3}\partial_{x}h^{\mu}_{b}-g^{2}\epsilon^{a1d}\epsilon^{d1b}(t^{2}C_{2}^{2}+D_{12}^{2})h^{\mu}_{b}$
$\displaystyle\quad\quad+g^{2}tC_{2}D_{3}h^{\mu}_{b}(\epsilon^{a1d}\epsilon^{d3b}+\epsilon^{a3d}\epsilon^{d1b})-g^{2}\epsilon^{a3d}\epsilon^{d3b}D_{3}^{2}h^{\mu}_{b}$
$\displaystyle\quad\quad+2g^{2}\epsilon^{a2b}(\delta^{\mu
3}D_{12}D_{3}h^{x}_{b}-\delta^{\mu 1}D_{3}D_{12}h^{z}_{b})=0.$ (441)

We may now Fourier transform in space, but not in time, and the linear
equation for the modes of the fluctuations reduce to a set of 12 coupled
linear differential equations which we can split into four $3\times 3$
matrices corresponding to the $\mu=0,1,2,3$ sectors introduced earlier, with
cross-terms coupling the 0 and 1 sectors. The diagonal $3\times 3$ matrix has
the form

$M^{t}_{ANA}=\begin{bmatrix}(\Box+g^{2}D_{3}^{2})&2igk_{x}D_{3}&g^{2}tD_{3}C_{2}\\\
-2igk_{x}D_{3}&\Box+g^{2}(t^{2}C_{2}^{2}+D_{12}^{2}+D_{3}^{2})&2ig(D_{12}k_{z}-tC_{2}k_{x})\\\
g^{2}tD_{3}C_{2}&-2ig(D_{12}k_{z}-tC_{2}k_{x})&\Box+g^{2}(t^{2}C_{2}^{2}+D_{12}^{2})\end{bmatrix}$
(442)

where $\Box$ $=$ $\frac{\partial^{2}}{\partial t^{2}}+\textbf{k}^{2}$,
$\textbf{k}^{2}$ $=$ $k_{x}^{2}+k_{y}^{2}+k_{z}^{2}$. Diagonalizing this set
of equations is non-trivial, but would result in a set of 12 non-linear
differential equations in $t$, which most likely require numerical evaluation.

## References

  * (1) S. Mrowczynski, B. Schenke and M. Strickland, Phys. Rept. 682 (2017), 1-97 doi:10.1016/j.physrep.2017.03.003 [arXiv:1603.08946 [hep-ph]].
  * (2) T. Epelbaum and F. Gelis, Phys. Rev. Lett. 111 (2013), 232301 doi:10.1103/PhysRevLett.111.232301 [arXiv:1307.2214 [hep-ph]].
  * (3) T. Epelbaum and F. Gelis, Phys. Rev. D 88 (2013), 085015 doi:10.1103/PhysRevD.88.085015 [arXiv:1307.1765 [hep-ph]].
  * (4) T. Lappi, Phys. Rev. C 67 (2003), 054903 doi:10.1103/PhysRevC.67.054903 [arXiv:hep-ph/0303076 [hep-ph]].
  * (5) T. Lappi and L. McLerran, Nucl. Phys. A 772 (2006), 200-212 doi:10.1016/j.nuclphysa.2006.04.001 [arXiv:hep-ph/0602189 [hep-ph]].
  * (6) F. Gelis, K. Kajantie and T. Lappi, Phys. Rev. Lett. 96 (2006), 032304 doi:10.1103/PhysRevLett.96.032304 [arXiv:hep-ph/0508229 [hep-ph]].
  * (7) R. Baier, A. H. Mueller, D. Schiff and D. T. Son, Phys. Lett. B 502 (2001), 51-58 doi:10.1016/S0370-2693(01)00191-5 [arXiv:hep-ph/0009237 [hep-ph]].
  * (8) S. Mrowczynski, Phys. Lett. B 314 (1993), 118-121 doi:10.1016/0370-2693(93)91330-P
  * (9) S. Mrowczynski, Acta Phys. Polon. B 37 (2006), 427-454 [arXiv:hep-ph/0511052 [hep-ph]].
  * (10) A. Kurkela and G. D. Moore, JHEP 12 (2011), 044 doi:10.1007/JHEP12(2011)044 [arXiv:1107.5050 [hep-ph]].
  * (11) A. Kurkela and G. D. Moore, JHEP 11 (2011), 120 doi:10.1007/JHEP11(2011)120 [arXiv:1108.4684 [hep-ph]].
  * (12) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, Phys. Rev. D 89 (2014) no.7, 074011 doi:10.1103/PhysRevD.89.074011 [arXiv:1303.5650 [hep-ph]].
  * (13) A. Kurkela, A. Mazeliauskas, J. F. Paquet, S. Schlichting and D. Teaney, Phys. Rev. Lett. 122 (2019) no.12, 122302 doi:10.1103/PhysRevLett.122.122302 [arXiv:1805.01604 [hep-ph]].
  * (14) A. Kurkela, A. Mazeliauskas, J. F. Paquet, S. Schlichting and D. Teaney, Phys. Rev. C 99 (2019) no.3, 034910 doi:10.1103/PhysRevC.99.034910 [arXiv:1805.00961 [hep-ph]].
  * (15) A. Kurkela and Y. Zhu, Phys. Rev. Lett. 115 (2015) no.18, 182301 doi:10.1103/PhysRevLett.115.182301 [arXiv:1506.06647 [hep-ph]].
  * (16) J. Berges and S. Schlichting, Phys. Rev. D 87 (2013) no.1, 014026 doi:10.1103/PhysRevD.87.014026 [arXiv:1209.0817 [hep-ph]].
  * (17) J. Berges, S. Scheffler, S. Schlichting and D. Sexty, Phys. Rev. D 85 (2012), 034507 doi:10.1103/PhysRevD.85.034507 [arXiv:1111.2751 [hep-ph]].
  * (18) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, Phys. Rev. D 89 (2014) no.7, 074011 doi:10.1103/PhysRevD.89.074011 [arXiv:1303.5650 [hep-ph]].
  * (19) J. Randrup and S. Mrowczynski, Phys. Rev. C 68 (2003), 034909 doi:10.1103/PhysRevC.68.034909 [arXiv:nucl-th/0303021 [nucl-th]].
  * (20) P. B. Arnold and J. Lenaghan, Phys. Rev. D 70 (2004), 114007 doi:10.1103/PhysRevD.70.114007 [arXiv:hep-ph/0408052 [hep-ph]].
  * (21) P. B. Arnold, J. Lenaghan, G. D. Moore and L. G. Yaffe, Phys. Rev. Lett. 94 (2005), 072302 doi:10.1103/PhysRevLett.94.072302 [arXiv:nucl-th/0409068 [nucl-th]].
  * (22) A. Rebhan, P. Romatschke and M. Strickland, Phys. Rev. Lett. 94 (2005), 102303 doi:10.1103/PhysRevLett.94.102303 [arXiv:hep-ph/0412016 [hep-ph]].
  * (23) P. B. Arnold, G. D. Moore and L. G. Yaffe, Phys. Rev. D 72 (2005), 054003 doi:10.1103/PhysRevD.72.054003 [arXiv:hep-ph/0505212 [hep-ph]].
  * (24) P. Romatschke and A. Rebhan, Phys. Rev. Lett. 97 (2006), 252301 doi:10.1103/PhysRevLett.97.252301 [arXiv:hep-ph/0605064 [hep-ph]].
  * (25) D. Bodeker and K. Rummukainen, JHEP 07 (2007), 022 doi:10.1088/1126-6708/2007/07/022 [arXiv:0705.0180 [hep-ph]].
  * (26) S. Bazak and S. Mrowczynski, Phys. Rev. D 105 (2022) no.3, 034023 doi:10.1103/PhysRevD.105.034023 [arXiv:2111.11396 [hep-ph]].
  * (27) S. Bazak and S. Mrowczynski, Phys. Rev. D 106 (2022) no.3, 034031 doi:10.1103/PhysRevD.106.034031 [arXiv:2205.08282 [hep-ph]].
  * (28) Z. G. Mou, P. M. Saffin and A. Tranberg, JHEP 02 (2021), 189 doi:10.1007/JHEP02(2021)189 [arXiv:2006.13620 [hep-th]].
  * (29) M. Hindmarsh, S. J. Huber, K. Rummukainen and D. J. Weir, Phys. Rev. Lett. 112 (2014), 041301 doi:10.1103/PhysRevLett.112.041301 [arXiv:1304.2433 [hep-ph]].
  * (30) C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No and K. Rummukainen, et al. JCAP 03 (2020), 024 doi:10.1088/1475-7516/2020/03/024 [arXiv:1910.13125 [astro-ph.CO]].