# Development of nonlinear two fluid interfacial structures by combined action
of Rayleigh-Taylor, Kelvin-Helmholtz and Richtmyer-Meshkov
instabilities:Oblique shock

M. R. Gupta, Labakanta Mandal, Sourav Roy , Rahul Banerjee,Manoranjan Khan
Dept. of Instrumentation Science & Centre for Plasma Studies
Jadavpur University, Kolkata-700032, India
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The nonlinear evolution of two fluid interfacial structures like bubbles and
spikes arising due to the combined action of Rayleigh-Taylor and Kelvin-
Helmholtz instability or due to that of Richtmyer-Meshkov and Kelvin-Helmholtz
instability resulting from oblique shock is investigated. Using Layzer’s model
analytic expressions for the asymptotic value of the combined growth rate are
obtained in both cases for spikes and bubbles. However, if the overlying fluid
is of lower density the interface perturbation behaves in different ways.
Depending on the magnitude of the velocity shear associated with Kelvin-
Helmholtz instability both the bubble and spike amplitude may simultaneously
grow monotonically (instability) or oscillate with time or it may so happen
that while this spike steepens the bubble tends to undulate. In case of an
oblique shock which causes combined action of Richtmyer-Meshkov instability
arising due to the normal component of the shock and Kelvin Helmholtz
instability through creation of velocity shear at the two fluid interface due
to its parallel component, the instability growth rate-instead of behaving as
$1/t$ as $t\rightarrow\infty$ for normal shock, tends asymptotically to a
spike peak height growth velocity
$\sim\sqrt{\frac{5(1+A_{T})}{16(1-A_{T})}(\Delta v)^{2}}$ where $\Delta v$ is
the velocity shear and $A_{T}$ is the Atwood number. Implication of such
result in connection with generation of spiky fluid jets in astrophysical
context is discussed.

## I. INTRODUCTION

Rayleigh-Taylor (RTI) and Kelvin-Helmholtz (KHI) instabilities are associated
with the perturbation of the interface of two fluids of different densities
subject to the action of continuously acting acceleration (with respect to
time) and under the action of velocity shear,respectively. The perturbation
and the consequent instability may also be induced by a shock generated
impulsive acceleration known as Richtmyer-Meshkov (RMI) instability. Such
interfacial hydrodynamic instabilities occur in a wide range of physical
phenomenon from those associated with problems on wave generation by wind
blowing over water surface to problems related to Inertial Confinement Fusion
(ICF) or astrophysical problems like that of supernova explosion remnant which
belong to the domain of high energy density (HED) physics
${}^{\cite[cite]{[\@@bibref{}{rd06}{}{}]}}$. In ICF experiment HED plasmas may
be created due to multi kilo Joule laser with a pressure $\sim$ Mbar. In
ICF,in addition to RT and RM instabilities nonspherical implosion generate
shear flows; the later is also formed when shocks pass through irregular fluid
interfaces. The KHI and shear flow effects in general are also of practical
importance in a number of HED system. They should be considered in a multi
shock implosion schemes for direct drive capsule for ICF, since KHI may
accelerate the growth of turbulent mixing layer at the interface between the
ablator and solid deuterium-tritium nuclear fuel. In HED and astrophysical
system, it has been seen that structures driven by shear flow appear on the
high density spikes produced by R-T and R-M
instabilities${}^{\cite[cite]{[\@@bibref{}{kk03}{}{}]}}$. They may develop in
course of evolution of these instabilities
${}^{\cite[cite]{[\@@bibref{}{br06}{}{}]}}-^{\cite[cite]{[\@@bibref{}{db07}{}{}]}}$
and cover enormous range of spatial scales from $10^{17}$cm for jets from
young stellar objects to $10^{24}$cm for jets from quasars or active galactic
nuclei${}^{\cite[cite]{[\@@bibref{}{br06}{}{}]}}$. Examples are suggested to
be provided by pillars (”elephant trunk”) of Eagle Nebula which are identified
with spikes of a heavy fluid penetrating a light
fluid${}^{\cite[cite]{[\@@bibref{}{ls54}{}{}]}},^{\cite[cite]{[\@@bibref{}{ef54}{}{}]}}$.
Another example in astrophysics is the Herbig-Haro (HH) object like HH34,
where jets are observed with knots. Buhrke et al
${}^{\cite[cite]{[\@@bibref{}{tb88}{}{}]}}$ explained that Kelvin-Helmholtz
instability is the reason for knots in the jets. The jet must be $\sim 10$
times denser than its surrounding medium having velocity $\sim$300 km/sec and
Mach no. 30\. Steady isolated jets may form structure through the growth of
K-H modes. The stability properties of super magnetosonic astrophysical jets
are subject of current interest.

The linear theory of the combined effects of RT,KH and RM instabilities have
been investigated earlier ${}^{\cite[cite]{[\@@bibref{}{m94}{}{}]}}$. Weakly
nonlinear theoretical results of Kelvin-Helmholtz and Rayleigh-Taylor
instability growth rates together with different aspects of density and shear
velocity gradients have also been
discussed${}^{\cite[cite]{[\@@bibref{}{lw09}{}{}]}}-^{\cite[cite]{[\@@bibref{}{lw10}{}{}]}}$.

In case of the temporal evolution of these instabilities nonlinear structures
develop at the two fluid interface. The structure is called a bubble if the
lighter fluid pushes across the unperturbed surface into the heavier fluid and
a spike if the opposite takes place. The dynamics of such RTI and RMI
generated nonlinear structures have been studied
${}^{\cite[cite]{[\@@bibref{}{jh94}{}{}]}}-^{\cite[cite]{[\@@bibref{}{ps03}{}{}]}}$
under different physical situation using an expression near the tip of the
bubble or spike up to second order in the transverse coordinate to unperturbed
surface following Layzer’s
${}^{\cite[cite]{[\@@bibref{}{mr09}{}{}]}}$approach.

In the present paper, we investigate the combined effect of Rayleigh-
Taylor,Richtmyer-Meshkov and Kelvin Helmholtz instabilities by extending the
above method so as to include the effect of velocity shear induced
contribution to the growth rate of the tip of the nonlinear mushroom like
structures generated by shock wave (normal or oblique) incident on the
unperturbed interface.

In the event of excitation of RM instability due to normal incidence of shock
in absence of velocity shear of the growth rate of the height of the finger
like structures decay as
$(1/t)$${}^{\cite[cite]{[\@@bibref{}{vn02}{}{}]}},^{\cite[cite]{[\@@bibref{}{ps03}{}{}]}}$.
It is however interesting to note that if the shock incidence is oblique (or
if it passes across an irregular surface) the growth rate of the tip of the
spiky structure does not decrease as $(1/t)$ but attains a saturation value
proportional to $\sqrt{k^{2}(\Delta v)^{2}/(1-A_{T})}$ where $\Delta
v$=difference is the tangential velocity of the fluids at the interface and
$A_{T}$ is the Atwood number. Thus the growth rate may be quite large if
$A_{T}\rightarrow 1$ which may be likely in astrophysical situation and thus
play an important role in formation of
jets${}^{\cite[cite]{[\@@bibref{}{br06}{}{}]}},^{\cite[cite]{[\@@bibref{}{dd02}{}{}]}}$.

The paper is organized in the following manner. In section II is developed the
basic equations describing the dynamics of nonlinear structures which evolve
in consequence of the combined effects of these different types of
hydrodynamical instabilities. In section III it is shown that the classical
results${}^{\cite[cite]{[\@@bibref{}{dl55}{}{}]}}$ follow on linearization of
the evolution equation describing the bubbles and spikes. Numerical as well as
some analytical results regarding the saturation growth rates are presented in
section IV. Finally section V presents a brief summary of this results.

## II. BASIC EQUATIONS OF EVOLUTION OF THE HYDRODYNAMIC INSTABILITIES

Let the $y=0$ plane denote the unperturbed surface of separation of two fluids
(the line $y=0$ in the two dimensional form of this problem). The fluid with
density $\rho_{a}$ is assumed to overlie the fluid with density $\rho_{b}$.
The gravity $\overrightarrow{g}$ is assumed to act along the negative y- axis.
Any perturbation of the horizontal interface or a shock driven impulse gives
rise to Rayleigh-Taylor instability($\rho_{a}>\rho_{b}$) or Richtmyer -Meshkov
instability which in course of temporal evolution gives rise to nonlinear
interfacial structures.

The two fluids separated by the horizontal boundary are further assumed to be
in relative horizontal motion and thus subjected to Kelvin-Helmholtz
instability arising due to horizontal velocity shear. Thus we are faced with
the problem of the combined action of Rayleigh-Taylor and Kelvin-Helmholtz
instabilities.We shall see the same formulation will be applicable to
Richtmyer-Meshkov instability associated with an oblique shock incident on the
two fluid interface.
After perturbation the finger shaped interface is assumed to take up a
parabolic shape given by

$\displaystyle y(x,t)=\eta(x,t)=\eta_{0}(t)+\eta_{2}(t)(x-\eta_{1}(t))^{2}$
(1)

For a bubble (here the lower fluid is pushing across the interface into the
upper fluid with density $\rho_{a}>$ density $\rho_{b}$) we have,

$\displaystyle\mbox{}\qquad\qquad\eta_{0}>0\quad\mbox{and}\quad\eta_{2}<0$ (2)

and for spike:

$\displaystyle\qquad\qquad\eta_{0}<0\quad\mbox{and}\quad\eta_{2}>0$ (3)

The height of the vertex of the parabola i.e, the height of the peak of the
bubble (or spike) above the x-axis is $|\eta_{0}(t)|$. The position of the
peak at time t is at $x=\eta_{1}(t)$ and because of the relative streaming
motion of the two fluids the peak moves parallel to the x-axis with velocity
$\dot{\eta}_{1}(t)$. The densities of both fluids are uniform and fluid motion
is supposed to be single mode potential flow.

For the upper fluid with density $\rho_{a}$ we take the velocity potential

$\displaystyle\phi_{a}(x,y,t)=\left[\alpha_{a}(t)\cos{(k(x-\eta_{1}(t)))}+\beta_{a}(t)\sin{(k(x-\eta_{1}(t)))}\right]e^{-k(y-\eta_{0}(t))}-xu_{a}(t);\quad
y>0$ (4)

and for the lower fluid (density $\rho_{b}$) the velocity potential

$\displaystyle\phi_{b}(x,y,t)=\left[\alpha_{b}(t)\cos{(k(x-\eta_{1}(t)))}+\beta_{b}(t)\sin{(k(x-\eta_{1}(t)))}\right]e^{k(y-\eta_{0}(t))}-xu_{b}(t)+yb_{0}(t);\quad
y<0$ (5)

Before proceeding with the analysis of the kinematic and boundary conditions
using the two fluid interface perturbation $y=\eta(x,t)$ we forward the
following justification for restricting the expansion to terms
O$(x-\eta_{1}(t))^{2}$. We are concerned only motion very close to the tip of
the bubble or spike i.e., only in the region $k\mid x-\eta_{1}(t)\mid<<1.$
Consequently one is justified in neglecting terms O $(\mid
x-\eta_{1}(t)\mid)^{4}$ unless the coefficients of such terms are sufficiently
large. Further it has been shown ${}^{\cite[cite]{[\@@bibref{}{mr10}{}{}]}}$
that even it terms
$\sim\eta_{4}(t)(x-\eta_{1}(t))^{4}+\eta_{6}(t)(x-\eta_{1}(t))^{6}$ are
retained the contribution from coefficients
$\mid\eta_{4}\mid,\mid\eta_{6}\mid<<$ that from $\mid\eta_{2}\mid$ at least in
the asymptotic state $\tau\rightarrow\infty$.

The kinematical boundary conditions satisfied at the interfacial surface
$y=\eta(x,t)$are

$\displaystyle\frac{\partial\eta}{\partial t}-\frac{\partial\phi_{a}}{\partial
x}\frac{\partial\eta}{\partial x}=-\frac{\partial\phi_{a}}{\partial y}$ (6)
$\displaystyle-\frac{\partial\phi_{a}}{\partial x}\frac{\partial\eta}{\partial
x}+\frac{\partial\phi_{a}}{\partial y}=-\frac{\partial\phi_{b}}{\partial
x}\frac{\partial\eta}{\partial x}+\frac{\partial\phi_{b}}{\partial y}$ (7)

The dynamical boundary conditions are next obtained from Bernoulli’s equation
for the two fluids

$\displaystyle-\frac{\partial\phi_{a}}{\partial
t}+\frac{1}{2}(\vec{\nabla}\phi_{a})^{2}+gy\rho_{a}=-p_{a}+f_{a}(t)$ (8)
$\displaystyle-\frac{\partial\phi_{b}}{\partial
t}+\frac{1}{2}(\vec{\nabla}\phi_{b})^{2}+gy\rho_{b}=-p_{b}+f_{b}(t)$ (9)

by using the surface pressure equality

$\displaystyle p_{a}=p_{b}$ (10)

leading to

$\displaystyle-(\frac{\partial\phi_{a}}{\partial
t}-\frac{\partial\phi_{b}}{\partial
t})+\frac{1}{2}(\vec{\nabla}\phi_{a})^{2}-\frac{1}{2}(\vec{\nabla}\phi_{b})^{2}+g(\rho_{a}-\rho_{b})y$
$\displaystyle=f_{a}(t)-f_{b}(t)$ (11)

satisfied at the interface $y=\eta(x,t)$ Now from Eq.(1)

$\displaystyle\frac{\partial\eta}{\partial
t}=\dot{\eta}_{0}(t)-2\dot{\eta}_{1}(t)\eta_{2}(t)(x-\eta_{1}(t))+\dot{\eta}_{2}(t)(x-\eta_{1}(t))^{2}$
(12)

Also utilizing the property that close to the tip of the bubble or spike,
$k|x-\eta_{1}(t)|<<1$, we express the velocity components in the following
form

$\displaystyle v_{ax}=-\frac{\partial\phi_{a}}{\partial
x}=(u_{a}-k\beta_{a})+k^{2}\alpha_{a}(x-\eta_{1})+\beta_{a}k^{2}(\eta_{2}+k/2)(x-\eta_{1})^{2}$
(13) $\displaystyle v_{ay}=-\frac{\partial\phi_{a}}{\partial
y}=k\alpha_{a}+k^{2}\beta_{a}(x-\eta_{1})-k^{2}\alpha_{a}(\eta_{2}+k/2)(x-\eta_{1})^{2}$
(14)

and similar expressions for $v_{bx}$ and $v_{by}$.

Following Layzer’s${}^{\cite[cite]{[\@@bibref{}{dl55}{}{}]}}$ model we
substitute for $\eta_{t},\eta_{x},(v_{a{(b)}})_{x},(v_{a{(b)}})_{y}$in the
kinematic and boundary conditions represented by Eqs.(6),(7)and (11)and equate
coefficients of $(x-\eta_{1}(t))^{i};(i=0,1,2)$ and neglect terms
$O((x-\eta_{1}(t))^{i});(i\geq 3)$.This yields the following three algebraic
equations for the three unknown $b_{0},\alpha_{b},\beta_{b}$ :

$\displaystyle b_{0}=-\frac{6\eta_{2}}{(3\eta_{2}-k/2)}k\alpha_{a}$ (15)
$\displaystyle\alpha_{b}=\frac{(3\eta_{2}+k/2)}{(3\xi_{2}-k/2)}\alpha_{a}$
(16)
$\displaystyle\beta_{b}=\frac{(\eta_{2}+k/2)k\beta_{a}-{\eta_{2}(u_{a}-u_{b})}}{k(\eta_{2}-k/2)}$
(17)

and regarding the five other unknowns,viz
$\eta_{0}(t),\eta_{1}(t),\eta_{2}(t),\alpha_{a}(t),\beta_{b}(t)$ the following
five nonlinear ODE’s [Eqs.(18)-(22)].

$\displaystyle\frac{d\xi_{1}}{d\tau}=\xi_{3}$ (18)
$\displaystyle\frac{d\xi_{2}}{d\tau}=-\frac{1}{2}(6\xi_{2}+1)\xi_{3}$ (19)
$\displaystyle\frac{d\xi_{3}}{d\tau}=\frac{N_{1}(\xi_{2},r)}{D_{1}(\xi_{2},r)}\frac{\xi_{3}^{2}}{(6\xi_{2}-1)}+\frac{2(1-r)\xi_{2}(6\xi_{2}-1)}{D_{1}(\xi_{2},r)}+\frac{N_{2}(\xi_{2},r)}{D_{1}(\xi_{2},r)}\frac{(6\xi_{2}-1)\xi_{4}^{2}}{2\xi_{2}(2\xi_{2}-1)^{2}}\hskip
70.0pt$
$\displaystyle+2\frac{(4\xi_{2}-1)(6\xi_{2}-1)}{D_{1}(\xi_{2},r)(2\xi_{2}-1)^{2}}\left[(V_{a}-V_{b})^{2}\xi_{2}-(V_{a}-V_{b})(2\xi_{2}+1)\xi_{4}\right]$
(20)
$\displaystyle\frac{d\xi_{4}}{d\tau}=\frac{(2\xi_{2}-1)}{D_{2}(\xi_{2},r)}\left[(f_{b}-rf_{a})-r\frac{\xi_{3}\xi_{4}}{2\xi_{2}}\right]+\frac{2(f_{a}-f_{b})}{D_{2}(\xi_{2},r)}\xi_{2}\hskip
200.0pt$
$\displaystyle+\frac{(6\xi_{2}+1)\xi_{3}}{2D_{2}(\xi_{2},r)(6\xi_{2}-1)(2\xi_{2}-1)}\left[4(V_{a}-V_{b})(4\xi_{2}-1)-\frac{\xi_{4}}{\xi_{2}}(28\xi_{2}^{2}-4\xi_{2}-1)\right]$
(21)
$\displaystyle\frac{d\xi_{5}}{d\tau}=V_{a}-\frac{\xi_{4}(2\xi_{2}+1)}{2\xi_{2}}$
(22)

where

$\displaystyle\xi_{1}=k\eta_{0};\xi_{2}=\eta_{2}/k;\xi_{5}=k\eta_{1}$ (23)
$\displaystyle\xi_{3}=k^{2}\alpha_{a}/\sqrt{kg};\xi_{4}=k^{2}\beta_{a}/\sqrt{kg},\tau=t\sqrt{(}kg)$
(24) $\displaystyle
V_{a}=u_{a}\sqrt{(}k/g);V_{b}=u_{b}\sqrt{(}k/g);f_{a}=\frac{dV_{a}}{d\tau};f_{b}=\frac{dV_{b}}{d\tau}.$
(25)

The functions $N_{1,2}(\xi_{2},r),D_{1,2}(\xi_{2},r)$ where
$r=\frac{\rho_{a}}{\rho_{b}}$ is the density ratio are given by

$\displaystyle N_{1}(\xi_{2},r)=36(1-r)\xi_{2}^{2}+12(4+r)\xi_{2}+(7-r)$ (26)
$\displaystyle D_{1}(\xi_{2},r)=12(r-1)\xi_{2}^{2}+4(r-1)\xi_{2}-(r+1)$ (27)
$\displaystyle N_{2}(\xi_{2},r)=16(1-r)\xi_{2}^{3}+12(1+r)\xi_{2}^{2}-(1+r)$
(28) $\displaystyle D_{2}(\xi_{2},r)=2(1-r)\xi_{2}+(1+r)$ (29)

The above set of five Eqs. (18)-(22) together with Eqs. (23)-(29)which define
the different variables and functions describe the combined effect of RT and
KH instabilities.

On the other hand the impingement of an oblique shock on the two fluid
interface causes the joint effect of Richtmyer-Meshkov and Kelvin-Helmholtz
instability. The impact gives rise to an instantaneous acceleration which will
change the normal velocity (y-component) by an amount $\Delta
v=v_{after}-v_{before}$ and transverse velocity (x-component) by $\Delta
u_{a(b)}=(u_{a(b)})_{after}-(u_{a(b)})_{before}$. Taking nonzero values only
for the post shock velocities we replace the acceleration by their impulsive
values. We set:

$\displaystyle\frac{du_{a(b)}}{dt}=u_{a}\delta(t)\rightarrow\Delta v(t)$ (30)

and replace $g\rightarrow\Delta v\delta(t)$ The dynamical variables are non
dimensionalized using normalization in terms of $(k\Delta v)$instead of
$\sqrt{kg}$.

The combined effect of RM-KH instability resulting from oblique incidence of
shock on the two fluid interface is then described by the same set of
equations as Eqs.(18)-(22) together with the following replacements:
$(i)$The second term on the RHS of Eq.(20)drops out.
(ii) $\xi_{3},\xi_{4}$ and $\tau$ to be replaced by
$\overline{\xi_{3}}=\alpha_{a}k^{2}/(k\Delta
v),\overline{\xi_{4}}=\beta_{a}k^{2}/(k\Delta v)$ and $\tau=t(k\Delta v)$
respectively.
(iii) $V_{a}$ and $V_{b}$ by $\overline{V}_{a}=u_{a}/\Delta
v,\overline{V}_{b}=u_{b}/\Delta v$.
(iv) $f_{a}$ by

$\displaystyle\overline{f}_{a}=\frac{d\overline{v}_{a}}{d\overline{\tau}}=\frac{u_{a}}{\Delta
v}\Delta(\overline{\tau})\quad\mbox{and}\quad
f_{b}\quad\mbox{by}\quad\overline{f_{b}}=\frac{d\overline{v}_{b}}{d\overline{\tau}}=\frac{u_{b}}{\Delta
v}\Delta(\overline{\tau})$ (31)

## III. LINEAR APPROXIMATION

We now show that the usual combined RT and KH instability growth rates
${}^{\cite[cite]{[\@@bibref{}{sc81}{}{}]}}$ are recovered on linearization of
Eqs. (18)-(22). Let us put

$\displaystyle\frac{d(k\eta_{1})}{d\tau}=\frac{d\xi_{5}}{d\tau}=\alpha_{a}V_{a}+\alpha_{b}V_{b}\,;\,\,\,\,\left(\alpha_{a,(b)}=\frac{\rho_{a,(b)}}{\rho_{a}+\rho_{b}}\right)$

in Eq.(22)giving

$\displaystyle\xi_{4}=2\alpha_{b}(V_{a}-V_{b})\frac{\xi_{2}}{2\xi_{2}+1}\approx
2\alpha_{b}(V_{a}-V_{b})\xi_{2}$ (32)

on linearization . In absence of velocity shear $V_{a}-V_{b}=0$,we get
$\xi_{4}=0$. Thus the problem reduces to that of RT instability alone with no
contribution from KH instability. Linearizing Eqs. (19),(20)and (21) we get

$\displaystyle\frac{d\xi_{2}}{d\tau}=-\frac{1}{2}\xi_{3}$ (33)
$\displaystyle\frac{d\xi_{3}}{d\tau}=-2\left[A_{T}+\alpha_{a}\alpha_{b}(V_{a}-V_{b})^{2}\right]\xi_{2}$
(34) $\displaystyle\frac{d\xi_{4}}{d\tau}=-\rho_{a}(V_{a}-V_{b})\xi_{3}$ (35)

$A_{T}=\frac{\rho_{a}-\rho_{b}}{\rho_{a}+\rho_{b}}$ is the Atwood number.
Eq.(32) connecting $\xi_{2}$ and $\xi_{4}$ provides the consistency condition.
The exponential growth rate due to combined effect of RT and KH instability
coincides with the classical linear theory result
${}^{\cite[cite]{[\@@bibref{}{sc81}{}{}]}}$

$\displaystyle\gamma(k)=\sqrt{kg\left[A_{T}+\alpha_{a}\alpha_{b}(V_{a}-V_{b})^{2}\right]}$
(36)

## IV. RESULTS AND DISCUSSIONS

(A) Combined effect of RT and KH instability:

The growth rate of the RT instability induced nonlinear interfacial structures
is further enhanced due to KH instability. Setting $\frac{du_{a}}{dt}=0$ and
$\frac{du_{a}}{dt}=0$ the growth rate of the peak height of the bubbles and
spikes are obtained by numerical integration of Eqs. (18)-(22) and the results
are shown in Fig.1. The dependence of the growth rate on $V_{a}$ and $V_{b}$
keeping $(V_{a}-V_{b})$ unchanged are also indicated in the same diagrams.It
is found that for $V_{b}>V_{a}$ the growth rate is greater than that for
$V_{a}>V_{b}(\mid V_{a}-V_{b}\mid$ is the same for both cases); the asymptotic
values is the two cases are however identical. Moreover for
$\frac{\rho_{a}}{\rho_{b}}>1$ Eqs.(18)-(22) show that as
$\tau\rightarrow\infty$ there occurs growth rate saturation given by

$\displaystyle(\xi_{3})_{bubble}^{asym}=\sqrt{\frac{2A_{T}}{3(1+A_{T})}+\frac{5(1-A_{T})}{16(1+A_{T})}(V_{a}-V_{b})^{2}}$
(37)

and

$\displaystyle(\xi_{3})_{spike}^{asym}=\sqrt{\frac{2A_{T}}{3(1-A_{T})}+\frac{5(1+A_{T})}{16(1-A_{T})}(V_{a}-V_{b})^{2}}$
(38)

while

$\displaystyle(\xi_{4})^{asym}=0$

for both bubble and spike respectively.

Thus both saturation growth rate are enhanced for due to further
destabilization caused by the velocity shear.

On the other hand if
$r=\rho_{a}/\rho_{b}<1(A_{T}=\frac{\rho_{a}-\rho_{b}}{\rho_{a}+\rho_{b}}<0)$
there is no RT instability but it follows from Eqs.(37)and (38) that
instability due to velocity shear (Kelvin -Helmholtz instability) persists on
both the wind ward side and leeward side (i.e; both for bubbles and spikes) if
(see Fig.2)

$\displaystyle\frac{32|A_{T}|}{15(1+A_{T})}=\frac{16(1-r)}{15r}<(V_{a}-V_{b})^{2}$
(39)

and stabilized on both sides if (see Fig.3 which shows oscillation of
$\xi_{1}$and $\xi_{3}$ with respect to $\tau$)

$\displaystyle(V_{a}-V_{b})^{2}<\frac{16}{15}(1-r)$ (40)

If however $(V_{a}-V_{b})^{2}$ lies in the interval specified by the above
inequalities,i.e;

$\displaystyle\frac{16}{15}(1-r)<(V_{a}-V_{b})^{2}<\frac{16(1-r)}{15r};\mbox{(}r<1)$
(41)

it follows from the same two Eqs.(37)-(38) that the peak of the spike
continues to steeper (instability) with $\tau$ as the heavier fluid (density
$\rho_{b}$)pushes across the interface into the lighter fluid (density
$\rho_{a}$) while the bubble height will execute low finite amplitude
undulations. The above observation is shown to be suppressed in Fig.4. At time
t, the peak height which of the spike or the bubble occurs at
$x=\eta_{1}(\tau)$ and thus moves to the right (x-increases) as
$\eta_{1}(\tau)$ increases with $\tau$. The spike peak height increase
monotonically with t while that of the bubble undulates with low amplitude.
The three dimension representation of the steepening of the peak of the spike
as it moves along x-direction with time is shown in Fig.5.In this respect
there exists approximate qualitative agreement exists with the results of the
weakly nonlinear analysis${}^{\cite[cite]{[\@@bibref{}{lw10}{}{}]}}$.

(B) Combined effect of Richtmyer-Meshkov and Kelvin-Helmholtz instability:
oblique shock

The time evolution of the two fluid interfacial structure resulting from the
combined effect of Richtmyer-Meshkov and Kelvin-Helmholtz instabilities
consequent to impingement of an shock is described by the set of
Eqs.(18)-(22),(26)-(29) with modifications as shown in the set of Eq.(31). If
the shock incidence is oblique then the normal component generates velocity
shear and causes KH instability.${}^{\cite[cite]{[\@@bibref{}{m94}{}{}]}}$ The
shock generated initial values of $\overline{\xi}_{3}$ and
$\overline{\xi}_{4}$ are obtained from the impulsive accelerations represented
by the $\delta-$ function terms in Eq.(30) giving

$\displaystyle(\overline{\xi}_{3})_{\tau=0}=\left[\frac{2(1-r)\xi_{2}(6\xi_{2}-1)}{D_{1}(\xi_{2},r)}\right]_{{(\xi_{2})}{\tau=0}}$
(42)
$\displaystyle(\overline{\xi}_{4})_{\tau=0}=\frac{1}{D_{2}(\xi_{2},r)}\left[\frac{(2\xi_{2}-1)(u_{b}-ru_{a})+2\xi_{2}(u_{a}-u_{b})}{\Delta
v}\right]_{(\xi_{2})_{\tau=0}}$ (43)

Results obtained from numerical solution of Eqs.(18)-(22) with modifications
given by Eq.(31) subject to initial conditions (42) and (43) are presented in
Fig.6. The growth rate contributed in absence of velocity shear,i.e; by
normally incident shock induced Richtmyer-Meshkov instability varies as
$t\rightarrow\infty$. However in presence of velocity shear the growth rate
due to combined influence of RM and KH instability the growth rate approaches
finite saturation value asymptotically.

For RM-KH instability induced spikes it is given by the following closed
expression

$\displaystyle(\overline{\xi}_{3})_{t\rightarrow\infty}^{spike}=(\frac{d\xi_{1}}{dt})_{t\rightarrow\infty}^{spike}=\sqrt{\frac{5(1+A_{T})}{16(1-A_{T})}(u_{a}-u_{b})^{2}/(\Delta
v)^{2}}$ (44)

which becomes large as the Atwood number $A_{T}\rightarrow$1
(equivalently$\rho_{a}/\rho_{b}>>1)$

The following discussions suggest a higher plausibility of the effectiveness
of the joint influence of RM and KH instability in the explanation of certain
astrophysical phenomena.

Corresponding to parameter values for Eagle
Nebula($\rho_{a}/\rho_{b}=0.5\times 10^{2}$ and $|u_{a}-u_{b}|=2\times 10^{6}$
cm
sec-1)${}^{\cite[cite]{[\@@bibref{}{br06}{}{}]}},^{\cite[cite]{[\@@bibref{}{dd02}{}{}]}}$
the velocity of rise of the spike peak height hspike(the height of the pillar)
according to Eq.(44)is $(\frac{dh}{dt})_{t\rightarrow\infty}^{spike}\approx
0.79\times 10^{7}$cm sec-1. Modification through inclusion of Rayleigh-Taylor
instability effect (see Eq.(38)) can only slightly increases this value to
$\approx 10^{7}$cm sec-1.This gives the time to reach the observed pillar
height of $3\times 10^{19}$ cm of the Eagle Nebula $\approx 10^{4}$ years.
There are different time scales involved in the problem of development of the
pillar of the Eagle Nebula. As pointed out by Pound,
${}^{\cite[cite]{[\@@bibref{}{mw98}{}{}]}}$ there is a characteristic time
scale for hydrodynamic motion $\tau_{dyn}\approx(\Delta v)^{-1}$ where $\Delta
v$ is the velocity shear inside the cloud. Corresponding to data given in
ref.(3) this turns out to be $\tau_{hydrodanmic}\approx 10^{5}$yrs which is
the upper time limit for development of the Eagle Nebula pillar(”elephant
trunk”). But it is at least two orders of magnitude greater than the time
scale $\tau_{cool}\sim 10^{2}-10^{3}$ yrs imposed due to radiative cooling of
the cloud
${}^{\cite[cite]{[\@@bibref{}{br06}{}{}]}},^{\cite[cite]{[\@@bibref{}{dd02}{}{}]}}$.
In comparison the time scale of the development of the pillar is found here
$\approx 10^{4}$yrs. Thus consequent to the hydrodynamic model based on the
combined influence of Richtmyer-Meshkov and Kelvin-Helmholtz instability the
gap between the two time scales $\tau_{cool}$ and $\tau_{hydrodynamic}$is
reduced by one order of magnitude.

A high Mach number, radiatevily cooled jet of astrophysical interest has been
produced in laboratory using intense laser irradiation of a gold
cone${}^{\cite[cite]{[\@@bibref{}{df99}{}{}]}}$. The evolution of the jet was
imaged in emission and radiography.

K-H instability growth rate has recently been observed in HED plasma
experiment using Omega laser ($\lambda$)=0.351$\mu m$ delivering 4.3 $\pm
0.1$kJ to the target overlapping 10 drive beams on to the ablator
${}^{\cite[cite]{[\@@bibref{}{eh09}{}{}]}}$.Incompressible K-H growth rate
peak to valley at Foam-Plastic interface has been compared with several
analytical modes.

## V. Summary

Finally we summarize the results:

(a)If the heavier fluid overlies the lighter fluid the growth rate of both the
bubble and spike peak heights due to RT instability are enhanced due to
concurrent presence of velocity shear, i.e, K.H instability Fig.1. The
asymptotic growth rates are given by Eqs. (37) and (38).

(b) In the opposite case,i.e, if the overlying fluid is lighter and lower one
is heavier ($r=\rho_{a}/\rho_{b}<1$) both the spike and bubble peak
displacement increases continuously with time if
$\frac{16(1-r)}{15r}<(V_{a}-V_{b})^{2}$,i.e, instability persists (Fig.3)
while stabilization occurs if $(V_{a}-V_{b})^{2}<\frac{(16(1-r)}{15})$ (Fig. 4
shows oscillation of peak heights of bubbles and spikes).

(c) For $\frac{16(1-r)}{15}<(V_{a}-V_{b})^{2}<\frac{16(1-r)}{15r}$ for r$<1$
the spike steepens with time (the peak height continuously increases with time
as indicated in Fig.5. gives a three dimensional graph of displacement y
against x and $\tau$). But the peak displacement of the bubble undulates
within a small range Fig.4.

(d)If the two fluid interface is subjected to an oblique shock Kelvin-
Helmholtz instability due to generation of velocity shear occurs
simultaneously with Richtmyer-Meshkov instability.The growth rates of bubbles
and spikes due to this joint action are shown in Fig.6. respectively. It is
important to note that the growth rate of the combined action tends
asymptotically to a saturation value given by Eq.(44); this is in contrast to
that due to generation of RM instability due to normal shock incidence for
which the growth rate behaves as $\frac{1}{t}$ as t$\rightarrow\infty$.
Moreover this growth rate as shown by Eq.(44) the rate of growth of this spike
height has sufficiently large magnitude if the Atwood number
$A_{T}\rightarrow$1($\rho_{a}<<\rho_{b}$).This may have interesting
implication in the hydrodynamic explanation of formation of sufficiently long
spiky jets in astrophysical situation, e.g, in case of the Eagle Nebula.

## ACKNOWLEDGEMENTS

This work is supported by the Department of Science & Technology, Government
of India under grant no. SR/S2/HEP-007/2008.

## References

  * [1] R.P.Drake, High Energy Density Physics, Spinger, (2006).
  * [2] K.Kifonidis,T.Plewa,H.T.Janka,E.Muller, Astron. and Astrophys. 408,621 (2003).
  * [3] B.A.Remington,R.P.Drake and D.D.Ryutov, Rev. Mod. Phys. 78, 755 (2006).
  * [4] D.D.Ryutov and B.A Remington, Plasma Phys. Control Fusion 44, B407 (2002).
  * [5] B.A. Remington,R.P.Drake,H.Takabe and D.Arnett, Phys. Plasmas 7, 1641 (2000).
  * [6] L.Spitzer,Jr., Astrophys. J. 120, 1 (1954).
  * [7] E.A.Frieman, Astrophys. J. 120, 18 (1954).
  * [8] T.Buhrke,R.Mundt and T.P.Ray Astron. and Astrophys. 200,99 (1988).
  * [9] K.O.Mikaelian, Phys. Fluids 6, 1943 (1994).
  * [10] L.F.Wang,W.H.Ye,Z.F.Fan,Y.J.Li,X.T.He and M.Y.Yu, Europhys. Lett. 86, 15002 (2009).
  * [11] L.F.Wang,W.H.Ye,Y.J.Li, Europhys. Lett. 87,54005 (2009).
  * [12] L.F.Wang,C.Xue,W.H.Ye and Y.J.Li, Phys. Plasmas 16, 112104 (2009).
  * [13] L.F.Wang,W.H.Ye and Y.J.Li, Phys. Plasmas 17,052305 (2010).
  * [14] J.Hecht, U.Alon and D.Shvarts, Phys. Fluids 6,4019 (1994).
  * [15] A.L.Velikovich and G.Dimonte, Phys.Rev.Lett.76,3112 (1996).
  * [16] G.Hazak, Phys.Rev.Lett. 76,4167 (1996).
  * [17] Q.Zhang, Phys.Rev.Lett. 81,3391 (1998).
  * [18] V.N.Goncharov, Phys.Rev.Lett. 88,134502 (2002).
  * [19] S.I.Sohn, Phys.Rev. E 67,026301 (2003).
  * [20] M.R.Gupta,S.Roy,M.Khan,H.C.Pant,S.Sarkar and M.K.Srivastava, Phys.Plasma. 16,032303 (2009).
  * [21] D.Layzer, Astrophys. J. 122,1 (1955).
  * [22] S.Chandrasekhar, Hydrodynamics and Hydrodynamic Stability, Dover,New York 1981.
  * [23] M.R.Gupta,L.K.Mandal,S.Roy and M.Khan, Phys.Plasma. 17,012306 (2010).
  * [24] M.W.Pound, Astrophys. J. 493,L113 (1998).
  * [25] D.R.Farely,K.G.Estabrook,S.G.Glendinning,S.H.Glenzer,B.A.Remington,K.Shigemori,
J.M.Stone,R.J.Wallace,G.B.Zimmerman and J.A.Harte, Phy.Rev.Lett. 83, 1982
(1999).

  * [26] E.C.Harding,J.F.Hansen,O.A.Hurricane,R.P.Drake,H.F.Robey,C.C.Kuranz,B.A.Remington,
M.J.Bono,M.J.Grosskopf and R.S.Gillespie, Phy.Rev.Lett. 103, 045005 (2009).

> Figure 1: Initial values
> $r=\frac{\rho_{a}}{\rho_{b}}=1.5,\xi_{1}=-\xi_{2}=\xi_{3}=\xi_{4}=\xi_{5}=0.1$
> for bubble;$-\xi_{1}=\xi_{2}=-\xi_{3}=\xi_{4}=\xi_{5}=0.1$ for spike.Plot
> showing variation of $\xi_{1},\xi_{2}$,growth rate $\xi_{3},\xi_{4}$ and
> transverse displacement $\xi_{5}$ of bubble and spike with $V_{a}=V_{b}=0.0$
> for solid black line-spike and broken black line for
> bubble.$V_{a}=0.1,V_{b}=0.5$ for broken blue line-bubble and solid blue for
> spike,$V_{a}=0.5,V_{b}=0.1$,broken red line for bubble and solid red line-
> spike.

> and for following relation$\frac{16(1-r)}{15r}<(V_{a}-V_{b})^{2}.$

Figure 2: r=$\frac{\rho_{a}}{\rho_{b}}=0.4$;Lower fluid denser.Dashed line for
spike (heavier fluid pushes into lighter fluid) and unbroken line for
bubble.$V_{a}=0.8,V_{b}=-0.6$.Initial condition as in Fig.1.

> Figure 3: r=0.4,$V_{a}=0.0,V_{b}=0.2.(V_{a}-V_{b})^{2}<\frac{16(1-r)}{15}$
> .Initial condition as in Fig.1. Unbroken line for bubble and dashed line for
> spike.

> Figure 4:
> r=0.4,$V_{a}=0.6,V_{b}=-0.4;\frac{16(1-r)}{15}<(V_{a}-V_{b})^{2}<\frac{16(1-r)}{15r}$.Initial
> condition as before (Fig.3.).Unbroken line for bubble and dashed line for
> spike;height of spike peak increases monotonically with time
> (steepening);bubble depth undulates.

> Figure 5: 3 dimensional plot of spike(Interface
> $Y=\eta_{0}(\tau)+\eta_{2}(\tau)(x-\eta_{1}(\tau))^{2}$)belonging to the
> plot given in fig.4.

> Figure 6: Oblique shock: RM and KH instability for spike (dashed line) and
> bubble (unbroken line).Initial values as in fig.1. and
> $V_{a}=0.1,V_{b}=0.5.$


# Simulating Coverage Path Planning with Roomba

Robert Chuchro
<EMAIL_ADDRESS>
###### Abstract

Coverage Path Planning involves visiting every unoccupied state in an
environment with obstacles. In this paper, we explore this problem in
environments which are initially unknown to the agent, for purposes of
simulating the task of a vacuum cleaning robot. A survey of prior work reveals
sparse effort in applying learning to solve this problem. In this paper, we
explore modeling a Cover Path Planning problem using Deep Reinforcement
Learning, and compare it with the performance of the built-in algorithm of the
Roomba, a popular vacuum cleaning robot.

##  1 Introduction

A Roomba is a robot vacuum cleaner built with sensors to automatically vacuum
a living space while avoiding obstacles and stair ledges. The Roomba contains
several sensors which allow it to receive information from the outside world.
In this paper, the relevant sensors are the bumper sensors in front of the
Roomba. This sensor will be triggered when the Roomba encounters an obstacle,
such as a wall or piece of furniture. The Roomba can move by turning itself or
driving forward.

##  2 Problem Statement

The environment consists of a bounded two-dimensional plane which may contain
obstacles. Any location which is not an obstacle is reachable from any other
location in the environment. The goal is for the Roomba to traverse the entire
open space in under an hour, avoiding collisions with obstacles.

Initially, the environment is unknown. The agent will need to adopt a policy
to explore the space around it.

###  2.1 Representation

To simplify the problem, we will discretize the environment as a grid world,
with each grid cell being a square with length equal to the diameter of the
Roomba, as shown in Figure (1).

The agent will initialize its world as only having knowledge of the starting
cell, with the remaining neighboring states as hidden states, represented with
’?’. The agent has three actions available to it: drive forward, turn left 90
degrees, and turn right 90 degrees. All three actions are assigned to take one
second, which is represented as one time step.

When the agent attempts to move onto a hidden state, it will observe either an
empty spaces or an obstacle, based on feedback from its bumper sensor. Here,
the goal is to explore the hidden cells, ensuring that every empty space was
visited.

Figure 1: Grid world representation of Roomba’s environment. The green cell is
the base / start location.

###  2.2 Success Measurement

Success of a particular strategy can be measured by time taken to fully
explore a given environment. In the context of simulating, this can be
measured by the number of actions taken before visiting all empty cells.

In order to avoid modeling a desired or complex behavior through the reward
function, the reward function remains kept simple:

$\begin{split}R(s,a,s^{\prime})=\begin{cases}1&s^{\prime}\text{ is an
unvisited cell},\\\ -1&s^{\prime}\text{ is an obstacle},\\\
0&\text{otherwise}\end{cases}\end{split}$ (1)

The environment also has a natural decay, since we are limited to one hour
(1800 time steps).

One of the goals is to come up with a strategy that can outperform Roomba’s
baseline algorithm, which will be outlined in a subsequent section.

Performances amongst the agents can be compared using the following metrics:
percentage of open space explored, number of collisions with obstacles, and
number of time steps to completion.

##  3 Related Work

###  3.1 Path Planning using Discretized States

The coverage algorithm described in [9] involves discretizing a state space
with obstacles, and planning safe paths betweens states. Based on simulation
results, it was shown that using a 4-connected grid representation of the
environment created paths that were no more than 2% less efficient than if it
were 8-connected. This means that restricting the turn angles of the agent to
90 degrees does not have a drastically negative effect on simulation results,
in exchange for reducing the size of the state and action space.

###  3.2 Boustrophedon Cellular Decomposition

Boustrophedon or S-shaped pathing, follows the same pathing that farmers use
on fields of crops, or lawn mowers on a field of grass. The algorithm involves
finding a perimeter, then traversing back and forth along the length of the
perimeter across the entire space [5].

This pathing concept was extended in [2] by attempting to divide the region
into segments, each of which can be completed by executing Boustrophedon
pathing.

###  3.3 Combination of Basic Policies

A combination of policies was shown to outperform any individual strategy in
[5]. The basic pathing strategies were:

  * -

Outward Spiral

  * -

Wall Walking

  * -

Random Walk

  * -

Boustrophedon Pathing

It was not stated how the combination policy was determined; however, a
possible approach would be to learn a scoring function to evaluate which
strategy to use based on the current set of observations.

###  3.4 Deep Reinforcement Learning

One of the main challenges of this task is for the agent to learn a policy
that generalizes well so that its policy is not specific to a particular
environment. A classic example of generalized reinforcement learning is
presented in [8]. Minh, et al. introduced a Deep Q-Network (DQN), a neural
network used for Q-Learning. They used a convolutional neural network to train
seven Atari 2600 games. They were able to achieve super-human results in six
of the seven games it was tested on, with no adjustment of the architecture or
hyperparameters between different games.

We will apply Deep Q-Learning to the simulated environment to produce an agent
that can explore well in unseen environments.

##  4 Methods

###  4.1 Baseline Roomba Algorithm

The basic Roomba strategy consists of exploring outward in a spiral fashion
until the nearest obstacle is reached. The Roomba then alternates between
moving along a discovered obstacle edge and turning then driving forward until
colliding with another obstacle [6]. The Roomba repeats this behavior for an
hour, then attempts to return to its base.

In the context of our grid world representation, the turn angles will need to
be restricted to cardinal directions. After colliding with an obstacle, the
agent will then do random walks for the duration of the episode – driving
forward until an obstacle is encountered, then turning. The algorithm is
outlined below:

⬇

1function baselineRoomba()

2 initialize s, agent

3 while not done

4 if agent.usingSpiral

5 a = spiral(agent, s)

6 else

7 a = randomWalk(agent, s)

8 s = observeNextState(s, a)

Figure 2:  Depiction of Roomba’s basic exploration algorithm [6].

###  4.2 Deep Q-Learning

Pseudo code for Q-Learning algorithm is as follows:

⬇

1initialize Q[num_states,num_actions]

2observe initial state s

3while condition

4 select and carry out an action a

5 observe reward r and new state s’

6␣␣$Q[s,a]\ +=\alpha(r+\gamma maxQ[s^{\prime},a^{\prime}]-Q[s,a])$

7␣␣s=s’

$\alpha$ is the learning rate that controls how much is difference between
previous Q-value and newly proposed Q-value. When $\alpha$ is equal to 1 then
update equation is similar to Bellman equation.

For Deep Q network given a transition $<s,a,r,s^{\prime}>$, the Q-table update
rule can be re written as:

  * Do a feedforward pass for the current state s to get predicted Q-values for all actions.

  * Do a feedforward pass for the next state s’ and calculate maximum overall network outputs $maxQ(s\textquoteright,a\textquoteright)$.

  * Set Q-value target for action to $r+\gamma*maxQ(s\textquoteright,a\textquoteright)$. For all other actions make output equal to 0.

  * Update the weights using back propagation (gradient descent).

  * Optimize with L2 loss: $Loss=(r+maxQ(s^{\prime},a^{\prime})-Q(s,a))^{2}$

Our implementation of Deep Q-Learning is implemented using Tensor Flow [1].
The programming language used was Julia-0.6, and we utilized a Tensor Flow
wrapper library [7]. We use Tensor Flow’s Adam optimization algorithm with an
initial learning rate of 2e-4.

###  4.3 Network Architecture

It is typical that a neural network will consist of several convolutional
layers, followed by fully connected layers, such as in [8]. However, in image
recognition problems, the input data is typically very large and the
convolutional layers are seen as feature identifiers, whereas the fully
connected layers take in the features from previous layers and use their
weights to identify correlations or non-linear combinations between the input
features [3].

The input to our network which represents the state consists of two
components: the first is the number of time steps lapsed, and the second is a
$9\times 9$ array consisting of the partially observed surrounding
environment. That is, the center is marked by the current location, while the
surrounding cells consist of the following values:

  * -1

: Unobserved cell

  * 0

: Observed empty cell

  * 1

: Observed obstacle

  * 2

: Base/Starting cell

By restricting our input to a fixed subset of the actual observed environment,
we are not allowing the model to focus on specific features of an environment,
but instead on learning to make decisions based on potentially reoccurring
subsets of environments.

Since the input to our network is already fairly structured, we omit the
convolutional layers and use only two fully connected layers, the first with
width 64, and the second with a width equal to the number of actions, which in
this simulation is 3. These outputs represent our Q-values. The fully
connected layers are separated by a ReLU non-linearity activation. [4].

###  4.4 Generating Random Environments

Each simulation episode generated a unique environment, so that agents could
not simply learn how to navigate a specific layout. The environment generation
algorithm first generated a 2D room layout with dimensions ranging between 10
and 20. Then, out of 17 ’furniture’ pieces, up to 6 were randomly placed
inside the layout with random orientation. A sample layout is depicted in
figure (3).

Figure 3:  Array representation of sample environment generation.

##  5 Experiments

###  5.1 Sinusoidal Exploration Decay

Typically, reinforcement learning algorithms begin with a high exploration
rate, and decrease it as training progresses. This models the behavior of
initially exploring, followed by exploiting as the model learns the
environment. The decay method is traditionally a linear or exponential decay.
In our early trials, we found that after a certain number of iterations,
training would progress would flatten, but then continue to improve when a new
epoch was started. When a new epoch begins, the most significant difference in
our implementation is the reset of the exploration hyperparameter, $\epsilon$.
We have introduced a new $\epsilon$ decaying function: one which exponentially
decays over episodes in a sinusoidal fashion.

$\begin{split}\epsilon=\epsilon_{0}\cdot\epsilon_{d}^{x}\cdot\frac{1}{2}(1+cos(\frac{2\pi
xn}{X}))\end{split}$ (2)

  * -

$\epsilon_{0}$ is initial epsilon

  * -

$\epsilon_{d}$ is decay rate

  * -

n is number of mini epochs

  * -

X is number of training episodes

  * -

x is current training episode number

Figure 4: Plot of equation (2) (purple) alongside with $\epsilon^{x}$ (blue).
The model can be saved at the troughs of the waves, serving as an epoch.
Parameters used: X=10000, n=5, $\epsilon_{0}$=1, $\epsilon_{d}$=.9997

The main motivation behind this function was for our model to be able to
escape local optima while training with an already reduced $\epsilon$. Since
it is difficult to model an ideal exploration rate for a particular
environment, a decaying sinusoidal function makes less assumptions about the
expected decay of $\epsilon$, hopefully finding the correct decay rate
periodically during the entire epoch. An analogy can be made about the
suggested strategy. It focuses on the idea how one should use car brakes on a
slippery road (without ABS): it is difficult for a driver to apply the brakes
an optimal amount while maintaining rolling friction, but if the driver pumps
the brakes repeatedly, they will achieve the correct amount of pressure
periodically.

###  5.2 Baseline vs DQN

Despite the simplicity and randomness of Roomba’s basline algorithm, a simple
combination of spiraling and straight walks with random turns demonstrated to
be fairly thorough in coverage; however, due to the randomness, it frequently
collided with obstacles. Since no learning occurred, results were consistent
across the 10000 iterations (5).

Figure 5:  Running averages of percent coverage (green) and number of
collisions (red) over 10000 episodes for the baseline algorithm

Preliminary results for only 800 episodes of training using Deep Q-Learning
indicate a strong ability of learning to avoid obstacles, however, there
appears to be a policy trade off between minimizing collisions and coverage
percentage (6).

Figure 6:  Running averages of percent coverage (green) and number of
collisions (red) over 800 episodes for the DQN

Table (7) compares coverage percentage and obstacle avoidance across several
preliminary epochs of the DQN as well as the Roomba baseline algorithm. At one
point, the DQN is able to get over 80% coverage, however, this decreases as
the agent learns to avoid obstacles. It is also worth noting that coverage is
better when $\epsilon$ is slightly greater than 0, indicating that some
randomness may assist with exploring an unknown environment.

Performance Comparison
---
Agent | % Coverage | # of Collisions
Baseline Roomba | 96 | 268
DQN-e1 | 55 | 64
DQN-e2 | 72 | 63
DQN-e3 | 74 | 49
DQN-e4 | 66 | 32
Figure 7:  An apparent trade off between coverage and collisions.

##  6 Conclusion

Preliminary results show that Deep Q-Learning is effective at learning to
avoid obstacles, and that the graphs had not yet leveled off after 800
iterations of training. This warrants for further exploration of this model,
and more training. One Limitation of this study is the insufficient number of
training iterations for the DQN, which should be well into the thousands,
instead of 800. As prior studies suggest, it is likely difficult to properly
model this problem with a pure learning algorithm, and a hybrid approach of
applying some variation of a search algorithm with learning will likely prove
to be more effective.

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] Howie Choset and Philippe Pignon.  Coverage Path Planning: The Boustrophedon Cellular Decomposition, pages 203–209.  Springer London, London, 1998.
  * [3] Adit Deshpande.  A beginner’s guide to understanding convolutional neural networks.  https://adeshpande3.github.io/adeshpande3.github.io/A-Beginner%27s-Guide-To-Understanding-Convolutional-Neural-Networks, 2016\.
  * [4] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In Geoffrey Gordon, David Dunson, and Miroslav Dudík, editors, Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, volume 15 of Proceedings of Machine Learning Research, pages 315–323, Fort Lauderdale, FL, USA, 11–13 Apr 2011\. PMLR.
  * [5] K. M. Hasan, Abdullah-Al-Nahid, and K. J. Reza.  Path planning algorithm development for autonomous vacuum cleaner robots.  In 2014 International Conference on Informatics, Electronics Vision (ICIEV), pages 1–6, May 2014.
  * [6] Julia Layton.  How robotic vacuums work.  https://electronics.howstuffworks.com/gadgets/home/robotic-vacuum.htm, 2005\.
  * [7] malmaud.  Julia wrapper for tensorflow.
  * [8] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin A. Riedmiller.  Playing atari with deep reinforcement learning.  CoRR, abs/1312.5602, 2013.
  * [9] A. Zelinsky, R.A. Jarvis, J. C. Byrne, and S. Yuta.  Planning paths of complete coverage of an unstructured environment by a mobile robot.  In In Proceedings of International Conference on Advanced Robotics, pages 533–538, 1993.


# Stability analysis of non-Abelian electric fields

Jude<EMAIL_ADDRESS>Tanmay<EMAIL_ADDRESS>Physics
Department, Arizona State University, Tempe, Arizona 85287, USA.

###### Abstract

We study the stability of fluctuations around a homogeneous non-Abelian
electric field background that is of a form that is protected from Schwinger
pair production. Our analysis identifies the unstable modes and we find a
limiting set of parameters for which there are no instabilities. We discuss
potential implications of our analysis for confining strings in non-Abelian
gauge theories.

##  I Introduction

In pure non-Abelian gauge theory, the gauge fields carry non-Abelian electric
charge. Hence a non-Abelian electric field is susceptible to decay via the
Schwinger pair production of gauge field quanta Matinyan and Savvidy (1978);
Brown and Weisberger (1979); Yildiz and Cox (1980); Ambjorn and Hughes (1982a,
b); Cooper and Nayak (2006); Nayak and van Nieuwenhuizen (2005); Cooper et al.
(2008); Nair and Yelnikov (2010); Ragsdale and Singleton (2017); Cardona and
Vachaspati (2021). However, confining electric field tube configurations do
not decay, leading to a quandary – how are electric flux tubes stable to
Schwinger pair production? This question was raised and investigated in
Vachaspati (2022). The essential idea is that there are many gauge
inequivalent ways of constructing an electric field in non-Abelian theory
Brown and Weisberger (1979). The straight-forward embedding of the Maxwell
gauge field in the non-Abelian theory is indeed unstable to Schwinger pair
production. However, other inequivalent gauge fields, that nonetheless produce
the same electric field, are protected against quantum dissipation. Such gauge
field configurations are candidates for describing confining electric flux
tubes.

In a quantum theory, there will be fluctuations about the electric field
background and these fluctuations will ultimately be quantized. It is
therefore of interest to determine the fluctuation eigenfrequencies and
eigenmodes, and especially to determine if there are any unstable
fluctuations. The question has been addressed in Refs. Bazak and Mrowczynski
(2022a, b) for a homogeneous electric field where a number of unstable modes
were found333Stability analyses of non-Abelian fields are also of interest in
heavy ion collisions (see for example Ref. Berges et al. (2021).). However
those analyses did not eliminate fluctuations that were inconsistent with
their adopted gauge conditions; neither did they account for extra conditions
imposed by the reality of the gauge fields. Indeed we shall see that these
conditions are critical for the stability analysis. The earlier analyses were
also limited to either a special point in the parameter space of the
background electric field Bazak and Mrowczynski (2022a) or to only the zero
momentum modes Bazak and Mrowczynski (2022b).

We start by describing the homogeneous electric field in SU(2) non-Abelian
gauge theory in Sec. II. Then we consider small fluctuations around the
background electric field in Sec. III, expand the fluctuations in modes in
Sec. IV. The modes get classified according to whether they are longitudinal
or transverse, and whether they are orthogonal to the electric field. The
transverse-orthogonal (TO) modes are discussed in Sec. V while the transverse-
nonorthogonal (TN) and longitudinal (L) modes are discussed together in Sec.
VI as they are coupled. Unstable TO modes are found to exist in the infrared
and depend on the parameters entering the background configuration, and an
interesting limit is found for which the unstable TO modes are absent. The
analysis for the TN and L modes is significantly more complicated and we limit
ourselves to some special cases, for example large or small wavenumbers, and
for wave vectors parallel and orthogonal to the electric field. Our results
again show some unstable modes in the infrared and once more, just as in the
case of TO modes, we find that there are no unstable modes in the special
limit of background parameters.

Our results are summarized in Sec. VII where we also discuss other ideas. The
reader not interested in the technicalities of the analysis can find details
about the electric field background in Sec. II and then proceed to the
conclusions in Sec. VII.

##  II Electric field background

Consider the $SU(2)$ pure gauge theory,

$L=-\frac{1}{4g^{2}}W_{\mu\nu}^{a}W^{\mu\nu a}+j_{\mu}^{a}W^{\mu a}$ (1)

where $g$ is the gauge coupling, $\mu,\nu=0,1,2,3$ are Lorentz indices and
$a=1,2,3$ is the color index. The current $j_{\mu}^{a}$ is an external current
which will be specified below. The field strength $W_{\mu\nu}^{a}$ is given in
terms of the gauge potential $W_{\mu}^{a}$ by,

$W_{\mu\nu}^{a}\equiv\partial_{\mu}W_{\nu}^{a}-\partial_{\mu}W_{\nu}^{a}+\epsilon^{abc}W_{\mu}^{b}W_{\nu}^{c}.$
(2)

The gauge field equations of motion are

${\cal D}_{\nu}W^{\mu\nu a}=j^{\mu a}$ (3)

where

${\cal D}_{\nu}W^{\mu\nu a}\equiv\partial_{\nu}W^{\mu\nu
a}+\epsilon^{abc}W_{\nu}^{b}W^{\mu\nu c}.$ (4)

We wish to consider the stability of a class of gauge fields that give rise to
homogeneous electric fields that we treat as a background. The gauge fields
are

$A_{\mu}^{\pm}\equiv A_{\mu}^{1}\pm iA_{\mu}^{2}=-\epsilon e^{\pm i\Omega
t}\partial_{\mu}z,\ \ A_{\mu}^{3}=0$ (5)

where $\epsilon$ and $\Omega$ are parameters that label members of the class,
and $z$ is the spatial $z$ coordinate. The electric field is gauge equivalent
to Vachaspati (2022),

$E^{a}_{i}=\epsilon\Omega\delta^{a3}\delta_{iz}$ (6)

and the amplitude of the electric field is

$E=\epsilon\Omega$ (7)

As shown in Ref. Brown and Weisberger (1979), gauge fields with distinct
values of $\Omega^{2}$, even for the same value of $E$, are gauge
inequivalent.

In the two dimensional parameter space $(\epsilon,\Omega)$, the electric field
is constant whenever $\epsilon\Omega$ is constant. We will find that the limit
$\epsilon\to 0$, $\Omega\to\infty$ but with $E=\epsilon\Omega$ held constant
to be of interest from the point of view of stability.

The external currents $j_{\mu}^{a}$ in (1) are chosen such that the background
is a solution of the classical equations of motion. Therefore

$j^{\mu a}={\cal D}^{(A)}_{\nu}A^{\mu\nu a}$ (8)

which gives

$j_{\mu}^{\pm}\equiv j_{\mu}^{1}\pm
ij_{\mu}^{2}=-\frac{\epsilon\Omega^{2}}{2}e^{\pm i\Omega t}\partial_{\mu}z,\ \
j_{\mu}^{3}=0.$ (9)

For the purposes of the stability analysis we simply assume that this is an
external current, though it is possible that the currents can arise
semiclassically as discussed in Ref. Vachaspati (2022) and summarized in Sec.
VII.

##  III Fluctuations

We now consider small perturbations around the background,

$W_{\mu}^{a}=A_{\mu}^{a}+q_{\mu}^{a}.$ (10)

Inserting this into the equations of motion, (3), and working to linear order
in the perturbations $q_{\mu}^{a}$ we get,

$\partial_{\nu}q^{\mu\nu a}+\epsilon^{abc}(A_{\nu}^{b}q^{\mu\nu
c}+q_{\nu}^{b}A^{\mu\nu c})=0$ (11)

where

$q_{\mu\nu}^{a}=\partial_{\mu}q_{\nu}^{a}-\partial_{\nu}q_{\mu}^{a}+\epsilon^{abc}(A_{\mu}^{b}q_{\nu}^{c}-A_{\nu}^{b}q_{\mu}^{c})$
(12)

and

$A_{\mu\nu}^{a}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}+\epsilon^{abc}A_{\mu}^{b}A_{\nu}^{c}$
(13)

which, with our chosen background,

$A^{\pm}_{\mu}=-\epsilon e^{\pm i\Omega t}\,\partial_{\mu}z,\ \
A^{3}_{\mu}=0,$ (14)

gives

$A_{\mu\nu}^{\pm}=\pm iEe^{\pm i\Omega
t}(\partial_{\mu}z\partial_{\nu}t-\partial_{\nu}z\partial_{\mu}t),\ \
A_{\mu\nu}^{3}=0$ (15)

We will be adopting temporal gauge ($W_{0}^{a}=0$), so $q_{0}^{a}=0$.

##  IV Mode Expansion

We first define fluctuations in a “rotating frame”, $Q^{a}_{i}$, as follows,

$q^{1}_{i}+iq^{2}_{i}\equiv e^{i\Omega t}(Q^{1}_{i}+iQ^{2}_{i}),\ \
q^{3}_{i}\equiv Q^{3}_{i}$ (16)

where $Q^{a}_{i}$ are real. Next we expand $Q_{i}^{a}$ in spatial and temporal
Fourier modes as follows,

$\displaystyle Q_{i}^{a}$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}e^{-i\omega_{\bf k}t}e^{i{\bf
k}\cdot{\bf x}}p_{i,{\bf k}}^{a}$ (17)

where $p_{i,{\bf k}}^{a}$ are the Fourier amplitudes.

While $\omega_{\bf k}$ can in general be complex, the reality of the fields
$Q_{i}^{a}$ constrain physical values of $\omega_{\bf k}$ and $p_{i}^{a}$ to
satisfy,

$\omega_{\bf k}^{*}=-\omega_{-{\bf k}},\ \ (p^{a}_{i,{\bf
k}})^{*}=p^{a}_{i,-{\bf k}}.$ (18)

In what follows we will consider a single ${\bf k}$ mode and drop the ${\bf
k}$ subscripts, e.g. we write $\omega_{\bf k}$ simply as $\omega$.

Inserting the Fourier expansion into (11) gives 3 constraint equations (the
Gauss constraints) and 9 equations of motion for the 9 components of
$p_{i}^{a}$. The constraints are,

$\displaystyle\omega{\bf k}\cdot{\bf p}^{1}-i\Omega{\bf k}\cdot{\bf
p}^{2}+\epsilon\Omega p_{z}^{3}$ $\displaystyle=$ $\displaystyle 0$ (19)
$\displaystyle\omega{\bf k}\cdot{\bf p}^{2}+i\Omega{\bf k}\cdot{\bf
p}^{1}-i\epsilon\omega p_{z}^{3}$ $\displaystyle=$ $\displaystyle 0$ (20)
$\displaystyle\omega{\bf k}\cdot{\bf p}^{3}+i\epsilon\omega
p_{z}^{2}-2\epsilon\Omega p_{z}^{1}$ $\displaystyle=$ $\displaystyle 0$ (21)

and the equations of motion are,

$\displaystyle\hskip-42.67912pt(-\omega^{2}+k^{2}-\Omega^{2}){\bf
p}^{1}+i2\omega\Omega{\bf p}^{2}-({\bf k}\cdot{\bf p}^{1}){\bf k}=0$ (22)
$\displaystyle\hskip-42.67912pt(-\omega^{2}+k^{2}-\Omega^{2}+\epsilon^{2}){\bf
p}^{2}-i2\omega\Omega{\bf p}^{1}-({\bf k}\cdot{\bf p}^{2}){\bf k}$
$\displaystyle\hskip-14.22636pt-\epsilon(-i{\bf k}\cdot{\bf p}^{3}+\epsilon
p_{z}^{2}){\hat{z}}+i\epsilon p_{z}^{3}{\bf k}-i2\epsilon k_{z}{\bf p}^{3}=0$
(23) $\displaystyle\hskip-28.45274pt(-\omega^{2}+k^{2}+\epsilon^{2}){\bf
p}^{3}-({\bf k}\cdot{\bf p}^{3}){\bf k}-i\epsilon({\bf k}\cdot{\bf
p}^{2}){\hat{z}}$ $\displaystyle\hskip 42.67912pt+i2\epsilon k_{z}{\bf
p}^{2}-i\epsilon p_{z}^{2}{\bf k}-\epsilon^{2}p_{z}^{3}{\hat{z}}=0$ (24)

where we have now employed the vector notation: ${\bf
p}^{a}=(p^{a}_{1},p^{a}_{2},p^{a}_{3})$ and ${\bf k}\cdot{\bf
p}^{a}=k_{i}p^{a}_{i}$. It is straightforward to check that any solution of
Eqs. (19)-(24) with $\omega_{\bf k}^{*}=-\omega_{-{\bf k}}$ will also satisfy
$(p^{a}_{i,{\bf k}})^{*}=p^{a}_{i,-{\bf k}}$, i.e. Eqs. (19)-(24) are
consistent with the reality conditions in (18).

Figure 1: The triad of orthonormal vectors
$\\{{\hat{k}},{\hat{n}},{\hat{\xi}}\\}$ and the direction of the electric
field along ${\hat{z}}$. The vectors ${\hat{k}}$, ${\hat{\xi}}$ and
${\hat{z}}$ are in a plane while ${\hat{n}}$ is normal to the plane.

The variables ${\bf p}^{a}$ have a natural decomposition in a basis of spatial
vectors $\\{{\hat{k}},{\hat{n}},{\hat{\xi}}\\}$ (see Fig. 1) where

$\hat{k}=\frac{{\bf k}}{k},\ \ \hat{n}=\hat{\xi}\times\hat{k},\ \
\hat{\xi}=\frac{\hat{z}-c\hat{k}}{s}$ (25)

For convenience, we have denoted $c\equiv\hat{k}\cdot\hat{z}=\cos\theta$ and
$s=|\hat{z}\times\hat{k}|=\sin\theta$ where $\theta$ is the angle between
$\hat{k}$ and $\hat{z}$. Then we have the useful relation

$\hat{z}=s\hat{\xi}+c\hat{k}.$ (26)

Next write

${\bf p}^{a}=\alpha_{a}\hat{k}+\beta_{a}\hat{n}+\gamma_{a}\hat{\xi}$ (27)

The $\\{\alpha_{a}\\}$ modes, with polarization in the $\hat{k}$ direction,
are longitudinal. The $\\{\beta_{a}\\}$ modes are polarized in the $\hat{n}$
direction and are transverse and also polarized orthogonal to the electric
field. The $\\{\gamma_{a}\\}$ modes are transverse and polarized in the
$\hat{\xi}$ direction which is at an angle $\pi/2-\theta$ to the electric
field. We shall call $\\{\alpha_{a}\\}$ the “longitudinal” (L) modes, the
$\\{\beta_{a}\\}$ as the “transverse-orthogonal” (TO) modes, and the
$\\{\gamma_{a}\\}$ as the “transverse-nonorthogonal” (TN) modes. Note that for
$\theta=0$, the TN modes are also orthogonal to the electric field and
coincide with TO modes.

The constraints and equations of motion in (19)-(21) and (22)-(24) can be
written in terms of the 9 functions $\\{\alpha_{a},\beta_{a},\gamma_{a}\\}$.
The constraints are,

$\displaystyle\omega k\alpha_{1}-i\Omega k\alpha_{2}+\epsilon\Omega
c\alpha_{3}+\epsilon\Omega s\gamma_{3}=0$ (28) $\displaystyle\omega
k\alpha_{2}+i\Omega k\alpha_{1}-i\epsilon\omega c\alpha_{3}-i\epsilon\omega
s\gamma_{3}=0$ (29) $\displaystyle\hskip-28.45274pt\omega
k\alpha_{3}+i\epsilon\omega c\alpha_{2}+i\epsilon\omega
s\gamma_{2}-2\epsilon\Omega c\alpha_{1}-2\epsilon\Omega s\gamma_{1}=0$ (30)

Note that the constraints do not involve the $\\{\beta_{a}\\}$ functions.

The equations of motion are

$\displaystyle(-\omega^{2}-\Omega^{2})\alpha_{1}+i2\omega\Omega\alpha_{2}=0$
(31)
$\displaystyle\hskip-35.56593pt(-\omega^{2}-\Omega^{2}+\epsilon^{2}s^{2})\alpha_{2}-i2\omega\Omega\alpha_{1}-\epsilon^{2}cs\gamma_{2}$
$\displaystyle\hskip 113.81102pt+i\epsilon ks\gamma_{3}=0$ (32)
$\displaystyle(-\omega^{2}+\epsilon^{2}s^{2})\alpha_{3}-i\epsilon
ks\gamma_{2}-\epsilon^{2}cs\gamma_{3}=0$ (33)
$\displaystyle(-\omega^{2}+k^{2}-\Omega^{2})\beta_{1}+i2\omega\Omega\beta_{2}=0$
(34)
$\displaystyle\hskip-21.33955pt(-\omega^{2}+k^{2}-\Omega^{2}+\epsilon^{2})\beta_{2}-i2\omega\Omega\beta_{1}-i2\epsilon
kc\beta_{3}=0$ (35)
$\displaystyle(-\omega^{2}+k^{2}+\epsilon^{2})\beta_{3}+i2\epsilon
kc\beta_{2}=0$ (36)
$\displaystyle(-\omega^{2}+k^{2}-\Omega^{2})\gamma_{1}+i2\omega\Omega\gamma_{2}=0$
(37)
$\displaystyle\hskip-21.33955pt(-\omega^{2}+k^{2}-\Omega^{2}+\epsilon^{2}c^{2})\gamma_{2}-i2\omega\Omega\gamma_{1}+i\epsilon
ks\alpha_{3}$ $\displaystyle\hskip
85.35826pt-\epsilon^{2}cs\alpha_{2}-i2\epsilon kc\gamma_{3}=0$ (38)
$\displaystyle\hskip-28.45274pt(-\omega^{2}+k^{2}+\epsilon^{2}c^{2})\gamma_{3}-i\epsilon
ks\alpha_{2}+i2\epsilon kc\gamma_{2}$ $\displaystyle\hskip
113.81102pt-\epsilon^{2}sc\alpha_{3}=0$ (39)

The $\\{\beta_{a}\\}$ functions do not appear in the constraint equations, nor
do they depend on the $\\{\alpha_{a},\gamma_{a}\\}$. Hence they can be treated
separately. In the next subsection we will first consider the
$\\{\beta_{a}\\}$ problem and in the following subsection come to the more
complicated $\\{\alpha_{a},\gamma_{a}\\}$ problem.

##  V TO modes ($\beta_{a}$)

The equations for $\beta_{a}$ can be written as a matrix equation: $MX=0$,

$M=\begin{pmatrix}-\kappa^{2}-\Omega^{2}&i2\omega\Omega&0\\\
-i2\omega\Omega&-\kappa^{2}-\Omega^{2}+\epsilon^{2}&-i2\epsilon ck\\\
0&+i2\epsilon ck&-\kappa^{2}+\epsilon^{2}\end{pmatrix}$ (40)

where $\kappa^{2}\equiv\omega^{2}-k^{2}$ and
$X^{T}=(\beta_{1},\beta_{2},\beta_{3})$.

To find the eigenvalues, we set the determinant of the matrix to zero. This
yields a cubic equation in $\lambda\equiv\omega^{2}$,

$\displaystyle\hskip-25.6073ptP(\lambda)\equiv\lambda^{3}-\lambda^{2}(3k^{2}+2\epsilon^{2}+2\Omega^{2})$
$\displaystyle\hskip
0.0pt+\lambda(3k^{4}+4\epsilon^{2}s^{2}k^{2}+\epsilon^{4}+\epsilon^{2}\Omega^{2}+\Omega^{4})$
$\displaystyle\hskip-25.6073pt-(k^{2}-\Omega^{2})[(k^{2}+\epsilon^{2}-\Omega^{2})(k^{2}+\epsilon^{2})-4\epsilon^{2}c^{2}k^{2}]=0$
(41)

The cubic equation can be solved explicitly to obtain the eigenvalues, however
the expressions are opaque. We get more insight by considering a different
approach.

The cubic equation in (41) will have three roots and can be written as

$(\lambda-\lambda_{1})(\lambda-\lambda_{2})(\lambda-\lambda_{3})=0$ (42)

Note that the roots $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ for ${\bf
k}$ and $-{\bf k}$ are identical since $k^{2}$, $c^{2}$ and $s^{2}$ are
unchanged due to the sign flip. Hence, for example, $\lambda_{1,{\bf
k}}=\lambda_{1,-{\bf k}}$. Together with the reality condition of (18) this
relation implies,

$\omega_{1,{\bf k}}^{2}=\omega_{1,-{\bf k}}^{2}=(\omega_{1,{\bf k}}^{*})^{2}.$
(43)

Hence eigenfrequencies of physical modes satisfy

$\omega_{\bf k}=\pm\omega_{\bf k}^{*}.$ (44)

i.e. physical eigenfrequencies are purely real or purely imaginary. In terms
of $\lambda$, only the real roots of (41) are of physical interest.

Next consider the polynomial as in (41) but without the $\lambda$ independent
term,

$\displaystyle{\tilde{P}}(\lambda)\equiv\lambda^{3}-\lambda^{2}(3k^{2}+2\epsilon^{2}+2\Omega^{2})$
$\displaystyle\hskip
56.9055pt+\lambda(3k^{4}+4\epsilon^{2}s^{2}k^{2}+\epsilon^{4}+\epsilon^{2}\Omega^{2}+\Omega^{4})$

Then,

${\tilde{P}}(\lambda)\equiv\lambda(\lambda-\lambda_{+})(\lambda-\lambda_{-})$
(45)

where $\lambda_{\pm}$ are obtained by solving a quadratic that involves ${\bf
k}$ (as $k^{2}$ and $s^{2}$) and the parameters $\epsilon$ and $\Omega$. We
can check that the real parts of all three roots of ${\tilde{P}}$ are non-
negative. Therefore ${\tilde{P}}$ has the shape shown in Fig. 2.

Figure 2: The cubic curve ${\tilde{P}}(\lambda)$ (the middle curve) and the
cubics $P(\lambda)$ for $C>0$ (upper curve) and for $C<0$ (lower curve). The
zero root of ${\tilde{P}}$ shifts to negative $\lambda$ for $C>0$ and to
positive $\lambda$ for $C<0$. The root can become complex for sufficiently
large and negative $C$.

Figure 3: Stability plots for the TO modes for $\epsilon^{2}=\Omega^{2}/2$
(top), $\epsilon^{2}=\Omega^{2}$ (middle) and $\epsilon^{2}=3\Omega^{2}/2$
(bottom). As $\epsilon$ is decreased the two triangular regions in the
$\epsilon^{2}<\Omega^{2}$ plot (top) shrink and approach the
$k^{2}=\Omega^{2}$ vertical line as $\epsilon^{2}\to 0$.

Next let us return to the cubic in (41) which can be written as,

$P(\lambda)={\tilde{P}}(\lambda)+C$ (46)

where

$\displaystyle C$ $\displaystyle=$
$\displaystyle-(k^{2}-\Omega^{2})[(k^{2}+\epsilon^{2}-\Omega^{2})(k^{2}+\epsilon^{2})-4\epsilon^{2}c^{2}k^{2}]$
(47) $\displaystyle=$
$\displaystyle-(\kappa-\Omega^{2})(\kappa-\kappa_{+})(\kappa-\kappa_{-})$

where $\kappa\equiv k^{2}$,

$\displaystyle\kappa_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\biggl{[}\Omega^{2}+2\epsilon^{2}c_{2}$ (48)
$\displaystyle\hskip
14.22636pt\pm\sqrt{\left(\Omega^{2}+2\epsilon^{2}c_{2}\right)^{2}-4\epsilon^{2}(\epsilon^{2}-\Omega^{2})}\
\biggr{]}$

and $c_{2}\equiv\cos(2\theta)$. If $C>0$, the ${\tilde{P}}$ curve shifts
upwards (see Fig. 2) and the $\lambda=0$ root of ${\tilde{P}}$ shifts to the
left, i.e. $\lambda<0$. This indicates an instability. On the other hand, if
$C<0$, there is no instability.

It is simplest to analyze the case with $\epsilon^{2}=\Omega^{2}$ for then

$\kappa_{\pm}=\frac{\Omega^{2}}{2}\left[(4c^{2}-1)\pm|4c^{2}-1|\right]$ (49)

For $c^{2}>1/4$, $\kappa_{-}=0$ and $\kappa_{+}=\Omega^{2}(4c^{2}-1)$. Then
for $c^{2}>1/2$ we obtain an instability ($C>0$) for
$\Omega^{2}<k^{2}<\Omega^{2}(4c^{2}-1)$ and for $1/4<c^{2}<1/2$ modes are
unstable if $\Omega^{2}(4c^{2}-1)<k^{2}<\Omega^{2}$. For the opposite case of
$c^{2}<1/4$, $\kappa_{-}=4c^{2}-1<0$ and $\kappa_{+}=0$. Then the instability
occurs for $0<k^{2}<\Omega^{2}$. This explains the instability domains shown
in Fig. 3 in the $\epsilon^{2}=\Omega^{2}$ case.

The analysis and results are similar for $\epsilon^{2}<\Omega^{2}$ as is clear
from Fig. 3. With $c=0$, we find that $C>0$ for
$\Omega^{2}-\epsilon^{2}<k^{2}<\Omega^{2}$. Essentially, as $\epsilon^{2}$ is
reduced the domain of instability shrinks towards a vertical line along
$k^{2}=\Omega^{2}$. An interesting limit is $\epsilon\to 0$ and
$\Omega\to\infty$ with $\epsilon\Omega=E$ held constant. In this limit, there
are no unstable TO modes.

The analysis is a bit more involved for $\epsilon^{2}>\Omega^{2}$. First
consider $c^{2}=0$ ($c_{2}=-1$). Then (48) gives
$\kappa_{\pm}=-\epsilon^{2},\,-\epsilon^{2}+\Omega^{2}$ and both roots are for
negative $\kappa$. From (47) we then see that $C>0$ for $0<k^{2}<\Omega^{2}$
as is also seen in Fig. 3. Next consider $c^{2}=1$ ($c_{2}=+1$). Then (48)
gives

$\kappa_{\pm}=\frac{1}{2}\biggl{[}\Omega^{2}+2\epsilon^{2}\pm\sqrt{\Omega^{4}+8\epsilon^{2}\Omega^{2}}\
\biggr{]}.$ (50)

Now $\kappa_{+}>\Omega^{2}$ but $\kappa_{-}$ may be larger or smaller than
$\Omega^{2}$. If $\kappa_{-}<\Omega^{2}$, then from (47) we see that $C>0$ for
$0<k^{2}<\kappa_{-}$ and also for $\Omega^{2}<k^{2}<\kappa_{+}$ but $C<0$ for
$\kappa_{-}<k^{2}<\Omega^{2}$. Hence we see the gap in the instability domain
in Fig. 3 at $c^{2}=1$. If, on the other hand, $\kappa_{-}>\Omega^{2}$, then
$C>0$ for $0<k^{2}<\Omega^{2}$ and then again for
$\kappa_{-}<k^{2}<\kappa_{+}$.

The shapes of the unstable regions can be understood in terms of the roots of
$C$ in (47), namely 1 and $\kappa_{\pm}$. For example, in the case
$\epsilon^{2}>\Omega^{2}$, for $c^{2}=0$, there is one real root
($\kappa=\Omega^{2}$) and $\kappa_{\pm}$ are negative. As $c^{2}$ increases,
$\kappa_{\pm}$ become complex. At some critical value of $c^{2}$ the imaginary
part of $\kappa_{\pm}$ vanishes. This is at the minimum of the parabolic shape
in the plot of Fig. 3 and can occur for $k^{2}<\Omega^{2}$ or
$k^{2}>\Omega^{2}$ depending on the value of $\epsilon^{2}/\Omega^{2}$. The
left edge of the parabola is given by $\kappa_{-}$, and the right edge is
given by $\kappa_{+}$. For yet larger $c^{2}$, $\kappa_{+}$ becomes larger
than $\Omega^{2}$. The unstable region is when two of the factors in (47) are
positive and one is negative.

##  VI TN and L modes ($\\{\alpha_{a},\gamma_{a}\\}$)

Eqs. 31-33 and 37-39 can be written in matrix form $MX=0$ with
$X^{T}=(\alpha_{1},\alpha_{2},\alpha_{3},\gamma_{1},\gamma_{2},\gamma_{3})$
and

$M=\begin{pmatrix}-\omega^{2}-\Omega^{2}&i2\omega\Omega&0&0&0&0\\\
-i2\omega\Omega&-\omega^{2}-\Omega^{2}+\epsilon^{2}s^{2}&0&0&-\epsilon^{2}cs&i\epsilon
sk\\\ 0&0&-\omega^{2}+\epsilon^{2}s^{2}&0&-i\epsilon sk&-\epsilon^{2}cs\\\
0&0&0&-\omega^{2}+k^{2}-\Omega^{2}&i2\omega\Omega&0\\\
0&-\epsilon^{2}cs&i\epsilon
sk&-i2\omega\Omega&-\omega^{2}+k^{2}-\Omega^{2}+\epsilon^{2}c^{2}&-i2\epsilon
ck\\\ 0&-i\epsilon sk&-\epsilon^{2}cs&0&i2\epsilon
ck&-\omega^{2}+k^{2}+\epsilon^{2}c^{2}\end{pmatrix}$ (51)

The eigenvectors are also required to satisfy the constraints in (28)-(30) as
we will discuss further after considering some special cases.

###  VI.1 Special case: $c=1$

For $c=1$, the matrix in (51) becomes block diagonal in 3 blocks, the first is
the $\\{\alpha_{1},\alpha_{2}\\}$ $2\times 2$ block with two degenerate
eigenvalues $\omega^{2}=\Omega^{2}$, the second block is the
$\\{\alpha_{3}\\}$ $1\times 1$ block with eigenvalue $\omega^{2}=0$, and the
third $\\{\gamma_{1},\gamma_{2},\gamma_{3}\\}$ block is given by the matrix in
(40) with $c=1$. Then the analysis in Sec. V for the TO modes applies
immediately (with $c=1$). This is expected since for $c=1$ there is no
distinction between TO and TN modes.

The constraint equations with $c=1$ read,

$\displaystyle\omega k\alpha_{1}-i\Omega
k\alpha_{2}+\epsilon\Omega\alpha_{3}=0$ (52) $\displaystyle\omega
k\alpha_{2}+i\Omega k\alpha_{1}-i\epsilon\omega\alpha_{3}=0$ (53)
$\displaystyle\omega
k\alpha_{3}+i\epsilon\omega\alpha_{2}-2\epsilon\Omega\alpha_{1}=0$ (54)

Note that the $\gamma_{a}$ are unconstrained. On the other hand, (52)-(54)
over-constrain the eigensolutions in the $\\{\alpha_{1},\alpha_{2}\\}$ sector
with $\omega=\pm\Omega$, and in the $\\{\alpha_{3}\\}$ sector with $\omega=0$,
and neither of these two eigensolutions are physically admissible.

###  VI.2 Special case: $c=0$

With $c=0$, we have $s=1$, and the matrix in (51) becomes block diagonal in
the $\\{\alpha_{1},\alpha_{2},\gamma_{3}\\}$ and the
$\\{\alpha_{3},\gamma_{1},\gamma_{2}\\}$ blocks. The $3\times 3$ matrix for
the first block is

$M_{1}=\begin{pmatrix}-\omega^{2}-\Omega^{2}&i2\omega\Omega&0\\\
-i2\omega\Omega&-\omega^{2}-\Omega^{2}+\epsilon^{2}&i\epsilon k\\\
0&-i\epsilon k&-\omega^{2}+k^{2}\end{pmatrix}$ (55)

with constraint,

$\displaystyle\omega k\alpha_{1}-i\Omega
k\alpha_{2}+\epsilon\Omega\gamma_{3}=0$ (56) $\displaystyle\omega
k\alpha_{2}+i\Omega k\alpha_{1}-i\epsilon\omega\gamma_{3}=0$ (57)

and the matrix for the second block is

$\hskip-9.95863ptM_{2}=\begin{pmatrix}-\omega^{2}+\epsilon^{2}&0&-i\epsilon
k\\\ 0&-\omega^{2}+k^{2}-\Omega^{2}&i2\omega\Omega\\\ i\epsilon
k&-i2\omega\Omega&-\omega^{2}+k^{2}-\Omega^{2}\end{pmatrix}$ (58)

with constraint,

$\omega k\alpha_{3}+i\epsilon\omega\gamma_{2}-2\epsilon\Omega\gamma_{1}=0.$
(59)

We now discuss these $3\times 3$ blocks separately.

####  VI.2.1 $\\{\alpha_{1},\alpha_{2},\gamma_{3}\\}$ block

In this block, gauge fields of the first two colors are oscillating in the
longitudinal direction, whereas the third has amplitude in the transverse
direction and orthogonal to the background electric field.

A straight-forward procedure would be to first solve the eigenproblem for
$M_{1}$ and then check for the eigenvectors that satisfy the constraints.
However we find it simpler to first solve the constraints (56)-(57) and then
deal with the eigenproblem.

The constrains in (56)-(57) can be used to eliminate two of the three
variables, say $\alpha_{2}$ and $\gamma_{3}$, while the third variable can be
absorbed in the normalization of the resulting eigenvector. Hence we seek an
eigenvector of the form,

$V_{1}\equiv\begin{pmatrix}\alpha_{1}\\\ \alpha_{2}\\\
\gamma_{3}\end{pmatrix}=\begin{pmatrix}2\epsilon\omega\Omega\\\
-i\epsilon(\omega^{2}+\Omega^{2})\\\ k(-\omega^{2}+\Omega^{2})\end{pmatrix}$
(60)

Insertion in $M_{1}V_{1}=0$ shows that there is no solution for $\omega$ for
$\epsilon\neq 0$, $k\neq 0$. Hence these modes are over-constrained and
absent. For $k=0$, $M_{1}V_{1}=0$ gives

$(\omega^{2}-\Omega^{2})^{2}-\epsilon^{2}(\omega^{2}+\Omega^{2})=0$ (61)

which has the roots

$\omega_{\pm}^{2}=\Omega^{2}+\frac{\epsilon^{2}}{2}\pm\sqrt{\left(\Omega^{2}+\frac{\epsilon^{2}}{2}\right)^{2}+\Omega^{2}(\epsilon^{2}-\Omega^{2})}$
(62)

Therefore $\omega_{-}^{2}<0$ if and only if $\Omega^{2}<\epsilon^{2}$ and
$k=0$.

####  VI.2.2 $\\{\alpha_{3},\gamma_{1},\gamma_{2}\\}$ block

In this block, the gauge field of the third color oscillates in the
longitudinal direction, whereas the first two colors oscillate transversely
and orthogonal to the background electric field.

Now the constraint (59) reduces the eigenvector to be of the form,

$V_{2}\equiv\begin{pmatrix}\alpha_{3}\\\ \gamma_{1}\\\
\gamma_{2}\end{pmatrix}=\begin{pmatrix}2\epsilon\Omega\alpha_{3}\\\
\omega(k\alpha_{3}+i\epsilon\gamma_{2})\\\
2\epsilon\Omega\gamma_{2}\end{pmatrix}.$ (63)

Imposing $M_{2}V_{2}=0$ we find the solution $\omega^{2}=\epsilon^{2}$
provided $k^{2}=\epsilon^{2}+\Omega^{2}$, and this mode is stable. For $k=0$
there is a solution with $\omega^{2}=\epsilon^{2}$.

An unusual feature of the first of these two modes is that it has non-trivial
spatial dependence but it exists only for a fixed value of
$k=\sqrt{\epsilon^{2}+\Omega^{2}}$ and the direction of ${\bf k}$ is
perpendicular to the electric field, while the mode is polarized along the
electric field and in the longitudinal direction. This mode represents
fluctuations in the homogeneity of the electric field but with a definite
wavelength.

###  VI.3 Special case: $k^{2}\gg\Omega^{2},\epsilon^{2}$

We now consider the ultraviolet limit $k^{2}\gg\Omega^{2},\epsilon^{2}$. The
constraints (28)-(30) now give

$\alpha_{1}=0=\alpha_{2}=\alpha_{3}$ (64)

and only the $\\{\gamma_{a}\\}$ represent physical modes with the dispersion
relation $\omega^{2}=k^{2}$. There are no unstable modes in this limit.

###  VI.4 Special case: $k\to 0$

The problem simplifies in the $k\to 0$ limit as the
$\\{\alpha_{3},\gamma_{3}\\}$ block decouples from the
$\\{\alpha_{1},\alpha_{2},\gamma_{1},\gamma_{2}\\}$ block.

The $2\times 2$ matrix for the $\\{\alpha_{3},\gamma_{3}\\}$ block is

$M_{3}=\begin{pmatrix}-\omega^{2}+\epsilon^{2}s^{2}&-\epsilon^{2}cs\\\
-\epsilon^{2}cs&-\omega^{2}+\epsilon^{2}c^{2}\end{pmatrix}$ (65)

and the constraint reduces to $c\alpha_{3}+s\gamma_{3}=0$. $M_{3}$ has
eigenvalues $\omega^{2}=0$ and $\omega^{2}=\epsilon^{2}$ but only the latter
is consistent with the constraint. Thus there are no unstable modes in the
$\\{\alpha_{3},\gamma_{3}\\}$ block.

The $4\times 4$ matrix for the
$\\{\alpha_{1},\alpha_{2},\gamma_{1},\gamma_{2}\\}$ block is

$M_{4}=\begin{pmatrix}-\omega^{2}-\Omega^{2}&i2\omega\Omega&0&0\\\
-i2\omega\Omega&-\omega^{2}-\Omega^{2}+\epsilon^{2}s^{2}&0&-\epsilon^{2}cs\\\
0&0&-\omega^{2}-\Omega^{2}&i2\omega\Omega\\\
0&-\epsilon^{2}cs&-i2\omega\Omega&-\omega^{2}-\Omega^{2}+\epsilon^{2}c^{2}\end{pmatrix}$
(66)

and the constraint is

$i\omega(c\alpha_{2}+s\gamma_{2})-2\Omega(c\alpha_{1}+s\gamma_{1})=0.$ (67)

We solve (31) and (37) with $k=0$ to get

$\alpha_{2}=-i\frac{(\omega^{2}+\Omega^{2})}{2\omega\Omega}\alpha_{1},\ \
\gamma_{2}=-i\frac{(\omega^{2}+\Omega^{2})}{2\omega\Omega}\gamma_{1},$ (68)

which, together with (67), gives

$(\omega^{2}-3\Omega^{2})(c\alpha_{1}+s\gamma_{1})=0.$ (69)

Therefore to satisfy the constraint we must either have
$\omega^{2}=3\Omega^{2}>0$ or $c\alpha_{1}+s\gamma_{1}=0$.

Evaluation of the determinant of $M_{4}$ on Mathematica gives,

${\rm
Det}(M_{4})=(\omega^{2}-\Omega^{2})^{2}[(\omega^{2}-\Omega^{2})^{2}-\epsilon^{2}(\omega^{2}+\Omega^{2})]$
(70)

This has the root $\omega^{2}=3\Omega^{2}$ but only if
$\Omega^{2}=\epsilon^{2}$. In any case, $\omega^{2}=3\Omega^{2}>0$ and implies
a stable mode. So we now focus on the other case, namely

$c\alpha_{1}+s\gamma_{1}=0.$ (71)

Combining (71) with (68), and ignoring an overall normalization factor, the
Gauss constraint forces us to only consider the eigenvector,

$\hskip-5.69046ptV_{4}^{T}=\left(2\omega\Omega
s,-i{(\omega^{2}+\Omega^{2})}s,-2\omega\Omega
c,i{(\omega^{2}+\Omega^{2})}c\right).$ (72)

Requiring $M_{4}V_{4}=0$ leads once again to (61) and to the roots in (62).
Therefore $\omega_{-}^{2}<0$ if and only if $\Omega^{2}<\epsilon^{2}$ and
$k=0$ and the unstable eigenmode can be found by setting $\omega=\omega_{-}$
in (72).

###  VI.5 Special case: $\Omega=E/\epsilon$, $\epsilon\to 0$

In Sec. V we have seen that there are no unstable TO modes with
$\Omega=E/\epsilon$ and $\epsilon\to 0$. Now we consider the TN and L modes in
this regime.

With $\epsilon\to 0$, the matrix $M$ in (51) takes on a simple block diagonal
form. The $\\{\alpha_{1},\alpha_{2}\\}$ block has two degenerate eigenvalues
$\omega^{2}=\Omega^{2}$; the $\\{\alpha_{3}\\}$ block has eigenvalue
$\omega^{2}=0$; the $\\{\gamma_{1},\gamma_{2}\\}$ block has eigenvalues
$\omega^{2}=(k\pm\Omega)^{2}$; and the $\\{\gamma_{3}\\}$ block has eigenvalue
$\omega^{2}=k^{2}$. The corresponding eigenvectors can be inserted into
equations (28)-(30) to check if the Gauss constraints are satisfied. However,
since none of the eigenvalues for $\omega^{2}$ are negative, it is clear that
there are no unstable TN and L modes for these limiting parameters.

##  VII Conclusions

We have considered the stability of a homogeneous electric field background in
pure SU(2) gauge theory. The gauge fields underlying the electric field are
taken to be of the form in (5) and not of the Maxwell type:
$A_{i}^{a}=-Et\delta^{a3}\delta_{iz}$. This is because gauge fields of the
Maxwell type are unstable to Schwinger pair production while the gauge fields
in (5) are protected from decay due to this process Vachaspati (2022).
However, the gauge fields in (5) are not solutions of the vacuum classical
equations of motion; instead non-vanishing currents are necessary. There are
two ways to explain these non-vanishing currents. The first is that they are
due to classical external sources in which case they are simply postulated.
The second way is that the classical equations of motion should be replaced by
equations that take quantum effects into account and these “effective
classical equations” can contain sources. For example, in the semiclassical
approximation quantum fluctuations provide current sources for the background
Vachaspati (2022),

$\displaystyle j^{\mu a}$ $\displaystyle\equiv$
$\displaystyle\epsilon^{abc}\langle\partial_{\nu}q^{\nu b}q^{\mu c}-q^{\nu
b}\partial^{\mu}q_{\nu}^{c}+2q^{\nu b}\partial_{\nu}q^{\mu c}\rangle_{R}$ (73)
$\displaystyle+A_{\nu}^{b}\langle q^{\nu a}q^{\mu b}-2q^{\nu b}q^{\mu
a}\rangle_{R}+A_{\nu}^{a}\langle q^{\nu b}q^{\mu b}\rangle_{R}$
$\displaystyle+A^{\mu b}\langle q^{b}_{\nu}q^{\nu a}\rangle_{R}-A^{\mu
a}\langle q_{\nu}^{b}q^{\nu b}\rangle_{R}$

where $q^{a}_{\mu}$ are the quantum fluctuations in the background
$A^{a}_{\mu}$ and $\langle\cdot\rangle_{R}$ denotes a renormalized expectation
value taken in the quantum state of $q^{a}_{\mu}$. For stable modes, the
quantum state might be given by simple harmonic oscillator states for each of
the eigenmodes of $q^{a}_{\mu}$. However the quantum state of unstable modes
will not be simple harmonic oscillator states which is why it is important to
perform a stability analysis. We will comment further on the unstable modes
after summarizing our results.

The gauge field background in (5) is described by two parameters: $\epsilon$
and $\Omega$. The electric field strength is given by $E=\epsilon\Omega$. The
results of the fluctuation analysis depend on whether
$\epsilon^{2}>\Omega^{2}$ or $\epsilon^{2}\leq\Omega^{2}$. The fluctuations
naturally split into “TO modes”” that are transverse to the wavevector ${\bf
k}$ and orthogonal to the background electric field, “TN modes” that are
transverse to ${\bf k}$ but not orthogonal to the electric field, and “L
modes” that are in the longitudinal direction.

The TO modes decouple from the TN and L modes. The stability analysis of Sec.
V shows that TO modes are stable except in a range of $k^{2}$ that depends on
the angle $\theta$ between the electric field and the wavevector ${\bf k}$.
The instability regions depend on the background parameters and are plotted in
Fig. 3. There are two important results emerging from our analysis. The first
is that the region of parameter space $(k^{2},c^{2})$ ($c=\cos\theta$) where
unstable modes exist depends on the relation between $\epsilon^{2}$ and
$\Omega^{2}$. The instability region is smaller when $\epsilon^{2}<\Omega^{2}$
and shrinks to zero as $\epsilon^{2}\to 0$. Note that the electric field
strength is given by $E=\epsilon\Omega$ and can be held fixed in the limit by
also taking $\Omega\to\infty$. The second is that unstable modes exist only
for small values of $k^{2}$. For example, for $\epsilon=\Omega$, there are no
unstable modes for $k^{2}>3\epsilon^{2}$ for any value of $c^{2}$.

The TN and L modes are coupled in general and the analysis is more involved
than for the TO modes. In Sec. VI we discuss the stability of these modes in
various parameter regimes. The special cases of $\theta=0$ and $\theta=\pi/2$
are considered. For $\theta=0$ the analysis is identical to that of TO modes,
while for $\theta=\pi/2$ there is an instability if $\epsilon^{2}>\Omega^{2}$
and $k=0$. There is also a special stable mode that corresponds to
oscillations of the background electric field orthogonal to its direction,
similar to a sound wave. We have also considered the special case of large
$k^{2}$ and here the modes are simply those of free massless waves with
dispersion $\omega^{2}=k^{2}$. Finally, we examine the $\epsilon\to 0$ limit
with $E=\epsilon\Omega$ held fixed and show that there are no unstable TN and
L modes, just as there are no unstable TO modes in this limit.

As mentioned in Sec. I, we were motivated to perform this stability analysis
because confining strings in QCD are expected to be stable. The electric
fields we have considered as backgrounds do not excite Schwinger pair
production but, as we have seen, have classical instabilities for certain
infrared modes. How do these classical instabilities impact the possibility
that the electric fields we have considered are responsible for confining
strings? The first point we note is that there are no instabilities in the
limit of $\epsilon\to 0$ and $E=\epsilon\Omega$ fixed. So it could be that the
electric field in a confining string corresponds to this set of parameters.
Then there are no unstable modes and the quantum state for each mode is that
of a simple harmonic oscillator. The second point is that the instabilities we
have found are for a homogeneous electric field and only occur for small
values of $k^{2}$, ($k^{2}\lesssim\epsilon^{2}$ for $\epsilon^{2}>\Omega^{2}$)
that is, on large length scales. In contrast, the electric flux in a string
only has support in a finite area – the string cross-section – and we do not
expect any unstable modes on length scales larger than the thickness of the
string. (Though there is still the question of the infinite extent of the
string along the electric field direction and whether the instabilities for
$\theta=0$ will survive.) It would be worthwhile performing an explicit
stability analysis for a flux tube configuration such as Vachaspati (2022),

$A^{\pm}_{\mu}=-\epsilon e^{\pm i\Omega t}\,f(r)\,\partial_{\mu}z,\ \
A^{3}_{\mu}=0,$ (74)

where $f(r)$ is a profile function for the string and
$r\equiv\sqrt{x^{2}+y^{2}}$ is the cylindrical radial coordinate. Another
interesting question is if the homogeneous electric field background we have
considered is unstable towards forming an Abrikosov lattice Abrikosov (1957)
of electric flux tubes. After all we have identified certain unstable modes
with spatial dependence that is orthogonal to the background homogeneous
electric field.

###### Acknowledgements.

We are grateful to Jan Ambjorn, Guy Moore, Stanislaw Mrówczyński, Igor
Shovkovy and Raju Venugopalan for comments. T. V. thanks the CCPP at NYU for
hospitality where some of this work was done. This work was supported by the
U.S. Department of Energy, Office of High Energy Physics, under Award No. DE-
SC0019470.

## References

  * Matinyan and Savvidy (1978) S. G. Matinyan and G. K. Savvidy, Nucl. Phys. B 134, 539 (1978).
  * Brown and Weisberger (1979) L. S. Brown and W. I. Weisberger, Nucl. Phys. B 157, 285 (1979), [Erratum: Nucl.Phys.B 172, 544 (1980)].
  * Yildiz and Cox (1980) A. Yildiz and P. H. Cox, Phys. Rev. D 21, 1095 (1980).
  * Ambjorn and Hughes (1982a) J. Ambjorn and R. J. Hughes, Nucl. Phys. B 197, 113 (1982a).
  * Ambjorn and Hughes (1982b) J. Ambjorn and R. J. Hughes, Phys. Lett. B 113, 305 (1982b).
  * Cooper and Nayak (2006) F. Cooper and G. C. Nayak, Phys. Rev. D 73, 065005 (2006), eprint hep-ph/0511053.
  * Nayak and van Nieuwenhuizen (2005) G. C. Nayak and P. van Nieuwenhuizen, Phys. Rev. D 71, 125001 (2005), eprint hep-ph/0504070.
  * Cooper et al. (2008) F. Cooper, J. F. Dawson, and B. Mihaila, in _Conference on Nonequilibrium Phenomena in Cosmology and Particle Physics_ (2008), eprint 0806.1249.
  * Nair and Yelnikov (2010) V. P. Nair and A. Yelnikov, Phys. Rev. D 82, 125005 (2010), eprint 1005.2582.
  * Ragsdale and Singleton (2017) M. Ragsdale and D. Singleton, J. Phys. Conf. Ser. 883, 012014 (2017), eprint 1708.09753.
  * Cardona and Vachaspati (2021) C. Cardona and T. Vachaspati, Phys. Rev. D 104, 045009 (2021), eprint 2105.08782.
  * Vachaspati (2022) T. Vachaspati, Phys. Rev. D 105, 105011 (2022), eprint 2204.01902.
  * Bazak and Mrowczynski (2022a) S. Bazak and S. Mrowczynski, Phys. Rev. D 105, 034023 (2022a), eprint 2111.11396.
  * Bazak and Mrowczynski (2022b) S. Bazak and S. Mrowczynski (2022b), eprint 2205.08282.
  * Berges et al. (2021) J. Berges, M. P. Heller, A. Mazeliauskas, and R. Venugopalan, Rev. Mod. Phys. 93, 035003 (2021), eprint 2005.12299.
  * Abrikosov (1957) A. A. Abrikosov, Sov. Phys. JETP 5, 1174 (1957).