Quasistatic transfer protocols for atomtronic superfluid circuits

Quasi-static protocols for systems that feature a mixed phase-space with both chaos and quasi-regular regions are beyond the standard paradigm of adiabatic processes. We focus on many-body system of atoms that are described by the Bose–Hubbard Hamiltonian, specifically a circuit that consists of bosonic sites. We consider a sweep process: slow variation of the rotation frequency of the device (time dependent Sagnac phase). The parametric variation of phase-space topology implies that the quasi-static limit is not compatible with linear response theory. Detailed analysis is essential in order to determine the outcome of such transfer protocol, and its efficiency.


Results
The model. We consider a system with N bosons in a 3-site ring. The system is described by the Bose-Hubbard Hamiltonian [Methods] with hopping frequency K and on-site interaction U. The sweep control-parameter is the Sagnac phase , which is proportional to the rotation frequency of the device: it can be regarded as the Aharonov-Bohm flux that is associated with Coriolis field in the rotating frame 59,60 . There are 3 momentum orbitals k = 0, ±2π/3 . Initially all the particles are condensed in k = 0 . A caricature for the preparation is provided in Fig. 1 (left panel).
Following 51 we define a depletion coordinate n and an imbalance coordinate M, such that the occupations of the orbitals are n 0 = N−2n , and n ± = n±M . The model Hamiltonian can be written in terms of (n, M), and the conjugate phases (ϕ, φ) . Namely [Methods]: The first term H (0) is an integrable piece of the Hamiltonian that has M as a constant of motion: while the additional terms induce resonances that spoil the integrability, and give rise to chaos: The hopping frequency K and the Sagnac phase hide in the expression for the energy of the condensate, and in the detuning parameters: We also note that the energies of the totally depleted states ( n = (N/2) ) are (1) H(ϕ, n; φ, M) = H (0) (ϕ, n; M) + H (+) + H (−) (2) Φ < Φ mts m ts Φ > Φ ε ε Figure 1. Orbital occupation. For the purpose of illustration we consider a ring with N = 4 particles. The orbitals are represented by horizontal lines (the horizontal shifts hints the sign of the momentum). Initially (left panel) the particles are condensed in the #0 momentum orbital. As is increased beyond mts this configuration becomes metastable (right). We ask what is the moment when the #0 orbital is depleted, and what is the final distribution of the particles. The N = 4 system has 15 energy levels that corresponds to the different possibilities to distribute the particles between the orbitals. In the presence of non-zero interaction those levels are partially mixed. www.nature.com/scientificreports/ Note that the latter expression has zero contribution from the H (±) terms. The chaos affects the pathway between the initial condensate at n = M = 0 , and the peripheral depleted states at n = (N/2) , but has only little effect on the gross features of the energy landscape.
Metastability. The central point in phase space n = M = 0 is a stationary point (SP) of the Hamiltonian for any , meaning that we have there ṅ = 0 . But this does not mean that this SP is stable. As implied by the caricature of Fig. 1, the condensate at n = 0 is no longer situate at the minimum of the energy landscape once Semiclassics. The classical (as opposed to semiclassical) treatment of the Hamiltonian is commonly termed Mean Field Theory (MFT). The evolving state is represented by a single point in phase space. We can scale the time such that t := Kt , and the occupations such that n := n/N . Then one finds that the dynamics is controlled by the dimensionless interaction parameter Upon quantization (aka "second quantization") the scaled value of the Planck constant is = 1/N , see e.g. 49,50,52 . Quantum states can be represented in phase space by their Wigner function. In particular the initial coherent state at n = 0 is represented in phase-space by a Gaussian-like distribution of radius R ∼ 1/N. What we call "semiclassical treatment" is far better and reliable compared to MFT, and is commonly called Truncated Wigner Approximation (TWA). Within the framework of TWA the Moyal brackets are approximated by Poisson brackets, which means that the Wigner function is propagated by the classical equations of motion.
The TWA is very accurate as long as quantum tunneling is neglected. The tunneling amplitude scales as exp[−Action/ ] , where = 1/N . Therefore it is much slower compared with any classical process. Discussion of tunneling in the BHH context can be found in 62 , and later we demonstrate numerically that it can be neglected for a simulation with N = 30 particles.

Simulations.
We describe the results of semiclassical simulations. Detailed analysis will follow after that.
The condensate preparation at = 0 is represented by a Gaussian cloud of points in phase space, at the central SP ( n = 0 ). The evolution of the cloud in a dynamical sweep simulation is demonstrated in Fig. 2. The colorcode shows the evolution of the depletion coordinate (n), and the vertical position of the cloud points indicate the population imbalance M (left panels), or the energy E = H (right panels), or the current I = −∂H/∂� as a function of time (inset). For the latter we use the following expression in terms of (n, M), Note that the cloud is a semiclassical representation of the evolving state. Accordingly, to get the expectation value of the energy or of the current, an average has to be taken over the ensemble of evolving trajectories. In Fig. 2 the average is not taken in order to provide an insight for the dispersion as well.
The cloud follows the ground state energy E GS only up to mts . Then it continues to follow the condensate energy E 0 during an additional time interval. The cloud starts spreading not before stb , and not later than dyn . The spreading is indicated by the departure of energy from E 0 . The depletion of the condensate is indicated by the color that changes abruptly from blue ( n = 0 ) to red ( n∼N/2 ). It takes place during a distinct short time interval when �(t) ∼ � swp . The depletion stage is also clearly reflected as a jump in the current-versus-time plot. Finally, the subsequent evolution after the depletion does not follow any of the adiabatic E n curves, as discussed further below.
We display in Fig. 2 three representative simulations: very slow sweep (top row), optimal sweep rate (middle row), and faster sweep (lower row). The results for many such simulations are gathered in Fig. 3, where the dependence of n and M on the sweep rate is demonstrated for different values of the interaction u. What we call optimal sweep rate provides the most coherent outcome (minimum dispersion). Contrary to the traditional dogma, it is not true that "slower is better". Adiabatic evolution. It is illuminating to discuss the dependence of the energy landscape using a quantum "energy level" language. The parametric evolution of the many body eigen-energies is presented in Fig. 4a. If the system were completely chaotic, then we could associate each E n with a micro-canonical energy surface that encloses a phase space volume  The interaction is u = 2.3 , and the associated vertical lines are from left to right � mts = π and � stb = 1.26π and � dyn = (3/2)π and � swp = 1.62π . The units of time have been chosen such that K = 1 . Each row is for a different sweep rate. From up to down we have � = 3π · 10 −4 (slow) and � = 5π · 10 −4 (optimal) and � = 3π · 10 −3 (faster).
Irrespective of chaos, a practical numerical procedure to find the phase-space volume is to invert the dependence E = E n where n = 1, 2, 3 . . . . The validity of this statement is implied by the Wigner-Wyle formalism. The representative E n curves in the background of Fig. 2 have been calculated using this procedure with N = 30.
For an adiabatic sweep, the phase-space volume Eq. (11) is the so-called adiabatic invariant 2-4 . This statement assumes a globally chaotic energy surface. In the classical context we say that during an adiabatic sweep the system stays in the same adiabatic energy surface. In the quantum context we say that the system stays in the same adiabatic energy level.
In a strictly quantum-adiabatic scenario, the system stays in its ground state with energy E GS (�) , and therefore the population is fully depleted from k = 0 to the other orbitals. Such quantum adiabaticity cannot be observed for a realistic sweep rate, because it requires many-body tunneling from a metastable minimum of the energy landscape 62 . Consequently, for large N, the semiclassical picture provides a sound approximation. In Fig. 4b we demonstrate that even a circuit with small number of particles ( N = 30 ) follows a semiclassical-like scenario.
The semiclassical adiabatic scenario excludes the possibility of tunneling, and therefore can start only when �(t) > � stb , namely, once the central SP becomes a saddle in the energy landscape. In order to determined stb we use the Bogolyubov procedure, which brings the Hamiltonian in the vicinity of the SP to a diagonalized form: Explicit results for the Bogolyubov frequencies are provided in the Methods section. The SP becomes a saddle once the ω q do not have the same sign. This happens for larger than The topography at the vicinity of the central SP, once it becomes a saddle is as follows: it is still a minimum in the M = 0 subspace, while away from M = 0 the energy floor is lower (see the Supplementary for plots of the (M, E) energy landscape). Nevertheless, we see from the simulation of Fig. 2 that spreading away from the central SP starts only at a later stage, whereafter the cloud departs the E 0 curve, neither follows any of the E n curves.
Quench-related spreading. Let us consider first the simpler scenario of preparing a cloud at n = 0 , which is the = 0 ground state, and then evolving it with H(� = 0) , aka a quench process. The SP for t > 0 (after the quench) is dynamically unstable if the Bogolyubov frequencies become complex. This happens (see Methods) for larger than After the quench the cloud spreads away from n = 0 in the landscape that is described by Fig. 5a, as illustrated in Fig. 5b. The Poincare section there shows that the stability island is taken-over by a chaotic strip. The points  Sweep-related spreading. We now consider again a quasi-static sweep process. Naively, we might expect that spreading will start once dyn is crossed. But a more careful inspections reveals that the QS limit is subtle. We see from the upper left panel of Fig. 2, and from Fig. 5c that for a slow sweep the cloud splits into two pieces. The dynamics is caricatured in Fig. 6. The reason for the splitting is related to the co-existence of two different mechanisms. One resembles the quench scenario. Namely, somewhere in the range [ stb , dyn ] spreading is initiated along the chaotic strip. But a different spreading mechanisms comes into play after dyn is crossed. This second mechanism dominates the "optimal sweep" of Fig. 2. For an optimal sweep the chaos-related spreading mechanism has no time to develop. The additional sweep-related mechanism is not related to chaos, but to the bifurcation of the stability island. It obeys the Kruskal-Neishtadt-Henrard theorem [13][14][15][16][17][18][19][20][21][22][23] , namely, the cloud is drained into the emerging stability island. The full optimal sweep scenario is displayed in the left panels of Fig. 7. The sweep rate is � = 3π · 10 −4 (slow). The snapshot is taken at ∼ stb . An inner piece of the cloud is still locked in the tiny n = 0 stability island, and therefore has energy close to E 0 . An outer piece of the cloud was formed due to very slow spreading in the chaotic corridor, and therefore has lower energy. The Poincare section at the background is adjusted to this lower energy. Lower inset (blue points): The further evolution of the same cloud after we stop the sweep at = stb and wait to see further ergodization in the chaotic strip. The upper inset would look like that if the sweep were much slower. Main panel: Zoom that displays the red and the blue clouds of the insets. www.nature.com/scientificreports/ Depletion process. As we already observed in Fig. 2, the spreading of the cloud starts before or latest at dyn . But looking at the color-code we see that the depletion happens at a distinct moment when �(t) ∼ � swp . This is the moment when a corridor connects the central SP n = 0 with the peripheral region n = N/2 . In the absence of chaos n = N/2 is formally an SP of the H (0) (ϕ, n; M = 0) Hamiltonian. Each SP has its own separatrix. For = swp the two SPs have the same energy, and therefore the two separatrices coalesce. From the equation E 0 = E ∞ (0) we get Once we add the H ± terms, this joint separatrix becomes a chaotic strip, what we call "corridor". The corridor is available for a small range of around ∼ swp . During the time interval that the corridor is opened, the central SP is depleted. Both the energy landscape and the evolution are demonstrated in Fig. 7.

Subsequent evolution.
We already pointed out that strict classical adiabaticity in the QS sense of Kubo does not hold for our scenario: for �(t) > � swp the system does not follow any of the E n curves. The reason for that is figured out by further inspection of the dynamics. For �(t) > � swp the chaotic strip decomposes into quasi regular tori. Consequently a different adiabatic scenario takes over, that of Einstein and Landau, where adiabatic invariants are the "actions" of the tori. Each piece of the cloud is locked in a different torus, and therefore we do not observe in Fig. 7 further ergodization in the M direction.
Quasi static average. Without any approximation we always have Ė = −�I� t˙ . In the Ott-Wilkinson-Kubo formulation of linear response theory 2-8 , it is assumed that for a QS process the instantaneous average can be replaced by an evolving microcanonical average I E due to quasi-ergodicity. But we are not dealing with a globally chaotic energy surface. Rather, the cloud occupies at any moment only a fraction of the energy shall, or a set tori that depart from the microcanonical shell. We use the notation I QS for the corresponding average. Accordingly For a system with 2 freedoms the QS average is well defined: at any moment the ergodic region that is accessible for the evolving cloud is bounded by KAM surfaces. This is not true if we had more than 2 freedoms: then the accessible region would likely exhibit a more complicated dependence on the rate of the sweep. Anyway, in the present context the current of Eq. (10) reflects the occupation of the orbitals, and therefore can be expressed in terms of (M, n). The expectation value of the current can be calculated for the evolving cloud of the simulation, see inset of Post-sweep ergodization. For a QS process it is expected to witness quasi-ergodic distribution at any moment. For faster sweep the cloud fails to follow the evolving energy landscape, and therefore a post-sweep ergodization stage is expected, as indeed observed in Fig. 2 for the "faster" sweep. But surprisingly post-sweep ergodization stage is also observed if the sweep rate is extremely slow, as observed in Fig. 2 for the "slow" sweep. The reason for that is explained by Fig. 7. Namely, in the case of a very slow dynamics, the cloud is split into Figure 6. The spreading mechanism. The dynamics of Fig. 5 is caricatured for an optimal sweep (left panels) and for a slow sweep (right panels). The panels are ordered by time (from top to bottom). For an optimal sweep, chaos has no time to induce spreading, therefore, even if the cloud is larger (not displayed) the spreading process looks the same. For a slow sweep the outer part of the cloud has the time to spread way from the center along the chaotic strip. This chaotic spreading is initiated in the range

Discussion
Disregarding the very well studied 2-site Bosonic Josephosn junction, the trimer is possibly the simplest building block for an atomtronic circuit. It is the smallest ring that possibly can be exploited as a SQUID-type Qubit device [29][30][31][32] . The first requirement is to have the possibility to witness a stable superflow [49][50][51] . The second requirement is to have the possibility to witness coherent operation. The latter is indicated by, say, coherent oscillations between clockwise and anti-clockwise superflow currents 32 . The third requirement is to have the possibility to execute protocols that do not spoil the coherence, meaning that the particles remain condensed in some evolving orbital 61 . In semiclassical perspective it means that an initial Gaussian cloud does not ergodize. One may say that ergodicity due to chaos, as opposed to stability, is the threat that looms over the condensation of bosons in optical lattices. Inspired by experiments with toroidal rings 60 , here we considered a lattice ring that undergoes a prototype sweep protocol: increasing from 0 to 3π such that the k = 0 orbital goes from the floor to the ceiling. During this process this orbital is depleted. The details of the process are as follows: As is increased beyond a value mts , the followed SP becomes a metastable minimum; For larger that stb it becomes a saddle in the energy landscape of the circuit; Depending on the sweep rate it can maintain dynamical stability up to some larger value dyn ; Beyond this value the SP becomes unstable, but this does not automatically implies that the coherent state is depleted; A fully developed depletion process requires a corridor that leads to ergodization within a chaotic sea; Such corridor is opened during a small interval around ∼ swp ; During the chaotic stage of the sweep we witness partial ergodization, and the final state of the system is in general not fully-coherent. An optimal sweep rate can be determined.
In a larger perspective we emphasize that the traditional view of adiabaticity is not enough in order to a address a QSTP for a system that has mixed integrable and chaotic dynamics. Some historical background is essential in order to appreciate this statement. On the one extreme we have the Einstein-Landau theory for adiabaticity for integrable systems 1 . On the other extreme we have the Kubo-Ott-Wilkinson picture of adiabaticity in chaotic systems [2][3][4][5][6][7][8] , which is associated with energy absorption in accordance with linear-response theory. But realistic systems are neither integrable nor chaotic, but rather have mixed phase space whose topological structure changes during the sweep process. The simplest scenario is separatrix crossing, that can be addressed using the Kruskal-Neishtadt-Henrard theorem [13][14][15][16][17][18][19][20][21][22][23] . More generally tori can merge into chaos, and new sets of tori can be formed later on. This leads to anomalous dissipation 9,10 and irreversibility in the QS limit 11,12 . With the same spirit we have explored in this work the mechanisms that are involved in QS transfer protocols, and also the non-trivial dependence of the outcome on the sweep rate.

Methods
The Hamiltonian. The BHH for an L-site rotating ring is where j mod(L) labels the sites of the ring, the a-s are the bosonic field operators, and is the Sagnac phase.
It is convenient to switch to momentum representation. For a clean ring the momentum orbitals have wavenumbers k = (2π/L) × integer . One defines annihilation and creation operators b k and b † k , such that b † k = 1 √ L j e ikj a † j creates bosons in the k-th momentum orbitals. Consequently the BHH takes the form where the constraint k 1 +k 2 +k 3 +k 4 = 0 mod(2π ) is indicate by the prime, and the single particle energies are Later we assume, without loss of generality, that the particles are initially condensed in the k = 0 orbital. This is not necessarily the ground-state orbital, because we keep as a free parameter. Note that we optionally use k as a dummy index to label the momentum orbitals.
Trimer hamiltonian. For the purpose of semiclassical treatment we express the Hamiltonian in terms of occupations and conjugate phases. For the 3-site ring ( L = 3 ) we get: We define q 1 = ϕ 1 − ϕ 0 and q 2 = ϕ 2 − ϕ 0 where the subscripts refers to k 1,2 = ±(2π/3) . Using the notation a † j a † j a j a j − K 2 e i(�/L) a † j+1 a j + h.c.
n 2 1 + n 2 2 + n 1 n 2 + 2U 3 (N−n 1 −n 2 ) √ n 1 n 2 cos q 1 + q 2 (23) H (+) = 2U 3 (N−n 1 −n 2 )n 1 n 2 cos q 1 − 2q 2 (24) φ[mod(4π)] = q 1 − q 2 = ϕ 1 − ϕ 2 ϕ[mod(2π)] = q 1 + q 2 = ϕ 1 + ϕ 2 − 2ϕ 0 www.nature.com/scientificreports/ For positive u and � < � stb the SP is the minimum of the energy landscape, and the Bogolyubov frequencies are positive, see Fig. 8. The SP becomes a saddle once ω − changes sign and becomes negative. The SP becomes dynamically unstable once the ω ± become complex. Note that the energy of the SP, once it becomes unstable, gets above the M = 0 floor, see Fig. S1 of the Supplementary. By inspection of Eq. (29) we can identify a critical value of the interaction u c = 9/4 . For large interaction ( u > u c ) the Bogolyubov frequencies remain complex up to the end of the sweep at � = 3π . This indicates that the SP in not at the maximum of the energy landscape, see Fig. S1. The upper most SPs in this region support self-trapped states. For weak interaction ( u < u c ) the Bogolyubov frequencies become real and negative once we cross � = 3 arccos (−(9/4)u) , indicating that the SP becomes a stable maximum.