# Chaos onset in large rings of Bose-Einstein condensates

Damian Wozniak  Institut für Physik, Universität Greifswald, 17487 Greifswald,
Germany  Department of Physics, Royal Holloway, University of London, Egham,
Surrey TW20 0EX, United Kingdom  Johann Kroha  Fachbereich Physik, Universität
Bonn, D-53115 Bonn, Germany  Anna Posazhennikova  Institut für Physik,
Universität Greifswald, 17487 Greifswald, Germany

###### Abstract

We consider large rings of weakly-coupled Bose-Einstein condensates, analyzing
their transition to chaotic dynamics and loss of coherence. Initially, a ring
is considered to be in an eigenstate, i.e. in a commensurate configuration
with equal site fillings and equal phase differences between neighboring
sites. Such a ring should exhibit a circulating current whose value will
depend on the initial, non-zero phase difference. The appearance of such
currents is a signature of an established coherence along the ring. If phase
difference falls between $\pi/2$ and $3\pi/2$ and interparticle interaction in
condensates exceeds a critical interaction value $u_{c}$, the coherence is
supposed to be quickly destroyed because the system enters a chaotic regime
due to inherent instabilities. This is, however, only a part of the story. It
turns out that chaotic dynamics and resulting averaging of circular current to
zero is generally offset by a critical time-scale $t_{c}$, which is almost two
orders of magnitude larger than the one expected from the linear stability
analysis. We study the critical time-scale in detail in a broad parameter
range.

##  I Introduction

Ring-coupled Bose-Einstein condensates (BECs) with an initially finite phase
difference between neighboring sites constitute a particularly interesting
system. They allow for circulating currents which results in a controlled
formation of topological defects such as vortices. Possible applications of
such systems range from interferometry Kasevich1997 ; Anderson1998  to quantum
computation, atomtronics and SQUIDS Hallwood2010 ; Amico2014 ; Amico2015 ;
Arwas2016 ; Ryu2013 .

A while ago, a ring of three condensates was studied in the whole range of
initial phase differences between neighboring sites Tsubota2000 ; Tsubota2002
with the goal to find the probability of vortex generation via the Kibble-
Zurek mechanism in nonuniform, domain-structured superfluids. Experimentally
the idea was tested by three 87Rb condensates merging, which indeed led to the
formation of vortices, whose number strongly depended on the merging velocity
Scherer2007 .

Although in a ring of three coupled condensates, circular current can be
nonzero when all three phase differences differ from each other, the maximum
value of the circular current is reached only for the commensurate case, i.e.
when all the three phase differences are the same Tsubota2000 . This is
because only a commensurate case corresponds to an eigenstate of the system.
Importantly, this circular current depends on the system parameters, in
particular on the interaction between condensed particles. Linear stability
analysis provides a critical interaction value $u_{c}$, above which some of
the eigenmodes become unstable and chaotic dynamics sets in for $u>u_{c}$. The
time-averaged circular current gradually tends to zero as the interaction
increases apart from the remaining sharp peaks associated with the eigenmodes
Tsubota2000 . It is therefore not clear why nonzero circular currents are
still present for $u\gg u_{c}$ in the numerical results of Ref. Tsubota2000 .

Further studies on condensate rings (with number of sites $N_{s}\geq 3$)
investigated different theoretical aspects including dynamical and
thermodynamical stability Paraoanu2003 ; dePassos2009 , chaos and ergodicity
Arwas2017 ; Arwas2019 , symmetry analysis and effects of quantum many-body
dynamics Nemoto2000  and quantum quenches Zurek2011 . However, in those works,
the values of the circular currents were not investigated in detail in the
chaotic regime, and time scales associated with the currents were not
discussed.

From an experimental point of view, it is rather challenging to realize
circulating currents in large rings and for large winding numbers due to their
quick decay to flows with lesser winding numbers (see for instance experiments
on 87Rb annular condensates in Ref. Hadzibabic2012 ). Recently, a substantial
progress in the creation of stable superflows has been achieved in six- and
seven-site rings of polaritonic condensates confined in microcavities
Cookson2021 . Particularly interesting is that persistent circular currents
with large winding numbers have been observed for nominally unstable initial
configurations for $N_{s}=7$ and winding numbers $k=2$ and $k=3$.Cookson2021

Motivated by these developments and open questions, we investigate the effect
of chaos on circular currents in detail, i.e. dependence on interaction,
initial conditions and system size. We show that, although chaos obstructs
circular flow, it does not set in immediately even if the system is tuned to
an unstable eigenmode. We identify a critical time scale associated with this
type of dynamics and demonstrate how the time scale depends on the systems
various parameters. We show that the time scale is much larger than the one
expected from the linear stability analysis, which explains numerical results
of Ref. Tsubota2000 , and can provide an insight into experimental results of
Ref. Cookson2021 .

The paper is organized as follows: In section II we formulate the model,
derive equations of motion and the expression for circular current. In section
III we analyze the circular current in the non-interacting system for a
special case of symmetric initial conditions and show the current becomes a
simple sine wave in the limit of large $N_{s}$. In section IV after a brief
discussion of unstable modes and characteristic interaction, we analyze time-
dependent circular current at unstable modes and identify a particular time
scale $t_{c}$ associated with system transition to a chaotic regime. We show
that the slide to the chaotic regime occurs exponentially and study the time-
averaged circular current showing how the current at unstable modes gradually
disappear upon increasing interaction. We conclude in section V.

##  II Model and equations of motion

We consider a system of $N$ condensed bosons trapped in a one-dimensional
periodic potential consisting of $N_{s}$ wells with periodic boundary
conditions. The average filling factor in the system

$\rho=\frac{N}{N_{s}}$ (1)

is assumed to be macroscopic so that the semiclassical Gross-Pitaevskii
approximation is applicable for the system description (for example, in
experiments on long arrays of 87Rb condensates $\rho\approx 1000$ according to
Ref. Cataliotti2001 ). The local condensates are considered to be weakly
linked in order for Josephson current to be induced between the sites
Smerzi1997 . Note that in all calculations we keep the number of particles per
site, $\rho$, constant when changing the system size $N_{s}$. The general
Gross-Pitaevskii equation reads

$\mathrm{i}\hbar\frac{\partial}{\partial t}\Psi({\bf
r},t)=\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{ext}({\bf r})+g|\Psi({\bf
r},t)|^{2}\right)\Psi({\bf r},t),$ (2)

where $\Psi({\bf r},t)$ is the mean-field averaged bosonic field operator
$\langle\hat{\Psi}({\bf r},t)\rangle$, $V_{ext}({\bf r})$ is the multi-well
external potential, $g$ is a repulsive contact interaction constant. By
expanding the semi-classical wave function in terms of a set of localized
basis functions (Wannier functions $\Phi_{i}({\bf r}$),

$\Psi({\bf r},t)=\sum_{i=1}^{N_{s}}\phi_{i}({\bf r})\psi_{i}(t),$ (3)

and integrating out the spatial degrees of freedom, we obtain the standard
discrete nonlinear Schrödinger equations (DNLSEs) for $\psi_{i}$-s
Trombettoni2001

$\displaystyle\mathrm{i}\hbar\frac{\partial}{\partial
t}\psi_{i}(t)=[E_{i}+U_{i}|\psi_{i}(t)|^{2}]\psi_{i}(t)-K_{i,i-1}\psi_{i-1}(t)$
$\displaystyle-K_{i,i+1}\psi_{i+1}(t),\quad i=1,2,...,N_{s}.$ (4)

The periodic boundary conditions imply ($N_{s}+1\rightarrow 1$). DNLSEs
adequately capture the dynamics of multiple coupled condensates, as was
verified in the experimental work Cataliotti2001 .

Figure 1:  Schematic setup of the ring system for five wells ($N_{s}=5$). Blue
circles represent condensates wave-functions, $K$ is the Josephson coupling
constant (5), $n_{i}$-s are condensate populations according to (6).

The model parameters in Eq. (4), the zero-point energies $E_{i}$, on-site
interaction $U_{i}$ and Josephson couplings between neighbouring wells
$K_{i,i\pm 1}$ are given by Smerzi1997

$\displaystyle E_{i}$ $\displaystyle=$ $\displaystyle\int d{\bf
r}\left[\frac{\hbar^{2}}{2m}|\nabla\phi_{i}({\bf r})|^{2}+|\phi_{i}({\bf
r})|^{2}V_{ext}({\bf r})\right],$ $\displaystyle U_{i}$ $\displaystyle=$
$\displaystyle g\int d{\bf r}|\phi_{i}({\bf r})|^{4},$ $\displaystyle
K_{i,i\pm 1}$ $\displaystyle=$ $\displaystyle-\int d{\bf
r}\left[\frac{\hbar^{2}}{2m}\nabla\phi_{i}({\bf r})\nabla\phi_{i\pm 1}({\bf
r})\right.$ (5) $\displaystyle+$ $\displaystyle\left.\phi_{i}({\bf
r})V_{ext}({\bf r})\phi_{i\pm 1}({\bf r})\right].$

We make use of the ansatz

$\psi_{i}(t)=\sqrt{n_{i}}(t)e^{\mathrm{i}\theta_{i}(t)}$ (6)

where we introduced the site populations $n_{i}$ normalized by the filling
factor $\rho$

$n_{i}(t)=\frac{N_{i}(t)}{\rho},$ (7)

where $N_{i}(t)$ is the number of particles in the condensate on site $i$ at
time $t$. Thus, an initially homogeneous distribution of atoms means
$n_{i}(0)=1$ for all $i$. With that Eqs. (4) can be rewritten as a set of
differential equations for $n_{i}$ and phase differences
$\theta_{i,i+1}=\theta_{i+1}-\theta_{i}$ as follows

$\displaystyle\dot{n}_{i}$ $\displaystyle=$
$\displaystyle-2\sqrt{n_{i}}(\sqrt{n_{i+1}}\sin\theta_{i,i+1}-\sqrt{n_{i-1}}\sin\theta_{i-1,i}),$
$\displaystyle\dot{\theta}_{i,i+1}$ $\displaystyle=$ $\displaystyle
u(n_{i}-n_{i+1})+\left(\sqrt{\frac{n_{i}}{n_{i+1}}}-\sqrt{\frac{n_{i+1}}{n_{i}}}\right)\cos\theta_{i,i+1}$
(8) $\displaystyle+$
$\displaystyle\sqrt{\frac{n_{i-1}}{n_{i}}}\cos\theta_{i-1,i}-\sqrt{\frac{n_{i+2}}{n_{i+1}}}\cos\theta_{i+1,i+2}.$

In deriving these equations we assumed the following simplifications:
$E_{i}=0$, $U_{i}\equiv U$ and $K_{i,i+1}\equiv K$ for all values of $i$. We
also expressed the time argument $t$ in units of $\hbar/K$, and introduced the
dimensionless interaction parameter

$u=\frac{U\rho}{K}.$ (9)

These equations conserve $\sum_{i}n_{i}=N_{s}$ and the total energy

$E=\frac{\rho
K}{\hbar}\left(\frac{u}{2}\sum_{i}n_{i}^{2}-2\sum_{i}\sqrt{n_{i}n_{i+1}}\cos\theta_{i,i+1}\right).$
(10)

We define the circular current as the average current in a clockwise loop
around the ring of condensates

$I=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}I_{i,i+1},$ (11)

where $I_{i,i+1}$ is the particle current from site $i$ to $i+1$ defined as

$I_{i,i+1}=2\operatorname{Im}(\psi_{i}^{*}(t)\psi_{i+1}(t))=2\sqrt{n_{i}n_{i+1}}\sin(\theta_{i,i+1}).$
(12)

Note that the current defined in this way has units $\rho K/\hbar$. In the
next section, we analyze the current analytically for $u=0$ and a special case
of initial conditions.

##  III Circular current in the non-interacting case

The non-interacting case with $U_{i}=0$ and coupling constants
$K_{i,i+1}\equiv K$ can be solved exactly for any number of sites $N_{s}$, as
shown in Appendix A. However, one does not need these solutions in order to
calculate circular current, as it can be straightforwardly derived from the
current conservation condition (one can show, for example, that
$\partial_{t}I=0$ with the help of Eqs.(4)). We thus get for the circular
current

$I=I_{0}=\frac{2}{N_{s}}\operatorname{Im}\sum_{i=1}^{N_{s}}\psi_{i}^{*}(0)\psi_{i+1}(0).$
(13)

We see that the current is constant and the value of this constant depends on
initial $n_{i}(0)$-s and $\theta_{i,i+1}(0)$-s. For the homogeneous condensate
distribution $n_{i}(0)=1$, the current only depends on the initial phase
differences.

Figure 2:  Average current of the noninteracting system $I_{0}$ (16) divided
by its maximum value (=2) versus initial phase difference $\theta_{0}$ for
initial conditions (15), plotted for condensate rings of different sizes (
$N_{s}=3;4;5;8;15$ and $20$).

We chose the following initial conditions for the phase differences

$\theta_{i,i+1}=\left\\{\begin{array}[]{ccl}\theta_{0}&\mbox{for}&i=1,...,N_{s}-1,\\\
-(N_{s}-1)\theta_{0}&\mbox{for}&i=N_{s},\end{array}\right.$ (14)

so that

$\psi_{j}(0)=e^{\mathrm{i}(j-1)\theta_{0}},\quad j=1,2,..,N_{s},$ (15)

We now get for the current in Eq. (13) the simple expression

$I_{0}(\theta_{0})=\frac{2}{N_{s}}\left\\{(N_{s}-1)\sin\theta_{0}-\sin[(N_{s}-1)\theta_{0}]\right\\}.$
(16)

Since the current is an odd function of the initial phase difference
$I_{av}^{0}(\theta_{0})=-I_{av}^{0}(-\theta_{0})$ and is $2\pi$-periodic, it
is sufficient to consider $\theta_{0}\in[0,\pi]$. In Fig. 2 we plotted $I_{0}$
normalized by its maximum value $(I_{0})_{max}=2$ versus $\theta_{0}$ for
various ring sizes. We see that the current has local maxima at discrete
values of $\theta_{0}$ given by

$\theta_{0}^{k}=\frac{2\pi k}{N_{s}},$ (17)

where the integer $k=1,...,N_{s}$ has the meaning of the winding number. This
result is not surprising and could be alternatively found from the extrema of
the total energy of the noninteracting system (10)

$E_{0}=-2\frac{\rho
K}{\hbar}\left[(N_{s}-1)\cos\theta_{0}+\cos[(N_{s}-1)\theta_{0}]\right].$ (18)

This is because $\theta_{0}^{k}$-s play the role of quantized components of
the effective quasi-momentum in Fourier space, and the group velocity and as a
consequence the circular current is proportional to $\partial
E/\partial\theta_{0}^{k}$. Naturally, the energy values at (17) coincide (up
to our normalization factor) with the eigenvalues $\lambda_{k}$ of
noninteracting Hamiltonian calculated in Appendix A and read

$E_{k}=-2\frac{NK}{\hbar}\cos\frac{2\pi k}{N_{s}}.$ (19)

Note that the discrete modes (17) correspond to the commensurate case with
homogeneous initial conditions. Values of $\theta_{0}$ in between the discrete
modes belong to a special case of incommensurate configuration when all but
one (the phase difference between first and the last sites) initial phase
differences are the same as shown in Eq.(14).

At the discrete modes (17) the current (16) is simply a sine-function

$I_{0}(\theta_{0}^{k})=2\sin\left(\theta_{0}^{k}\right).$ (20)

##  IV Circular current for ring-coupled, interacting condensates

###  IV.1 Chaotic behavior and the critical interaction

Once the interaction is switched on, the system dynamics become more
complicated, and analytical solutions are no longer possible. It is important
that, contrary to the non-interacting case, the dynamics become chaotic for a
specific parameter range. To identify the parameter range, we start from the
linear stability analysis of the coupled real Eqs. (8). The stability is
decided by the eigenvalues of the corresponding Jacobian matrix (see Appendix
B for details). These $2N_{s}$ eigenvalues can be derived analytically due to
the Jacobian matrix’s special, block-wise circulant structure. The eigenvalues
are

$\begin{split}\lambda_{j}(k)&=2\mathrm{i}\left\\{\sin\theta_{0}^{k}\sin{\left(\frac{2\pi
j}{N_{s}}\right)}\pm\right.\\\ &\left.\sin{\left(\frac{\pi
j}{N_{s}}\right)}\sqrt{2\cos\theta_{0}^{k}\left[2\cos\theta_{0}^{k}\sin^{2}\left(\frac{\pi
j}{N_{s}}\right)+u\right]}\right\\},\end{split}$ (21)

where $j=1,2,..,N_{s}$. The eigenvalues depend on the eigenmodes
$\theta_{0}^{k}$, and the dimensionless interaction $u$ from Eq.(9). Note that
$\epsilon_{j}(k)=-\mathrm{i}\lambda_{j}(k)$ correspond to the discrete
Bogoliubov excitation spectrum of the system, discussed previously in the
context of condensate arrays Smerzi2002 , and the Bogoliubov-de Gennes
description of a circular array of Bose-Einstein condensates Paraoanu2003 .
The eigenvalues $\lambda_{j}$ are zero or purely imaginary, and both Eqs. (8)
are stable unless the expression under the square root turns negative. Since
we consider only repulsive interactions, $u>0$, the condition for at least one
eigenvalue to acquire a real part is

$\displaystyle\cos\theta_{0}^{k}<0,$ (22) $\displaystyle
2\cos\theta_{0}^{k}\sin^{2}\left(\frac{\pi j}{N_{s}}\right)+u>0.$ (23)

The analysis of these inequalities provides the expression for the upper bound
of interaction below which the stationary points of the nonlinear equations
are still stable. The first eigenvalue acquiring a real part turns the system
unstable, which first occurs for $j=1$, giving us the critical interaction

$u_{c}=-2\cos\theta_{0}^{k}\,\sin^{2}\left(\frac{\pi}{N_{s}}\right).$ (24)

Note that $u_{c}>0$ because of the condition (22). We see that $u_{c}$ depends
on the initial conditions through the modes $\theta_{0}^{k}$, with the mode
selection criterion Eq. (22). Accordingly, the range of unstable modes is
defined by the interval

$\frac{\pi}{2}<\theta_{0}^{k}<\frac{3\pi}{2}\ \ \text{or}\ \ N_{s}<4k<3N_{s}.$
(25)

We will refer to such modes as "unstable discrete modes", keeping in mind that
they become unstable for $u>u_{c}$. Conditions similar to Eqs. (24), (25) were
derived also in Ref. [Paraoanu2003, ] from the Bogoliubov-de Gennes equations.
Conversely, modes in the complementary range,

$\theta_{0}^{k}\in[-\pi/2,\pi/2]$

are stable and will be called "stable discrete modes".

It follows from Eqs.(24) and (25) that the maximum possible value of the
critical interaction $u_{c}=1$, and it is reached for the $\pi-$mode of a
four-site ring. Since $u_{c}$ is just proportional to the corresponding cosine
of the corresponding discrete mode, one can plot a universal, mode-independent
graph for normalized characteristic interaction $u_{c}/(-\cos\theta_{0}^{k})$,
which we display in Fig. 3. One can see that $u_{c}$ tends to zero relatively
fast with the increasing number of sites $N_{s}$. The inset of the graph shows
the stability diagram, where modes denoted by open circles represent stable
solutions independent of the value of $u$ (given that only non-negative
interactions are considered), whereas solid circles represent modes that
become unstable for $u>u_{c}$.

Figure 3: Dependence of characteristic interaction $u_{c}$ (divided by
$\cos(\theta_{0}^{k})$) on the number of sites $N_{s}$, ranging between 3 and
40. Inset: stability diagram showing discrete modes, i.e. values of initial
phase difference $\theta_{0}$ from (17) versus $N_{s}$. The modes denoted by
open circles represent modes that are stable independently of the interaction
$u$-value. Discrete modes denoted by solid circles showcase the modes which
are stable for $0<u\leq u_{c}$ but become unstable for $u>u_{c}$.

The linear stability analysis introduces the instability exponent
$\alpha_{0}$, given by the real part of the Jacobian eigenvalue $\lambda_{1}$,

$\alpha_{0}\equiv
Re(\lambda_{1})=2\sin{\left(\frac{\pi}{N_{s}}\right)}\sqrt{2|\cos\theta_{0}^{k}|\left[u-u_{c}\right]}.$
(26)

This $\alpha_{0}$ is the rate by which the circular current $I(t)$ of a given
mode $\theta_{0}^{k}$ is expected to deviate exponentially in time from its
stationary value $I_{0}(\theta_{0}^{k})$. To be more specific, it should be
$2\alpha_{0}$, since the current is proportional to the product of two
condensate wave functions. The equation (26) hence establishes the deviation
time scale $\sim 1/\alpha_{0}$ which would diverge at the stability
boundaries, i.e., for $\theta_{0}\to\pi/2$ or $3\pi/2$, or $u\to u_{c}$. We
will show in the following section that, interestingly, chaotic behavior sets
in abruptly at a much larger, critical time $t_{c}\gg 1/\alpha_{0}$, due to
the nonlinear character of the system, not captured by the linear stability
analysis.

###  IV.2 Temporal onset of chaotic behavior

In this section, we fully characterize the decay of the unstable, discrete
modes by studying the circulating current numerically beyond linear stability
analysis. Initially, the circular current has a nonzero value equal to its
noninteracting value $I_{0}$, derived in section III.

In Figs. 4 and 5 we show examples of the time-dependent circular current of
the interacting system defined in Eqs. (11) and (12) for a ring of $N_{s}=5$
and $N_{s}=20$ condensates, respectively. For $N_{2}=5$, according to Eq. (25)
there are only two unstable modes $k=2$ and $k=3$. These modes are symmetric
with respect to inflection at $\theta_{0}=\pi$ (see the inset in Fig. 3) and,
therefore, the critical value of interaction $u_{c}$ is the same for both of
the modes, $u_{c}\approx 0.56$. In Fig. 4 we consider three values of
interaction, which are all greater than $u_{c}$: $u=0.6$, $u=1.2$ and $u=2.4$,
and we, therefore, expect the current to behave chaotically in all three
cases, which is indeed observed. However, chaos sets in at a certain time
$t_{c}>0$ which strongly depends on the interaction $u>u_{c}$. One can see
that $t_{c}$ decreases as the interaction increases. For example, for $u$
close to $u_{c}$ (upper panel), $t_{c}$ is about 112, for $u=2.4$ it is about
15.

Figure 4:  Time-dependent circular current (11) for homogeneous initial conditions: $N_{s}=5$, $\theta_{0}=\theta_{0}^{k}$, $k=2$ and $k=3$. The current is calculated numerically for three different values of dimensionless interaction $u$ shown in the panels. Time $t$ is in the units of $1/K$.  Table 1: Values of characteristic interaction $u_{c}$ for different modes $\theta_{0}^{k}$ for the case of $N_{s}=20$.  $k$ | 6; 14 | 7; 13 | 8; 12 | 9; 11 | 10
---|---|---|---|---|---
$u_{c}$ | 0.015 | 0.029 | 0.040 | 0.047 | 0.049

To find out whether $t_{c}$ has a mode dependence, we consider in Fig. 5 a
larger ring of $N_{s}=20$. It has 9 discrete, unstable modes whose critica
linteraction strengths $u_{c}$ are listed in table 1. The circular currents
are shown for different values of $u$, all greater than the respective
$u_{c}$. We observe that $t_{c}$ does not only depend on $u$, but also on the
mode. For example, for $u=3$ values of $t_{c}$ of the "outer" modes ($k=6$ and
$k=14$) are clearly greater than $t_{c}$ for other modes.

Figure 5: Time-dependent circular current $I(\theta_{0}^{k})$ of the unstable
modes $k=6,7,8,...,14$ in the case of $N_{s}=20$. The current is calculated
numerically for three different values of $u$ shown in the panels. Time $t$ is
in the units of $1/K$.

In order to quantify the onset of chaos, we show in Fig. 6 the time evolution
of the deviation of the circular current $I(t)$ from its noninteracting value
$I_{0}(\theta_{0}^{k})\equiv I_{0}$ on a logarithmic scale. As expected from
the linear stability analysis, this deviation is initially exponential in time
and, thus, initially the current deviation remains exponentially small, too
small to be resolved in the linear plots of Figs. 4 and 5.

The exponential behaviour of the deviation is governed by two parameters
$\alpha$ and $b$

$\Delta I(t)\equiv|I_{0}-I(t)|/2=\ e^{\alpha t+b}\,,$ (27)

which can be determined from linear fits to the logarithmic plot within the
exponential time range (see dashed lines in Fig. 6). The comprehensive
analysis for a wide range of system paramters $u\in[0.5,10]$ and initial
conditions $\theta_{0}^{k}\in[0,\pi]$ shows that the coefficient $b$ has no
systematic dependence on system parameters, with an average value of $\langle
b\rangle\approx-72.34$ and a standard deviation of $\Delta b\approx 1.51$.
Furthermore, from Fig. 6 and similar plots for the wide parameter range
mentioned above (not shown), we find that chaotic behavior, i.e. deviation
from the exponential time evolution, sets in abruptly when the normalized
current deviation $\Delta I(t)$ reaches a universal value of about $\ln[\Delta
I(t_{c})]\equiv\lambda=-2$. This defines, via Eq. (27), the critical time for
onset of chaos as,

$t_{c}\approx\frac{\lambda-b}{\alpha}\approx\frac{70}{\alpha}.$ (28)

It shows that both, the initial exponential deviation and the onset of chaotic
evolution, are controlled by the instability exponent $\alpha$ alone, where,
however, the critical time $t_{c}$ is a factor of 70 larger than the
exponential time scale $1/\alpha$. We note in passing that similar long-time
coherent evolution with abrupt chaos onset was found in single Bose-Josephson
junctions Trujillo2009 . Thus, the observed dependence of $t_{c}$ (see Figs.
4, 5) on system parameters enters through the dependence of $\alpha$, which we
explore next.

Figure 6:  Logarithmic plot of the deviation of the current from its initial
value $I_{0}$ versus time for $N_{s}=20$, $k=6$ and three different $u$-s.
Time $t$ is in the units of $1/K$. Dashed lines are linear fits of linear
parts of the graphs with intercept $b$ and slope $\alpha$: $\alpha\approx
0.5$, $b\approx-71.62$ (for $u=0.25$); $\alpha\approx 0.98$, $b\approx-70.55$
($u=0.5$); $\alpha\approx 1.98$, $b\approx-71.62$ ($u=1$).

In Fig. 7 we show results of the numerical evaluation of $\alpha$ depending on
initial mode $\theta_{0}^{k}$ and interaction $u$. The plot contains data for
various $N_{s}$, $k$ and interaction $u$. Short vertical lines of data at
$\theta_{0}^{k}=2\pi/3\approx 0.666\,\pi$ correspond to rings with
$N_{s}=3,6,9...$, because they all have such a mode. However, the dependence
of a particular mode on $N_{s}$ is very weak, so that we neglect it in the
discussion and analysis. As expected, $\alpha$ tends to zero for
$\theta_{0}^{k}\rightarrow\pi/2$, since this value of $\theta_{0}^{k}$ marks
the stability boundary (see the stability diagram in Fig. 3). Away from the
stability boundary, i.e. close to $\pi$-modes, $\alpha$ has a weaker
dependence on $\theta_{0}^{k}$ and is significantly dependent on $u$. This
dependence is approximately linear, as shown in the inset: $\alpha(\pi)\approx
2u$ on the interval $u\in(u_{c},5)$, with $u_{c}\ll 1$ in this case. For small
$u$, $\alpha$ tends to zero, i.e., the modes become stable.

Figure 7:  The exponent $\alpha$ versus initial $\theta_{0}^{k}$ for different
values of $u$. The inset shows the dependence of the maximum value of $\alpha$
equal to $\alpha(\theta_{0}^{k}=\pi)$ versus $u$.

At this point, it is interesting to compare the numerically estimated $\alpha$
with the analytical value from (26). It turns out that the agreement holds
only for small rings ($N_{s}=3$, $N_{s}=4$), whereas for larger rings the
agreement with this relation quickly degrades and holds only in the close
vicinity of $u_{c}$. This defines the applicability of linear stability
analysis in our systems - essentially only in the vicinity of $u_{c}$.

Coming back to the chaos onset time $t_{c}$, as expected, it becomes
arbitrarily large in two cases: close to stability boundaries in terms of
initial phase difference, i.e. close to $\theta_{0}=\pi/2$ and
$\theta_{0}=3\pi/2$, and for interactions close to $u_{c}$. This can be seen
in Fig. 8, where we present numerically calculated $t_{c}$ as a function of
$u$ for a broad range of initial conditions $\theta_{0}^{k}$, in particular,
for all unstable modes $\theta_{0}^{k}\in(\pi/2,\pi]$ of all rings between
$N_{s}=3$ and $N_{s}=20$. To not overload the graph with information, we
color-coded curves according to their initial conditions split into intervals
of $\theta_{0}$: "black" corresponds to $\theta_{0}^{k}\in[0.8\pi,\pi]$,
closest to the anti-phase mode and farthest from the stable modes. One can see
that most of them (in particular those with small $u_{c}$) bunch around the
curve marked by the black arrows. The marked curve is easily fit with a two-
parameter fitting function resulting in $t_{c}\approx 28/(u-u_{c})^{0.69}$.
"Red" curves represent the next interval $\theta_{0}^{k}\in[0.6\pi,0.78\pi)$.
They start to deviate from the black curves, especially for large $u$. This
can be better seen in Fig. 9, where we plot $t_{c}^{-1}$ versus $u-u_{c}$.
Last, the blue curves span the interval of $\theta_{0}^{k}\in[0,0.59\pi]$. One
curve is calculated for $N_{s}=47$ and $k=12$, being an example of a curve
relatively close to $\pi/2$-mode with $\theta_{0}^{k}\approx 0.511$ (this is
the uppermost curve in the graph). The fitting of this curve results in
$t_{c}\approx 62/(u-u_{c})^{0.53}$, leading to a sizable $t_{c}$ even for
large values of $u$.

Figure 8:  Dependence of the onset of chaotic behavior $t_{c}$ on interaction
$u$ measured from $u_{c}$ for all unstable discrete modes in the interval
$(\pi/2,\pi]$ for $3\leq N_{s}\leq 20$. The dashed curve is for $N_{s}=47$ and
$k=12$ (this amounts to $\theta_{0}\approx 0.511\pi$). Initial conditions are
specified by $\theta_{0}^{k}$. We split the initial conditions into color-
coded intervals.  Figure 9:  Dependence of the inverse onset of chaotic
behavior $t_{c}^{-1}$ on interaction $u-u_{c}$ for all unstable modes in the
interval $(\pi/2,\pi]$ for $3\leq N_{s}\leq 20$. The lowest curve is for
$N_{s}=47$ and $k=12$. The curves are color-coded according to initial
conditions in the same way as in Fig. 8.

Note that $t_{c}$ can be determined differently, for example, from the
deviation from its initial value by $0.0001$ of any site occupation. Finding
$t_{c}$ from site occupation would lower the $t_{c}$ of entire Fig.8 about an
average of $10\%$, which does not change any conclusions or analysis.

To get an idea of how large $t_{c}$ can be in terms of experimental values, we
translate the value $t_{c}=100$ to ms from experimental data on arrays of
condensates Cataliotti2001 . For example, from experiments on long arrays of
cold atoms ($N_{s}=200$), we take the value of Josephson coupling $K\sim
0.07E_{R}$, where $E_{R}$ is the recoil energy of 87Rb atom of mass $m$
absorbing one of the lattice photons

$E_{R}=\frac{h^{2}}{2m\lambda^{2}}\sim 1.5\cdot 10^{-11}eV$ (29)

with $\lambda=795\cdot 10^{-9}$m, Cataliotti2001 . This gives us the value of
$K\sim 1.05\cdot 10^{-12}$eV. From this analysis it follows that the
dimensionless $t_{c}=100$ will correspond to $62.86$ ms. Such values of
$t_{c}$ could be achievable for very weak interactions as follows from Fig. 8.

Note, that when initial conditions fall in between the eigenmodes, one expects
time-dependent current, as shown in Fig. 10 for $\theta_{0}=0.75\pi$ and
$N_{s}=20$. We see that with increasing $u$ and therefore non-linearity,
current oscillations become larger, more chaotic and eventually average to
zero over time for large $u$, as we will see in the next section. The notion
of $t_{c}$ does not make sense in this case, since there appears no additional
energy scale associated with real parts of eigenvalues.

Figure 10: Current oscillations with time for initial conditions
$\theta_{0}=0.75\pi$ and $N_{s}=20$, i.e. for a non-eigenmode. The results are
shown for three different values of dimensionless interaction $u$ shown in the
inset.

###  IV.3 Time-averaged circular current and coherence

We now evaluate the time-averaged circular current numerically

$\langle I(\theta_{0})\rangle=\lim_{T\rightarrow\infty}\frac{1}{T}\int
dt\,I(\theta_{0})$ (30)

for initial conditions (15) with $\theta_{0}$ ranging from $0$ to $2\pi$
covering stable and unstable regions of the stability diagram in full. We
chose a ring with $N_{s}=20$ sites and four values of $u$ greater than the
maximum $u_{c}$ listed in Table I. Given the discussion about $t_{c}$ in the
previous section, it is clear that the resulting time-averaged current at the
discrete unstable modes will depend on the numerically available time
interval, over which we can let our program run, providing reliable results.
In our case $t\in[0,200]$. The results for the averaged current are presented
in Fig. 11. The four panels of Fig. 11 correspond to four different values of
$u$ as indicated. The dashed curve is the current for $u=0$ and is shown on
all the panels for reference. For small values of the dimensionless
interaction $u=0.1$ the deviations of the current from $u=0$ values in the
stable part of the plot are hardly visible (i.e. for $\theta_{0}\in[0,\pi/2]$
and $[3\pi/2,2\pi]$). For the unstable part, $\theta_{0}\in(\pi/2,3\pi/2)$, it
is striking that $t_{c}$ is, in fact, greater than $200$, otherwise the
averaged current at the unstable discrete modes stays unaffected by the
chaotic regime even though $u>(u_{c})_{max}$. Since there are no additional
time scales for $\theta_{0}$ in between the discrete modes, the average
current values begin to deviate from their $u=0$ values due to chaotic
dynamics (see Fig. 10).

Figure 11:  Time-averaged circular current normalized by its maximum value
(=2) versus initial phase difference $\theta_{0}$. The current is calculated
for $N_{s}=20$ ring and four different values of $u$ indicated in the panels.
The red dashed line corresponds to the $u=0$ non-interacting case and is put
as a reference, the solid lines correspond to the interacting systems. The
current is averaged over the time interval $t\in[0,200]$.

Upon increasing interaction, $t_{c}$ decreases, so that for $u=0.2$ the
$t_{c}$ becomes of the order of $180$, which is comparable with the maximum
value of $t$ in our numerics. As a result, the current values at discrete
unstable modes remain at most unaffected, whereas the in-between
$\theta_{0}$-s correspond to circular currents quickly averaging out to zero.
This effect is also visible around the stability boundaries $\theta_{0}=\pi/2$
and $\theta_{0}=3\pi/2$. Although current values at stable discrete modes
remain constant, the values in between begin to feel the effect of the
increasing interaction. These tendencies become more pronounced as we increase
the interaction further. As $t_{c}$ rapidly decreases ($t_{c}\sim 75$ for
$u=0.5$ and $t_{c}\sim 35$ for $u=1$), so does the time-averaged current in
the unstable region. For $u=1$, this current is practically zero everywhere in
the unstable region. The stable region is less affected by interaction, and
the effect is most visible for values of $\theta_{0}$ falling in between the
stable discrete modes, close to the stability boundaries. In the limit of
large $N_{s}$, the unstable modes fill the interval $(\pi/2,3\pi/2)$, while
the stable modes fill the intervals $[0,\pi/2]$ and $[3\pi/2,2\pi]$. Large
$N_{s}$ approaches the continuous limit where all $\theta_{0}$ are eigenmodes
(because $\theta_{0}^{k}$ merge together). The current in the limit of large
$N_{s}$ is described by $2\sin\theta_{0}$, but in the unstable region the
current averages to zero above $u_{c}$, provided averaging time period is
greater than $t_{c}$.

Figure 12:  The real part of the first-order correlation function
Re$\,g_{ij}^{(1)}$ averaged over time on the interval $t\in[0,200]$. The
correlation function is plotted as a function of distance w.r.t. to the first
site (the distance is normalized by a lattice constant $a\equiv 1$).

The circular current averaging to zero is equivalent to a loss of coherence in
the system. The coherence can also be quantified by the first-order coherence
function

$g_{ij}^{(1)}=\frac{\langle\psi_{i}^{*}\psi_{j}\rangle}{\sqrt{\langle|\psi_{i}|^{2}\rangle\langle|\psi_{i}|^{2}\rangle}}.$
(31)

It can be experimentally measured by interferometry as a function of a
distance from a fixed site. In Fig. 12 we show the real part of the coherence
function versus distance from site $i=1$. We chose two different initial
values of the phase difference: $\theta_{0}^{(3)}$ for $N_{s}=6$, i.e. the
$\pi$-mode of a 6-site ring, and a mode $\theta_{0}^{(3)}$ close to $\pi$ of a
7-site ring. Although both of these modes are unstable, experiments on
polariton condensates Cookson2021  demonstrate nonzero circulating currents
and coherence corresponding to small values of $u$, e.g. $uN_{s}=1$ in Fig.
12. It means that either the interaction was rather small in the experiment,
and $t_{c}$ relatively large, or there were other stabilizing factors in the
experimental system. For example, the experimental system is an open system,
whereas our system is closed.

In Fig. 12 we further demonstrate how the coherence is destroyed by
interaction and averages to zero for increasing $u$ (see data for $u\cdot
N_{s}=10$ and $u\cdot N_{s}=50$).

##  V Conclusions and discussion

We analyzed circular currents and their stability in rings of condensates
under specific initial conditions of equal filling and homogeneous phase
differences. We found a set of discrete eigenmodes (phase differences)
differentiated by their winding numbers. When such a mode falls into the
interval $(\pi/2,3\pi/2)$, it is stable until the interaction exceeds a
certain value $u_{c}$. This critical interaction depends on the mode and on
the ring size. When $u_{c}$ is exceeded, the system dynamics and the circular
current become chaotic, with the current quickly averaging to zero over time.
This marks the effective loss of coherence in the system. We showed that this
dephasing occurs not immediately upon entering the unstable regime but rather
after a chaos onset period $t_{c}$. This period $t_{c}$ can be arbitrary large
when the system is close to the stability boundary, i.e. when the unstable
mode under consideration is close to $\pi/2$ or $3\pi/2$, or when $u$ is close
to $u_{c}$. For modes close to the antiphase $\pi$ mode in the instability
region, i.e. for modes falling onto interval $[\pi,\pi\pm 0.2\pi]$ one can
even establish a universal behavior of $t_{c}\approx 28(u-u_{c})^{-0.69}$.

We also established that the critical time for chaos onset, that is, the
dephasing time, is proportional to but about two orders of magnitude larger
than the time scale of exponential deviations, $1/\alpha$, where $\alpha$ is
the current instability exponent. This may be relevant for technological
applications of quantum coherent dynamics of rings of interacting Bose-
Einstein condensates. The presence of the large chaos onset period $t_{c}$
explains why in the previous works about three-site condensate rings, the
circular current was non-vanishing in the chaotic regime. It may also shed
light on recent experimental observations of circulating currents in loops of
polaritonic condensates with large winding numbers, although there can be
other stabilizing factors.

In the future, it would be interesting to explore how fluctuations would
affect the present description and extension of the work to polaritonic
condensates would be of interest.

## Acknowledgments

We acknowledge S. Ray and especially M. Eschrig for many fruitful discussions
regarding the project. We are also grateful to an anonymous Referee for their
excellent comments which helped us improve the manuscript. This work was
financially supported in part (J.K.) by the Deutsche Forschungsgenmeinschaft
(DFG) via SFB/TR 185 (277625399) and the Cluster of Excellence ML4Q (EXC
$2004/1-390534769$).

## Appendix A: Exact solution of the non-interacting case.

The non interacting case when $u=0$ in Eq.(4) can be solved analytically by
rewriting these simultaneous equations in a matrix form. Albeit it cannot be
written in a succinct final solution for $\psi$ due to a large number of
unknown coefficients. However with the simplifications of $\hbar\equiv 1$,
$E_{i}\equiv E_{0}$, $U_{i}\equiv U$, $K_{i,i+1}\equiv K$ for $i=1,..,N_{s}$
and expressing the time argument $t$ in the units of $1/K$ we get

$\displaystyle\mathrm{i}\frac{\partial}{\partial
t}\vec{\psi}(t)=\hat{H}\vec{\psi}(t)$ (32)

where $\epsilon=E_{0}/K$, vector
$\vec{\psi}=(\psi_{0},\psi_{1},..,\psi_{N_{s}-1})$ contains $N_{s}$ condensate
functions which are numbered from $0$ to $N_{s}-1$ for convenience, and
$\hat{H}$ is a $N_{s}\times N_{s}$ matrix

$\hat{H}=\begin{pmatrix}\epsilon&-1&0&0&\dots&0&-1\\\
-1&\epsilon&-1&0&\dots&0&0\\\ 0&-1&\epsilon&-1&\dots&0&0\\\ &&&\ddots&&&\\\
0&0&0&0&\dots&\epsilon&-1\\\ -1&0&0&0&\dots&-1&\epsilon\end{pmatrix}.$ (33)

This matrix is circulant and real-symmetric with real eigenvalues (in our
case, coinciding with energy eigenvalues). The eigenvectors $\vec{v}$ of
circulant matrices are well-known and do not depend on circulant matrix
entries Circulant_matrix . The $k$-component of $\vec{v}$, corresponding to an
eigenvalue $\lambda_{j}$ reads

$v_{k}(\lambda_{j})=\exp{\mathrm{i}\frac{2\pi jk}{N_{s}}}.$ (34)

Here indices $k,j=0,1,...,N_{s}-1$. The eigenvalues of $\hat{H}$ are also
readily found, as they are the discrete Fourier transforms of the first row of
the matrix $\hat{H}$

$\lambda_{j}=\sum^{N_{s}-1}_{k=0}h_{0k}\exp{\mathrm{i}\frac{2\pi jk}{N_{s}}},$
(35)

where $h_{00},h_{01},...,h_{0,N_{s}-1}$ are the entries of the first row of
the circulant matrix (33). In our case the eigenvalues acquire a simple form

$\lambda_{j}=\epsilon-2\cos{\left(\frac{2\pi j}{N_{s}}\right)}.$ (36)

This gives us a general solution

$\psi_{n}(t)=\sum_{j=0}^{N_{s}-1}c_{j}v_{n}(\lambda_{j})\,e^{-\mathrm{i}\lambda_{j}t},$
(37)

where the coefficients $c_{j}$ are defined by initial values of the condensate
wave-functions $\vec{\psi}(0)$

$\vec{\psi}(0)=\hat{F}\vec{c}$ (38)

where the discrete Fourier transform matrix $\hat{F}$ contains eigenvectors
$\vec{v}(\lambda_{j})$ as columns

$F_{kj}=v_{k}(\lambda_{j})$ (39)

with $k,j=0,1,...,N_{s}-1$. $\vec{c}$ is hence determined by the inverse
Fourier transform

$c_{j}=\frac{1}{N_{s}}\sum_{k=0}^{N_{s}-1}e^{-\mathrm{i}\frac{2\pi}{N_{s}}k\cdot
j}\psi_{k}(0).$ (40)

Finally,

$\displaystyle\psi_{n}(t)=\frac{1}{N_{s}}\sum_{j,k=0}^{N_{s}-1}e^{\mathrm{i}\frac{2\pi}{N_{s}}(n-k)j}\psi_{k}(0)e^{-\mathrm{i}\lambda_{j}t}.$
(41)

For the particle number per site $N_{n}=\psi_{n}\psi_{n}^{*}$ we get

$\displaystyle N_{n}$ $\displaystyle=$
$\displaystyle\sum_{j,k=0}^{N_{s}-1}c_{n,j,k}\
e^{-\mathrm{i}t(\lambda_{j}-\lambda_{k})},$ (42) $\displaystyle c_{n,j,k}$
$\displaystyle=$
$\displaystyle\frac{1}{N_{s}^{2}}\sum_{l,m=0}^{N_{s}-1}e^{\mathrm{i}\frac{2\pi}{N_{s}}(n-k-l+m)j}\psi_{l}(0)\psi^{*}_{m}(0).$
(43)

As an example, we present solutions for $N_{s}=3$. The eigenvalues are
$\lambda_{1}=\lambda_{2}=\epsilon+1$, $\lambda_{3}=\epsilon-2$ and the three
wavefunction have the form

$\psi_{i}(t)=\frac{S}{3}e^{\mathrm{i}(2-\epsilon)t}+\left(\psi_{i}(0)-\frac{S}{3}\right)e^{-\mathrm{i}(\epsilon+1)t},\quad
S=\sum_{i=1}^{3}\psi_{i}(0).$ (44)

The occupation number per site is then

$n_{i}=n_{i}(0)+2\real\left[\frac{S}{3}\left(\psi^{*}_{i}(0)-\frac{S^{*}}{3}\right)\left(e^{3\mathrm{i}t}-1\right)\right].$
(45)

For initial conditions (15) and $\theta_{0}$ coinciding with the eigenmodes
(17), the $S=0$ giving the simple expressions

$\displaystyle\psi_{i}(t)$ $\displaystyle=$
$\displaystyle\psi_{i}(0)e^{-\mathrm{i}(\epsilon+1)t},$ $\displaystyle
n_{i}(t)$ $\displaystyle=$ $\displaystyle n_{i}(0)=1.$ (46)

## Appendix B: Linear stability analysis of the interacting system.

Eqs. (8) are real equations describing the dynamics of interacting system. In
order to analyze their linear stability we construct the Jacobian matrix for
$2N_{s}$ variables
$n_{1},n_{2},...,n_{N_{s}},\varphi_{1},\varphi_{2},...,\varphi_{N_{s}}$

$J=\begin{pmatrix}\derivative{\dot{n}_{1}}{n_{1}}&\dots&\derivative{\dot{\varphi}_{N_{s}}}{n_{1}}\\\
\vdots&\ddots&\vdots\\\
\derivative{\dot{n}_{1}}{\varphi_{N_{s}}}&\dots&\derivative{\dot{\varphi}_{N_{s}}}{\varphi_{N_{s}}}\end{pmatrix}$
(47)

where $\varphi_{i}=\theta_{i,i+1}$.

Fixed points are determined from the steady state condition of Eqs.(8)

$x=n_{i}=1\ ,\ \ y=\varphi_{i}=\frac{2\pi}{N_{s}}k,\quad i,k=1,...,N_{s}.$
(48)

The Jacobian matrix $J$ at the fixed points can be written as a two by two
block matrix

$J=\begin{pmatrix}S&C\\\ D&S\\\ \end{pmatrix}.$ (49)

Here the matrices $S$, $C$ and $D$ are circulant $N_{s}$x$N_{s}$ matrices,
whose elements can be written as

$\displaystyle S_{ij}$ $\displaystyle=$
$\displaystyle\sin(y)(\delta_{i,j+1}-\delta_{i,j-1}),$ $\displaystyle D_{ij}$
$\displaystyle=$ $\displaystyle 2x\cos(y)(-\delta_{i,j}+\delta_{i,j+1}),$
$\displaystyle C_{ij}$ $\displaystyle=$
$\displaystyle\left(u+\frac{3\cos(y)}{2x}\right)(\delta_{i,j}-\delta_{i,j-1})$
(50) $\displaystyle+$ $\displaystyle\frac{\cos(y)}{2x}\
(\delta_{i,j-2}-\delta_{i,j+1}).$

We see that only $C$-matrix depends on interaction $u$. All circulant matrices
of the same size have the same eigenvectors, thus all circular matrices of the
same size can be diagonalized by $J_{S}=U^{-1}SU$ where $U$ matrix columns are
the circulant matrix eigenvectors. Using the rules for inverses of block
matrices, we perform a simple transformation to diagonalize the circulant
matrices in the Jacobian matrix

$\begin{pmatrix}U^{-1}&0\\\ 0&U^{-1}\\\ \end{pmatrix}\begin{pmatrix}S&C\\\
D&S\\\ \end{pmatrix}\begin{pmatrix}U&0\\\ 0&U\\\
\end{pmatrix}=\begin{pmatrix}J_{S}&J_{C}\\\ J_{D}&J_{S}\\\ \end{pmatrix}$ (51)

The similar transformation preserves the eigenvalues, the $J_{S}$ , $J_{C}$ ,
$J_{D}$ matrices are now diagonal matrices. As all of the new diagonal
matrices commute, we can take advantage of a block matrix determinant rule to
find the eigenvalues of $J$

$Det\begin{pmatrix}J_{S}-\lambda I&J_{C}\\\ J_{D}&J_{S}-\lambda I\\\
\end{pmatrix}=Det((J_{S}-\lambda I)^{2}-J_{D}J_{C})=0$ (52)

This equation can be rewritten as

$(\lambda_{j}^{S}-\lambda_{j})^{2}-\lambda_{j}^{D}\lambda_{j}^{C}=0,\quad
j=1,...,N_{s}.$ (53)

where $\lambda_{j}^{S},\lambda_{j}^{D},\lambda_{j}^{C}$ are the eigenvalues of
$S,D,C$ respectively and $\lambda_{j}$ is the eigenvalue of the Jacobian $J$.
We get

$\lambda_{j}=\lambda_{j}^{S}\pm\sqrt{\lambda_{j}^{D}\lambda_{j}^{C}},\quad
j=1,...,N_{s}.$ (54)

Now we collect the eigenvalues of the circulant matrices $S,D$ and $C$, which
can be found following the method explained in Appendix A:

$\displaystyle\lambda_{j}^{S}$ $\displaystyle=$ $\displaystyle
2\mathrm{i}\sin(y)\sin{\left(\frac{2\pi j}{N_{s}}\right)},$
$\displaystyle\lambda_{j}^{D}$ $\displaystyle=$ $\displaystyle
2x\cos{y}\left(\exp{-\mathrm{i}\frac{2\pi j}{N_{s}}}-1\right),$
$\displaystyle\lambda_{j}^{C}$ $\displaystyle=$
$\displaystyle\left(u+\frac{3\cos(y)}{2x}\right)\left(1-\exp{\mathrm{i}\frac{2\pi
j}{N_{s}}}\right)$ (55) $\displaystyle+$
$\displaystyle\frac{\cos(y)}{2x}\left(\exp{\mathrm{i}\frac{4\pi
j}{N_{s}}}-\exp{-\mathrm{i}\frac{2\pi j}{N_{s}}}\right).$

As a result we get for $\lambda_{j}=\lambda_{j}(k)$

$\begin{split}\lambda_{j}(k)&=2\mathrm{i}\left\\{\sin\theta_{0}^{k}\sin{\left(\frac{2\pi
j}{N_{s}}\right)}\pm\right.\\\ &\left.\sin{\left(\frac{\pi
j}{N_{s}}\right)}\sqrt{2\cos\theta_{0}^{k}\left[2\cos\theta_{0}^{k}\sin^{2}\left(\frac{\pi
j}{N_{s}}\right)+u\right]}\right\\},\end{split}$ (56)

The eigenvalues are purely imaginary if the expression under the square root
is non-negative for all $j$. This would correspond to a neutral center and
stable system. If at least one of the eigenvalues acquires a real part, this
will correspond to an exponential instability in the system. We discuss this
in more detail in the main text after Eq.(21).

## References

  * (1) T. L. Gustavson, P. Bouyer, and M. A. Kasevich, Precision Rotation Measurements with an Atom Interferometer Gyroscope, Phys. Rev. Lett.78, 2046 (1997).
  * (2) B. P. Anderson and M. A. Kasevich, Macroscopic Quantum Interference from Atomic Tunnel Array, Science 282, 1686 (1998).
  * (3) D. W. Hallwood, T. Ernst, and J. Brand, Robust mesoscopic superposition of strongly correlated ultracold atoms, Phys. Rev. A 82, 063623 (2010).
  * (4) L. Amico1, D. Aghamalyan, F. Auksztol, H. Crepaz, R. Dumke and L. C. Kwek, Superfluid qubit systems with ring shaped optical lattices, Sci. Rep. 4, 04298 (2014).
  * (5) D. Aghamalyan, M. Cominotti, M. Rizzi, D. Rossini, F. Hekking, A. Minguzzi, L. C. Kwek and L. Amico, Coherent superposition of current flows in an atomtronic quantum interference device, New J. Phys. 17, 045023 (2015).
  * (6) G. Arwas and D. Cohen, Chaos and two-level dynamics of the atomtronic quantum interference device, New. J. Phys. 18, 015007 (2016).
  * (7) C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, Experimental Realization of Josephson Junctions for an Atom SQUID, Phys. Rev. Lett. 111, 205301 (2013).
  * (8) M. Tsubota and K. Kasamatsu, Josephson Current Flowing in Cyclically Coupled Bose-Einstein Condensates, J. Phys. Soc. Jpn. 69, 1942 (2000).
  * (9) K. Kasamatsu and M. Tsubota, Vortex generation in cyclically coupled superfluids and the Kibble-Zurek machnism, J. Low Temp. Phys. 126, 315 (2002).
  * (10) D. R. Scherer, C. N. Weiler, T. W. Neely, and B. P. Anderson, Vortex formation by Merging of Multiple Trapped Bose-Einstein Condensates, Phys. Rev. Lett. 98, 110402 (2007).
  * (11) Gh.-S. Paraoanu, Persistent currents in a circular array of Bose-Einstein condensates, Phys. Rev. A 67, 023607 (2003).
  * (12) E. T. D. Matsushita and E. J. V. de Passos, Stability of Bose-Einstein condensates in a circular array, arXiv:0909.0920 (cond-mat) (2009).
  * (13) C. Arwas and D. Cohen, Chaos, metastability and ergodicity in Bose-Hubbard superfluid circuits, AIP Conference Proceedings 1912, 0200001 (2017).
  * (14) C. Arwas and D. Cohen, it Monodromy and chaos for condensed bosons in optical lattices, Phys. Rev. A 99, 023625 (2019).
  * (15) K. Nemoto, C. A. Holmes, G. J. Milburn, and W. J. Munro, Quantum dynamics of three coupled atomic Bose-Einstein condensates, Phys. Rev. A 63, 013604 (2000).
  * (16) J. Dziarmaga, M. Tylutki, and W. H. Zurek, Ring of BEC pools as a trap for persistent flow, Phys. Rev. B 84, 094528 (2011).
  * (17) S. Moulder, S. Beattie, R. P. Smith, N. Tammuz, and Z. Hadzibabic, Quantized supercurrent decay in an annular Bose-Einstein condensate, Phys. Rev. A 86, 013629 (2012).
  * (18) T. Cookson, K. Kalinin, H. Sigurdsson, J. Töpfer, S. Alyatkin, M. Silva, W. Langbein, N. G. Berloff and P. G. Lagoudakis, Geometric frustration in polygons of polariton condensates creating vortices of varying topological charge, Nat. Comm. 12, 1 (2021).
  * (19) F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi, M. Inguscio, Josephson Junction Arrays with Bose-Einstein Condensates, Science 293, 843 (2001).
  * (20) 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).
  * (21) A. Trombettoni and A. Smerzi, Discrete Solitons and Breathers with Dilute Bose-Einstein Condensates, Phys. Rev. Lett. 86, 2353 (2001).
  * (22) A. Smerzi, A. Trombettoni, P. G. Kevrekidis, and A. R. Bishop, Dynamical Superfluid-Insulator Transition in a Chain of Weakly Coupled Bose-Einstein Condensates, Phys. Rev. Lett. 89, 170402 (2002).
  * (23) M. Trujillo-Martinez, A. Posazhennikova, and J. Kroha, Nonequilibrium Josephson oscillations in Bose-Einstein condensates without dissipation, Phys. Rev. Lett. 103, 105302 (2009).
  * (24) P. J. Davis, Circulant Matrices, 2nd Edition, Chelsea Publishing, (1994).