# Interference induced anisotropy in a two-dimensional dark state optical
lattice

E. Gvozdiovas  Institute of Theoretical Physics and Astronomy, Vilnius
University, Saulėtekio Ave. 3, LT-10257 Vilnius, Lithuania  I. B. Spielman
<EMAIL_ADDRESS>http://ultracold.jqi.umd.edu Joint Quantum Institute,
University of Maryland, College Park, Maryland 20742-4111, USA  National
Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA  G.
Juzeliūnas<EMAIL_ADDRESS>Institute of Theoretical Physics
and Astronomy, Vilnius University, Saulėtekio Ave. 3, LT-10257 Vilnius,
Lithuania

###### Abstract

We describe a two-dimensional optical lattice for ultracold atoms with spatial
features below the diffraction limit created by a bichromatic optical standing
wave. At every point in space these fields couple the internal atomic states
in a three-level Lambda coupling configuration. Adiabatically following the
local wavefunction of the resulting dark state yields a spatially uniform
Born-Oppenheimer potential augmented by geometric scalar and vector potentials
appearing due to spatially rapid changes of the wavefunction. Depending on
system parameters, we find that the geometric scalar potential can interpolate
from a 2D analogue of the Kronig-Penney lattice, to an array of tubes with a
zig-zag shaped barrier. The geometric vector potential induces a spatially
periodic effective magnetic field (the Berry’s curvature) that can be tuned to
cause destructive interference between neighboring tubes, thereby decoupling
them at a critical point in parameter space. We numerically investigate the
energy spectrum including decay from the excited state, and find that the
adiabatic approximation is sound for strong coupling strengths, leading to
negligible loss in the dark state manifold. Furthermore, the spectrum is well-
described by a non-Hermitian tight binding model with on-site losses, and
hopping characterized by both loss and, surprisingly, gain.

Realizing long-lived strongly correlated quantum matter with ultracold atoms
in optical lattices is an ongoing challenge. Despite the now decades old
realization of the superfluid to Mott insulator transition in 1D, 2D and 3D
Greiner  _et al._ (2002); Stöferle  _et al._ (2004); Spielman  _et al._
(2007), there has been little progress in realizing strongly correlated
systems such as fractional quantum Hall states. In both cases, interactions
are enhanced by reducing the contribution of the kinetic energy: inhibiting
tunneling in a deep optical lattice in the case of a Mott insulator, or
quenching the kinetic energy with a magnetic field in the case of fractional
quantum Hall states. The former case simply localizes particles to lattice
sites, producing an uncorrelated insulator. Here we describe a new technique
for creating nearly flat bands, even in the presence of strong tunneling,
using Aharonov-Bohm like quantum interference between sites.

We consider a 2D extension to existing 1D dark state optical lattices studied
theoretically Łącki  _et al._ (2016); Jendrzejewski  _et al._ (2016); Ge and
Zubairy (2020); Kubala  _et al._ (2021); Gvozdiovas  _et al._ (2021), and
realized experimentally Wang  _et al._ (2018); Tsui  _et al._ (2020), enabling
interference phenomena that are not possible in 1D. As indicated in Figs.
1(a,b), our lattice is created from a pair of orthogonal standing waves and a
transverse running wave coupling three internal atomic internal states in a
Lambda configuration. The local dark state of this scheme has zero energy and
no excited state contribution. Atomic motion introduces geometric scalar and
vector potentials Dum and Olshanii (1996); Juzeliūnas  _et al._ (2005);
Goldman  _et al._ (2014), as well as non-adiabatic mixing to the excited
state. The geometric potentials are maximal at the nodes of the optical
standing wave where the atomic dark state changes rapidly. For atoms
adiabatically following the dark state in 1D, this gave rise to a Kronig-
Penney like lattice with barriers far narrower than the optical wavelength
Łącki  _et al._ (2016); Jendrzejewski  _et al._ (2016); Wang  _et al._ (2018).
For specific parameters we find the natural 2D analog of this lattice
consisting of square tiles spaced by narrow barriers. However, generically the
geometric scalar and vector potentials—the latter quantified by the Berry
curvature—can form a lattice of Dirac $\delta$-function like needles, or can
take on a serpentine appearance, creating an array of undulating tubes.
Unexpectedly, we observe that these tubes abruptly decouple at critical points
in parameter space where Aharonov-Bohm interference from the geometric vector
potential inhibits tunneling. This produces nearly completely flat bands
transverse to the tubes with barriers that would otherwise allow substantial
tunneling. The band flattening is analogous to the formation of dispersionless
Landau levels with the application of a uniform magnetic field.

Figure 1(b) schematically illustrates our proposed experimental geometry. A
laser beam traveling along $\boldsymbol{e}_{z}$ drives one arm of a
$\Lambda$-scheme that is then completed by a second arm consisting of four
mutually interfering laser beams in the
$\boldsymbol{e}_{x}$-$\boldsymbol{e}_{y}$ plane. This geometry adds two
additional degrees of freedom as compared to the existing 1D dark-state
lattices: the relative intensity intensity between the in-plane lasers as well
as their relative phase (controlled by displacing a retro-reflection mirror).
Changing the relative intensity converts linear barriers in to serpentine
ones. Tuning the relative phase morphs the lattice from a 2D Kronig-Penney
like lattice of linear barriers to one with needle-like potential maxima. In
addition, the phase difference breaks time reversal symmetry, introducing a
non-zero Berry curvature.

Our manuscript is organized as follows. We introduce the basic formulation of
our 2D $\Lambda$-lattice in Sec. I, and identify the associated symmetries in
Sec. II. Section III describes our numerical method, and presents our main
results. Lastly, in Sec. IV we conclude with a discussion and outlook.

##  I Formulation

###  I.1 Hamiltonian for the 2D Lambda scheme

Figure 1: a), Lambda($\Lambda$) coupling scheme. Two atomic ground states
$\ket{1}$ and $\ket{2}$ are laser-coupled with detuning $\Delta$ to an excited
state $\ket{e}$ (with spontaneous decay rate $\Gamma$) with strengths
$\Omega_{1,2}(\boldsymbol{r})$. b), Schematic. $\Omega_{1}(\boldsymbol{r})$ is
a plane wave traveling along $\boldsymbol{e}_{z}$ and
$\Omega_{2}(\boldsymbol{r})$ consists of two orthogonal standing waves created
by the pictured interfering beams. A movable mirror imparts a controllable
phase shift on the retro-reflected beam traveling along
$\boldsymbol{e}_{x}-\boldsymbol{e}_{y}$. c), Geometric scalar $U$ and d),
Berry curvature ${\mathbf{B}}$ for motion in the dark state computed for
$\epsilon=0.4$, $\epsilon_{c}=1$ and $\chi=\pi/4$.

We consider ultracold atoms subject to a 2D atom-light interaction

$\frac{\hat{V}(\boldsymbol{r})}{\hbar}=\left(-\Delta-\frac{i}{2}\Gamma\right)\ket{e}\bra{e}+\sum_{j=1}^{2}\left[\;\frac{\Omega_{j}(\boldsymbol{r})}{2}\ket{e}\bra{j}+{\rm
H.c.}\right]\,,$ (1)

describing the $\Lambda$-type coupling scheme Dum and Olshanii (1996);
Juzeliūnas  _et al._ (2005); Goldman  _et al._ (2014) shown in Fig. 1(a). Here
the atomic ground states $\ket{1}$ and $\ket{2}$ are coupled with strength
$\Omega_{j}(\boldsymbol{r})$ to the excited state $\ket{e}$ with detuning
$\Delta$ 111Here we adopt the rather unnatural “spectroscopy” convention where
positive $\Delta$ corresponds to negative energy difference for transitions to
the excited state.. The excited state has a spontaneous decay rate $\Gamma$.
Altogether this gives the atomic Hamiltonian

$\hat{H}=\frac{\hat{\mathbf{p}}^{2}}{2M}+\hat{V}(\boldsymbol{r})$ (2)

in terms of the position $\boldsymbol{r}=(x,y)$, momentum
$\hat{\mathbf{p}}=-i\hbar\boldsymbol{\nabla}$, and atomic mass $M$. The above
Hamiltonian is non-Hermitian due to the imaginary contribution $i\Gamma/2$
arising from time-irreversible decay from the excited state $\ket{e}$ in Eq.
(1). In Sec. III.2 we numerically demonstrate that losses due to $i\Gamma/2$
are minimal in the so-called dark state and thus Hermitian dynamics are
maintained. Additionally, even the Hermitian contribution to (1) can break
time reversal symmetry when any of the $\Omega_{j}$ coefficients are complex.

The Hamiltonian acts in the space of state-vectors

$\ket{\psi\left(\boldsymbol{r}\right)}=\sum_{j=1,2,e}\psi_{j}\left(\boldsymbol{r}\right)\ket{j}\,,$
(3)

containing the atomic internal states $\ket{j}$ and the associated wave-
functions $\psi_{j}\left(\boldsymbol{r}\right)$ for atomic center of mass
motion. The corresponding full abstract state vector would be given by
$\ket{\psi}=\int
d\boldsymbol{r}\ket{\psi\left(\boldsymbol{r}\right)}\otimes\ket{\boldsymbol{r}}$,
with
$\ket{\psi\left(\boldsymbol{r}\right)}=\innerproduct{\boldsymbol{r}}{\psi}$.

###  I.2 New basis with dark state

We now re-express $\ket{\psi\left(\boldsymbol{r}\right)}$ in a basis
containing a long-lived dark state in which geometric potentials with sub-
wavelength features can emerge. A dark state is a (generally position-
dependent) superposition of atomic ground states for which
$\hat{V}(\boldsymbol{r})\ket{D}=0$. Therefore, in such a basis,
$\hat{V}(\boldsymbol{r})$ contributes no potential energy, no coupling terms,
and no spontaneous decay for $\ket{D}$. This allows the two coupling arms in
Fig. 1(a) to be driven on resonance with $\ket{e}$ without loss from
$\ket{D}$.

Here we consider orthogonal dark

$\ket{D(\boldsymbol{r})}=\frac{1}{\Omega}\left[\;\Omega_{2}(\boldsymbol{r})\ket{1}-\Omega_{1}(\boldsymbol{r})\ket{2}\;\right]$
(4)

and bright

$\ket{B(\boldsymbol{r})}=\frac{1}{\Omega}\left[\;\Omega_{1}(\boldsymbol{r})\ket{1}+\Omega_{2}(\boldsymbol{r})\ket{2}\;\right]$
(5)

state superpositions, where
$\Omega=\sqrt{\left|\Omega_{1}\right|^{2}+\left|\Omega_{2}\right|^{2}}$ is an
averaged coupling strength. Unlike $\ket{D(\boldsymbol{r})}$, the bright state
couples to the excited state $\ket{e}$.

In the basis of dark, bright and excited states, the state vector is

$\ket{\psi\left(\boldsymbol{r}\right)}=\psi_{\rm
D}\left(\boldsymbol{r}\right)\ket{D(\boldsymbol{r})}+\psi_{\rm
B}\left(\boldsymbol{r}\right)\ket{B(\boldsymbol{r})}+\psi_{\rm
E}\left(\boldsymbol{r}\right)\ket{e}\,,$ (6)

where $\psi_{\rm D}\left(\boldsymbol{r}\right)$, $\psi_{\rm
B}\left(\boldsymbol{r}\right)$ and $\psi_{\rm E}\left(\boldsymbol{r}\right)$
are wave-functions for the atomic center of motion in the corresponding
internal states.

The atom-light coupling operator $\hat{V}(\boldsymbol{r})$ is diagonalized by
the trio of dressed states $\ket{D}$ and $\ket{\pm}$, where $\ket{\pm}$ are
superpositions of $\ket{B}$ and $\ket{e}$ only. When these states depend on
position, they are not eigenstates of the full Hamiltonian $\hat{H}$ due to
the kinetic energy $\hat{\mathbf{p}}^{2}/(2M)$; this leads to geometric
potentials for the projected dynamics in each dressed state Goldman  _et al._
(2014); Mead (1992). In the present case, geometric potentials are introduced
by the spatially varying coupling strengths $\Omega_{1,2}(\boldsymbol{r})$.

###  I.3 Effective potentials for adiabatic dark state

When the total Rabi frequency $\Omega$ at every point in space greatly exceeds
the characteristic energy of the atomic center of mass motion, the atoms will
adiabatically follow their initial dressed state with negligible transitions
to the other dressed states. For dark-state atoms, the state vector (6) can be
approximated as

$|\psi\left(\boldsymbol{r}\right)\rangle\approx\psi_{\rm
D}\left(\boldsymbol{r}\right)\ket{D(\boldsymbol{r})}\,.$ (7)

The validity of this approximation for 1D dark state lattices has been
extensively studied Łącki  _et al._ (2016); Ge and Zubairy (2020);
Jendrzejewski  _et al._ (2016). We correspondingly arrive at the 2D adiabatic
Hamiltonian in the dark state manifold Juzeliūnas  _et al._ (2005); Goldman
_et al._ (2014)

$\hat{H}_{\rm D}=\frac{1}{2M}(-i\hbar\boldsymbol{\nabla}-{\mathbf{A}}_{\rm
D})^{2}+U_{\rm D}\,,$ (8)

where $U_{\rm D}$ and ${\mathbf{A}}_{\rm D}$ are the geometric scalar and
vector potentials. Because our focus is on the dark state manifold, we
suppress the subscript ${\rm D}$ in what follows. The scalar potential

$U_{\rm D}\equiv
U(\boldsymbol{r})=\frac{\hbar^{2}}{2M}\frac{(\boldsymbol{\nabla}\xi^{\ast})\cdot(\boldsymbol{\nabla}\xi)}{\left(1+|\xi|^{2}\right)^{2}}\,$
(9)

is plotted in Fig. 1(c) and Fig. 2. Here we introduce
$\xi(\boldsymbol{r})\equiv\Omega_{2}(\boldsymbol{r})/\Omega_{1}(\boldsymbol{r})$,
the complex valued ratio of coupling strengths in Eq. (4). The geometric
vector potential

${\mathbf{A}}_{\rm
D}\equiv{\mathbf{A}}(\boldsymbol{r})=i\hbar\frac{\xi^{\ast}\boldsymbol{\nabla}\xi-\xi\boldsymbol{\nabla}\xi^{\ast}}{2\left(1+|\xi|^{2}\right)}\,$
(10)

is non-zero only when $\xi(\boldsymbol{r})$ has an imaginary component.
Lastly, the geometric magnetic field shown in Fig. 1(d) is the curl of the
vector potential

${\mathbf{B}}(\boldsymbol{r})=\boldsymbol{\nabla}\times{\mathbf{A}}=i\hbar\frac{(\boldsymbol{\nabla}\xi^{\ast})\times(\boldsymbol{\nabla}\xi)}{\left(1+|\xi|^{2}\right)^{2}}\,.$
(11)

We confirm the adiabatic assumption for our 2D lattice in Sec. III.2 by
comparing numerical results from the full non-Hermitian Hamiltonian (2) to the
adiabatic approximation (Fig. 4).

###  I.4 Sub-wavelength effective potentials

Figure 2: Geometric scalar potential $U$ for various $\chi$ and $\epsilon_{c}$
with $\epsilon=0.1$.

We now describe a configuration of couplings $\Omega_{1,2}(\boldsymbol{r})$
shown in Fig. 1(b) which yield rapid changes in the dark state wavefunction
(4), parameterized by $\boldsymbol{\nabla}\xi$ [entering Eqs. (9)–(11)],
resulting in 2D geometric potentials with features below the optical
diffraction limit. The first laser field
$\Omega_{1}(\boldsymbol{r})=\Omega_{p}\,\exp(ik_{\mathrm{R}}z)$ coupling
$\ket{1}$ to $\ket{e}$ is a plane wave traveling along $\boldsymbol{e}_{z}$
with amplitude $\Omega_{p}$ and wavevector $k_{\mathrm{R}}=2\pi/a$. We
consider atoms that are tightly confined in the $z=0$ plane, such that
$\exp(ik_{\mathrm{R}}z)\approx 1$. The second coupling field
$\Omega_{2}(\boldsymbol{r})$ results from a crossed pair of standing waves
with amplitudes $\Omega^{(\pm)}_{c}$ in the
$\boldsymbol{e}_{x}$-$\boldsymbol{e}_{y}$ plane Hemmerich  _et al._ (1991,
1992); Juzeliūnas and Spielman (2012); a controllable path-length-difference
$d$ for the retro-reflected field along
$(\boldsymbol{e}_{x}-\boldsymbol{e}_{y})/\sqrt{2}$ introduces a phase
$\chi=d/a$ which breaks time-reversible symmetry and allows for non-zero Berry
curvature.

After making the rotating wave approximation (RWA) and phase shifting
$\ket{2}$ by $e^{-i\chi/2}$, the Rabi frequencies of the coupling fields are

$\displaystyle\Omega_{1}(\boldsymbol{r})=\Omega_{p}\,,$ (12)
$\displaystyle\Omega_{2}(\boldsymbol{r})=\sum_{\pm}\pm\Omega^{(\pm)}_{c}e^{\mp
i\chi/2}\cos(k_{\mathrm{R}}x\pm k_{\mathrm{R}}y)\,.$

Inserting (12) into (9)–(11), the explicit forms of the geometric potentials
for this choice of coupling fields are

$\displaystyle\frac{U(\boldsymbol{r})}{E_{\mathrm{R}}}$
$\displaystyle=\frac{\beta^{2}_{+}+\epsilon^{2}_{c}\beta^{2}_{-}}{\alpha^{2}}2\epsilon^{2}(1+\epsilon_{c}^{2})\,,$
$\displaystyle\frac{{\mathbf{A}}(\boldsymbol{r})}{\hbar k_{\mathrm{R}}}$
$\displaystyle=\frac{\sin{(2k_{\mathrm{R}}y)}\mathbf{e}_{x}+\sin{(2k_{\mathrm{R}}x)\mathbf{e}_{y}}}{\alpha}\epsilon_{c}\sin{\chi}\,,$
(13) $\displaystyle\frac{B_{z}(\boldsymbol{r})}{\hbar k_{\mathrm{R}}^{2}}$
$\displaystyle=\frac{\cos{(2k_{\mathrm{R}}x)}-\cos{(2k_{\mathrm{R}}y)}}{\alpha^{2}}2\epsilon^{2}(1+\epsilon_{c}^{2})\epsilon_{c}\sin{\chi}\,,$

where $E_{\mathrm{R}}=\hbar^{2}k_{\mathrm{R}}^{2}/(2M)$ is the single photon
recoil energy, and
${\mathbf{B}}(\boldsymbol{r})=B_{z}(\boldsymbol{r})\,\boldsymbol{e}_{z}$
implying that the magnetic field is orthogonal to the
$\boldsymbol{e}_{x}$-$\boldsymbol{e}_{y}$ plane. Here we defined a factor

$\displaystyle\alpha{(\boldsymbol{r})}$
$\displaystyle=\epsilon^{2}(1+\epsilon_{c}^{2})+\eta^{2}_{+}+\epsilon_{c}^{2}\eta^{2}_{-}-2\epsilon_{c}\eta_{+}\eta_{-}\cos{\chi}\,$

present in (13), and

$\displaystyle\eta_{\pm}(\boldsymbol{r})$
$\displaystyle=\cos{(k_{\mathrm{R}}x\pm k_{\mathrm{R}}y)}\,,$
$\displaystyle\beta_{\pm}(\boldsymbol{r})$
$\displaystyle=\sin{(k_{\mathrm{R}}x\pm k_{\mathrm{R}}y)}\,,$

as well as the ratios of the laser field amplitudes

$\displaystyle\epsilon$
$\displaystyle=\frac{\Omega_{p}}{\sqrt{\Omega_{c}^{(+)2}+\Omega_{c}^{(-)2}}}$
and $\displaystyle\epsilon_{c}$
$\displaystyle=\frac{\Omega^{(-)}_{c}}{\Omega^{(+)}_{c}}.$ (14)

The ratios $\epsilon$ and $\epsilon_{c}$ determine how rapidly the internal
structure of the dark state changes near the zeros of
$\Omega_{2}(\boldsymbol{r})$, allowing control of both the height and spatial
extent of the effective potentials in Eqs. (13).

The scalar potential $U$ is plotted Fig. 2 for different values of $\chi$ and
$\epsilon_{c}$ ($B_{z}$, not shown, is graphically very similar to $U$). The
peak values of the scalar and magnetic fields

$\displaystyle U_{\text{max}}$
$\displaystyle=\frac{2E_{\mathrm{R}}}{\epsilon^{2}}$ and $\displaystyle
B_{\text{max}}(\boldsymbol{r})$ $\displaystyle=\frac{4\hbar
k_{\mathrm{R}}^{2}\epsilon_{c}\sin{\chi}}{\epsilon^{2}(1+\epsilon_{c}^{2})}\,$

are proportional to $2/\epsilon^{2}$, and thus diverge as $\epsilon\rightarrow
0$. Additionally, $B_{\text{max}}$ depends on $\epsilon_{c}$ and $\chi$,
reaching a maximum value with $\epsilon_{c}=1$ and $\chi=\pi/2$. We also
numerically computed the full width at half maximum of the maxima of $U$ and
${\mathbf{B}}$ along their thinnest direction (as seen in Fig. 2(d), this
direction has no particular association with $\boldsymbol{e}_{x}$ or
$\boldsymbol{e}_{y}$). We find that when $\epsilon\ll 1$, i.e.,
$\Omega_{2}\gg\Omega_{1}$, the lattice has tall sub-wavelength barriers that
can be further tuned by adjusting $\epsilon_{c}$ and $\chi$.

The ratio $\epsilon_{c}$ determines the degree of serpentine bending in the
geometric potential. This leads to an effective 1D to 2D transition; for
example (with $\chi=0$), the potential transitions from a brick-like structure
(with holes at the crossing points) at $\epsilon_{c}=1$ [Fig. 2(a)], to an an
array of modulated walls [$\epsilon_{c}=0.5$ in Fig. 2(d)], finally arriving
at straight 1D walls [$\epsilon_{c}=0$]. Additionally, the lattices shown in
Fig. 2(d)–(f) can be rotated by 90 degrees by replacing
$\epsilon_{c}\rightarrow 1/\epsilon_{c}$.

When $\chi=\pi/2$ and $\epsilon_{c}=1$ [Fig. 2(c)], the magnetic field and
scalar potential reduce to a 2D array of needle-like peaks. Moreover, with
$\epsilon\rightarrow 0$ the 2D integral of $U$ over the peak converges to
$E_{\mathrm{R}}a^{2}/2\pi=\hbar^{2}\pi/M$. As such, even for
$\epsilon\rightarrow 0$ the surface integral of the scalar potential does not
diverge, leading to a 2D array of Dirac $\delta$ function potentials—a 2D
Dirac comb—with strength $E_{\mathrm{R}}a^{2}/2\pi$. The same applies to the
magnetic field ${\mathbf{B}}$ with strength $2\pi\hbar$.

Just as in the 1D case, intensity imbalances between the different arms of the
$\Omega_{2}({\bf r})$ field ultimately limit the minimum width of the barriers
Wang  _et al._ (2018). Since our primary focus is on interference effects
rather than minimizing the barrier widths, this is not a significant
consideration in this work.

In the next Section we consider the symmetries of the atom-light coupling
which impose requirements on the eigensolutions of both the full and dark
state Hamiltonians. These symmetries will be later utilized in the numerical
treatment to unfold the energy bands.

##  II Symmetries of the Hamiltonian

Including the couplings $\Omega_{1,2}$ in Eq. (12), the full Hamiltonian (2)
is invariant with respect to spatial shifts along $\boldsymbol{e}_{x}$ and
$\boldsymbol{e}_{y}$ by the lattice constant $a$, i.e.,
$\hat{H}(x+a,y)=\hat{H}(x,y+a)=\hat{H}(x,y)$, so that

$\left[\hat{H},\exp\left(-\frac{i\,\mathbf{a}_{l}\cdot\hat{\mathbf{p}}}{\hbar}\right)\right]=0\,,\quad\mathrm{with}\quad
l=1,2\,,$ (15)

with elementary unit vectors

$\displaystyle\mathbf{a}_{1}$ $\displaystyle=a\,\boldsymbol{e}_{x},$
$\displaystyle{\rm and}$ $\displaystyle\mathbf{a}_{2}$
$\displaystyle=a\,\boldsymbol{e}_{y}\,.$ (16)

Interestingly, the dark state geometric potentials are symmetric with regards
to translations by $a/2$, as evident in Figs. 1(c,d) and Fig. 2. By contrast,
the full Hamiltonian (2) does not obey this symmetry. The couplings
$\Omega_{1}$ and $\Omega_{2}$ are symmetric and anti-symmetric, respectively,
with the $a/2$ spatial shifts

$\displaystyle\Omega_{1}(x+a/2,y)=\Omega_{1}(x,y+a/2)=\Omega_{1}(x,y)\,,$ (17)
$\displaystyle\Omega_{2}(x+a/2,y)=\Omega_{2}(x,y+a/2)=-\Omega_{2}(x,y)\,.$

Thus the Hamiltonian $\hat{H}$ commutes with two combined shift operators

$\hat{T}_{\mathbf{a}_{l}/2}=\hat{U}\exp\left(-\frac{i\,\mathbf{a}_{l}\cdot\hat{\mathbf{p}}}{2\hbar}\right)\,,\quad\mathrm{where}\quad
l=1,2\,,$ (18)

and

$\displaystyle\hat{U}$ $\displaystyle=\left|2\right\rangle\left\langle
2\right|-\left|e\right\rangle\left\langle
e\right|-\left|1\right\rangle\left\langle 1\right|,$ $\displaystyle{\rm with}$
$\displaystyle\hat{U}^{2}=\hat{I}\,.$ (19)

The operator (18) combines a spatial translation by $\mathbf{a}_{l}/2$ with a
$\pi$ phase-flip of the states $\ket{e}$ and $\ket{1}$. Thus the square of the
combined operator
$\hat{T}_{\mathbf{a}_{l}/2}^{2}=e^{-i\,\mathbf{a}_{l}\cdot\hat{\mathbf{p}}/\hbar}$
returns to a state-independent spatial shift by $a$. The Hamiltonian $\hat{H}$
and the combined shift operator $\hat{T}_{\mathbf{a}_{l}/2}$ therefore share a
set of eigenstates following the Bloch ansatz

$\ket{\psi_{s}^{(\boldsymbol{q})}(\boldsymbol{r})}=e^{i\boldsymbol{q}\cdot\boldsymbol{r}}\ket{g_{s}^{(\boldsymbol{q})}(\boldsymbol{r})}\,,$
(20)

with

$\hat{H}\left|\psi_{s}^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle=E_{s}(\boldsymbol{q})\left|\psi_{s}^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle\,,$
(21)

and

$\hat{T}_{\mathbf{a}_{l}/2}\left|\psi_{s}^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle=e^{i\boldsymbol{q}\cdot\mathbf{a}_{l}/2}\left|\psi_{s}^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle\,,$
(22)

with eigenenergy $E_{s}(\boldsymbol{q})$, crystal momentum $\boldsymbol{q}$,
and dark state band index $s=1,2,3,...$. In what follows we focus on the
lowest band with $s=1$ and therefore omit the band index 222In practice, we
diagonalize the full Hamiltonian and identify the dark state manifold as those
bands with the lowest minimum imaginary contribution to the energy..

The periodic part of the Bloch solution (20) satisfies

$\hat{T}_{\mathbf{a}_{l}/2}\left|g^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle=\hat{U}\left|g^{(\boldsymbol{q})}(\boldsymbol{r}+\mathbf{a}_{l}/2)\right\rangle=\left|g^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle\,$
(23)

for a spatial shift of a half of the lattice constant. Expanding in terms of
atomic internal states gives

$|g^{(\boldsymbol{q})}(\boldsymbol{r})\rangle=\sum_{j=e,1,2}g_{j}^{(\boldsymbol{q})}\left(\boldsymbol{r}\right)\ket{j}\,,$
(24)

subject to the conditions

$g_{j}^{(\boldsymbol{q})}(\boldsymbol{r}+\mathbf{a}_{l}/2)=-g_{j}^{(\boldsymbol{q})}(\boldsymbol{r})\quad\mathrm{for}\,\,j=e,1\,,$
(25)

and

$g_{j}^{(\boldsymbol{q})}(\boldsymbol{r}+\mathbf{a}_{l}/2)=g_{j}^{(\boldsymbol{q})}(\boldsymbol{r})\quad\mathrm{for}\,\,j=2\,.$
(26)

The Bloch ansatz given by Eq. (20) is characterized by a 2D crystal momentum
$\boldsymbol{q}=q_{x}\mathbf{e}_{x}+q_{y}\mathbf{e}_{y}$ covering an extended
Brillouin zone (BZ) with $q_{x,y}\in[-k_{\mathrm{R}},k_{\mathrm{R}})$, a four
fold increase in area compared to the BZ of a square lattice with period $a$.
Correspondingly, the area of the unit cell is reduced by a factor of four
333One could choose $\hat{U}=-\hat{U}$ in (19) for the combined translation
operator, as this too would satisfy $\hat{U}^{2}=\hat{I}$. However, the form
given by (19) unfolds the energy bands such that the band corresponding to the
ground state in the dark state manifold typically has minimum energy and
minimum losses at the center of the BZ with $\boldsymbol{q}=\mathbf{0}$, e.g.,
see Fig. 3. Furthermore, there are in total four $\hat{U}$ operators one could
choose from: two each for independent shifts along $\boldsymbol{e}_{x}$ and
$\boldsymbol{e}_{y}$.

We note that when $\chi=\pi/2$, the Hamiltonian supports an additional
symmetry with respect to the spatial shifts by
$\mathbf{a}_{1}/4\pm\mathbf{a}_{2}/4$, as can be seen in Fig. 2(c,f). In this
particular case the BZ can be further unfolded into a rhombus.

For $\chi=0$ (and neglecting the decay rate $\Gamma$), the Hamiltonian
$\hat{H}$ obeys time reversal symmetry. In that case a simultaneous complex
conjugation and inversion of the quasi-momentum leaves the eigenvalue equation
(21) unchanged, giving

$E^{(-\boldsymbol{q})}=E^{(\boldsymbol{q})}\quad\mathrm{and}\quad
g_{j}^{(-\boldsymbol{q})}(\boldsymbol{r})=\left[g_{j}^{(\boldsymbol{q})}(\boldsymbol{r})\right]^{*}\,.$
(27)

We numerically confirmed that this condition is well maintained for atomic
dynamics in the dark state manifold where atomic decay is suppressed. Note
that even for $\chi\neq 0$ the condition
$E^{(-\boldsymbol{q})}=E^{(\boldsymbol{q})}$ holds because complex conjugation
does not change $|\Omega_{2}|$ and thus does not alter the energy spectra
plotted in Figs. 3, 4.

##  III Numerical results

We now describe our numerical method for obtaining the energy spectra and
present our findings 444See Supplementary Material at
http://web.vu.lt/ff/g.juzeliunas/2D_Lambda_Supplementary.zip for the Python3
codes used to obtain the results featured in this publication. The material
contains instructions on installing and operating the software..

###  III.1 Numerical method

The band structure is most easily solved by factoring out the plane wave
component $e^{i\boldsymbol{q}\cdot\boldsymbol{r}}$ in the Bloch eigen-function
(20)

$\hat{H}^{(\boldsymbol{q})}\left|g^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle=E^{(\boldsymbol{q})}\left|g^{(\boldsymbol{q})}(\boldsymbol{r})\right\rangle\,,$
(28)

transforming the Hamiltonian from Eq. (21) to

$\hat{H}^{(\boldsymbol{q})}(\boldsymbol{r})=\frac{1}{2M}\left(-i\hbar\boldsymbol{\nabla}+\hbar\boldsymbol{q}\right)^{2}+\hat{V}(\boldsymbol{r})\,.$
(29)

We obtain the energy spectrum numerically using the Fourier representation of
the eigen-value equation (28) and of the periodic Bloch functions

$\ket{g^{(\boldsymbol{q})}(\boldsymbol{r})}=\sum_{n_{x},n_{y}=-N}^{N}\>e^{ik_{R}(xn_{x}+yn_{y})}\ket{g^{(\boldsymbol{q})}(n_{x},\,n_{y})}\,,$
(30)

with $(2N+1)^{2}$ Fourier components. The resulting right-handed matrix
eigenvalue problem

$E^{(\boldsymbol{q})}\,g^{(\boldsymbol{q})}_{n_{x}n_{y}j}=\left[H^{(\boldsymbol{q})}\right]^{n^{\prime}_{x}n^{\prime}_{y}j^{\prime}}_{n_{x}n_{y}j}\,g^{(\boldsymbol{q})}_{n_{x}n_{y}j}\,$
(31)

is encoded with the combined set of indices $(n_{x},n_{y},j)$ including both
the Fourier components ($n_{x}$, $n_{y}$) and the atomic internal states $j$.
The Hamiltonian-matrix $[H^{(\boldsymbol{q})}]$ is sparsely populated with
$3^{2}(2N+1)^{4}$ elements and a typical filling ratio $\lesssim 10^{-5}$ for
$N\approx 100$. We use shift inversion to amplify solutions near the bottom of
the dark state manifold using libraries optimized for sparse matrix
diagonalization Lehoucq  _et al._ (1998); Virtanen  _et al._ (2020). Section
II implies that certain Fourier components of the periodic Bloch function
$g_{j}^{(\boldsymbol{q})}(n_{x},\,n_{y})$ must be zero to unfold the BZ; we
strictly enforce this condition by zeroing out some of matrix elements as
described in Appendix A.

Diagonalization of the adiabatic dark state Hamiltonian (8) is much more
challenging numerically due to many non-zero Fourier components associated
with the effective potentials $U$ and ${\mathbf{A}}$. Although the adiabatic
dark state Hamiltonian-matrix $[H^{(\boldsymbol{q})}_{D}]$ is 9 times smaller
with $(2N+1)^{4}$ elements, for $N\approx 100$ it has a filling ratio of
$\lesssim 0.04$, making it significantly more dense. By removing Fourier
components with negligible amplitudes, we reduce the filling ratio to
$\lesssim 0.01$.

###  III.2 Energy dispersions

Figure 3: Real (Hermitian) and imaginary (anti-Hermitian) parts of the energy
dispersion (represented by a non-linear colormap and blue contour lines,
respectively) describing the lowest dark Bloch band. The parameters for
(a,b,c,d) match that of Figs. 2(a,d,e,b). Other parameters are $\epsilon=0.1$,
$\Delta=0$, $\hbar\Omega_{p}=2000\,E_{\mathrm{R}}$,
$\hbar\Gamma=1000\,E_{\mathrm{R}}$ and $N=250$. The maximum absolute value of
the anti-Hermitian part of the energy is $\approx 0.022E_{\mathrm{R}}$ in part
c) at the corners of the BZ; the minimum is $\approx 8.6\cdot
10^{-6}E_{\mathrm{R}}$ in part b) for $\boldsymbol{q}=\mathbf{0}$; the average
over parts a)–d) is $\approx 0.01E_{\mathrm{R}}$.

The real and imaginary parts of the energy dispersion describing the lowest
Bloch band in the dark manifold are plotted in Fig. 3 for four combinations of
$\epsilon_{c}$ and $\chi$ (we avoid $\chi=\pi/2$, corresponding to Figs.
2(c,f), which results in a gapless energy dispersion). The real part is
qualitatively different for each combination of parameters in Fig. 3: in (a)
the dispersion is reminiscent of that of a 2D square lattice; in (b) the
curvature near ${\boldsymbol{q}=\mathbf{0}}$ has become anisotropic and tiny
local minima have appeared at the corners of the BZ; in (c) the dispersion at
${\boldsymbol{q}=\mathbf{0}}$ has become a saddle point and the energies at
the corners of the BZ continue to fall; the trend is completed in (d) where
the dispersion has a global maximum at ${\boldsymbol{q}=\mathbf{0}}$.

The imaginary part (blue contours) in Fig. 3 results from a small admixture of
the excited state and is everywhere negative; it contains no contribution from
the gauge field ${\mathbf{A}}$, and thus quantifies only anti-Hermitian
losses. From the perspective of the dark state adiabatic potentials, this
admixture results from non-adiabatic coupling to the bright states. However,
despite the large value of $\hbar\Gamma/E_{\mathrm{R}}=1000$ used in Fig. 3,
we observe a relatively small population transfer into the excited state: for
our parameters the imaginary energy can be as low as $8.6\times
10^{-6}E_{\mathrm{R}}$ at the center of the BZ [Fig. 3(b)]. For a wide range
of $\epsilon_{c}$ and $\chi$, the imaginary part of the energy averaged over
the BZ is $\approx-0.01E_{\mathrm{R}}$. The excited state occupation
probability (and therefore losses) is further reduced at blue-detuning
$\Delta>0$, by reducing sharpness in the potential peaks with a larger value
of $\epsilon$, or by increasing the $\Omega$’s, as observed for 1D dark state
lattices Łącki  _et al._ (2016); Wang  _et al._ (2018); Jendrzejewski  _et
al._ (2016); Ge and Zubairy (2020); Gvozdiovas  _et al._ (2021); Kubala  _et
al._ (2021).

Figure 4: Real part of the energy dispersions from the full Hamiltonian [Eq.
(2)] (a,c), and adiabatic dark state Hamiltonian [Eq. (8)] (b,d). (a,b) were
computed for $\chi=0$, $\epsilon_{c}=1$; (c,d) for $\chi=1.4$ rad and
$\epsilon_{c}=0.09$. The other parameters are $N=100$, $\epsilon=0.1$; for
(a,c), $\hbar\Omega_{p}=2000E_{\mathrm{R}}$, $\hbar\Gamma=1000E_{R}$, and
$\Delta=0$.

Figure 4 compares the real part of the ground dark band dispersion computed
using the full Hamiltonian [Eq. (2)] in (a,c), with that computed using the
adiabatic approximation [Eq. (8)] in (b,d). In each case we computed the band
structure for $\chi=0$, $\epsilon_{c}=1$ and $\chi=1.4\ {\rm rad}$,
$\epsilon_{c}=0.09$ to highlight regimes where the adiabatic approximation is
at its best [in (a,b)] and worst [in (c,d)], respectively. We find that these
dispersions are visually indistinguishable even in the presence of needle-like
barriers [Fig. 2(c)]. Quantitatively, the largest discrepancy is at the edge
of the BZ, where losses (absent in the adiabatic approximation) are maximal.
Conversely, states with crystal momentum $\boldsymbol{q}\approx\mathbf{0}$
retain near-perfect adiabaticity even in the worst case scenario [Fig.
4(c,d)]. The regions of validity of the adiabatic approximation are analogous
to those in 1D systems Wang  _et al._ (2018); Jendrzejewski  _et al._ (2016);
Ge and Zubairy (2020); Gvozdiovas  _et al._ (2021); Kubala  _et al._ (2021):
$\Delta=0$; large $\Omega$’s, non-infinitesimal $\epsilon$, and small, but
non-zero $\Gamma$ to avoid bright-dark resonances Łącki  _et al._ (2016), all
tend to reduce leakage from the dark state. The influence on adiabaticity of
$\chi$ and $\epsilon_{c}$ is non-trivial and the parameters used in (c)
reflect the global maximum in the $\chi$ and $\epsilon_{c}$ parameter space.

Generally, the relationship between losses and the parameters $\epsilon$,
$\chi$ and $\epsilon_{c}$ is complicated: all of these parameters have a
significant influence on the non-adiabatic corrections. Overall, losses
averaged over the BZ are minimal when $\epsilon_{c}=1$ and $\chi=0$; in this
limit the effective magnetic field vanishes, leaving only the scalar potential
which takes the form of an array of 2D cages with gaps as shown in Fig. 2(a)
with the resulting energy dispersion in Fig. 3(a) and Fig. 4(a,b). This
scenario also yields a highly flat ground band with a large band gap ($\approx
3E_{\mathrm{R}}$ with $\epsilon\ll 1$). In this limit, the dark state
Hamiltonian (8) can effectively be approximated as a sum of two orthogonal 1D
potentials with sub-wavelength barriers:

$\hat{H}_{D}(\boldsymbol{r})\approx\hat{H}^{\rm(sep)}_{D}(x)+\hat{H}_{D}^{\rm(sep)}(y)\,,$
(32)

giving energy bands

$E(\boldsymbol{q})=E^{\rm(sep)}(q_{x})+E^{\rm(sep)}(q_{y})\,.$ (33)

The 1D Hamiltonian

$\hat{H}_{D}^{\rm(sep)}(x)=\frac{\hat{p}_{x}^{2}}{2M}+U_{D}^{\rm(sep)}(x)$
(34)

contains a Kronig-Penney like potential

$U^{(\text{sep})}_{D}(x)=\frac{E_{\mathrm{R}}\epsilon_{\text{sep}}^{2}\,\cos^{2}{(k_{\mathrm{R}}x)}}{\left[\epsilon_{\text{sep}}^{2}+\sin^{2}{(k_{\mathrm{R}}x)}\right]^{2}}\,,$
(35)

which appears in 1D analogues of our setup Łącki  _et al._ (2016); Wang  _et
al._ (2018), with $\Omega_{2}(x)=\Omega_{c}\sin{(k_{\mathrm{R}}x)}$ and
$\epsilon_{\text{sep}}=\Omega_{p}/\Omega_{c}$. The real part of the energy
scales approximately as $s^{2}_{x}+s^{2}_{y}$, with positive integers $s_{x}$
and $s_{y}$. For a sufficiently deep lattice with $\epsilon\lesssim 0.2$, we
numerically confirmed that the separable (32) and non-separable (8)
Hamiltonians give similar energy dispersions when
$\epsilon\approx\epsilon_{\text{sep}}$. Such an effective scenario depicts an
array of gapless 2D cages with a constant wall height equal to half of the
maximum height of the true scalar potential (13). Finally, we note that the
full Hamiltonian (2) cannot be treated this way due to its internal structure.

###  III.3 Non-Hermitian tight binding model

Figure 5: a), Hopping directions in a periodic array of tubes. The left(right)
half illustrates odd(even) order tunneling between lattice sites. b),
Hermitian (real) part of hopping amplitudes $J_{m_{x},m_{y}}$ corresponding to
the directions in part a) and describing the ground state in the dark manifold
vs $\epsilon_{c}$ with $\chi=0.2$ rad. Dashed(solid) lines mark
negative(positive) values. c), Real energy dispersion for the thin vertical
line in b) with energy represented by the colormap in Fig. 3. d), Combinations
of $\epsilon_{c}$ and $\chi$ for which the real part of $J_{0,1}$ vanishes.
The orange point coincides with c) and the orange line marks $\chi=0.2$ rad
used in b)–c). Other parameters for b)–d) are $\epsilon=0.1$, $\Delta=0$,
$\hbar\Omega_{p}=2000\,E_{\mathrm{R}}$, $\hbar\Gamma=1000\,E_{\mathrm{R}}$ and
$N=260$.

We express the band structure, such as shown in Fig. 3, as a Fourier transform
Goringe  _et al._ (1997); Delerue  _et al._ (2001); Bellec  _et al._ (2013):

$E^{(\boldsymbol{q})}=\sum_{m_{x}=-\infty}^{\infty}\,\sum_{m_{y}=-\infty}^{\infty}J_{m_{x},\,m_{y}}\exp(-i\boldsymbol{q}\,\cdot\,\delta\mathbf{R}_{m_{x},\,m_{y}})\,,$
(36)

where $J_{\mathbf{m}}$ describing hopping with range
${\mathbf{m}\equiv(m_{x},m_{y})}$ has both real and imaginary parts. The
element $J_{0,0}$ is the on-site energy. These complex tight binding
parameters describe conventional tunneling, and as is well established in
photonic systems, can incorporate both gain and loss Rudner and Levitov
(2009); Schomerus (2013); El-Ganainy  _et al._ (2018); Lieu (2018); Ashida
_et al._ (2020); Han  _et al._ (2021). We obtain the tight binding parameters
from the Fourier transform

$J_{m_{x},\,m_{y}}=k_{\mathrm{R}}^{-2}\int_{\rm
BZ}E^{(\boldsymbol{q})}\exp(i\boldsymbol{q}\,\cdot\,\delta\textbf{R}_{m_{x},\,m_{y}})\,\mathrm{d}\boldsymbol{q}\,$
(37)

of the numerically obtained band structure $E^{(\boldsymbol{q})}$. For any
value of $\chi$ and $\epsilon_{c}$, we find that
$\text{Re}\,J_{m_{x},\,m_{y}}$ fully describes the band structure of the
Hermitian part of the Hamiltonian, and therefore
$\text{Im}\,J_{m_{x},\,m_{y}}$ fully accounts for the anti-Hermitian
contribution of $i\Gamma/2$.

We begin by commenting on the impact of the available physical parameters
$\epsilon$, $\epsilon_{c}$, $\chi$. Starting with the square lattice scenario
depicted in Fig. 2(a), such that $\epsilon_{c}=1$ and $\chi=0$; $\epsilon$
determines the lattice depth, simultaneously modifying all of the hopping
parameters while moving between shallow ($\epsilon\rightarrow 1$) and deep
lattice ($\epsilon\rightarrow 0$) regimes. Next, tuning $\chi\rightarrow\pi/2$
delocalizes atoms as the scalar potential walls shrink to point-like barriers
as shown in Fig. 2(c); this can be counteracted by moving away from
$\epsilon_{c}=1$ towards $\epsilon_{c}\rightarrow 0$ or
$\epsilon_{c}\rightarrow\infty$, restoring the longitudinal extent of barriers
as they approach 1D walls, as can be seen by comparing Figs. 2(b,e).
Alternatively, the spatial extent of barriers can be restored by increasing
$\epsilon$ slightly, consequently reducing their sharpness.

The dependence of the real part of the hopping parameters on $\epsilon_{c}$ is
shown in Fig. 5(b) (with $\epsilon=0.1$ and $\chi=0.2$), and the corresponding
imaginary part is plotted in Fig. 6. Many of the hopping parameters are
identical due to the symmetries discussed in Sec. II—leading to
$J_{m_{x},m_{y}}=J_{m_{y},m_{x}}=J_{-m_{x},-m_{y}}=J_{-m_{y},-m_{x}}$—and are
thus omitted (this includes the non-Hermitian part). Furthermore, when
$\epsilon_{c}=1$ the lattice becomes symmetrical with respect to $x$ and $y$
leading to $J_{\pm m_{x},\pm m_{y}}=J_{\pm m_{y},\pm m_{x}}$.

###  III.4 Band flattening

Figure 5(b) shows our main finding: in a narrow region of $\epsilon_{c}$
($\epsilon_{c}\approx 0.498$, thin vertical line) the real part of the
nearest-neighbor (NN) hopping $J_{1,0}$ vanishes, and the remaining odd
hopping terms such as $J_{2,-1}$, $J_{2,1}$, $J_{3,-2}$ approach zero. The
remaining even-order tunneling processes become dominant, with a leading
contribution from $J_{1,-1}=J_{-1,1}$ describing diagonally oriented tunneling
within tubes (magenta); the next leading contributions are $J_{1,1}$ and
$J_{2,-2}$, giving coupling between next nearest neighboring tubes, and longer
range tunneling within tubes, respectively. This effectively describes an
array of nearly decoupled tubes represented by the highly anisotropic energy
dispersion in Fig. 5(c). The odd terms $J_{2,1}$ and $J_{2,-1}$ couple
neighboring tubes, but are weaker than $J_{1,-1}$ by up to 4 orders of
magnitude. As we describe in Sec. III.5, this results from Aharonov-Bohm like
quantum interference from the geometric vector potential.

Similarly, one can observe the decoupling point by fixing $\epsilon_{c}$ and
tuning $\chi$. In fact, for every $\chi\gtrsim 0.05\pi$, two values of
$\epsilon_{c}$ give the decoupled tube scenario [Fig. 5(d)]; these are related
by $\epsilon_{c}\rightarrow 1/\epsilon_{c}$ and result from identical lattices
rotated by 90 degrees. These two branches merge at $\epsilon_{c}=1$ where the
lattice is symmetric with respect to 90 degree rotations. A more detailed
discussion on these issues is presented in Sec. III.5.

Figure 6: Anti-Hermitian contribution to the tight binding parameters
$J_{m_{x},m_{y}}$ shown in Fig. 5(a,b). Dashed(solid) lines mark
negative(positive) values.

We also examined the energy band gap at the special points. For $\chi=0.2$
rad, $\epsilon=0.1$ and $\epsilon_{c}\approx 0.498$, both the indirect and
direct energy gaps are $\approx 1.25E_{\mathrm{R}}$. They can be made wider
(with an upper limit of $3E_{\mathrm{R}}$) while maintaining the weakly
coupled tube scenario by approaching $\epsilon_{c}\rightarrow 1$, and by
reducing $\chi$ and $\epsilon$.

We now turn to the non-Hermitian part of the energy where our findings are no
less interesting. The imaginary contribution to the on-site energy
$\text{Im}\,J_{0,0}\approx-0.01E_{\mathrm{R}}$ is always negative, describing
on-site atom loss. However, the imaginary part of the energy is nearly zero
for crystal momentum $\boldsymbol{q}=\mathbf{0}$, implying that the sum of the
imaginary tight binding parameters $J_{m_{x},\,m_{y}}$ entering Eq. (36) is
nearly zero. In Fig. 6 we demonstrate that this results from hopping matrix
elements with imaginary components of both signs – an example of tight-binding
gain-loss balance also observed in 1D dark state lattices with decay
Gvozdiovas  _et al._ (2021), which appears despite a strictly lossy
Hamiltonian (2).

###  III.5 Quantum interference

Figure 7:  a), b), Branching paths (blue, red) representing the nearest-
neighbor tunneling $J_{0,1}$ drawn on top of the effective potentials $U$ and
${\mathbf{B}}$ for parameters matching Fig. 5(c). c), Areas of solid and
dashed parallelograms in a)–b). d), Barrier height $U_{p}$ (blue), and width
$d_{p}$ (red), both for the solid contour. c)–d) are plotted for the special
points in Fig. 5(d). The orange points mark parameters describing a)–b).

Here we qualitatively explain the suppression of intertube tunneling by
considering trajectories linking neighboring tubes as sketched in Fig. 7(a,b).
A quantitatively complete path integral description involving the sum over all
paths is not needed to understand the basic origin of the suppression.

To this end, we consider simple ray-like paths connecting the centers of
neighboring lattice sites that undergo Snell’s law type refraction at the
potential barriers (for this argument we do not consider the reduction in
transmission amplitude due to reflections). We compute the phase difference

$\phi_{B}=\frac{1}{\hbar}\left(\int_{C_{1}}{\mathbf{A}}\cdot\mathrm{d}\boldsymbol{r}-\int_{C_{2}}{\mathbf{A}}\cdot\mathrm{d}\boldsymbol{r}\right)=\frac{1}{\hbar}\oint_{C}{\mathbf{A}}\,\mathrm{d}\boldsymbol{r}\,$
(38)

associated with paths $C_{1}$ (red) and $C_{2}$ (blue) that combine to
encircle the tall barrier. As illustrated in Fig. 7(a,b), these together form
a closed contour $C$. The accumulated phase $\phi_{B}$ is thus the line
integral of ${\mathbf{A}}$ along $C$ (equal to the integral of ${\mathbf{B}}$
within $C$ by Stokes’ theorem); when $\phi_{B}=\pi+2\pi n$, for integer $n$,
these two paths destructively interfere, suppressing tunneling.

We investigated two families of contours.

  1. A

Dashed: these cross the scalar potential barrier at its minimum [Fig. 7(a)]
and fully enclose the magnetic field peak [Fig. 7(b)].

  2. B

Solid: these paths are derived from the dashed contour by symmetrically moving
the left and right corners along the scalar potential wall until
$\phi_{B}=\pi$.

Contour (A) was selected to minimize the potential energy cost of the path at
the expense of increased length and kinetic energy owing to the larger corner
angle. For the parameters used in Fig. 5(c), this contour results in
$\phi_{B}\approx 1.3\pi$: larger than needed for destructive interference.
Indeed, Fig. 7(c) plots the area enclosed by these trajectories as a function
of $\epsilon_{c}$, and shows that the optimal trajectory (B) is always reduced
in size. This indicates that the representative (i.e. saddle point) trajectory
minimizes a combination of potential and kinetic energy. Furthermore, at
$\epsilon_{c}=1$ (the point where the lattice has 90 degree rotational
symmetry) the trajectories’ areas are maximized, and as shown in Fig. 7(d) the
barrier height (blue) and width (red) pertaining to contour (B) are minimized.
This corresponds to paths with the most extreme trade-off: minimal potential
energy and maximal kinetic.

Our argument qualitatively describes first order tunneling such as $J_{1,0}$.
More generally, this description also explains suppression of only odd-order
tunneling processes. As an example, consider the even-order hopping parameters
$J_{1,1}$ and $J_{2,0}$. $J_{1,1}$ tunneling is achieved by a single classical
path that cuts through the scalar potential minima (thus without any option
for interference effects). We can explain $J_{2,0}$ in terms of a stacked pair
of solid trajectories, but in this case $\phi_{B}=2\pi$, leading to
constructive interference. In general, even order tunneling terms are
associated with even integer multiples of $\pi$ (either constructive
interference, or none at all) and odd-order trajectories have odd-integer
multiples of $\pi$ (destructive interference). Lastly, $J_{1,-1}$ dominates
because its path is completely unobstructed.

Moreover, we observe destructive interference of NN tunneling in the higher
energy bands of the dark manifold for similar parameter values, supporting the
generic applicability of our classical ray model.

##  IV Discussion and Outlook

The 2D lattice featuring sub-wavelength structures considered here yields
highly tunable geometric scalar and vector potentials with minimal spontaneous
emission. The scalar potential can yield: a 2D square lattice with sub-
wavelength barriers, an array of Delta function-like peaks (a 2D Dirac comb),
or a lattice of interacting zigzag tubes. Furthermore, the band structure is
greatly affected by the geometric vector potential where tunneling between
tubes can be suppressed due to Aharonov-Bohm type destructive interference.

These lattices can be used to realize novel many body phases. When tunneling
is suppressed in conventional deep lattices, the associated maximally
localized Wannier orbitals are very strongly confined to individual lattice
sites. In the present case, both intra- and inter-tube interactions are
enhanced even at near-zero inter-tube tunneling, owing to the relatively
shallow barriers and concomitantly extended Wannier orbitals.

From a broader perspective, this technique can create lattices with features
well below the optical diffraction limit wherever the interfering laser beams
in the $\boldsymbol{e}_{x}$-$\boldsymbol{e}_{y}$ plane approach zero. Changing
the number and intersection angles of these in-plane beams therefore allows
for a range of lattice geometries, including quasi-crystalline. In addition,
the dark state lattice discussed here can also be extended to disordered
configurations by using an optical speckle field for $\Omega_{2}$ in Fig. 1
rather than a standing wave potential. The resulting dark states feature
disordered geometric potentials, including a disordered magnetic field. From
our observation that the synthetic magnetic field can be used to destroy many
of the hopping parameters, one can expect that a disordered magnetic field
could create non-trivial tunneling paths. In the broader context of
localization in 2D disordered systems Abrahams  _et al._ (1979), the
localization properties of such a time-reversal symmetry breaking disorder
potential is unclear Galitski (2005).

###### Acknowledgements.

The authors thank E. Gutierrez, S. Subhankar and E. Benck for carefully
reading the manuscript. This work was supported by the Research Council of
Lithuania (Grant No. S-MIP-20-36). IBS acknowledges support by the National
Institute of Standards and Technology, and the National Science Foundation
through the Quantum Leap Challenge Institute for Robust Quantum Simulation
(Grant No. OMA-2120757). GJ and IBS conceptualized the work; EG carried out
all numerical simulations and analytical derivations, and created all the
figures. All authors contributed equally to writing the manuscript.

##  Appendix A Unfolding the BZ

Here we explain our numerical recipe for obtaining the Hamiltonian-matrix
$[H^{(\boldsymbol{q})}]$ that describes the unfolded BZ. Using the symmetries
discussed in Sec. II, we modify the Hamiltonian-matrix by zeroing out some of
the matrix elements.

The conditions for non-zero matrix elements are determined from Eqs. (19),
(25) and (26). We first define $a/2$ symmetry constants $\mathcal{M}_{j}$ for
each internal state. After choosing $\hat{U}$ according to Eq. (19),
$\mathcal{M}_{j}$ are given by:

$\mathcal{M}_{j}=\begin{cases}0&\text{for $j=2$}\\\ 1&\text{for
$j=e,1$}\end{cases}\,.$ (39)

Thus $\mathcal{M}_{j}$ is even for $j=2$ and odd for $j=e,1$ following Eqs.
(25)–(26). It then follows, for example in the case of $j=2$, that the even
Fourier components $g^{(\boldsymbol{q})}_{2}(2n_{x},2n_{y})$ describing
$g^{(\boldsymbol{q})}_{2}(\boldsymbol{r})$ must be non-zero since it is an
even function with regards to shifts by $a/2$ in $\boldsymbol{e}_{x}$ and
$\boldsymbol{e}_{y}$, giving $\mathcal{M}_{2}=0$. The opposite is true for
$j=e,1$. This can be written mathematically as

$\displaystyle\begin{aligned} &g^{(\boldsymbol{q})}_{j}(n_{x},n_{y})\neq
0\,,\quad\text{if}\\\
&(-1)^{\mathcal{M}_{j}+n_{x}}=1\quad\text{and}\quad(-1)^{\mathcal{M}_{j}+n_{y}}=1\,.\end{aligned}$
(40)

The non-zero Hamiltonian-matrix elements are then

$\displaystyle\begin{aligned}
&\left[H^{(\boldsymbol{q})}\right]^{n^{\prime}_{x}n^{\prime}_{y}j^{\prime}}_{n_{x}n_{y}j}\neq
0\,,\quad\text{if}\\\
&(-1)^{\mathcal{M}_{j}+n_{x}}=1\quad\text{and}\quad(-1)^{\mathcal{M}_{j}+n_{y}}=1\quad\text{and}\\\
&(-1)^{\mathcal{M}_{j^{\prime}}+n^{\prime}_{x}}=1\quad\text{and}\quad(-1)^{\mathcal{M}_{j^{\prime}}+n^{\prime}_{y}}=1\,.\end{aligned}$
(41)

We thus arrive at a new Hamiltonian-matrix with eigensolutions characterized
by the extended BZ $q_{x,y}\in[-k_{\mathrm{R}},k_{\mathrm{R}})$. The Fourier-
space eigenvectors (30) diagonalizing this matrix are not truncated: all of
the $3(2N+1)^{2}$ Fourier components (including the ones equal to zero)
describe the solution, and so the numerically obtained eigenvectors for each
$\boldsymbol{q}$-value describe a real space unit cell of area $a^{2}$ (as
though the BZ were not unfolded). We note that this method of zeroing out
matrix elements is sub-optimal – a better solution would be to truncate the
Fourier space, reducing the size of the Hamiltonian-matrix.

The same filtering procedure is valid for the adiabatic Hamiltonian-matrix
$[H^{(\boldsymbol{q})}_{D}]$. It is described by one internal state $j=j_{D}$
– the dark state (4) – which is invariant with respect to the combined shift
operator $\hat{T}_{\mathbf{a}_{l}/2}$ defined in (18):

$\hat{T}_{\mathbf{a}_{l}/2}\left|D(\boldsymbol{r})\right\rangle=\left|D(\boldsymbol{r})\right\rangle\,,$
(42)

therefore one has $\mathcal{M}_{D}=0$ ($g_{D}^{(\boldsymbol{q})}$ is an even
function, see also Fig. 1(c,d)).

## References

  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. Hänsch, and I. Bloch, Nature 415, 39 (2002).
  * Stöferle  _et al._ (2004) T. Stöferle, H. Moritz, C. Schori, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004).
  * Spielman  _et al._ (2007) I. B. Spielman, W. D. Phillips, and J. V. Porto, Phys. Rev. Lett. 98, 080404 (2007).
  * Łącki  _et al._ (2016) M. Łącki, M. A. Baranov, H. Pichler, and P. Zoller, Phys. Rev. Lett. 117, 233001 (2016).
  * Jendrzejewski  _et al._ (2016) F. Jendrzejewski, S. Eckel, T. G. Tiecke, G. Juzeliūnas, G. K. Campbell, L. Jiang, and A. V. Gorshkov, Phys. Rev. A 94, 063422 (2016).
  * Ge and Zubairy (2020) W. Ge and M. S. Zubairy, Phys. Rev. A 101, 023403 (2020).
  * Kubala  _et al._ (2021) P. Kubala, J. Zakrzewski, and M. Łącki, Phys. Rev. A 104, 053312 (2021).
  * Gvozdiovas  _et al._ (2021) E. Gvozdiovas, P. Račkauskas, and G. Juzeliūnas, SciPost Phys. 11, 100 (2021).
  * Wang  _et al._ (2018) Y. Wang, S. Subhankar, P. Bienias, M. Łącki, T.-C. Tsui, M. A. Baranov, A. V. Gorshkov, P. Zoller, J. V. Porto, and S. L. Rolston, Phys. Rev. Lett. 120, 083601 (2018).
  * Tsui  _et al._ (2020) T.-C. Tsui, Y. Wang, S. Subhankar, J. V. Porto, and S. L. Rolston, Phys. Rev. A 101, 041603(R) (2020).
  * Dum and Olshanii (1996) R. Dum and M. Olshanii, Phys. Rev. Lett. 76, 1788 (1996).
  * Juzeliūnas  _et al._ (2005) G. Juzeliūnas, J. Ruseckas, and P. Öhberg, J. Phys. B: At. Mol. Opt. Phys. 38, 4171 (2005).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Rep. Prog. Phys. 77, 126401 (2014).
  * Note (1) Here we adopt the rather unnatural “spectroscopy” convention where positive $\Delta$ corresponds to negative energy difference for transitions to the excited state.
  * Mead (1992) C. A. Mead, Rev. Mod. Phys. 64, 51 (1992).
  * Hemmerich  _et al._ (1991) A. Hemmerich, D. Schropp, and T. W. Hänsch, Phys. Rev. A 44, 1910 (1991).
  * Hemmerich  _et al._ (1992) A. Hemmerich, D. Schropp, T. Esslinger, and T. W. Hänsch, Europhys. Lett. 18, 391 (1992).
  * Juzeliūnas and Spielman (2012) G. Juzeliūnas and I. B. Spielman, New. J. Phys. 14, 123022 (2012).
  * Note (2) In practice, we diagonalize the full Hamiltonian and identify the dark state manifold as those bands with the lowest minimum imaginary contribution to the energy.
  * Note (3) One could choose $\hat{U}=-\hat{U}$ in (19) for the combined translation operator, as this too would satisfy $\hat{U}^{2}=\hat{I}$. However, the form given by (19) unfolds the energy bands such that the band corresponding to the ground state in the dark state manifold typically has minimum energy and minimum losses at the center of the BZ with $\boldsymbol{q}=\mathbf{0}$, e.g., see Fig. 3. Furthermore, there are in total four $\hat{U}$ operators one could choose from: two each for independent shifts along $\boldsymbol{e}_{x}$ and $\boldsymbol{e}_{y}$.
  * Note (4) See Supplementary Material at http://web.vu.lt/ff/g.juzeliunas/2D_Lambda_Supplementary.zip for the Python3 codes used to obtain the results featured in this publication. The material contains instructions on installing and operating the software.
  * Lehoucq  _et al._ (1998) R. B. Lehoucq, D. C. Sorensen, and C. Yang,  _ARPACK Users’ Guide_ (Society for Industrial and Applied Mathematics, 1998) p. 152.
  * Virtanen  _et al._ (2020) P. Virtanen _et al._ , Nat. Methods 17, 261 (2020).
  * Goringe  _et al._ (1997) C. M. Goringe, D. R. Bowler, and E. Hernández, Rep. Prog. Phys. 60, 1447 (1997).
  * Delerue  _et al._ (2001) C. Delerue, M. Lannoo, and G. Allan, Phys. Stat. Sol. B 227, 115 (2001).
  * Bellec  _et al._ (2013) M. Bellec, U. Kuhl, G. Montambaux, and F. Mortessagne, Phys. Rev. B 88, 115437 (2013).
  * Rudner and Levitov (2009) M. S. Rudner and L. S. Levitov, Phys. Rev. Lett. 102, 065703 (2009).
  * Schomerus (2013) H. Schomerus, Opt. Lett. 38, 1912 (2013).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nat. Phys. 14, 11 (2018).
  * Lieu (2018) S. Lieu, Phys. Rev. B 97, 045106 (2018).
  * Ashida  _et al._ (2020) Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020).
  * Han  _et al._ (2021) Y. Z. Han, J. S. Liu, and C. S. Liu, New J. Phys. 23, 123029 (2021).
  * Abrahams  _et al._ (1979) E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
  * Galitski (2005) V. M. Galitski, Phys. Rev. B 72, 214201 (2005).


# FSER: Deep Convolutional Neural Networks for Speech Emotion Recognition

Bonaventure F. P. Dossou
Jacobs University, Germany
<EMAIL_ADDRESS>Yeno K. S. Gbenou
Drexel University, USA
<EMAIL_ADDRESS>
###### Abstract

Using mel-spectrograms over conventional MFCCs features, we assess the
abilities of convolutional neural networks to accurately recognize and
classify emotions from speech data. We introduce FSER, a speech emotion
recognition model trained on four valid speech databases, achieving a high-
classification accuracy of 95,05%, over 8 different emotion classes: anger,
anxiety, calm, disgust, happiness, neutral, sadness, surprise. On each
benchmark dataset, FSER outperforms the best models introduced so far,
achieving a state-of-the-art performance. We show that FSER stays reliable,
independently of the language, sex identity, and any other external factor.
Additionally, we describe how FSER could potentially be used to improve mental
and emotional health care and how our analysis and findings serve as
guidelines and benchmarks for further works in the same direction.

##  1 Introduction

Emotions are integral parts of daily communication between humans. Whether
oral or written, properly understanding the interlocutor is key for a good
communication. Even though emotions can be easily concealed in written
communication, several psychological studies have shown that it is more
difficult for humans to hide their feelings physically or vocally. More
generally, facial expressions and voice tones are very good indicators of
one’s emotional state. Many psychological and physiological studies [2, 6, 8,
10, 12, 16] (not limited to the ones cited here) have proved, that emotions
make us feel and act; stimulating and influencing both our facial expressions
and voice tone. For instance, adrenaline is released in fearful situations to
help us run away from danger, as excitement or joy can be expressed while we
are talking with friends, family, cuddling our pets or practicing risky sports
or activities such as mountaineering or skydiving.

The recognition of emotional states in speech, so called Speech Emotion
Recognition (SER) is not a new concept in the fields of Artificial
Intelligence and Machine Learning, but still is a very challenging task, that
deserves more attention.

Convolutional neural networks (CNNs) are types of artificial neural networks,
one of the most popular that have helped make breakthroughs in face
recognition, segmentation and object recognition, handwriting recognition, and
of course speech recognition. CNNs get their name from a well-known
mathematical operation called convolution, and in most cases are used for
efficient patterns recognition in data, notably including accurate images
classification.

#### Structure of the paper:

In sections 2 and 3, we respectively describe the speech corpora used in this
study, and the preprocessing pipeline. In section 4, after introducing our
model architecture, we describe our experiments, results and provide a use
case of FSER. The last section is dedicated to the conclusion.

##  2 Datasets, and Emotion Classes Creation

For the current study, we collected, united (i.e. speech samples
representatives of the same emotion were combined) and used speeches from 4
existing datasets: Emodb111http://emodb.bilderbar.info/docu/,
Emovo222http://voice.fub.it/activities/corpora/emovo/index.html,
Savee333http://personal.ee.surrey.ac.uk/Personal/P.Jackson/SAVEE/, RAVDESS
[9]. The speech files across all datasets are sampled at 48 kHz.

Emodb is a publicly available German speech database containing speech
recordings of seven emotions: sadness, anger, happiness, fear, disgust,
boredom and neutrality. The recordings were made by five men and five women,
all actors, creating 10 statements for each emotion that was tested and
evaluated. The entire dataset contains 535 speeches.

Emovo is a publicly available Italian speech database with seven emotions:
happiness, sadness, anger, fear, disgust, surprise and neutral. Six Italian
actors (three women and three men) generated 14 samples for each emotion. The
entire dataset contains 588 speeches.

Savee is a British English language database of public speaking, made of
speech recordings with seven emotions: happiness, sadness, anger, fear,
disgust, surprise and neutral. The recordings were made by four English male
actors, generating 15 samples for each emotion. The entire dataset contains
480 speeches.

Ravdess has been created by 24 authors: 12 men and 12 women with eight
emotions: neutral, calm, happy, sad, angry, scared, disgusted, surprised. The
entire dataset contains 1440 speeches.

Despite the fact that speeches across the datasets are sampled at the same
frequency, they are structured in different and specific ways. Hence, we had
to classify them appropriately and combine the same emotions from all
datasets. The average length of speech samples is 3-4 seconds.

After the collection and classification of all emotions from the four
different datasets, we got 10 classes of emotions: anger, anxiety, fear,
boredom, calm, disgust, happiness, neutral, sadness, surprise. Due to the
small number of $fear$ and $boredom$ emotions samples, we decided to add them
respectively to $anxiety$ and $calm$ emotion classes. This was also motivated
by the fact that, those speech samples have similar pitches and amplitudes.
This results in 8 classes of emotion.

We chose to disregard gender information and to focus only on emotions. The
final distribution of emotion classes with the numbers of speech samples of
each class, are stated in Table 1.

Emotion class | Amount of data (%)
---|---
Anger | 463 (15.22%)
Anxiety | 405 (13.31%)
Calm | 273 (8.97%)
Disgust | 382 (12.55%)
Happiness | 407 (13.37%)
Neutral | 379 (12.45%)
Sadness | 398 (13.09%)
Surprise | 336 (11.04%)
Total | 3043 (100%)
Table 1:  Distribution of emotion classes with the numbers of speech samples
of each class.

##  3 Data Processing and Visualisation

Figure 1: Speech sample visualisation using the $Librosa$ library

Figure 1 shows the amplitudes of a random speech sample along the time
measured in seconds (s). Speech signals are made up of amplitudes and
frequencies. To get more information from our speech samples, we decided to
map them into the frequency domain using the Fast Fourier Transformation (FFT)
[3, 4].

Figure 2: Speech sample characteristics using $Librosa$, and scipy.fft
libraries. Figure 3: Speech sample converted into mel-spectrogram, with time,
frequency and amplitudes information.

(a)

(b)

Figure 4:  Mel-Spectrograms of speakers speaking a sentence with anger (a) and
anxiety (b).

Applying the FFT to a sample using the scipy.fft and $librosa$ packages, gives
Figure 2, which shows the initial speech plot, the double-side (negative and
positive) FFT spectrum, and the positive FFT spectrum.

However, this is still lacking of time information. To remedy this, and make
sure we preserve frequencies, time and amplitudes information about the speech
samples, in reasonable and adequate range, we decided then to use mel-
spectrograms.

In a mel-spectrogram444http://www.glottopedia.org/index.php/Spectrogram the
abscissa represents time, the ordinate axis represents frequency, and
amplitudes are showed by the darkness of a precise frequency at a particular
time: low amplitudes are represented with a light-blue color, and very high
amplitudes are represented by dark red (see Figure 3).

There are two types of spectrograms: broad-band spectrograms and narrow-band
spectrograms. In our study, we used narrow-band spectrograms because they have
higher frequency resolution, and larger time interval for every spectrum than
broad-band spectrograms: this allows the detection of very small differences
in frequencies. Moreover, they show individual harmonic structures, which are
vibration frequency folds of the speech, as horizontal striations.

In Figure 4 are provided examples of mel-spectrograms of speakers speaking a
sentence with anger 4(a) and anxiety 4(b). We converted all original speech
signals to mel-spectrograms using the librosa
melspectrogram555https://librosa.org/doc/main/generated/librosa.feature.melspectrogram.html
module of the $librosa$ library, and transformed the SER task into a pattern
recognition and image classification problem.

We used 512 as length of the FFT window, 512 as the hop-length (number of
samples between successive frames) and a hanning windows size is set to the
length of FFT window.

#### Dataset Splitting:

2434 images (80%) have been allocated for the training phase. The remaining
609 (20%) have been used as testing set, to evaluate the model performance. As
a standard in the machine learning field, we then split the initial training
into two subsets: training set (80% - 1947 images) and validation set (20% -
487 images). To increase the size of the final training set, we applied a very
common technique that has shown much success in recent years with CNNs: data
augmentation [13, 18].

#### Data Augmentation:

Using the algorithm 1, we implemented the data augmentation (DA) method,
generating 20 additional pictures from every original input picture of our
training set. To each input image (mel-spectrogram) object, we applied width
and height shifting, zooming and horizontal flipping, using the
keras.preprocessing.image module from the $Keras$ deep learning
framework666https://keras.io/api/preprocessing/image/. After the DA, our
training set size increased to 40887 images (mel-spectrograms).

Figure 5: FSER Model Architecture with: 1) CL: Convolutional Layer + Relu
Activation function, 2) MPL: Maximum 2D Pooling Layer, 3) Dr: Dropout, and 4)
FCL: Fully-Connected (Dense) Layer. Figure 6: FSER ROC-AUC on 8-emotions
Recognition and Classification

| Precision | Recall | F1-Score | Support | AUC-Score
---|---|---|---|---|---
Anger | 96 | 96 | 96 | 91 | 99.84
Anxiety | 94 | 97 | 95 | 80 | 99.85
Calm | 93 | 97 | 95 | 54 | 99.92
Disgust | 97 | 92 | 94 | 77 | 99.72
Happiness | 95 | 93 | 94 | 83 | 99.75
Neutral | 94 | 94 | 94 | 77 | 99.74
Sadness | 95 | 94 | 94 | 78 | 99.76
Surprise | 94 | 95 | 94 | 69 | 99.82
accuracy | - |  | 95 | 609 | -
macro avg | 95 | 95 | 95 | 609 | 99.79
weighted avg | 95 | 95 | 95 | 609 | 99.79

Table 2: Classification Report of the Predictions of FSER on the Testing set

| Predicted Emotions
---|---
Expected Emotions $\downarrow$ | Anger | Anxiety | Calm | Disgust | Happiness | Neutral | Sadness | Surprise
Anger | 96.92 | 0.31 | 0.16 | 0.42 | 1.15 | 0.42 | 0.37 | 0.26
Anxiety | 0.53 | 95.46 | 0.47 | 0.29 | 1.12 | 0.53 | 0.94 | 0.65
Calm | 0.00 | 0.18 | 96.59 | 0.00 | 0.18 | 1.97 | 0.81 | 0.27
Disgust | 0.88 | 1.19 | 1.00 | 91.49 | 0.56 | 2.06 | 1.44 | 1.38
Happiness | 1.72 | 0.98 | 0.40 | 0.40 | 93.16 | 1.03 | 0.86 | 1.44
Neutral | 0.37 | 0.50 | 1.43 | 0.50 | 0.50 | 95.45 | 1.06 | 0.19
Sadness | 0.41 | 0.94 | 1.00 | 0.35 | 0.23 | 0.59 | 95.84 | 0.64
Surprise | 0.43 | 1.58 | 0.14 | 0.50 | 0.65 | 0.14 | 1.08 | 95.46

Table 3: Normalized Confusion Matrix of FSER on the SER Task

##  4 Model Architecture, Training and Results

###  4.1 Model Architecture and Training

Figure 5 presents the architecture of FSER. The input is an RGB image of shape
(64, 64, 3). Our model is made of four blocks of convolutional layers (CLs)
with reLU activation [11], maximum 2D pooling layers (MPLs) [19] and Dropouts
[5] interposed between them.

| Speech Dataset | Features | Accuracies (%)
---|---|---|---
Relevant related works $\downarrow$ | Emodb | Emovo | Savee | Ravdess | Mel-Spectrogram (*) | MFCC (**) | * | **
Qayyum _et al_. [1] | - | - | + | - | - | + | - | 83.61
Kannan, V. and Haresh, R. [17] | - | - | - | + | + | + | 86 | 53
Zhao _et al_. [20] | + | - | - | - | + | - | 95.33 | -
Our work (Emodb) | ✓ | - | - | - | ✓ | - | 97.74 | ✗
Our work (Savee) | - | - | ✓ | - | ✓ | - | 99.10 | ✗
Our work (Ravdess) | - | - | - | ✓ | ✓ | - | 98.67 | ✗

Table 4: Comparison of our work to other relevant related works on Emodb,
Savee and Ravdess datasets.

The output of the four CLs is flattened and fed into a series of three fully-
connected layers, which is finally fed into a softmax activation unit [7] to
output the probability of the input to belong to each class of emotion.

#### Hyper-parameters and Training:

We used a batch size of 64, a learning rate of 0.001 and Stochastic Gradient
Descent [14, 15] as optimizer. The loss function is the categorical cross-
entropy.

Our hyper-parameters (stated above), including the number of layers, the
number of filters in each layer, the dropout probability, the size of the
kernel, padding and stride, have been selected after many manual fine-tuning
trials.

FSER has been trained for 400 epochs, during 8 days on a 16GB Tesla K80 GPU,
using the platform PaperSpace.

###  4.2 Results and potential use case

Tables 2 and 3 present respectively the classification report of FSER on the
SER task and the normalized confusion matrix. As shown in Tables 2 and 3, and
in Figure 6, FSER has a nearly-perfect score for every emotion class.

Additionally, Table 4 shows the comparison of our work to other works done on
the SER task, with the datasets of interest. We can interpret Table 4 as
follow:

  * •

On Emodb, our FSER (97.44) outperformed the mel-spectrograms-based approach
proposed by Zhao _et al_. [20] (95.53).

  * •

On Savee, our FSER (99.10) outperformed the traditional MFCCs-based approach
proposed by Qayyum _et al_. [1] (83.61)

  * •

On Ravdess, our FSER (98.67) outperformed both mel-spectrograms and
traditional MFCCs approaches proposed by Kannan, V. and Haresh, R. [17] (86
and 53 respectively).

To test the robustness of FSER, we collected real-life audio samples from
French, Fon, and Igbo speakers. Each of them was asked to imitate as much as
they can those emotions, without providing any guidelines. In total, we had 24
audios (3 samples per emotion class), out of which FSER predicted correctly 20
audios. The 4 misclassed audios were found to belong to the classes $fear$ and
$boredom$ that we had to respectively add and consider as $anxiety$ and $calm$
classes. Recalling that FSER has been trained in English, German, and Italian,
its performance on these real-time audios, shows that FSER stays reliable,
independently of the language, sex identity, and any other external factor.

#### Potential use case:

For humans, emotions constitute an important factor to personal and global
health. Being emotionally unstable can affect not only mental wellness but
also physical health. During this pandemic, many health reports demonstrated a
significant increase in feelings of loneliness and emotional instability
leading to and increased rate of suicide for example. Being emotionally
healthy is hence crucial for our well-being. As AI-based systems are nowadays
helping doctors to quickly diagnose diseases such as cancers, FSER-like
systems could also be useful to doctors, particularly psychologists, in
identifying the emotional state of their patient. It is no news that some
people have trouble expressing their emotions with words, so FSER brings a
potential solution to help physicians to better understand their patients, and
provide them more efficient treatments. They could also be useful to
organizations providing mental care and support to families, clients, or
communities at large.

##  5 Conclusion

In the current study, we evaluated the abilities of CNNs to effectively
recognize and classify emotions from speech data. We introduced FSER, which
outperformed all models that have been introduced for SER task. We showed that
FSER is reliable, with no regards to the language spoken, the sex identity of
the individual, or any other external factor. FSER-like systems could be
beneficial to the healthcare system, as they can help to better understand
patients, to provide them better and quicker treatments. It is worth
mentioning the limitations of this work, notably the data augmentation we had
to perform to deal with the limited number of samples.

## References

  * [1] A. B. Abdul Qayyum, A. Arefeen, and C. Shahnaz.  Convolutional neural network (cnn) based speech-emotion recognition.  In 2019 IEEE International Conference on Signal Processing, Information, Communication Systems (SPICSCON), pages 122–125, 2019.
  * [2] Daniel C. Albas, Ken W. McCluskey, and Cheryl A. Albas.  Perception of the emotional content of speech: A comparison of two canadian groups.  Journal of Cross-Cultural Psychology, 7(4):481–490, 1976.
  * [3] B. Boashash.  Time-Frequency Signal Analysis and Processing: A Comprehensive Reference.  Oxford: Elsevier Science, 2003.
  * [4] R. N Bracewell.  The Fourier Transform and Its Applications.  Boston: McGraw-Hill, 2000.
  * [5] François Chollet.  Deep Learning with Python.  Manning, Nov. 2017.
  * [6] Pablo Fossa, Raymond Madrigal Pérez, and Camila Muñoz Marcotti.  The relationship between the inner speech and emotions: Revisiting the study of passions in psychology.  Human Arenas, 2020.
  * [7] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [8] Tom Johnstone.  The effect of emotion on voice production and speech acoustics, Aug 2017\.
  * [9] Steven R. Livingstone and Frank A. Russo.  The Ryerson Audio-Visual Database of Emotional Speech and Song (RAVDESS), Apr. 2018.  Funding Information Natural Sciences and Engineering Research Council of Canada: 2012-341583 Hear the world research chair in music and emotional speech from Phonak.
  * [10] Sydney Lolli, Ari Lewenstein, Julian Basurto, Sean Winnik, and Psyche Loui.  Sound frequency affects speech emotion perception: results from congenital amusia.  Frontiers in Psychology, 6:1340, 2015.
  * [11] Vinod Nair and Geoffrey E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In Johannes Fürnkranz and Thorsten Joachims, editors, ICML, pages 807–814. Omnipress, 2010.
  * [12] Kislova O.O. and Rusalova M.N.  Perception of emotions in speech. a review of psychological and physiological research.  Usp Fiziol Nauk, 2013.
  * [13] Luis Perez and Jason Wang.  The effectiveness of data augmentation in image classification using deep learning, 2017.
  * [14] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  Ann. Math. Statist., 22(3):400–407, 09 1951.
  * [15] Sebastian Ruder.  An overview of gradient descent optimization algorithms, 2016.
  * [16] K.R. Scherer.  Personality, emotion, psychopathology and speech.  In HOWARD GILES, W PETER ROBINSON, and PHILIP M SMITH, editors, Language, pages 233–235. Pergamon, Amsterdam, 1980.
  * [17] Kannan Venkataramanan and Haresh Rengaraj Rajamohan.  Emotion recognition from speech, 2019.
  * [18] Sebastien C. Wong, Adam Gatt, Victor Stamatescu, and Mark D. McDonnell.  Understanding data augmentation for classification: when to warp?, 2016\.
  * [19] Haibing Wu and Xiaodong Gu.  Max-pooling dropout for regularization of convolutional neural networks, 2015.
  * [20] Jianfeng Zhao, Xia Mao, and Lijiang Chen.  Speech emotion recognition using deep 1d and 2d cnn lstm networks.  Biomedical Signal Processing and Control, 47:312–323, 2019.

⬇

1from keras.preprocessing.image import *

2import numpy as np

3

4for image in images:

5 # images is the list of all images

6 image = load_img(image)

7 image = img_to_array(image)

8 image = np.expand_dims(image, 0)

9

10 # construct the image data

11 # generator for data augmentation

12 aug = ImageDataGenerator(*args)

13

14 # generating images

15 imageGen = aug.flow(image, *args)

16 total = 0

17 # loop over image data generator

18 for image_generated in imageGen:

19 total += 1

20 if total == 20:

21 break

Algorithm 1: Algorithm for the generation of new samples for Data Augmentation