Blocked populations in ring-shaped optical lattices

We study a special dynamical regime of a Bose-Einstein condensate in a ring-shaped lattice where the populations in each site remain constant during the time evolution. The states in this regime are characterized by equal occupation numbers in alternate wells and non-trivial phases, while the phase differences between neighboring sites evolve in time yielding persistent currents that oscillate around the lattice. We show that the velocity circulation around the ring lattice alternates between two values determined by the number of wells and with a specific time period that is only driven by the onsite interaction energy parameter. In contrast to the self-trapping regime present in optical lattices, the occupation number at each site does not show any oscillation and the particle imbalance does not possess a lower bound for the phenomenon to occur. These findings are predicted with a multimode model and confirmed by full three-dimensional Gross-Pitaevskii simulations using an effective onsite interaction energy parameter.

We study a special dynamical regime of a Bose-Einstein condensate in a ring-shaped lattice where the populations in each site remain constant during the time evolution. The states in this regime are characterized by equal occupation numbers in alternate wells and non-trivial phases, while the phase differences between neighboring sites evolve in time yielding persistent currents that oscillate around the lattice. We show that the velocity circulation around the ring lattice alternates between two values determined by the number of wells and with a specific time period that is only driven by the onsite interaction energy parameter. In contrast to the self-trapping regime present in optical lattices, the occupation number at each site does not show any oscillation and the particle imbalance does not possess a lower bound for the phenomenon to occur. These findings are predicted with a multimode model and confirmed by full three-dimensional Gross-Pitaevskii simulations using an effective onsite interaction energy parameter.

I. INTRODUCTION
The self-trapping phenomenon has been extensively studied in double-well systems by means of a two-mode model [1][2][3][4][5][6][7][8], and experimentally observed by Albiez et al. [9][10][11]. In this regime the population in one site remains higher than the one in the other well over all the evolution. This imbalance of particles performs oscillations around the non vanishing mean value, whereas the phase difference between the sites exhibits a running phase behavior. Theoretical studies of this phenomenon has been also carried out by several authors in extended regular lattices [12][13][14][15]. More recently, the study of self-trapping has also been addressed in ringshaped optical lattices [16][17][18]. Such works treat threeand four-well systems. In Refs. [16,17] the dynamics has been investigated through a multimode (M) model that utilized ad-hoc values for the hopping and onsite energy parameters. Whereas in Ref. [18], such parameters have been extracted from a mean-field approach using threedimensional localized "Wannier-like" (WL) onsite functions and including an effective onsite interaction energy parameter [19,20]. For large filling numbers, the inclusion of such a realistic interaction parameter has been shown to be crucial for the accurate description of the dynamics, yielding a sizable change on the time periods respect to those obtained by the standard model. In contrast, for filling number around unity mean-field approaches are not applicable and hence other microscopic methods have to be used [21,22]. In Ref. [23], it has been demonstrated that an effective interaction can also be extracted from the Bogoliubov excitations in the case of the Josephson regime. A systematic study of the selftrapping regime and the crossover to the Josephson oscillations in four-well systems including non-symmetric configurations has been developed in Ref. [18]. It is worthwhile noticing that the dynamics in multiple well condensates constitutes a promising area provided that successful efforts have been performed to experimentally construct ring-shaped optical lattices [24].
In this work we demonstrate theoretically the existence of a dynamical regime that exhibits a novel behavior. If the number of wells of the lattice is a multiple of four, there exists a family of nonstationary states with constant site populations and special non-trivial phases. These states could be regarded as a special variation of a ST regime where, in contrast to that observed in twoand multiple-well condensates [1,2,7,18], the population imbalance between neighboring sites can be arbitrarily low and do not exhibit any oscillation in time. For such states the M model order parameter can be expressed as a linear combination of particular degenerate Gross-Pitaevskii (GP) stationary states. However, due to the nonlinear nature of the GP equation, the states are nonstationary. The dynamics of these states is governed only by the onsite interaction energy parameter. We explicitly show that the angular momentum exhibits a simple oscillating behavior and that the velocity circulation around the ring alternates periodically between values −N c /4 and N c /4, being N c the number of weakly linked condensates. A goal of this work is to obtain an analytical expression for such a time period which involves only the imbalance and the effective interaction parameter. By comparing the evolution of the phase differences obtained through GP simulations for a four-well system and with the M model, we can establish the accuracy of such a parameter. The existence of these states is confirmed numerically by means of full three-dimensional Gross-Pitaevskii (GP) simulations showing a perfect accordance to the M model predictions for several population imbalances. Furthermore, a Floquet stability analysis confirms that for the imbalances studied here the dynamics turns out to be regular.
The paper is organized as follows. In Sec. II we briefly review the main concepts of the multimode model. In particular, we rewrite the equations of motion which include an effective onsite interaction parameter [18] and we outline the construction of the localized states in terms of the GP stationary ones. In Sec. III we describe the specific four-well system used in the numerical simulations. Section IV is devoted to study the properties of these states with blocked occupation numbers. As a first step we introduce a continuous family of states corresponding to fixed points of the M model in the phase diagram defined by the populations and phase differences. On the other hand, we demonstrate that they turn out to be quasi-stationary solutions of the GP equation. Such states are defined with a particular combination of phases which give rise the non-stationary blocked-occupation-number (BON) states. Secondly, we show that these BON states describe closed orbits in the phase diagram whose time period is solely determined by the onsite interaction energy. By performing GP numerical simulations with a four-well potential we analyze the hidden dynamics which includes variations of density in the interwell regions, oscillations of the velocity field circulation, and an active vortex dynamics. We end this section with a study of the Floquet stability of the BON states and a proposed experimental test. In Sec. V we show how to generalize the previous results for systems with larger number of sites. To conclude, a summary of our work is presented in Sec. VI and the definition of the parameters employed in the equations of motion are gathered in the Appendix.

II. MULTIMODE MODEL
The equations of motion of the multimode model has been previously studied both for multiple-well systems in general [20,25] and also in the case of a four-well system [16,18]. Here, we only review its main ingredients, focusing in the definition of their localized states extracted from the stationary solutions of the GP equations.

A. Multimode model equations of motion including interaction-driven corrections
Using the multimode model order parameter, written in terms of three-dimensional WL functions localized at the k-site, w k (r) [18], one obtains the equations of motion for the time dependent coefficients b k (t) = e iφ k |b k |, by replacing the order parameter in the time dependent GP equation. The 2N c real equations, written in terms of the populations n k = |b k | 2 = N k /N and phase differences ϕ k = φ k − φ k−1 including effective onsite interaction corrections [18], are where U eff = f 3D U . The definitions of the tunneling parameters J and F , and of the onsite interaction energy parameter U are given in the Appendix. The coefficient f 3D = 1 − α is obtained from the slope of the onsite interaction energy as function of ∆N k −N/N c . As shown in Refs. [18,20], the introduction of f 3D is crucial for obtaining an accurate dynamics. From this system of equations only 2N c −2 are independent since the variables must fulfill k n k = 1 and k ϕ k = 0.

B. Localized states
In previous works, we have described in detail the method for obtaining the localized states in terms of GP stationary states [18,20,25]. Summarizing, first the stationary states ψ n (r, θ, z) are obtained as the numerical solutions of the three-dimensional GP equation [26] with different winding numbers n, with n restricted to the values −[(N c − 1)/2] ≤ n ≤ [N c /2] [27] for large barrier heights [25]. Since the ψ n are orthogonal for different n [20,25], one can define orthogonal WL functions localized on the k-site by the following expression: A discussion of how to choose the global phases of ψ n (r, θ, z) in order to achieve the maximum localization of w k is given in Ref. [18].
In its turn, the stationary wavefunctions can be written in terms of the localized WL wavefunctions in Eq. (4) as For the four-well problem, N c = 4, the states with n = ±1 are degenerate and can be regarded as vortexantivortex states since ψ ±1 = 1 2 k w k e ±i π 2 k have opposite circulation. It is worthwhile remarking that as the GP equation is nonlinear, linear combinations of degenerate stationary states, as , e.g., the vortex-antivortex states, are in general nonstationary.

III. THE SYSTEM, TRAPPING POTENTIAL AND PARAMETERS
Although the states investigated in this work also exist for larger number of wells, in our numerical simulations we will consider a four-well ring-shaped trapping potential given by where r 2 = x 2 + y 2 and m is the atom mass. The harmonic frequencies are given by ω r = 2π × 70 Hz and ω z = 2π×90 Hz, and the lattice parameter is q 0 = 5.1µm. Hereafter, time and energy will be given in units of ω −1 r and ω r , respectively. The length will be given in units of the radial oscillator length l r = /(mω r ) 1.3 µm. We also fix the barrier height parameter at V b = 25 ω r and the number of particles to N = 10 4 .
For a system of Rubidium atoms in the above configuration we have obtained the following multimode parameters, the hopping J = −6.60 × 10 −4 ω r , the interaction driven hopping parameter F = 2.08×10 −3 ω r , the onsite interaction energy U = 3.16 × 10 −3 ω r , and the effective onsite interaction energy U eff = 2.27 × 10 −3 ω r , being α = 0.28. We will numerically solve the GP equation on a grid of up to 512 × 512 × 256 points and using a second-order split-step Fourier method for the dynamics with a time step of ∆t = 10 −4 ω −1 r . For more details see Ref. [18].

IV. THE STATES
In this section we will first analyze a set of stationary points of the M model with equally populated sites whose associated order parameters are in general not exact GP stationary states. These states shall be called peculiar. In a second step, we shall show that for states with conveniently chosen initial occupation numbers and the same distribution of initial phases as the peculiar states, the populations remain blocked during all the evolution. The properties of such BON states shall be studied next. A.
Peculiar stationary states The GP stationary states used for constructing the multimode model give rise to stationary points in the M model. However, in addition to these standard points, we found a peculiar set of stationary points in a condensate with N c = 4l sites. These states are defined by |b k | = 1/ √ N c and the following local phases: φ 0 = 0, φ 1 = f 0 − π, φ 2 = π, and φ −1 = f 0 for a four-well trap (l = 1). Whereas for larger l, the sequence of phases is repeated l times along the ring. We refer to these states as peculiar because f 0 could take any value, so that instead of having isolated points in the phase diagram we have a continuous family of stationary points parametrized by f 0 . This family contains the two stationary points f 0 = ±π/2 which correspond to singly-quantized vortex states, namely, GP stationary states with winding numbers ±1. In Fig. 1 (a) a scheme of the trap and the condensate is depicted qualitatively showing states with different populations, and in Fig. 1 (b) the localized WL function in the z = 0 plane are shown together with the peculiar initial phases.
We further investigate if the peculiar order parameter could also be another stationary solution of the GP equation [26], where µ is the chemical potential and g = 4π 2 a/m is the interaction strength among atoms with a being their s-wave scattering length.
The normalized-to-unity order parameter associated to the peculiar points reads, (8) which in terms of GP stationary states can be written as The peculiar states are therefore a superposition of vortex states with opposite circulation. Since the states ψ 1 (r) and ψ −1 (r) have the same chemical potential µ 1 = µ −1 and verify ψ 1 (r) = ψ * −1 (r), applying the GP equation (7) to ψ M we obtain, where in addition, we have that (11) is almost vanishing if the WL functions are well localized as it is in the present case. Therefore, these peculiar stationary points can be regarded as quasi-stationary solutions of the GP equation.
For the particular case of f 0 = ±π/2, the second term of the right hand side of Eq. (10) vanishes and thus the order parameter is an exact solution of the GP equation. Otherwise, a general value of f 0 generates an entire continuous family of states that shows a collective motion independent of time. The ψ ±1 stationary solutions can be regarded as particular cases of Eq. (9) with maximum angular momentum. On the other hand, for which is real, hence its angular momentum is zero. Nevertheless an active vortex dynamics is present, due to the nonzero circulation of ψ ±1 . The same holds for f 0 = π.
We want to remark that such a family of stationary solutions of the M model does not necessarily exist in ring lattices with an arbitrary number of wells as it can be straightforwardly deduced from the dynamical equations (2) and (3). For example, for N c = 3 even though the degenerate stationary states with winding numbers n = ±1 are also present, the corresponding stationary points in the phase diagram only exist as isolated points.

B. Nonstationary BON states
When the numbers of particles of alternate sites are equal and the phases maintain their peculiar relation: φ 0 = 0, φ 1 = f 0 − π, φ 2 = π, and φ −1 = f 0 , the site populations do not evolve. This condition gives rise to very special dynamical states where f 0 becomes time dependent. This is shown in Fig. 2 where we compare the evolution of the occupation numbers and phase differences using full three-dimensional GP simulations with the dynamics arising from the M model. f 0 (0) = −π. From Fig. 2 it may be seen that in all cases the number of particles in each site remains fixed in time, whereas their phase differences evolve faster for larger imbalances. We have further investigated the GPE dynamics for imbalances up to ∆N = 3 × 10 3 and verified that the populations remain constant within 0.1% accuracy.
These family of states bear some resemblance to selftrapped states; however, there are many important differences compared to the well-known ST dynamics in double well potentials. First of all, the population imbalance can be arbitrarily small. Instead, to reach these states, it is only necessary to achieve the peculiar phases described above, being f 0 (0) an arbitrary value (given that n k = n k+2 ). Second, the hopping parameters J and F play no role in the dynamics, hence we cannot associate the emergence of the BON states to the small enough tunneling energy splitting like in the ST regime in two wells [9].
The BON state normalized to unity reads, which written in terms of GP stationary states yields, In order to obtain f 0 (t) one can rewrite Eqs. (3) and extract the evolution of f 0 (t) from φ k = (n k−1 −n k )N U eff . This yields which in turns completely defines all the phase differences at any time within the M model. In particular, we have ϕ 0 (t) = −f 0 (t) and ϕ 1 (t) = f 0 (t) − π. We note that for n 0 = n 1 the state (13) coincides with the peculiar state (9). It is important to point out here, that the perfect agreement between the results of GP equation and the M model observed in Fig. 2 is due to the proper definition of the onsite interaction energy parameter U eff = 2.27 × 10 −3 ω r [19,20], instead of using the bare value U = 3.16 × 10 −3 ω r . To illustrate such a difference, we included in Fig. 2 the evolution of ϕ 1 (t) for ∆N = 100 using the bare parameter U . Hence, we confirm the accuracy on the calculation of the effective onsite interaction energy parameter also in this dynamical regime.
To conclude we note that using Eq. (14) one can obtain the time period for the phase differences, which turns to be also the time period of the persistent and collective oscillation around the ring.

Angular momentum
An additional evidence of this dynamical regime is reflected in the time evolution of other observables. In particular, we shall show that within the M model the angular momentum exhibits a sinusoidal behavior as a function of time with a period T M . The expectation value per particle of a general observable O, considering an arbitrary state in the M model is given by Since each well is equivalent to all others unless a discrete rotation, we have w k | O|w k = w 0 | O|w 0 and w k+1 | O|w k = w 1 | O|w 0 for all k. Hence, the expectation value becomes For the z-component of the angular momentum we have O = L z = −i ∂ ∂θ and then taking into account that the localized states can be chosen as real functions, one obtains the expectation value of angular momentum In a BON state with initial condition f 0 (0) = −π, as sin ϕ k = − sin f 0 (t) for every k, we can write where the bracket involving the localized states is a negative number. As expected, the period of this sinusoidal function is T M . Furthermore, one can see that the stationary state Eq. (8), corresponding to ∆N = 0, yields a constant angular momentum proportional to sin(f 0 ).

Underlying dynamics
Although the population in each well remains completely fixed, the order parameter evolves in time and exhibit spatial oscillations. In order to analyze such a dynamics we first investigate the evolution of the density profile. Using the BON state expression given by Eq. (12), the evolution of the density ρ M (r, t) = |ψ M (r, t)| 2 within the M model is given by with ḟ 0 (t) = U eff ∆N (cf. Eq. (14)). Equation (20) implies that ρ M (r, t) is approximately stationary within each well, where the overlap between the WL functions of neighboring sites is negligible. Whereas the density variations are confined to the inter-well regions or junctions where the localized states do overlap. Moreover, one can infer the change of sign of ∂ρ M (r, t)/∂t at the junctions by analyzing Eq. (20). One can thus conclude that particles oscillate across both junctions of a given site without changing its net population. Furthermore, the maximum and minimum density variations during the evolution occur at times t M and t m when f 0 (t M ) = 0 and f 0 (t m ) = π, respectively. The sense of the particles flow across the junctions can be read off from the spatial profiles of the phases φ(r, t) of the wavefunction. In Fig. 3 we show snapshots of these phases at several times in the z = 0 plane obtained from both ψ M (r, t) and ψ GP (r, t) in the left and right panels, respectively. In both cases we have subtracted a global phase φ 0 (t) from φ(r, t) in order to better observe the dynamics. The initial condition is ∆N = 200 and f 0 = −π, which yield a multimode period T M 13.8ω −1 r , in sharp contrast to 2π /(U ∆N ) = 9.95ω −1 r that would be obtained with the bare onsite interaction.
In the left column, from top to bottom, we show the phase φ M (r, t) obtained from the order parameter ψ M (r, t) for the aforementioned configurations at several times: a) at t = T M /4 (f 0 = −π/2), there is a π/2 difference between neighboring sites and the velocity field corresponds to that of a vortex with a phase gradient in the counterclockwise direction, b) at t = T M /2 (f 0 = 0), the phase difference between the right and left sites is π, which corresponds to a vanishing velocity field, c) at t = 3T M /4 (f 0 = π/2), there is a −π/2 difference between neighboring sites, being the circulation clockwise as for an antivortex. Finally, d) at t = T M (f 0 = π), there is a π phase difference between the top and bottom sites. In the figure we have marked with a plus symbol and with an empty circle the presence of a vortex and an antivortex, respectively. For the M model, it can be seen that, in the left column of Fig. 3, there exists a vortex and an antivortex at the origin for the configurations: a) t = T M /4 and c) t = 3/4T M , in agreement with distributions described above. In the model, the vortex (antivortex) remains fixed at x = 0, y = 0 during the interval 0 < t < T M /2 (T M /2 < t < T M ).
On the other hand, on the right column of Fig. 3 we show phase snapshots obtained from full 3D GP simulations for times near the four different situations previously discussed. We have observed that the GP evolutions incorporate additional fluctuations and hence the velocity circulation does not change exactly at quarters of the period T M . Moreover, the velocity field never vanishes, as the change of its circulation is associated with a passage of vortices instead of with the appearance of a nodal surface [6]. Nevertheless, as shown together in Fig. 3, the order parameter from the M model is able to capture rather accurately the spatial distribution of phases present in the exact GP dynamics.
In particular, from top to bottom in Fig. 3, we show the results for the GP times: t a = 3.4ω −1 r , t b = 6.9ω −1 r , t c = 10.3ω −1 r , and t d = 13.7ω −1 r . In each site we indicate the value of the local phase φ GP (r, t) evaluated at the center of the corresponding well to be compared with that obtained in the M model. It may be confirmed that at every time the phase difference between alternated sites is always π as predicted by the model. It becomes clear from the change of sign in the phase differences that the velocity field is inverted near each half period, when the extreme variations in the density at the junctions are achieved. Except for some fluctuations around such a transition, in the intermediate times the total topological charge is conserved, whereas the number and the position of the vortices may change. In particular, in the third row of the right column of Fig. 3, one vortex and two antivortices are observed with a total negative charge of −1 instead of the single fixed antivortex predicted by the M model.
It is worthwhile to recall that the velocity circulation is quantized along any closed curved inside the superfluid and, as established in the celebrated Helmholtz-Kelvin theorem [28], it is conserved during the evolution if the superfluid condition is not broken [29]. As a consequence, the value of the circulation can only change when a vortex passes through the curve (phase slip) or when the density goes to zero.
Although both the GP equation and the M model must obey the Helmholtz-Kelvin theorem, the order parameter given by the multimode model cannot predict the motion of vortices or the generation of vortex-antivortex pairs, hence the change of the velocity field circulation could be only provided through the appearance of nodal surfaces. The nodal surfaces arise when the minimum in the local density is achieved, i.e., at f 0 = 0, π. For example, at f 0 = 0 the order parameter in Eq. (12) reduces to which corresponds to the configuration b) on the left column of Fig. 3. If all the populations were equal this condition would lead to the x = 0 plane. In our case, the deviation from a plane is due to the difference in the populations. The intersection of the nodal surface with the plane z = 0 can be viewed in the graph by the sharp π change of the phase where the density goes to zero. Similarly, one can obtain the nodal surfaces for f 0 = π, which corresponds to the configuration d). In this case the curve where the density goes to zero is around y = 0.
In contrast with the M model, the change of the velocity circulation in the GP frame is produced by the dynamics of vortices passing through the potential barriers and may include generation of vortex-antivortex pairs. In fact, we have observed that several vortex-antivortex pairs may be spontaneously generated along the barriers, thus simulating a density closer to that of the M model nodal surface. This active dynamics of vortices around the transitions is produced in a timescale much smaller than T M and hence it is not possible to access the details of the vortex motion within the present numerical precision. As an illustration, we note that the last time of the depicted GP snapshots is slightly smaller than the T M period and there still exists an antivortex around the center of the system.

Velocity field circulation
Taking into account the previous findings for the multimode model one can conclude that in one T M period the system passes through a sequence of phases that yields an alternating velocity field circulation between values 1 and −1 along a curve that connects the four wells. The transition between these two values occurs at f 0 = 0 and f 0 = π when the order parameter develops a nodal surface. In Fig. 4 we show the velocity field circulation C = v · dr, as a function of time using the M model and GP simulations. It may be seen that the same behavior is observed with both approaches. The M model is thus able to reproduce the behavior of the circulation although the details of the internal vortex dynamics is lost. In the GP dynamics the change of circulation is caused by the motion of vortices together with the creation or annihilation of vortex-antivortex pairs. Signatures of such a vortex dynamics could be observed in Fig. 3 where we have shown the phases around the transition. Another evidence of a vortex dynamics can also be visualized in the middle panel of Fig. 4, where an additional change of sign is produced near the transition.

Stability analysis
In this section we investigate the stability of the BON states by means of a Floquet analysis [18,30] of the multimode dynamical equations. This analysis is based on the characterization of the linear dynamics around its periodic orbits. In our case, the BON states are periodic solutions with constant populations n i (t) = n i and linear phase differences ϕ i (t) = ϕ 0 i + (−1) i+1 2π t/T M where ϕ 0 i fulfill the relation: The linearization of the dynamics around these states yield the non-autonomous where δ is a vector comprising both density and phasedifference fluctuations. As the BON states correspond to symmetric initial populations with peculiar phases, it is natural to consider as variables δ the departures from a symmetric case, namely, we define Given that the matrix A has a period T M , the linearized dynamics can be characterized by the so-called Monodromy matrix M which contains the change of δ after one period, i.e., M · δ(0) = δ(T M ). The matrix is built from the solutions of Eq. (22) with canonical initial conditions evaluated at T M [18,30]. In Fig. 5 we depict elements of M showing the effect of an initial population fluctuation. The orbits are regular if the perturbed system remains near the initial one after a period. This happens when M ii 1 and M i =j 0. On the contrary, when the fluctuations are enhanced (|M ij | 1), the orbits are unstable. This may be observed in Fig. 5 for low imbalances.
For the peculiar stationary states (∆N = 0), the linear system is time-independent and the problem reduces to a straightforward diagonalization of A to obtain the excitation frequenciesω corresponding to the Bogoliubov collectives modes in the case of the full GPE. The four frequencies are found to bẽ where K = 2J + F . As N U eff K, F , the most stable frequency for a given system is attained for peculiar states with f 0 = ±π/2 which in turn yield an imaginary frequencyω 2 = −(K + F ) 2 . Therefore all f 0 give rise to dynamically unstable peculiar states. The stability of stationary vortex states (f 0 = ±π/2) has been previously investigated in [31] for circular arrays of BECs, finding that only states with circulation below N c /4 are stable.

Proposed experimental test
The correct preparation of BON states requires a special sequence of phases and symmetric initial populations (n k = n k+1 ). While the common approach to experimentally measure both of them is by means of TOF and absorption images, the simple dynamics of BON states offers an alternative way to confirm its correct realization using TOF images only. Given that the relative phases among neighboring sites are revealed in the interference patterns during the TOF expansion [9], it could be verified that they obey the peculiar sequence of phases at all times. According to Eq. (14), in this case f 0 (t) must be a linear function whose slope γ relates to the particle imbalance as which might probe to be a more accurate measure than the direct estimate from absorption images. On the other hand, since Eq. (26) requires the use of U eff instead of the bare U , it may also serve to confirm its numerical value. By using U the relative error on the imbalances could be as large as of order 20-30%, depending on the number of particles [7]. Due to the experimental uncertainty, absorption images may not be able to reveal a slightly broken symmetry of the population configuration. However, the evolution of the phases will depart from linearity and will not be determined by the single function f 0 .

V. EXTENSION TO LARGER NUMBER OF WELLS
It is possible to extend the peculiar and BON states to larger number of wells provided the sequence of phases ...0, f 0 − π, π, f 0 , 0... is repeated l = N c /4 times around the ring lattice, and the populations alternate between two values, with n 2k = n 0 and n 2k+1 = n 1 . This is only possible when the number of wells are multiples of 4. Taking into account these conditions in Eq. (1) and using Eq. (4) to eliminate the WL functions w k , we can write the following BON order parameter in terms of GP stationary states, It is straightforward to show that f 0 (t) still obeys Eq. (14), and thus the corresponding time period is also given by Eq. (15). Therefore, the analysis performed in the previous section can be repeated using the same procedure, including the Floquet theory. However, for these configurations the velocity field circulation alternates between ±N c /4 and the number of nodal surfaces at each half period is equal to l. Equation (27) shows that an arbitrary linear combination of ψ ± Nc 4 leads to the BON dynamics. For example, even though for N c = 8 it is not possible to generate the BON dynamics with a linear combination of the degenerate ψ ±1 states; any linear combination of ψ ±2 will indeed give rise to a BON dynamics.
Using Eq. (18) the mean value of the z-component of the angular momentum is given by, If we let n 0 = n 1 then we obtain the most general quasi-stationary states described in section IV.A: that satisfies Since Im ψ 2 the states Eq. (29) can be regarded as quasi-stationary solutions of the GP equation when the w k are well localized functions.

VI. SUMMARY AND CONCLUDING REMARKS
We have studied a particular dynamical regime of a Bose-Einstein condensate in a ring-shaped lattice which possesses a set of states with fixed number of particles in each site and a simple dynamics in their phases. The same distribution of phases along the sites that gives rise to such nonstationary states has been shown to generate a continuous family of stationary points in the phase space of the multimode model. Such peculiar states have constant nonzero angular momentum, when all the populations are equal, and include two states that correspond to exact GP stationary solutions.
We have shown that the nonlinearity of the GP equation governs the dynamics within this regime and that it is responsible for the population blocking in the nonstationary states. In contrast to the self-trapping phenomenon this effect does not possess a lower bound for the population imbalance.
We have studied the time evolution of BON states using both the multimode model and the three-dimensional GP equation finding an excellent agreement in the populations in each site and in their phase differences. This accuracy was possible due to the inclusion of the effective interaction energy parameter instead of the bare one. Even though the multimode model was unable to account for the motion of individual vortices and the creation/annihilation of vortex-antivortex pairs, it was demonstrated that it correctly predicts the evolution of the velocity circulation and angular momentum, characterizing this regime as a persistent current oscillating around the lattice.
By performing a Floquet stability analysis of the blocked populations states, we have verified that their dynamics is regular for the particle imbalances here considered. In a four-well system these states could, in principle, be experimentally achieved by initially manipulating the position of the potential barriers in order to have different populations or by using an elliptic trap in the (x, y) plane with their axis forming a π/4 angle during a short time and then reverting the potential to a circular harmonic trap. A simple way to produce the initial distribution of phases would be to start with the same state as we have used in our numerical calculations. This could be achieved by illuminating half of the condensate (e.g., x > 0 ) with an additional laser for a period of time until it develops a π phase difference between the half spaces x > 0 and x < 0. However, any other initial distribution of phases seems feasible using a Spatial Light Modulator (SLM) [32][33][34], and hence also the whole family of stationary M-model states could be directly generated. Furthermore, given that all the phases loose their dependence on a single linear function as soon as the symmetric condition on the site populations is lifted, these states could be first tested to adjust the population in alternate wells with arbitrary imbalances. A second phase imprinting application could then be used to generate the desired state. Since the BON states present a simple analytical form for the phase difference between neighboring sites, it could also allow to measure the initial population imbalance by means on interference patterns in TOF images, rather than absorption images. Together with the calculation of these parameters by the preceding definitions we have followed the alternative method outlined in Ref. [20] which involves directly the energies of the GP stationary states. Both approaches have proven to yield values equal in less than one percent. We note that we have disregarded the parameter that involves products of neighboring densities because, for the present system, its contribution turned out to be negligible.
This work was supported by CONICET and Universidad de Buenos Aires through grants PIP 11220150100442CO and UBACyT 20020150100157BA, respectively.