Superfluidity and Chaos in low dimensional circuits

The hallmark of superfluidity is the appearance of “vortex states” carrying a quantized metastable circulating current. Considering a unidirectional flow of particles in a ring, at first it appears that any amount of scattering will randomize the velocity, as in the Drude model, and eventually the ergodic steady state will be characterized by a vanishingly small fluctuating current. However, Landau and followers have shown that this is not always the case. If elementary excitations (e.g. phonons) have higher velocity than that of the flow, simple kinematic considerations imply metastability of the vortex state: the energy of the motion cannot dissipate into phonons. On the other hand if this Landau criterion is violated the circulating current can decay. Below we show that the standard Landau and Bogoliubov superfluidity criteria fail in low-dimensional circuits. Proper determination of the superfluidity regime-diagram must account for the crucial role of chaos, an ingredient missing from the conventional stability analysis. Accordingly, we find novel types of superfluidity, associated with irregular or chaotic or breathing vortex states.

Scientific RepoRts | 5:13433 | DOi: 10.1038/srep13433 d occupation differences and d relative phases). Accordingly the effective number of degrees of freedom becomes d.
Intensive studies have focused on the integrable Bosonic Josephson Junction (M = 2 hence d = 1). In few-site Bose-Hubbard systems with M > 2, such as the M = 3 trimer [33][34][35][36][37][38][39][40][41][42][43][44] , an utterly new perspective is essential due to the emergence of chaos. So far, the implications of chaos on superfluidity have not been illuminated, neither for the trimer ring 45 nor for M > 3 rings. We would like to highlight a crucial difference between the two cases: For M = 3 ring the d = 2 dimensional KAM tori in phase-space can divide the 2d − 1 energy-shell into separate regions, while for M > 3 (say M = 4, hence d = 3) it is not the case (a 3D surface cannot divide a 5D energy-shell into disjoint territories). Consequently for M > 3 all phase-space regions are interconnected via "Arnold diffusion".

Beyond the traditional view
The traditional criterion for superfluidity associates the metastability of the current with the existence of a stationary stable fixed point of the Hamiltonian flow, that supports a coherent vortex state. Accordingly a Bogoliubov de Gennes (BdG) linear-stability-analysis is assumed valid for determination of the regimes where superfluidity is anticipated. Low dimensional circuits have phase-space with both chaotic and quasi-regular motion. Consequently the traditional BdG paradigm is challenged: (i) Dynamical instability of a vortex state does not necessarily mean that superfluidity is diminished, because its collapse may be topologically arrested by KAM structures; (ii) Linear BdG stability of a vortex state does not always imply actual stability, because Arnold diffusion can provide detour paths out of seemingly elliptical regions; (iii) Due to quantum fluctuation, i.e. finite uncertainty width of a vortex-state, stability is required within a Plank cell around the fixed-point.
The result of these three observations is a novel superfluidity regime-diagram, quite distinct from the one that would be obtained using standard criteria. Considering that d = M − 1, and that Arnold diffusion can only take place when d > 2, there should be a dramatic difference between trimers (M = 3) and larger rings (M > 3): For the trimer, item (i) implies that superfluidity can persist even if the motion becomes chaotic; For larger rings, item (ii) implies that BdG (linear) dynamical stability is not a sufficient condition; while item (iii) implies that global analysis of phase-space topography is essential. In the extreme limit of M → ∞ one should remember that the dynamics become integrable due to rotational symmetry. Below we demonstrate how the above ideas affect the regime diagram of a few site ring. For the trimer (M = 3) superfluidity manifests itself beyond the regime of dynamical stability, while for M > 3 we find a much more intricate situation.

The Model and the superfluidity regime diagram
Consider N mass m Bosons in an M site ring that has a radius R. If the ring is rotated with frequency Ω , one may transform into the rotating frame, where the potential is time-independent. In this frame we have Coriolis force, which is formally like having a magnetic flux Φ = (2πmR 2 /ħ)Ω through the ring. Accordingly the system is described by the Bose-Hubbard Hamiltonian 29,30 ( ) Here j mod (M) labels the sites of the ring, a i and a i † are destruction and creation operators. K is the hopping frequency, and U is the on-site interaction. The Peierls phase factor includes the so-called Sagnac phase Φ . Optionally it can be realized by introducing a spatially-adiabatic variation of the atomic magnetic dipole orientation 46 , or gauge field as in 47 , see 48 . Without loss of generality Φ ∈ [0, π], and K > 0, and U > 0. Negative U is the same as positive U with a flipped energy landscape ( Negative Φ is related to positive Φ by time reversal. The Hamiltonian H commutes with the total particle number = ∑ † N a a i i i , hence the operator N is a constant of motion. In a semi-classical context one defines phase-space action-angle coordinates as follows: and the other is the phase Φ . Note that the re-scaling of the canonical variables implies that n is replaced by n/N. Hence upon quantization ϕ and n are conjugate with dimensionless Plank constant ħ = 1/N. The superfluidity regime-diagram. In Fig. 1 we plot the numerically determined (Φ , u) regime diagram for the superfluidity of rings with M = 3, 4, 5 sites. Image colors depict the current = (−∂ /∂Φ) H I for the eigenstate that carries maximal current. The solid and dashed lines indicate the energetic and the dynamical stability borders, as determined from the BdG analysis (see below). The regime diagrams do not agree with the traditional analysis: For the M = 3 ring superfluidity persists beyond the border of dynamical stability, while for M > 3 the dynamical stability condition is not sufficient.
In the absence of interaction (U = 0), time reversal is broken by having a non-zero Φ , and an eigenstate can carry a non-zero "persistent current". If we add some weak disorder W (random on-site potential) this current becomes smaller, and it diminishes in the limit Φ → 0. In contrast superfluidity features The I of the state that carries maximal current is imaged at the background. We observe that the stability regions (large current) are not as expected from the linear stability analysis: the solid line indicates the energetic-stability border; the dashed lines indicate the dynamical stability borders. For clarity we also include a negative u region which is in fact a duplication of the upper sheet. In (a) the dotted line indicates the "swap" transition (see text); and the dots labeled (a-f) mark (Φ , u) coordinates that are used in Fig. 3 to demonstrate the different regimes.
a macroscopically large metastable current that is achieved due to having a non-zero interaction (U ≠ 0).
Superfluidity is feasible if a middle vortex state (see below), or some irregular variant of it, maintains stability; otherwise it would mix with all the other eigenstates that reside in the same energy-shell, resulting in a micro-canonically small current. Accordingly, the stability of the current should be verified with respect to an added disorder W. In the model under study superfluidity is indeed maintained as long as W < U, and accordingly a finite strength of interaction U is required.

Metastability of vortex states
The stationary orbitals of a single particle in a clean ring are the momentum states with wavenumber k = (2π/M)m, where m is an integer modulo M. Coherent vortex states have N particles condensed into the same momentum orbital. From a semiclassical perspective, a coherent-state is represented by a Gaussian-like phase-space distribution. Such state is stable if it is supported by a region where the classical motion is "locked". Figure 2 is a cartoon that summarizes various possibilities. We first consider the traditional possibility of having a regular vortex-state that is supported by a stable fixed-point. In later sections we discuss the possibility of having irregular vortex-states that are supported by chaotic regions. tions. If the ω q are real, the linearized motion is "elliptic" around the fixed-point. But if some of the ω q acquire an imaginary part, the motion becomes unstable ("hyperbolic"), meaning that trajectories depart exponentially from the fixed point. In the latter case a chaotic motion with a non-zero Lyapunov exponent is implied. In the present context there is a cyclic degree of freedom that corresponds to N. Accordingly there is a pair of zero eigenvalues that should be excluded. Effectively we deal with 2d × 2d matrix where d = M − 1). It is also useful to notice that at the energetic-stability border, where the fixed-point becomes a saddle, we have

Stability analysis.
Our linear-stability analysis of the vortex states is analogous to that of 7 , and will be presented in the following sections, with additional details in the Methods section, leading to the solid and dashed lines of Fig. 1. Stability with respect to added disorder. Before we proceed we would like to make a comment regarding the implications of having weak disorder. The essence of superfluidity is the meta-stability of the current. This makes it distinguished from having merely a persistent current. Consider a clean ring with non-interacting particles and no disorder. We can prepare there a vortex state with macroscopically large current. However, any small disorder will randomize the velocity (as in Drude model) and the current will diminish. Adding interaction U between the particles change the picture dramatically. Now the energy landscape E = H(z) has a non-trivial topography. It may have (say) a valley with a local minimum, that can support a metastable state. If we add disorder it merely deforms the valley, while the fixed-point remains stable. Disregarding dimensionless prefactors this is true as long as W < U. It follows that finite U is essential in order to have physical metastability. We would like to emphasize that the same conclusion holds for any "energy landscape". Adding small disorder will affect only "degenerated" regions, say the separatrix region, but will not affect the overall phase-space structure. This is called "structural stability". To be on the safe-side, we have verified that the numerically determined regime diagrams are not affected by adding weak disorder.
Quantum fluctuations. As shown in Fig. 1 the observed superfludity regimes are not in accordance with the traditional BdG analysis. Irrespective of M, as we go higher in u, superfluidity is diminished even in the energetically-stable region. This is conspicuous for low N, as in Fig. 1c where we have ~2 particles per site, and can be explained as the consequence of having a finite uncertainty width. Namely, as u is increased the radius of the stability island (if exists) becomes smaller, until eventually it cannot support a stable vortex state. Equivalently, as N becomes smaller the uncertainty width of a vortex state becomes larger, until it 'spills' out of the stability island. This type of reasoning resembles the semiclassical view of the Mott transition (see below). Taking a closer look at the regime diagrams one observes that the above "quantum fluctuations" perspective is not enough in order to explain the observed differences. We therefore turn to provide a more detailed phase-space picture.
The ground state. The lowest energy fixed-point (m = 0 vortex state) is stable for any positive u.
However it is located in an island which is surrounded by a chaotic sea. As u is increased the island's size decreases until at u > N 2 /M, it becomes smaller than a single Planck cell. At this point it can no longer accommodate a vortex state and one observes a quantum Mott transition. See Fig. 2d for pedagogical illustration. Other studies of rotating ring lattices [20][21][22][23][24][25] have addressed additional quantum issues, such as the appearance of "cat states": In the trimer with Φ = π the ground state might be a macroscopic superposition of the degenerate vortex states m = 0 and m = 1.

Solitons.
While our main focus is on the stability of vortex states we briefly discuss other fixed points that were of interest in past work. The trimer model without rotation (Φ = 0) has been the subject of intense study [33][34][35][36][37][38][39][40][41][42][43][44] . In particular it has been noted that the Hamiltonian H(z) has additional non trivial fixed-points away from the symmetry axis. The simplest of which are self-trapped "bright solitons" obtained via bifurcation of a vortex state, in which the particles are localized in a single site. See Fig. 2e for pedagogical illustration. Other notable fixed-points correspond to single-depleted-well states in which one site is empty, while the remaining two sites are equally occupied by the particles. For M > 3 rings there are off-axis fixed points that support spatially modulated vortex states 49 .

The M = 3 ring
We focus our attention on the central region of energies, where the middle vortex state (m = 1) is located. Examples for phase-space trajectories at this energy are displayed in Fig. 3, and will be further discussed below. First we examine the linear stability of the pertinent vortex fixed-point.    The fixed-point becomes dynamically unstable if the eigenvalues acquire a real part, which is the so-called Lyapunov exponent. This happen when b 2 − 4c < 0, leading to dynamical instability in the region In principle for b 2 − 4c > 0 the condition b < 0 would imply an additional dynamical instability regime. But here b < 0 occurs inside the region of b 2 − 4c < 0. The stability borders are demonstrated in Fig. 1a.
Beyond BdG. In Fig. 1a we observe macroscopically large currents beyond the expected region. This has been noted in 45 without explanation. In particular we see that supefluidity survives in the limit Φ → 0, contrary to the expectation 50,51 that is based on the traditional stability argument. We can plot the wavefunction r rE 2 2 Ψ( ) = α of an eigenstate that supports superfluidity. Here the components of r are the coordinates r n n N 2 1 = ( − )/ and r ┴ = n 3 /N. The wavefunction of a standard regular vortex is merely a hump at the symmetry point n 1 = n 2 = n 3 = N/3. In Fig. 4 we display examples for the wavefunctions of non-standard vortex states: a chaotic vortex state, and a breathing vortex state. The terms "chaotic" and "breathing" is related to the underlying classical dynamics [options (ii) and (iii) below].
Launching trajectories at the vicinity of the vortex fixed-point it has been demonstrated that both stable and unstable oscillations can emerge 37 . Specifically we encounter the following possibilities: (i) the trajectories are locked in the vicinity of the vortex fixed point; (ii) the trajectories are quasi-periodic in phase-space; (iii) the trajectories are chaotic but uni-directional. The above possibilities are pedagogically illustrated by the cartoon of Fig. 2. Poincare sections of the trajectories (see below) reveal that a regular vortex state is supported by a regular island around the fixed point (case i); a breathing vortex is supported by a secondary island that has been created via bifurcation (case ii); while a chaotic vortex state is supported by a 'chaotic pond' of clockwise motion that does not mix with the anti-clockwise motion (case iii). Consequently the motion may become chaotic, but stay uni-directional, and superfluidity persists contrary to the common expectation. The route to chaos. As observed in Fig. 1a superfluidity is quite robust. The current diminishes only in the vicinity of what we call "swap transition", which is indicated in the figure by dotted line. In the Methods section we derive an explicit expression for the transition line: At the transition the two separatrixes that dominate the structure of phase-space coalesce, and consequently a global chaotic-sea is formed. The details of this transition are provided below. Such type of separatrix overlap bifurcation has been once encountered in molecular physics studies 52 , but not in chaotic context.
Phase-space tomography. In Fig. 3 we plot the spectrum of the trimer Hamiltonian for representative values of (Φ , u) that are indicted in Fig. 1a. We also plot in each case the Poincare sections at the energy of the m = 1 vortex fixed-point (see Methods section for details). From the quantum spectrum we can easily deduce the phase-space structure at any other energy. One can call it "quantum phase-space tomography". Consider for example Fig. 3c. We can easily correlate the largest current states with the red (upper) island; the secondary group of large current states with the yellow (left) island; and the small current states with the green chaotic sea. Additional information can be extracted from the purity of the states. For precise definition see 45 .
Here it is enough to say that S = 1 means that all the particle are condensed in a single orbital, while S < 1 means that the state is fragmented into 1/S orbitals. For ergodic state 1/S ~ M. Points in the spectrum are color-coded from black (S ~ 1) to purple (S ~ 1/M). One observes that the non-standard vortex states have high but not perfect purity. By inspection of Fig. 3 we observe the following regimes in the diagram of Fig. 1a: Regime (S) stands for simple phase-space structure with energetically-stable clockwise ("red") and anti-clockwise ("blue") islands that are separated by a forbidden region. In regime (B) we have two regular regions of clockwise motion, and "blue separatrix" that supports anti-clockwise motion. As we go up in u the blue separatrix becomes a chaotic sea. In regime (D) the middle vortex bifurcates, while the other clockwise island remains regular. In regime (D') the "vortex separatrix" swaps with the "blue separatrix". This swap is clearly demonstrated as we go from Fig. 3d,e. The border between regimes (D) and (D') is shown as a dotted line in Fig. 1a, see equation (8). Along this line the two separatrices coalesce. Crossing to regime (B') the bifurcation that is responsible to the "blue separatrix" is undone, and eventually we can go back to the (S) regime via a simply-connected (A) regime that has a simple structure with no separatrix.
In region (B) the vortex is not energetically stable: it is located on a saddle point in phase-space. Nevertheless dynamical stability is maintained. In region (D) the vortex is no longer dynamically stable, and the trajectories at the vicinity of the vortex are chaotic. Still the motion is confined by KAM tori within a "chaotic pond", and therefore remains uni-directional. In the vicinity of the swap, as we go up in u, the chaotic pond becomes a chaotic sea, and the superfluid current is diminished.
Upon quantization the chaotic pond can support a "chaotic vortex state", which has been illustrated in Fig. 4a. A second class of large current states are supported by stable periodic-orbits (POs) that has been bifurcated from the stationary vortex fixed-point, once the latter lost stability. These POs are elliptic fixed-points of the Poincare section, see Fig. 3d. Upon quantization the associated islands can support a "breathing vortex state", see Fig. 4b.

The M > 3 rings
Considering a no-rotating device, the traditional stability argument 50,51 implies marginally stable superfluidity for an M = 4 device, and stability if u is large enough for an M > 4 device. These observations are implied by the stability borders that are plotted in Fig. 1bc. The explicit expressions for the stability borders are provided in the following paragraph.

BdG analysis.
In the Methods section we derive an explicit expression equation (19) for the Bogoliubov frequencies ω q,± of the k m vortex. The frequencies are indexed by q = (2π/M) × integer. As we move in the (Φ , u) regime-diagram the ω q,± go through zero as we cross the following lines The border of the stable regime is determined by the first line that is encountered, which corresponds to the minimal value q = 2π/M. The dynamical stability is lost once one of the frequencies become complex. This happens when either of the following two expressions equal zero Again we have to substitute q = 2π/M, which is associated with the first encounter. The explicit results for M = 3 are easily recovered.
Beyond BdG. Looking on the numerically determined current (Fig. 1b) one observes that superfluidity can persist slightly beyond the dynamical stability border. But much more conspicuous is the diminishing of superfluidity within a large region where the BdG analysis predicts dynamical stability. We find (see below) that the latter effect is related to Arnold diffusion. Namely, if d > 2, the d dimensional KAM tori in phase space are not effective in blocking the transport on the 2d − 1 energy shell (as discussed in the Introduction). As u becomes larger this non-linear leakage effect is enhanced, stability of the motion is deteriorated, and the current is diminished. At this point it is helpful to distinguish between strict dynamical stability and linear dynamical stability 53 . For M > 3 the latter does not imply the former (all regions of phase space are interconnected via "Arnold diffusion" irrespective of their linear stability). As we go up in u the chaos becomes predominant, and consequently energetic-stability rather than dynamical-stability becomes the relevant criterion. It follows from this distinction that for M > 3 one has to distinguish between regular and irregular vortex states. This distinction is demonstrated in Fig. 5. A regular vortex is represented by a simple hump at the central (n i = N/M) fixed-point, whereas an irregular vortex has a richer structure that reflects the underlying fragmented phase-space structure.
In order to verify the above semiclassical reasoning, we try in Fig. 6 to reconstruct the quantum regime-diagram via classical simulations. This reconstruction provides a qualitative proof for the semiclassical reasoning, and furthermore demonstrates the N dependence of the the superfluity regime diagram. Namely, we launch a Gaussian cloud of trajectories that have an uncertainty width that corresponds to N. The fraction of trajectories that escape is used as a measure for the stability. The practical criterion for escape is having the average current I(t) getting below some threshold I ∞ within some time t ∞ . In principle t ∞ should be the Heisenberg time (proportional to N d ), and I ∞ can be (say) half I m . In practice the result is not sensitive to these thresholds. Regions where the current diminishes in-spite of BdG dynamical stability, and does not recover even if N is increased (smaller uncertainty width), establish the relevance of Arnold diffusion.

Discussion
A recent experiment with a toroidal ring 16 has demonstrated hysteresis in a quantized superfluid 'atomtronic' circuit. The procedure was to prepare a stable vortex state at rotation frequency that corresponds to Φ 1 , and then to check its stability after changing to Φ 2 . The theory there could be regarded as an extension of the traditional "energetic-stability" reasoning that underlays the Landau criterion. Our work has been motivated by the following question: what would be the results of the same type of experiment,  Fig. 4 with one extra dimension, and the color code is the same. The parameters (Φ , u) are (0.2π, 0) and (0.2π, 3) respectively. Note that a regular vortex is represented by a simple hump at the central (n 1 = n 2 = n 3 = n 4 ) fixed-point, whereas an irregular vortex has a richer structure that reflects the fragmented phase-space structure.  Fig. 1b. The lines are the BdG stability borders. Given (Φ , u) we launch a Gaussian cloud of trajectories that have an uncertainty width that corresponds to N. The fraction of trajectories (blue ~ 100% to red ~ 0%) that escape is used as a measure for stability (see text for details). Results are displayed for clouds that have uncertainty width Δ ϕ ~ π/2 (left) and Δ ϕ ~ π/4 (right).
if the toroidal ring had several weak links, or if it were replaced by a discrete M site ring. The realization of such ring that is described by equation (1) has been discussed in 21,32,46 . We argue that the chaotic nature of such circuit requires to go beyond the traditional analysis; else one would not be able to predict the borders of the stability-regime in the hysteresis loop. On the theoretical level we have highlighted a novel type of superfluidity that is supported by irregular or chaotic or breathing vortex states. Such states are supported by fragmented regions in phase-space (M > 3), or by chaotic ponds (M = 3), or by periodic-orbits respectively, hence they are missed by the traditional BdG analysis. Furthermore we have highlighted the limitations of the linear stability analysis for high dimensional chaos (M > 3).
As a secondary message, we would like to emphasize that the gross features of the classical phase-space can be easily extracted from the spectrum of the quantized Hamiltonian. To get the same information via classical analysis would be an extremely heavy task that would require generation of many trajectories in numerous phase space regions, on each possible energy shell, as opposed to our "quantum phase space tomography" which requires a single diagonalization of a finite matrix. If Nature were classical, Quantum Mechanics still would be invented as a valuable tool, just for the purpose of analysing mixed complex dynamics.

Methods
Poincare sections. Starting with the Hamiltonian equation (1) with M = 3, written in terms of action-angle variables, the classical dynamics is generated by the equation with scaled units such that K = N = 1. We solve this equation numerically. For plotting of trajectories it is convenient to use the coordinates (n 1 − n 3 , ϕ 1 − ϕ 3 ) and (n 3 − n 2 , ϕ 3 − ϕ 2 ). The section chosen is n 3 − n 2 = 0, at the energy of the m = 1 vortex E = (u/6) − cos[(2π/3) − (Φ /3)]. Given a phase space section point (n 1 − n 3 , ϕ 1 − ϕ 3 ), the equation H(z) = E has either zero or two solutions for the remaining coordinate ϕ 3 − ϕ 2 . This implies that the Poincare section has two sheets. For presentation purpose we pick the sheet where velocity ∂ t (n 3 − n 2 ) has a larger value. In Fig. 3 we plot Poincare sections at the energy of the m = 1 vortex fixed-point. The section is n 3 − n 2 = 0. For convenience we use the canonical coordinates (n 1 − n 3 , ϕ 1 − ϕ 3 ) and a scaled particle number n n N /  . Note that the m = 1 vortex fixed-point is always located at (0, 2π/3). The boundary of the allowed region is marked by a black line. Each trajectory in the Poincare section is a set of points that share the same color. The color reflects the average current, which is calculated by taking a time average over I.
We have here a d = 2 system, so the energy shell is 3D, while the Poincare section is 2D. According to KAM, some regular motion survives: these are the KAM tori. In the PS they look line lines that divide the motion into "territories". For example in Fig. 3b we see that most of the motion is regular. We have there some bounded chaotic motion in the vicinity of the blue separatrix, but near the vortex fixed-point the motion is predominantly regular (dense set of KAM tori). In Fig. 3c most of the KAM tori have been destroyed. Still the remaining KAM tori separate the red chaotic pond from the green chaotic sea. If these tori were destroyed, the motion would become globally chaotic (which is not the case).

Stability analysis for general M.
It is convenient to do the general stability analysis using the common many-body "quantum" notations. The canonical variables are a j and its conjugate instead of (ϕ, n). Then we transform to b k that are defined via the substitution