# Wannier solitons in spin-orbit-coupled Bose-Einstein condensates in optical
lattices with a flat-band.

Chenhui Wang<EMAIL_ADDRESS>International Center of Quantum Artificial
Intelligence for Science and Technology (QuArtist) and Department of Physics,
Shanghai University, Shanghai 200444, China  Centro de Física Teórica e
Computacional, Faculdade de Ciências, Universidade de Lisboa, Campo Grande,
Ed. C8, Lisboa 1749-016, Portugal  Yongping Zhang<EMAIL_ADDRESS>International Center of Quantum Artificial Intelligence for Science and
Technology (QuArtist) and Department of Physics, Shanghai University, Shanghai
200444, China  V. V. Konotop<EMAIL_ADDRESS>Centro de Física
Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa, Campo
Grande, Ed. C8, Lisboa 1749-016, Portugal  Departamento de Física, Faculdade
de Ciências, Universidade de Lisboa, Campo Grande, Ed. C8, Lisboa 1749-016,
Portugal

###### Abstract

We investigate families of soliton solutions in a spin-orbit coupled Bose-
Einstein condensate embedded in an optical lattice, which bifurcate from the
nearly flat lowest band. Unlike the conventional gap solitons the obtained
solutions have the shape well approximated by a Wannier function (or a few
Wannier functions) of the underlying linear Hamiltonian with amplitudes
varying along the family and with nearly constant widths. The Wannier solitons
(WSs) sharing all symmetries of the system Hamiltonian are found to be stable.
Such solutions allow for the construction of Wannier breathers, that can be
viewed as nonlinearly coupled one-hump solitons. The breathers are well
described by a few-mode model and manifest stable behavior either in an
oscillatory regime with balanced average populations or in a self-trapping
regime characterized by unbalanced atomic populations of the local potential
minima (similarly to the conventional boson Josephson junction), with the
frequencies controlled by the inter-atomic interactions.

##  I Introduction

Periodic modulation of parameters of a medium, where a wave propagates,
introduces artificial dispersion. If the medium is nonlinear the existence of
gap solitons becomes possible. This is a well-known phenomenon thoroughly
studied in diverse physical settings including Bose-Einstein condensates
(BECs) [1, 2, 3] and nonlinear optics [4]. Such soliton-bearing systems are
described by one or several coupled Gross-Pitaevskii (GP) or nonlinear
Schrödinger (NLS) equations with periodic potentials. For BECs loaded in
optical lattices (OLs), the attractive nonlinearity originated by a negative
scattering length of the inter-atomic interactions enables solitons
characterized by chemical potentials, $\mu$, located in either a semi-infinite
or in a finite gap of the spectrum of the underlying linear Hamiltonian.
Otherwise, if the nonlinearity is repulsive (i.e., the scattering length is
positive), solitons may exist with chemical potentials belonging only to
finite gaps. In either of the cases, gap solitons are localized solutions
belonging to families which are usually characterized by the dependence of the
number of atoms $N$ (alias norm) on the chemical potential. In a common
situation, families of solitons bifurcate from one of the edges of a linear
band $\mu_{\nu k}$ ($\nu$ is the index of a band and $k$ is the Bloch
wavenumber). In the small-amplitude limit, $N\to 0$, such gap solitons are
governed by an NLS equation for their envelopes, whose effective dispersion is
determined by $\partial_{k}^{2}\mu_{\nu k}$ (see e.g. [5]). Respectively, when
a soliton amplitude tends to zero its width infinitely increases approaching
the respective Bloch state it bifurcates from.

The described situation, however, is not applicable anymore if
$\partial_{k}^{2}\mu_{\nu k}=0$. Such a situation is encountered in optical
applications, where in the presence of non-zero dispersion,
$\partial_{k}^{4}\mu_{\nu k}\neq 0$, and the nonlinearity support so-called
quartic solitons [6], which can also exist in SO-BECs without lattices (as it
is discussed below).

It turns out that physical systems can have strongly or even completely
suppressed linear dispersion, when several or all derivatives
$\partial_{k}^{n}\mu_{\nu k}=0$ (for $n=1,2,...$) vanish. The last situation
occurs, for instance, in the case of flat bands existing in the spectra of
certain discrete systems. In that case, the spatial localization of waves is
possible even in the purely linear limit, while in the nonlinear case, it
allows for the existence of discrete compactons [7], unconventional families
of nonlinear modes [8], formation of solitons with abrupt edges [9], and
unusual expansion regimes [10].

Being ubiquitous in discrete systems [11, 12], flat bands in continuous
periodic systems were studied mainly in the tight-binding limit allowing for
approximating such systems by discrete lattices. The nonlinear effects of
genuinely continuous systems featuring one or a few flat bands well separated
from the rest of the spectrum did not receive attention yet. Meantime there is
an essential subtle difference between the concept of flat bands in discrete
systems and continuous systems with periodic potentials. A discrete system
bearing an ideally flat band, which is obtained as a tight-binding limit of
the respective continuous model features nonzero inter-site coupling which is
determined by the Fourier coefficients of the spectrum (see e.g. [13]). In
contrast, nearly perfect flat bands of a continuous system mean zero
coefficients of the Fourier expansion, and thus completely decoupled sites, in
the tight-binding limit (see the discussion below). Furthermore, the existence
of perfectly flat bands in continuous systems is impossible, i.e., the concept
of flatness itself must be specified.

The existence of nearly flat bands is a characteristic feature of continuous
spin-orbit-coupled BECs (SO-BECs) with judiciously chosen parameters [14, 15,
16, 17]. Discrete models describing a SO-BEC within the framework of a tight-
binding limit were considered in [18, 19, 20]. Being created in atomic gases
[21], SO-BECs remain in the focus of current experimental studies [22, 23].
However, the existence and properties of gap soliton bifurcating from a flat
band in continuous models have not been addressed.

In this work, we are interested in a minimal continuous SO-BEC model with
periodic coefficients (described in Sec. II) which features a single nearly
flat band that is separated by a finite gap from the rest of the spectrum. We
obtain numerically and describe analytically (in Sec. III) soliton families
which in the limit of small amplitudes are states which are well approximated
by the Wannier functions (WFs), differing from families bifurcating from the
Bloch states of non-flat bands, as it happens for the conventional gap
solitons. Such matter wavepackets remain localized even when their norm
becomes negligibly small and are termed [13] as Wannier solitons (WSs).
Another peculiarity of the WSs (addressed in Sec. IV) is that their enhanced
stability allows for the construction of stable multi-hump Wannier breathers.
Such excitations alternatively can be viewed as analogs of bosonic Josephson
junctions (BJJ) based on two, three, or even more BEC clouds, in which,
however, the atomic exchange among different clouds is enabled by inter-atomic
interactions, and hence is not observed in the linear regime. The outcomes are
summarized in the Conclusion.

##  II Flat bands enabled by spin-orbit coupling

We consider a quasi-one-dimensional SO-BEC loaded in an OL [24] which in the
mean-field approximation is described by the dimensionless Gross-Pitaevskii
(GP) equation for the spinor
${{\bm{\Psi}}}=\left({\Psi}_{1},{\Psi}_{2}\right)^{T}$:

$\displaystyle
i\partial_{t}{\bm{\Psi}}=H{\bm{\Psi}}+g\left({\bm{\Psi}}^{\dagger}{\bm{\Psi}}\right){\bm{\Psi}},$
(1)

where

$\displaystyle
H=-\frac{1}{2}\partial_{x}^{2}-i\gamma\sigma_{z}\partial_{x}+\frac{\Omega}{2}\sigma_{x}+V_{0}\sin^{2}x,$
(2)

is the linear Hamiltonian, $\gamma$ is the SOC strength, $\Omega$ is the Rabi
frequency, $\sigma_{x,y,z}$ are the Pauli matrices, $V_{0}$ is the amplitude
of the OL and $g$ is the effective atomic interaction. The units of the energy
and the spatial length are chosen as $2E_{L}$ and $1/k_{L}$, respectively.
Here, $E_{L}=\hbar^{2}k_{L}^{2}/2m$ is the recoil energy of the OL, $m$ is the
atom mass, and $k_{L}$ is the wave number of the lattice beams. The
dimensionless order parameter is normalized
$N=\int_{-\infty}^{\infty}{\bm{\Psi}}^{\dagger}{\bm{\Psi}}dx$ with the norm
$N$ related to the physical number of atoms $\mathcal{N}$:
$\mathcal{N}=NN_{0}$, where $N_{0}=gE_{L}/(\hbar\omega_{r}k_{L}a_{s})$,
$\omega_{r}$ is the trap frequency along the transverse directions, and
$a_{s}$ is the s-wave scattering lengths.

Considering the typical experimental data for a 87Rb BEC from Ref. [24]:
$\omega_{r}=2\pi\times 150\,$Hz, $a_{s}=5.29\,$nm, $k_{L}=2.88\mu\,{\rm
m}^{-1}$, $E_{L}\approx 3.2\times 10^{-31}$J, $\Omega=7.62E_{L}$ and
$V_{0}=-1.4E_{L}$, we obtain the dimensionless parameters
$\gamma=k_{R}/k_{L}=1.96$, with $k_{R}=5.63\,\mu{\rm m}^{-1}$ being the wave
number of the Raman beams, $\Omega=3.81$, $V_{0}=-0.7$ and $N_{0}\approx 217$.
In experiments, $\gamma$ can be varied either by changing the angle between
the incident Raman beams or through fast modulation of the intensities of the
Raman lasers [25, 26]. The Rabi frequency $\Omega$ and the amplitude of the OL
$V_{0}$ are free parameters that can be tuned by changing the intensity of
Raman and lattice beams, respectively.

It is known [14, 15, 16, 17] that $\gamma$ and $\Omega$ considered as control
parameters can be chosen to ensure the extreme flatness of the lowest band of
the linear spectrum $\mu_{\nu}(k)$ of the Hamiltonian $H$, defined by the
eigenvalue problem $H{\bm{\varphi}}_{\nu k}=\mu_{\nu}(k){\bm{\varphi}}_{\nu
k}$, where ${\bm{\varphi}}_{\nu k}(x)$ $\nu$ is the band index, and
$k\in[-1,1)$ is the Bloch wavenumber. Indeed, consider first the case where
$\Omega=2\gamma^{2}\gg 1$. Neglecting the OL in the leading approximation
obtain the linear dispersion relation for the lower branch $\mu_{-}(k)$ of the
spectrum:
$\mu_{-}=-\Omega/2\left[1-k^{4}/(2\Omega^{2})+\mathcal{O}(k^{6}/(2\Omega^{3}))\right]$.
Thus, for any given finite interval of wavenumbers, $k$, one can choose
$\Omega$ and $\gamma$ large enough to make as many orders of the dispersion
negligible as necessary. In particular, leaving only the leading term $\propto
k^{4}$ in the mentioned expansion one can look for quartic SO-BEC solitons as
approximate solutions. Adding now an OL that restricts the wavenumbers to the
first Brillouin zone, should result in the nearly flat lowest band.

There is however a limitation to such a straightforward approach. While
formally soliton families can bifurcate from such a flat band towards the
semi-infinite gap (what happens in the case of attractive condensates), the
finite gaps created by a very shallow OL are also negligible, i.e.,
practically no soliton families in finite gaps can be found. Therefore, we
start by briefly describing a situation where a nearly flat band is separated
from the rest of the spectrum by a non-vanishing finite gap. While the very
existence of such SO-BECs was illustrated in [14, 15], Fig. 1 we show that
flat bands can be encountered for a sufficiently wide range of physically
accessible parameters.

Figure 1: In the Rabi frequency $\Omega$ (a) and the gap width $D$ between the
flat lowest and first excited bands (black) and the maximum lowest band
flatness $\Delta_{0}$ (red) in (b) versus the SOC strength $\gamma$ with
$V_{0}=2$. The linear Bloch spectrum is illustrated for the points marked by
red dots and the respective letters for arbitrarily chosen parameters
$\gamma=2.01$, $\Omega=11.4844$ in (c) and $\gamma=4$, $\Omega=41.4649$ in
(d).

To characterize the width of the lowest band $\mu_{0}(k)$ we represent

$\displaystyle\mu_{0}(k)=\mu_{0}+\tilde{\mu}(k)$ (3)

where $\mu_{0}$ without argument stands for $\mu_{0}(0)$ (recall that
$k\in[-1,1)$). In a continuous model, exactly flat bands, i.e., bands with
$\tilde{\mu}(k)\equiv 0$, are impossible. Thus, the flatness of the lowest
band ($\nu=0$) implies a smallness of $\tilde{\mu}(k)$, which in this work,
will be quantified by the parameter

$\displaystyle\Delta_{0}=\max\mu_{0}(k)-\min\mu_{0}(k).$ (4)

We say that the band is nearly flat (or flat, for brevity) if
$\Delta_{0}\lesssim 10^{-4}$. The so-defined criterion is not strict and can
be modified. Our choice, in addition to ensuring the smallness of the
bandwidth, also aims to optimize several requirements, including the
possibility of operating within a sufficiently wide range of parameters (see
Fig. 1 below), to enlarge the region of validity of the estimate
$N\propto\mu$, and to squeeze the transition region between this region and
the edge of the linear band, holding at the same time a sufficient fidelity of
the numerical results. On the other hand, to characterize the size of the
lowest finite gap we introduce the parameter
$D=\min\mu_{1}(k)-\max\mu_{0}(k)$.

In Fig. 1(a) we illustrate the parameter dependence $\Omega(\gamma)$ in the
range of SOC strengths where the flat band coexists with a sufficiently large
gap [see $\Delta$ (red) and $D$ (black) curves in Fig. 1(b)]. The examples of
the two lowest bands, panels (c) and (d), illustrate the features of the SO-
BEC in a finite-depth lattice (in all calculations, we fix $V_{0}=2$)
featuring the only flat band, all other hands being non-flat. The flatness of
the shown lowest bands in the absence of SOC, i.e., at $\gamma=0$, is $0.17$.
Such band-gap structures significantly differ from the typical models of BEC
in deep optical lattices with several flat lower bands (see e.g. [13]). Note
also, that the flat lowest band of a SO-BEC without periodic potential would
require much higher values of the Rabi frequency (according to the above
discussion).

The conventional gap solitons, described by the NLS equation, require either
$|\mu_{0}^{\prime\prime}(k_{0})|\gg|\mu_{0}^{(m)}(k_{0})|$, where $k_{0}$ is a
Bloch wavenumber at which the edge of a finite or semi-infinite gap is
achieved. This allows one to neglect (or to account perturbatively) the higher
dispersion. When $\mu_{0}^{\prime\prime}(k_{0})$ is comparable with the forth-
order dispersion $\mu_{0}^{(4)}(k_{0})$ artificially induced by the OL, one
still can obtain gap solitons (such solitons were considered for optical
applications [27]). A real-world BEC, available in an experiment has a finite
extent. Hence below a certain amplitude, a conventional gap soliton with a
wide (exceeding the real size of the condensate) envelope of the respective
Bloch state cannot be observed experimentally. The case considered here is the
opposite one: in the leading approximation all orders of the dispersion of the
lowest band are negligible, i.e., the approximation of an ideally flat band is
valid. As we will show below, in the leading order widths of such solitons
remain finite until very small amplitudes, i.e., very close to the linear
limit.

##  III Wannier solitons

###  III.1 Small-amplitude expansion.

Having the Bloch basis given, one can define [28] the two-component WFs. For
the lowest flat band, one has

$\displaystyle{\bm{w}}_{n}(x)=\frac{1}{2}\int_{-1}^{1}{\bm{\varphi}}_{0k}(x)e^{-i\pi
kn}dk.$ (5)

where $n$ numbers the minima of the OL. Applying $H$ to the WF we obtain

$\displaystyle H{\bm{w}}_{n}=\mu_{0}{\bm{w}}_{n}(x)+\sum_{m\neq
n}f_{n-m}{\bm{w}}_{m}(x),$ (6)

where

$\displaystyle f_{m}$
$\displaystyle=\frac{1}{2}\int_{-1}^{1}\tilde{\mu}(k)e^{ikm\pi}dk$ (7)

and it is used that due to $\pi$ periodicity $\tilde{\mu}(k)=\tilde{\mu}(k+2)$
[see (3)]. It is relevant to mention that the last term in Eq. (6) while
having a small amplitude determined by the flatness of the band, is
delocalized in space, unlike the first term involving ${\bm{w}}_{n}(x)$, which
is localized in the vicinity of the points $x_{n}=n\pi$.

To describe the soliton families bifurcating from the flat band, we recall
that all $f_{m}$ are nearly zero (zero would correspond to an ideally flat
band). Thus, assuming (without loss of generality) that a soliton family
bifurcates from the linear mode centered at $x=0$, we look for a solution of
Eq. (1) in the form of the "mixed" (Wannier– for the lowest band and Bloch–
for all other bands) expansion

$\displaystyle\Psi\left(x,t\right)=\sqrt{\epsilon}e^{-i\mu
t}\left[a\left(t\right){\bm{w}}_{0}\left(x\right)+\epsilon{\bm{\psi}}_{1}+\mathcal{O}\left(\epsilon^{2}\right)\right],$
(8)

where $\epsilon\ll 1$ is a formal small parameter, $\mu=\mu_{0}+f_{0}$,

$\displaystyle{\bm{\psi}}_{1}$ $\displaystyle=\sum_{n\neq
0}b_{n}(t){\bm{w}}_{n}(x)+\sum_{\nu=1}^{\infty}\int_{-1}^{1}b_{\nu
k}(t){\bm{\varphi}}_{\nu k}(x)dk,$ (9)

and only the lowest-order terms are accounted for. Substituting expansions (8)
and (9) in Eq. (1), and projecting on the WFs ${\bm{w}}_{0}$ and
${\bm{w}}_{n}$ (with $n\neq 0$) withing the accepted accuracy we obtain
respectively:

$\displaystyle i\partial_{t}a=\epsilon\chi_{0}|a|^{2}a+\epsilon\sum_{n\neq
0}b_{n}f_{n},$ (10) $\displaystyle\epsilon
i\partial_{t}b_{n}-\epsilon\sum_{m\neq
0}b_{m}f_{m-n}=af_{-n}+\epsilon\chi_{n}|a|^{2}a,$ (11) $\displaystyle
i\partial_{t}b_{\nu k}+[\mu-\mu_{\nu}(k)]b_{\nu k}=\chi_{\nu k}|a|^{2}a,$ (12)

where

$\displaystyle\chi_{n}=$ $\displaystyle
g\langle{\bm{w}}_{n},({\bm{w}}_{0}^{\dagger},{\bm{w}}_{0}){\bm{w}}_{0}\rangle,$
(13) $\displaystyle\chi_{\nu k}=$ $\displaystyle g\langle{\bm{\varphi}}_{\nu
k},({\bm{w}}_{0}^{\dagger},{\bm{w}}_{0}){\bm{w}}_{0}\rangle,$ (14)

are the effective nonlinear coefficients, and we use the notation
$\langle{\bm{f}},{\bm{g}}\rangle=\int_{-\infty}^{\infty}{\bm{f}}^{\dagger}(x){\bm{g}}(x)dx$.
The solutions described by (8) belong to the family of WSs.

The expansion (8) does not allow for the formal limit $\epsilon\to 0$ because
${\bm{w}}_{0}(x)$ is not an eigenstate of the linear Hamiltonian. However, one
can show that for an interval of amplitudes, parametrized by a small parameter
$\epsilon$:

$\displaystyle\Delta_{0}\lesssim\epsilon\ll 1,$ (15)

formula (8) accurately describes the family of solutions. Indeed, subject to
the assumption $|b_{n}f_{n}|\sim|b_{n}|\Delta_{0}\ll|a(t)|$ [it is considered
that $a(t)\sim 1$], Eq. (10) is readily solved

$\displaystyle a(t)=a_{0}e^{-i\epsilon\chi_{0}|a_{0}|^{2}t},$ (16)

where $a_{0}$ is the initial amplitude. Then, in the l.h.s. of Eq. (11) one
can neglect the hopping terms and obtain

$\displaystyle\epsilon
b_{n}\approx\left(\frac{f_{-n}}{\epsilon\chi_{0}|a_{0}|^{2}}+\frac{\chi_{n}}{\chi_{0}}\right)\left(e^{-i\epsilon\chi_{0}|a_{0}|^{2}t}-1\right).$
(17)

Since $|\chi_{n}|\ll\chi_{0}$ due to the localization of the WFs, one
concludes that limiting inequality $\epsilon|b_{n}|\ll|a_{0}|$ remains valid
for times $t\ll 1/\Delta_{0}$ provided (15) is satisfied. Imposing the
condition $\Delta_{0}\ll\epsilon$ [cf. (15)] one obtains that the required
smallness of $\epsilon|b_{n}|$ holds for any time. Finally, we notice that it
follows from (12) that $b_{\nu k}$ terms do not display secular growth, and
thus give a small correction to the leading term in the expansion (8).

In the limit $\epsilon\to 0$, or more precisely when $\epsilon\ll\Delta_{0}$,
expansion (8) is not valid, because of the secular term in (11). To find the
soliton family in this limit one should consider the conventional multiple-
scale expansion [5, 1] characterized by the scaling
$N\sim\sqrt{|\mu-\mu_{0}|}\sim\epsilon$. It turns out that such scaling also
implies that the width of the soliton is of the order of $1/\epsilon$, which
in the case at hand corresponds to $10^{4}$ (recall that according to our
definition, the flatness of a band requires $\Delta_{0}\lesssim 10^{-4}$) or
to even more lattice periods in the dimensionless units. Meantime, lengths of
experimentally feasible condensates typically are of the order of $500\,\mu m$
[24], which is about $10^{3}$ in the dimensionless units. Thus the extremely
tiny "boundary" domain in the vicinity of the flat-band edges, which is
excluded from our consideration here, is also beyond the physical
accessibility in a real system.

###  III.2 Wannier soliton families

The condition (15) corresponds to a relatively wide interval of the system
parameters. In particular, it is valid for the whole range of SOC strength
shown in Fig. 1 (a) and (b). In the domain specified by (15) the family is
characterized by the linear dependence $N=\chi_{0}(\mu-\tilde{\mu}_{0})$. (As
it is customary, since $\epsilon$ is a formal small parameter, in the final
simulations it is set equal to one while $a_{0}$ is considered to be small).
This is verified in Fig. 2 using the numerical solution of Eq. (1) by means of
the Newton relaxation as well as the Difference methods. In the numerical
calculation, we used the step $dx=\pi/2^{11}$ (verifying the convergence for
the smaller steps).

Figure 2: Families of WSs are shown for $g=-0.1$ (in the interval
$\mu<\mu_{0}$) and $g=0.1$ (in the interval $\mu>\mu_{0}$). The position of
the flat band with $\mu_{0}\approx-5.1608$ and $\Delta_{0}\approx 2.8\cdot
10^{-6}$, which is not distinguishable on the scale of the plot is indicated
by the vertical grey line. Solid (dashed) lines are used for stable (unstable)
families bifurcating at $(\mu,N)=(\mu_{0},0)$ from the Bloch mode. Notice that
the b and c families are nearly indistinguishable on the scale of the figure.
In both panels, lines are for numerical results of Eq. (1), while the lines
with circles are for the results obtained from the approximation Eq. (10).
Different colors, also marked by different letters, correspond to the soliton
families whose profiles are illustrated in Fig. 3. The ln-ln plot in the inset
zooms in the region in the vicinity of the flat band $\mu_{0}$ for one-soliton
families bifurcating from the flat band in the limit $N\to 0$ for attractive
(the red line) and repulsive (the black line) condensates (the linear band in
the inset is located at $-\infty$ of the abscissa). Here and in all figures
below, unless specified otherwise, we set $V_{0}=2,\gamma=2.01$,
$\Omega=11.4844$ corresponding to point c in Fig. 1(a).

In Fig. 2, the one-soliton [an example is illustrated in Fig. 3(a)] families
for the attractive $g=-0.1$ ($\mu<\mu_{0}$) and repulsive $g=0.1$
($\mu>\mu_{0}$) nonlinearities are shown by blue solid lines. Each family
bifurcates from the linear band shown by a thin vertical line in Fig. 2. This,
however, occurs in the anomalously narrow vicinity of the band, which is a
peculiarity of the flatness of the band. Since the mentioned domain is not
visible on the scale of the main panel of the figure, in the inset, we zoom
this region showing it on the ln-ln plot. Clearly, the linear band on such a
plot is located at $-\infty$ of the abscissa, i.e., the bifurcation point
itself cannot be shown graphically. Beyond the tiny transition region near
$N=0$, these families are characterized by the predicted linear dependence
$N=\chi_{0}(\mu-\mu_{0})$ with the slopes $\chi_{0}>0$ and $\chi_{0}<0$ for
the repulsive and attractive BECs. In the region close to $\mu=\mu_{0}$
(zoomed in the inset), the family behaves as $N\sim\sqrt{|\mu-\mu_{0}|}$ and
the solitons become extremely wide at $\mu\to\mu_{0}$, thus breaking the
acceptable numerical accuracy. The transition from the gap-soliton interval of
the family to the WS linear dependence is confirmed numerically. The observed
localization of WSs indicates on the possibility of their experimental
observation in low-density condensates, where conventional gap solitons would
have extents much larger than the sizes of the typically used trapping
potentials.

Figure 3: Soliton profiles illustrated for the families indicated in Fig. 2 by
the respective letters, obtained for $g=-0.1$ arbitrarily chosen value
$\mu=-5.2217$ (solid lines), compared with the linear WFs (dots). In the left
(right) column we show the real and imaginary (absolute values) of the first
component $\Psi_{1}$. Panels (a,b) show fully symmetric solitons while panels
(c,d) show solitons with broken symmetry (see the text).

A one-peak soliton is confined to a unit cell provided the band is flat. A
numerical illustration of this feature is given in Fig. 4 where we show the
localization of WSs characterized by the inverse participation ratio (IPR):

$\displaystyle\text{IPR}=\frac{1}{N^{2}}\int_{-\infty}^{\infty}({\bm{\Psi}}^{\dagger}{\bm{\Psi}})^{2}dx\approx\frac{|\chi_{0}|}{gN^{2}}$
(18)

as a function of the flatness parameter $\Delta_{0}$ (notice the logarithmic
scale of the abscissa). To this end, we fix the norm $N=1$ and $\gamma=2.01$
and vary $\Omega$ between 0 and 14.48, i.e., from a conventional non-flat band
to (nearly) ideally flat band limits. One can clearly indicate the nonlinear
mobility edges separating areas of the localized and delocalized nonlinear
solutions (the conventional linear mobility edge does not exist since our
system is periodic). Delocalized gap solitons, IPR$\ll 1$, near a non-flat
band, become well-localized, IPR$\approx 0.6$, above the nonlinear mobility
edge even for weak nonlinearity $g=10^{-3}$ (the earthy yellow lines with
rotated triangles).

The widths of WSs remain practically unchanged in a wide interval of the
chemical potentials or nonlinearity strengths. In Fig. 5 we plot the
projection of a WS on the WF ${\bm{w}}_{0}$:
$P_{W}=\langle{\bm{w}}_{0},{\bm{\psi}}\rangle$ versus the nonlinearity
coefficient $g$ (for the fixed norm $N=1$). Since the states are strongly
localized in the vicinity of a local minimum of the potential, the quantity
$1-P_{W}$ characterizes the population of the upper non-flat bands. We observe
that in a sufficiently large interval of magnitude $g$, the projection $P_{W}$
is nearly one [Fig. 5 (a)] and it rapidly decreases in the close vicinity of
the flat band [Fig. 5 (b)] where the expansion (8) is not applicable.

Figure 4: IPR of the fundamental WS versus ${\ln}\Delta_{0}$ with
$\gamma=2.01$ for different non-linearity strengths $g$. Light cyan shadow
area corresponds to $\ln\Delta_{0}\approx-1.758$ corresponding to $\Omega=0$.

We have studied both, the linear stability of WSs by solving the standard
Bogoliubov–de Gennes equations, as well as the (in)stability obtained by
simulating their long-time behavior within Eq. (1) with the addition of
initial Gaussian noise of the order of 10% of the soliton profile. Both
approaches manifested fully consistent results. The results of this analysis
are summarized in Fig. 2 where stable (unstable) families are shown by solid
(dashed lines). The family of one-hump WSs is found stable in the whole region
of the parameters where the approximation $N\approx\chi_{0}(\mu-\mu_{0})$
holds.

Figure 5: The projection $P_{W}$ of the fundamental WS on the WF
${\bm{w}}_{0}$ vs. $g$ (a) and its zoom (shown using the logarithmic scale) in
the proximity of the flat band (b) for $N=1$.

In addition to the discrete translation symmetry the Hamiltonian (2), and thus
the GPE (1), has a Klein four–group of symmetries
$\\{1,\hat{\alpha}_{1},\hat{\alpha}_{2},\hat{\alpha}_{3}\\}$, where
$\hat{\alpha}_{1}=-\mathcal{P}\sigma_{x}$,
$\hat{\alpha}_{2}=-\mathcal{K}\sigma_{x}$,
$\hat{\alpha}_{3}=\mathcal{P}\mathcal{K}$, $\mathcal{P}$ is the inversion
operator, and $\mathcal{K}$ is the complex conjugation operator. The shown
one-hump WSs obey all mentioned symmetries. Meantime, the strong localization
of WSs suggests the existence of multi-hump solitons, which loosely can be
understood as bound states of the fundamental WSs and may obey either all or
only one of the symmetries $\hat{\alpha}_{j}$. Three examples are shown in
Fig. 3(b-d). A two-hump WS belonging to the fully-symmetric soliton family,
i.e., obeying the mentioned Klein four-group symmetries, is shown in Fig.
3(b). The soliton shown in Fig. 3(c) obeys only $\hat{\alpha}_{1}$ symmetry,
thus representing nonlinear symmetry breaking (solitons belonging to this
family are anti–$\hat{\alpha}_{2}$ and anti–$\hat{\alpha}_{3}$ symmetric).
Respectively, these families are characterized by nearly coinciding
dependencies $N(\mu)$ (cf. red solid line b and green dashed line c in Fig. 2
where they are practically indistinguishable). Every two-hump WS carries
nearly double the number of atoms compared with the one-hump soliton at the
same value of the chemical potential. We have found that the fully-symmetric
branch b, characterized by the linear dependence $N\sim\mu$, is stable.
Meantime, the $\hat{\alpha}_{1}$-symmetric branch c is unstable over the whole
domain of chemical potentials shown in Fig. 2 (except the very narrow region
in the vicinity of $\mu_{0}$ shown in the inset, which was not investigated
because of numerical constraints).

Finally, in Fig. 3(d) we illustrate a three-hump WS obeying only
$\hat{\alpha}_{3}$ symmetry ($\hat{\alpha}_{2}$ and $\hat{\alpha}_{3}$
involving inversion $\mathcal{P}$ being broken). Such solitons existing at
relatively weak nonlinearities are enabled by flat bands, but we found them
all to be unstable (see the respective branch in Fig. 2). Generally, the
obtained WSs emerging as strongly loclaized extra objects on the top of the
periodic OL contrast with gap solitons in the non-flat band whose internal
structure inherits the lattice periodicity [15].

##  IV Wannier breathers

###  IV.1 A few mode approximation

The enhanced stability of the WSs as well as their existence in the "quasi-
linear" limit determined by the relations (15), suggests a way of constructing
oscillatory solutions. Indeed, "trains" of WFs, similar to those shown in Fig.
3 (b)-(d), can be viewed as superpositions of single WSs. If the amplitudes of
the humps in such a train are weakly modulated, i.e., the train is not an
exact nonlinear solution, it is no longer a stationary solution too, and the
atomic transfer between the localized clouds can occur. In the case of two
coupled WSs such a setting is alike to a BJJ [29, 30, 31] and especially to
the realization of a BJJ in a quasi-periodic potential reported recently in
[32]. There is however an essential difference between the results reported
below with the previous studies of BJJ: due to the orthogonality of the WFs,
the coupling in our model is exclusively nonlinear (it vanishes in the linear
limit), i.e., the reported oscillating objects are rather Wannier breathers
than linearly coupled atomic clouds.

Nevertheless, the analogy with BJJ suggests that the theoretical approach for
the description of the Wannier breathers can use replication of the general
approach elaborated in the mean-field theory of BJJ [33]. Limiting the
consideration to the cases of two- and three-hump breathers with humps in the
nearest neighbor potential minima (the generalizations are straightforward),
we look for the leading order of a solution of the GPE (1) in the form of the
three-mode ansatz

$\displaystyle{\bm{\Psi}}=e^{-i\tilde{\mu}t}\sum_{n=-1}^{1}A_{n}(t){\bm{w}}_{n}(x).$
(19)

Here $A_{n}$ is time-dependent amplitudes of the WF ${\bm{w}}_{n}$ and
excitation of upper bands is neglected. This implies the conservation of the
norm
$\left|A_{-1}\right|^{2}+\left|A_{0}\right|^{2}+\left|A_{1}\right|^{2}=N$. The
dynamical equations for the time-dependent amplitudes are obtained in a
standard way (similar calculations for four coupled modes can be found say in
[32]). In terms of canonically conjugate variables, which are the phase
differences $\phi_{\pm}=\arg A_{0}-\arg A_{\pm 1}$ and the fractions of the
total number of atoms $N_{\pm}={\left|A_{\pm 1}\right|^{2}}/N$, these
equations are obtained from the Hamiltonian

$\displaystyle H=$
$\displaystyle-\frac{\chi_{0}}{2}\left[N_{+}^{2}+N_{-}^{2}+\left(1-N_{+}-N_{-}\right)^{2}\right]$
$\displaystyle-2\chi_{1}\sqrt{1-N_{+}-N_{-}}\left[\left(1-N_{-}\right)\sqrt{N_{+}}\cos\phi_{+}\right.$
$\displaystyle\left.+\left(1-N_{+}\right)\sqrt{N_{-}}\cos\phi_{-}\right].$
(20)

and read

$\displaystyle\frac{dN_{\pm}}{d\tau}=$ $\displaystyle
2\chi_{1}\left(1-N_{\mp}\right)\sqrt{N_{0}N_{\pm}}\sin\phi_{\pm}\,,$ (21)
$\displaystyle\frac{d\phi_{\pm}}{d\tau}=$
$\displaystyle\left(N_{\mp}-N_{0}\right)\left(\chi_{0}-\chi_{1}\frac{1-N_{\mp}}{\sqrt{N_{0}N_{\pm}}}\cos\phi_{\pm}\right)$
$\displaystyle-\chi_{1}\left(3\sqrt{N_{0}N_{\mp}}+\frac{N_{\mp}^{2}}{\sqrt{N_{0}N_{\pm}}}\right)\cos\phi_{\mp},$
(22)

with $\tau=Nt$ and $N_{0}=1-N_{-}-N_{+}$. We emphasize that the Hamiltonian
(IV.1), having apparent similarity with the three-level models governing atoms
in three-well traps [34, 35], describes essentially different dynamics,
because the coupling of BEC clouds in our case is non-linear.

###  IV.2 Evolution of a two-hump breather.

First, we consider the dynamics of a two-hump breather, which allows for
direct comparison with the standard mean-field model for the BJJ [33, 31]. To
this end, we set $A_{1}=0$. Now the population imbalance $z$ and the relative
phase $\phi$ are given by
$z=(\left|A_{-1}\right|^{2}-\left|A_{0}\right|^{2})/(\left|A_{-1}\right|^{2}+\left|A_{0}\right|^{2})$
and $\phi=\arg A_{0}-\arg A_{-1}$. One can obtain the reduced Hamiltonian from
Eq. (IV.1):

$\displaystyle H=-2\chi_{1}\sqrt{1-z^{2}}\cos\phi-\frac{1}{2}\chi_{0}z^{2},$
(23)

and dynamical equations

$\displaystyle\frac{dz}{d\tau}$
$\displaystyle=2\chi_{1}\sqrt{1-z^{2}}\sin\phi,\quad$ (24)
$\displaystyle\frac{d\phi}{d\tau}$
$\displaystyle=\chi_{0}z-2\chi_{1}\frac{z\cos\phi}{\sqrt{1-z^{2}}}.$ (25)

These equations are akin to the usual two-mode models, but the coupling and
hence the dynamics disappear in the linear limit (recall that
$\chi_{0,1}\propto g$). Furthermore, from the explicit expression for the
stationary points ($0,0$), ($0,\pm\pi$), and
($\pm\sqrt{1-\Lambda^{2}},\pm\pi$) where $\Lambda=2\chi_{1}/\chi_{0}$, one can
see that a self-trapping solution, i.e., a stationary point with imbalanced
populations $|z|=\sqrt{1-\Lambda^{2}}$, is determined by the properties of the
linear system and is independent of the nonlinearity since $\Lambda$ is
computed using the linear Wannier functions [see (13)].

Figure 6: Periods, $T$, of two-hump breathers versus the initial population
imbalance $z_{0}$ (a) and versus interaction strength $g$ (b) for
$\phi_{-}=0$, ${\chi_{0}/g\approx 0.6057}$ and ${\chi_{1}/g\approx-1.064\cdot
10^{-3}}$. Black lines are analytical results from Eqs. (24), (25), red
(repulsive) and blue (attractive) lines are numerical results from the direct
numerical solutions of the GPE (1). The population imbalance corresponding to
the separatrix is ${\color[rgb]{0,0,0}{z_{*}\approx 0.118}}$. In panel (b),
the grey lines are for transition points, and solid and dashed lines are for
${\color[rgb]{0,0,0}{z_{1}\approx 0.113<z_{*}}}$ and
${\color[rgb]{0,0,0}{z_{2}\approx 0.124}}>z_{*}$, respectively, shown in panel
(a) with light cyan lines. Two types of evolution at the points indicated in
panel (a) are explored in panels (c) (oscillations of the BJJ type) and (d)
(oscillations in the self-trapping regime). In both panels
$g=0.1,A_{1}=0,A_{0}=1$, while $A_{-1}=0.95$ in (c) and $A_{-1}=0.85$ in (d).

The dynamical regimes governed by the three-mode model (IV.1)–(IV.1) are
insensitive to the sign of interactions (up to the time reversal): the phase
portraits governed by Hamiltonian (23) coincide for positive and negative $g$.
Each phase portrait has a separatrix, with the imbalance denoted by $z_{*}$,
at which the period of oscillations diverges [see the black solid line in Fig.
6(a)]. However, the realistic dynamics, while preserving qualitative
separation between "Josephson oscillations", i.e., oscillations around $z=0$
[Fig. 6(c)], and oscillations in the vicinity of the "self-trapping" point
[Fig. 6(d)] manifests differences in the vicinity of the separatrix. These
differences are illustrated by blue and red dashed-dotted lines in Fig. 6(a).

Since the oscillations shown in Figs. 6(c) and (d) are due to nonlinear
coupling, they disappear in the linear limit. This is expressed in Fig. 6(b)
by the divergences of all the curves $T(|g|)$ in the limit $g\to 0$.
Furthermore, direct numerical simulations show that not all features of the
dynamics in the vicinity of the separatrix can be captured by the two-mode
model. In Fig. 6(b) one observes divergence of the periods of oscillations in
the points $g=-g_{1}$ with initial imbalance $z_{1}<z_{*}$ (solid blue line)
and $g=g_{2}$ with initial imbalance $z_{2}>z_{*}$ (red dashed line). These
singularities of the real dynamics can be explained by the "shift of the
separatrix" towards smaller and bigger imbalances $z_{1}$ and $z_{2}$,
respectively, as shown in Fig. 6(a). Thus, when $|g|$ increases at a fixed
population imbalance, the separatrix can be crossed either only by the
attractive condensate with the imbalance $z_{1}$ or only by the repulsive
condensate with the imbalance $z_{2}$. From the physical point of view, the
sensitivity of the system in the vicinity of the separatrix can be explained
by deviations of the localization, and thus of effective interactions, of the
WSs upon the effect of either attractive or repulsive nonlinearity.

Benefited from the enhanced stability of two-hump WSs, two-hump breathers
appear to be stable, verified by directly simulating their long-time (up to
$t=40000$) behavior with the addition of 5% Gaussian noise to the initial two-
hump profile. We also examined the stability for $|g|=1$ and found that such
breathers remain stable having relatively large amplitudes.

###  IV.3 On the evolution of the three-hump breather

The dynamics of a three-hump breather even within the framework of the three-
mode approximation described by Eqs. (IV.1)–(IV.1), where the central mode
${\bm{w}}_{0}$ is coupled to the two nearest-neighbor modes ${\bm{w}}_{\pm
1}$, can be remarkably complex. Therefore, to illustrate the existence of a
stable three-hump breather we address only the simplest case where
$N_{+}=N_{-}$, $\phi_{+}=\phi_{-}=\phi$. Within the framework of the mode
expansion, this assumption leads to a reduced Hamiltonian that is obtained
from Eq. (IV.1)

$\displaystyle
H=-\chi_{0}(3z^{2}-2z)-4\chi_{1}\sqrt{2(1-z^{2})}\left(3-z\right)\cos\phi,$
(26)

where
$z=(2\left|A_{-1}\right|^{2}-\left|A_{0}\right|^{2})/(2\left|A_{-1}\right|^{2}+\left|A_{0}\right|^{2})$
[cf. (23)]. The equations of motion read

$\displaystyle\frac{dz}{d\tau}$
$\displaystyle=4\chi_{1}\sqrt{2(1-z^{2})}(3-z)\sin\phi,$ (27)
$\displaystyle\frac{d\phi}{d\tau}$
$\displaystyle=2\chi_{0}(3z-1)+\frac{8\chi_{1}(2z^{2}-3z-1)}{\sqrt{2(1-z^{2})}}\cos\phi.$
(28)

Unlike in the case of a two-mode breather, now the fixed points which are
found as the roots of the polynomial

$\displaystyle\sqrt{1-z^{2}}(3z-1)\pm\sqrt{2}\Lambda(2z^{2}-3z-1)=0,$ (29)

where "$+$" and "$-$" stand for $\phi=0$ and $\phi=\pm\pi$, do not correspond
to equal populations of all three potential minima. In particular, for the
choice of the parameters used in Fig. 7 the four stationary points are found
numerically to be approximately $(0.33,0)$, $(0.34,\pm\pi)$, and $(\pm
1,\pm\pi)$.

In spite of these differences in the Hamiltonians (23) and (26), the latter
showing much richer dynamics, as it is seen from the set of the stationary
points, in Fig. 7 (a) and (b) we observe that the three-mode breather share
the same peculiarities of evolution as those discussed above in detail for the
two-mode case [cf. Fig. 6 (a) and (b)]. Examples of the numerical study of the
respective three-mode dynamics governed by the original model (1) are shown in
Fig. 7. Now, we again observe enhanced stability of the three-hump breathers
in different dynamical regimes, which was verified by considering the
evolution for a long time (up to $t=40000$).

Figure 7: Periods, $T$, of three-hump breathers versus the initial population
imbalance $z_{0}$ (a) and versus interaction strength $g$ (b) for
$\phi_{\pm}=0$, ${\chi_{0}/g=0.6067}$ and ${\chi_{1}/g=-1.064\cdot 10^{-3}}$.
Black lines are analytical results from Eq. (26). Red (repulsive) and blue
(attractive) lines are numerical results from direct numerical solutions of
GPE (1). The population imbalance corresponding to the separatrix is
${\color[rgb]{0,0,0}{z_{*}\approx 0.471}}$. In panel (b), the grey lines are
for transition points, and solid and dashed lines are for
${\color[rgb]{0,0,0}{z_{3}\approx 0.459}}<z_{*}$ and
${\color[rgb]{0,0,0}{z_{4}\approx 0.472}}>z_{*}$, respectively [they are shown
in panel (a) with light cyan lines]. In panels (c) and (d) two types of
evolution corresponding to the points indicated in panel (a) are shown for
$g=0.1,A_{0}=1$, and $A_{\pm 1}=0.95$ in (c) and $A_{\pm 1}=0.8$ (d).

##  V Conclusion

In this work, we described the families of Wanneir solitons that emerge in a
mean-field model of the spin-orbit coupled BEC with a nearly flat band. Such
families are characterized by the linear dependence of the number of condensed
atoms on the detuning of the chemical potential for both attractive and
repulsive condensates. Unlike the conventional gap solitons (existing in
systems with non-flat bands) Wannier solitons have shapes described by the
Wannier functions and widths that remain practically unchanged over the
significant domain of variation for the chemical potential. When solitons
share the whole symmetries of the Hamiltonian they were found to be stable.
This allows for the construction of Wannier breathers which represent
nonlinearly-coupled condensate clouds. The evolution of the breather is well
captured by few-mode models. Since the coupling in our model is nonlinear, the
type of evolution of the breathers is determined by the number of atoms (or
alternatively, by the scattering length). Two– and three–hump breathers,
considered as examples were also found stable. The considered system also
supports families of Wannire solitons, which obey only one of the symmetries
of the nonlinear Hamiltonian, i.e., manifests the nonlinear symmetry breaking.
The respective families, however, were found to be unstable. Finally, similar
Wannier solitons can be found in other continuous systems, chiefly in optical
nonlinear stratified media, having flat bands which are well separated from
the rest of the spectrum.

###### Acknowledgements.

The work was supported by the Portuguese Foundation for Science and Technology
(FCT) under Contracts PTDC/FIS-OUT/3882/2020 and UIDB/00618/2020, the China
Scholarship Council (CSC) under the Grant CSC N.202206890002. This work was
also supported by the National Natural Science Foundation of China with Grants
No.11974235 and 11774219\.

## References

  * [1] V. A. Brazhnyi and V. V. Konotop, Theory of nonlinear matter waves in optical lattices. Mod. Phys. Lett. 18, 627 (2004).
  * [2] O. Morsch and M. Oberthaler, Dynamics of Bose-Einstein condensates in optical lattices, Rev. Mod. Phys. 78, 179 (2006).
  * [3] Y. Zhang, Z. Liang, and B. Wu, Gap solitons and Bloch waves in nonlinear periodic systems, Phys. Rev. A  80, 063815 (2009).
  * [4] Y. V. Kartashov, V. A. Vysloukh, and L. Torner, Soliton Shape and Mobility Control in Optical Lattices, Prog. Opt. 52, 63 (2009).
  * [5] V. V. Konotop and M. Salerno, Modulational instability in Bose-Einstein condensates in optical lattices, Phys. Rev. A 65, 021602 (2002).
  * [6] C. M. de Sterke, A. F. J. Runge, D. D. Hudson, and A. Blanco-Redondo, Pure-quartic solitons and their generalizations—Theory and experiments, APL Photon. 6, 091101 (2021).
  * [7] A. V. Yulin and V. V. Konotop, Conservative and $\mathcal{PT}$-symmetric compactons in waveguide networks, Opt. Lett. 17, 4880 (2013).
  * [8] R. A. Vicencio and M. Johansson, Discrete flat-band solitons in the kagome lattice, Phys. Rev. A  87, 061803 (2013).
  * [9] V. Goblot, B. Rauer, F. Vicentini, A. Le Boité, E. Galopin, A. Lemaître, L. Le Gratiet, A. Harouri, I. Sagnes, S. Ravets, C. Ciuti, A. Amo, and J. Bloch, Nonlinear Polariton Fluids in a Flatband Reveal Discrete Gap Solitons, Phys. Rev. Lett.  123, 113901 (2019).
  * [10] D. Leykam, S. Flach, O. Bahat-Treidel, and A. S. Desyatnikov, Flat band states: Disorder and nonlinearity, Phys. Rev. B  88, 224203 (2013).
  * [11] D. Leykam, A. Andreanov, and S. Flach, Artificial flat band systems: from lattice models to experiments, Adv. Phys. X  3, 1473052 (2018).
  * [12] R. A. V. Poblete, Photonic flat band dynamics, Adv. Phys. X  6, 1878057 (2021).
  * [13] G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, Wannier functions analysis of the nonlinear Schroödinger equation with a periodic potential, Phys. Rev. E 66, 046608 (2002).
  * [14] Y. Zhang and C. Zhang, Bose-Einstein condensates in spin-orbit-coupled optical lattices: Flat bands and superfluidity, Phys. Rev. A  87, 023611 (2013).
  * [15] Y. Zhang, Y. Xu, and T. Busch, Gap solitons in spin-orbit-coupled Bose-Einstein condensates in optical lattices, Phys. Rev. A  91, 043629 (2015).
  * [16] Y. V. Kartashov, V. V. Konotop, D. A. Zezyulin, and L. Torner, Bloch Oscillations in Optical and Zeeman Lattices in the Presence of Spin-Orbit Coupling, Phys. Rev. Lett. 117, 215301 (2016).
  * [17] Y. V. Kartashov, V. V. Konotop, D. A. Zezyulin, and L. Torner, Dynamic localization in optical and Zeeman lattices in the presence of spin-orbit coupling, Phys. Rev. A 94, 063606 (2016).
  * [18] Z. Yan and S. Wan, Topological phases, topological flat bands, and topological excitations in a one-dimensional dimerized lattice with spin-orbit coupling, EPL 107, 47007 (2014).
  * [19] F. K. Abdullaev and M. Salerno, Flat bands and dynamical localization of spin-orbit-coupled Bose-Einstein condensates, Phys. Rev. A  98, 053606 (2018).
  * [20] G. Gligorić, A. Maluckov, Lj. Hadžievski, S. Flach, and B. A. Malomed, Nonlinear localized flatband modes with spin-orbit coupling, Phys. Rev. B  94, 144302 (2016).
  * [21] Y.-J. Lin, K. Jiménez-Garcí, and I. B. Spielman, Spin–orbit-coupled Bose–Einstein condensates, Nature 471, 83 (2011).
  * [22] C.-H. Li, C. Qu, R. J. Niffenegger, Su-Ju Wang, M. He, D. B. Blasing, A. J. Olson, C. H. Greene, Y. Lyanda-Geller, Qi Zhou, C. Zhang, and Y. P. Chen, Spin current generation and relaxation in a quenched spin-orbit-coupled Bose-Einstein condensate,  Nat. Comm. 10, 375 (2019).
  * [23] A. Putra, F. Salces-Cárcoba, Y. Yue, S. Sugawa, and I. B. Spielman, Spatial Coherence of Spin-Orbit-Coupled Bose Gases, Phys. Rev. Lett 124, 053605 (2020).
  * [24] C. Hamner, Y. Zhang, M. A. Khamehchi, M. J. Davis, and P. Engels, Spin-Orbit-Coupled Bose-Einstein Condensates in a One-Dimensional Optical Lattice, Phys. Rev. Lett.  114, 070401 (2015).
  * [25] Tunable Spin-Orbit Coupling via Strong Driving in Ultracold-Atom Systems, K. Jiménez-García, and L. J. LeBlanc, R. A. Williams, M. C. Beeler, C. Qu and M. Gong, and C. Zhang, and I. B. Spielman, Phys. Rev. Lett. 114, 125301 (2015).
  * [26] Y. Zhang, G. Chen, and C. Zhang, Tunable spin–orbit coupling and quantum phase transition in a trapped Bose-Einstein condensate, Sci. Rep. 3, 1937 (2013).
  * [27] J. T. Cole and Z. H. Musslimani, Band gaps and lattice solitons for the higher-order nonlinear Schrödinger equation with a periodic potential, Phys. Rev. A 90, 013815 (2014).
  * [28] W. Kohn, Analytic properties of Bloch waves and Wannier functions, Phys. Rev. 115, 809 (1959).
  * [29] J. Javanainen, Oscillatory exchange of atoms between traps containing Bose condensates, Phys. Rev. Lett. 57, 3164 (1986).
  * [30] A. J. Leggett, Bose-Einstein condensation in the alkali gases: Some fundamental concepts, Rev. Mod. Phys. 73, 307(2001).
  * [31] R. Gati and M. K. Oberthaler, A bosonic Josephson junction, J. Phys. B: At. Mol. Opt. Phys. 40, R61 (2007).
  * [32] H. C. Prates, D. A. Zezyulin, and V. V. Konotop, Bose-Einstein condensates in quasiperiodic lattices: Bosonic Josephson junction, self-trapping, and multimode dynamics, Phys. Rev. Research  4, 033219 (2022).
  * [33] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Quantum Coherent Atomic Tunneling between Two Trapped Bose-Einstein Condensates, Phys. Rev. Lett.  79, 4950 (1997).
  * [34] K. Eckert, J. Mompart, R. Corbalán, M. Lewenstein, and G. Birkl, Three level atom optics in dipole traps and waveguides, Opt. Comm. 264, 264 (2006).
  * [35] B. Liu, L. Fu, S. Yang, and J. Liu, Josephson oscillation and transition to self-trapping for Bose-Einstein condensates in a triple-well trap, Phys. Rev. A  75, 033601 (2007).