Persistent current of atoms in a ring optical lattice

We consider a small ensemble of Bose atoms in a ring optical lattice with weak disorder. The atoms are assumed to be initially prepared in a superfluid state with nonzero quasimomentum and, hence, may carry matter current. It is found that the atomic current persists in time for a low value of the quasimomentum but decays exponentially for a high (around one quarter of the Brillouin zone) quasimomentum. The explanation is given in terms of low- and high-energy spectra of the Bose–Hubbard model, which we describe using the Bogoliubov and random matrix theories, respectively.


2
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Introduction
Ultracold atoms in optical lattices constitute an intense research activity both in experimental and theoretical physics. Up to now this system has mostly been used for modelling the fundamental Hamiltonians of solid state theory (see, [1,2], for example) where the number of particles is macroscopically large. However, the recent progress with manipulating a countable number of atoms [3,4] makes it possible to build a system of arbitrary size, ranging from microscopic to macroscopic. In this border region between microscopic and macroscopic, one has to deal with a finite number of atoms which, on the one hand, is too large to use the single-particle approach but, on the other hand, is too small to justify the thermodynamic limit. In the present study, we theoretically analyse one of these problem related to superfluidity of a few (N ∼ 10) Bose atoms in a ring optical lattice [5] with a few (L ∼ 10) sites.
It should be stressed in the very beginning that, currently, there are two different definitions of superfluidity in the physics literature. One definition is based on the system's response to a phase twist. With respect to Bose atoms in a lattice this approach is discussed, in particular, in [6], and a method of how one can realize the twisted boundary conditions in a laboratory experiment is suggested in [5]. The other definition originates in the Landau criterion of superfluidity and involves a response of a superfluid flow to 'wall roughness' [7,8]. In this study, we try to reconcile both approaches. Specifically, we address the following problem. Assume that, we have N Bose atoms in a ring lattice with L sites in a superfluid state with given quasimomentum κ = 2πk/L: We are interested in the time evolution of this state (which we shall also refer to as the supercurrent state) in the presence of a weak scattering potential and atom-atom interactions. We note that for a Bose-Einstein condensate (BEC) of atoms (N 1) the problem of superfluid atomic current has been considered in a large number of papers (see [9]- [17], to cite few of them). The starting point of all these studies is the mean-field approach, which is sometimes rectified by taking into account the quantum fluctuations [15]- [17]. The mean-field theory predicts a destruction of the supercurrent as soon as the quasimomentum exceeds one quarter of the reciprocal lattice constant (κ > π/2 in the notations used). In order to justify the mean-field approach in a 1D lattice the mean number of atoms per one site should be much larger than unity. As stated above, in the present study, we focus on the opposite limit N/L ∼ 1, where the mean-field approach is not applicable. For this reason, we treat cold atoms in an optical lattice from a different viewpoint, in a sense closer to quantum optics than to condensed matter physics.
The paper essentially consists of two parts; in the first part (section 2), after a brief preliminary analysis, we report the results of numerical simulations of the system dynamics, and in the second part (sections 3 and 4), we explain the observed regimes in terms of the energy spectrum of the system. The main results are summarized in the concluding section 5.

Supercurrent dynamics
Before proceeding with numerical simulations, we shall briefly discuss possible regimes for the atomic current.
where |l are the Wannier functions, J the hopping matrix element, and V l the random scattering potential. In what follows, to be concrete, we shall consider 0 V l with J. The operator V = l V l |l l| couples the degenerate states with opposite quasimomentum, resulting in new eigenstates |κ c,s = (|κ is the Fourier transform of V l . Thus, in the course of time, an atom in a ring will periodically change its momentum to the opposite one with the frequency ∼ /hL. It is of worth stressing that this periodic dynamics is exclusively due to the finiteness of L and the assumed condition J, which means thatV couples only the degenerate states of the unperturbed system. 3 Next, we consider the multi-particle case, whereâ † l andâ l are the bosonic creation and annihilation operator,n l =â † lâ l , and U is the on-site interaction energy. Using the canonical transformation,b k = (1/ √ L) l exp(i2πkl/L)â l , it is convenient to present the Hamiltonian (3) in the form whereδ(k) = 1 if k is a multiple of L, andδ(k) = 0, otherwise. For U = 0 and = 0 the multi-particle eigenstates of the system (4) are the quasimomentum Fock states |n = |n 0 , n 1 , . . . , n L−1 , where k n k = N. Our state of interest corresponds to |κ = | . . . , 0, N k , 0, . . . , where all atoms have one and the same quasimomentum. Similar to the single-particle case, the random potential couples this state to the supercurrent state with the opposite quasimomentum |−κ = |. representation. Thus the time evolution of the state |κ is defined by the following (N + 1) × (N + 1) matrix, where E m = E κ ≡ −JN cos κ are the degenerate energies of the states |κ(m) and the next terms the transition matrix elements κ(m)|V |κ(m ) . The spectrum of the matrix (5) is equidistant with the level spacing 2|V(2k)|. Thus we have reproduced the result of the single-particle analysis, where the time evolution of the system is periodic with the frequency ∼ /hL. Now we switch on the interaction. Then the intermediate states |κ(m) acquire energy shifts E m = E m (U), which appear to be m-dependent. Using the first order perturbation theory, we obtain Due to the mismatch of the energy levels E m , the supercurrent states |κ ≡ |κ(0) and |−κ ≡ |κ(N) become effectively decoupled and, hence, the supercurrent should persist in time.
At this point, we would like to note the analogy of the problem discussed with that for a BEC in double well potential. 4 Drawing this analogy further, we can estimate the minimal U min required for stabilization of the supercurrent as For the sake of completeness, we present a derivation of the estimate (7) in the next subsection (which can be safely skipped if a reader is familiar with the subject).

Semiclassical approach
The standard method of treating the system ((5) and (6)) consists of mapping it on to an effective classical system (terms proportional to the identity matrix are omitted), followed by a semiclassical quantization, where 1/N plays the role of Planck's constant. The phase portrait of the system (8) is shown in figure 1 for g/|V | = 1 and g/|V | = 10. It is seen that when g is increased the phase portrait becomes similar to that of classical pendulum, with a separatrix separating the librational and rotational regimes. The maximal and minimal values of the classical action I along the separatrix are given by I * ≈ 1/2 ± √ |V |/2g. The quantum states, associated with I, are decoupled only if they lie above the separatrix. Then, by requiring |I * − 1/2| 1/2 and noting that |V | ∼ /L, we come to the estimate (7). Needless to say, the semiclassical approach described above requires 1/N 1 and is not accurate for small N. Nevertheless, even for N ∼ 10 the spectrum of the matrix A can be well understood in terms of the effective system (8). For the purpose of future reference, the right panel in figure 2 shows the numerical solution of the matrix eigenvalue problem for N = 7, L = 9, and |V| = 0.0168. In particular, at U = 0.2J one can identify the first four top levels with the phase trajectories below the separatrix, the next two levels with trajectories around the separatrix, and the last two almost degenerate levels with trajectories well above the separatrix.

Numerical results
The above conclusion about the persistent current relies on the applicability of a perturbative approach. Formally this means that the supercurrent state |κ , as well as the intermediate states |κ(m) , have to be approximate eigenstates of the system at = 0. This imposes the upper boundary U max on the interaction constant, which appears to depend on the quasimomentum κ. Indeed, the state |κ with the energy E κ ≈ −JN cos κ is an approximate eigenstate of the system only if U is smaller than the characteristic energy gap separating it from the other energy states, coupled to |κ by interaction. As the first guess one can set this gap to the mean level spacing, given by the inverse density of state E = 1/f(E). It is easy to show that for U/J 1 the density of states of (3) is given by the Gaussian distribution (see figure 4 below) Thus the characteristic gap for a supercurrent state, which belongs to the central part of the spectrum (i.e., for κ ∼ π/2), is essentially smaller than that for a supercurrent state with low quasimomentum κ π/2. As a consequence, U max for the supercurrent state with κ ∼ π/2 may be smaller than U min . In other words, the perturbative approach of subsection 2.1 (where we used first order perturbation theory to find corrections to the eigenenergies of states |κ(m) ) breaks down before the stabilization of the supercurrent is achieved. Figure 3 compares the dynamics of N = 7 atoms in a lattice with L = 9 sites, which were initially prepared in the supercurrent state (1) with high, κ = 6π/L, and low, κ = 2π/L, quasimomentum. The normalized mean momentum of the atoms, for κ = 2π/L but does not work for κ > 2π/L. We shall come back to this point later on in subsection 4.3.

High-energy spectrum
To get a better insight into the physics of the discussed phenomena, we shall discuss the results displayed in figure 3 in terms of the energy spectrum of the system (3). We begin with the case of a high quasimomentum which, as mentioned above, refers to the central part of the spectrum.

Spectral statistics
We have found that in the case of high initial quasimomentum a transition from oscillatory dynamics to irreversible decay is associated with the transition to chaos in the Bose-Hubbard model. Following [19], we shall monitor this transition by analysing the distribution of distances between the neighbouring levels, normalized to the mean level spacing: It should be stressed that the presence of random potential in the Hamiltonian (3) alone does not yet induce chaos in the system. The only consequence of a weak disorder (relevant to the spectral statistics) is that it breaks the translational symmetry and, hence, we need not worry about decomposition of the energy spectrum into the independent subsets (labelled, in the absence of a random potential, by total quasimomentum of the atoms [19]).
The results of the statistical analysis of the high-energy spectrum are presented in figure 4.
which is typical for a generic integrable system, and the Wigner-Dyson statistics, typical for non-integrable systems. These distributions reflect the different character of the parametric dependence of the energy levels E n = E n (λ) on some parameter in the Hamiltonian (λ = U in our case). Namely, in the integrable case the energy levels may cross and, hence, one finds an arbitrary small s. On the contrary, if the system is non-integrable, the energy levels show avoided crossings and the probability of finding small s tends to zero.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Panel (a) in figure 4 shows the density of states f(E) for U = 0.02, where only the data from the central part of the spectrum (marked by the inverse parabola) were used for the statistical analysis. It is seen in the lower panel that for U = 0.02 the level spacing distribution follows the Poisson statistics. Thus for this value of the interaction constant the system should be classified as integrable, which is consistent with the periodic dynamics of the mean momentum in figure 3(a). Panel (b) in figure 4 shows the density of states for U = 0.2. Apart from an uniform shift of the spectrum to positive values, no qualitative change in f(E) is observed. However, we do observe a qualitative change in the level spacing distribution. Now it reliably follows the Wigner-Dyson statistics which, as mentioned above, is a hallmark of quantum chaos.

Local density of states
The spectral statistics is only one (and, in fact, a rather poor) characteristic of the system. In particular, the level spacing distribution remains unchanged (Wigner-Dyson) in the interval 0.2 U 1, although the decay rate of the supercurrent changes with U. One gets more information about the system by studying its eigenfunctions. To this end we introduce a quantity R(m, n), closely related to the so-called local density of states. 5 In equation (12), | n (U) are the eigenfunctions of the Hamiltonian (3) calculated for a given U and ordered according to their energies. In what follows we shall fix U = 0.02, while U will be scanned in the interval 0.2 U 1. Since for U = 0.02 the system is integrable, the matrix (12) can be alternatively viewed as the matrix of the expansion coefficients of the chaotic states | n (U) over 'regular basis' |m = | m (U = 0.02) . The characteristic structure of the matrix (12) is shown in figure 5 for U = 0.2. It is seen that that R is a banded matrix with strongly fluctuating matrix elements. The mean values of the elements across the main diagonal, are shown in figure 6 on linear and logarithmic scales. (Here, as in the spectrum analysis, we consider an energy window of the order of unity in the central part of the spectrum.) It is seen thatR( m) converges to the Lorentzian, 6 We note, in passing, that a similar result is reported in the recent paper [21] devoted to the spectral properties of the three-site Bose-Hubbard model. The distribution (14), also known as the Breit-Wigner formula, implies the exponential decay of the supercurrent state. Indeed, considering the overlap integral κ|κ(t) , one has κ| exp − ihĤ t |κ = m,m ,n κ|m m| n exp − ih E n t n |m m |κ where we substitute the exact energy levels E n by their approximate positions, E n ≈ E κ + E m. (This approximation obviously holds till time th∼h/ E =hf(E k ),which increases exponentially with the system size.) Substituting hereR( m) from equation (14) we have κ|κ(t) = exp(− Et/h). We found that the width grows approximately quadratically with U in the interval 0.2 U 1.

Low-energy spectrum
We turn to the case of low quasimomentum. For small κ the energy of the supercurrent state (1) falls into the low-energy tail of the density of states (9), where the random matrix approach is not applicable. On the other hand, the low-energy spectrum of the interacting Bose atoms is believed to be described by the Bogoliubov theory. For this reason, we review the Bogoliubov approach for a finite size system. Through the section, if not stated otherwise, we assume the homogeneous case = 0.

Bogoliubov approach
As an intermediate step, let us show that the Bogoliubov approach amounts to the following two assumptions. (i) The low energy eigenstates of interacting Bose atoms are given by a linear superposition of the quasimomentum Fock states, where n atoms have quasimomentum κ, n quasimomentum −κ, and the rest N − 2n have zero quasimomentum [18], i.e., The analysis goes as follows. Substituting the wavefunction (15) in the eigenvalue equation with the Hamiltonian (4), we get a system of linear equations for the coefficients c n , 2E k n + U L a n c n + U L (b n c n−1 + b n+1 c n+1 ) = Ec n ,
Next, introducing the generating function, we present the system of linear equations (17) as a differential equation on the function (θ), 7 wheren = −i∂/∂θ and ε = E k /g. The general solution of (18) reads Finally, requiring (θ + 2π) = (θ) and calculating the relevant integral, 1 2π 2π 0 dθ 2 cos θ + 2 + ε = 1 we get the equidistant spectrum with the transition frequency The result (20) reproduces the famous Bogoliubov equation for the quasiparticle excitations of the Bogoliubov vacuum.

Bogoliubov spectrum
In the previous subsection, we have considered an excitation of the given quasimomentum state, with the single-particle excitation energy E k = J (1 − cos κ). To include the other quasimomentum states, the ansatz (15) should be generalized to | = n c n |N − 2 k n k , n 1 , n 2 , . . . , where n = (n 1 , . . . , n L/2 ). Substituting (23) in the stationary Schrödinger equation with the Hamiltonian (4), we obtain a system of rather complex equations on the coefficients c n , which can be solved analytically only in the thermodynamic limit. In this limit, as it is easy to show, the whole eigenvalue problem factorizes to L/2 eigenvalue problems of the form (17) and, hence, the whole spectrum is given by the direct sum of L/2 linear spectra. A remark about the total quasimomentum, which is a global symmetry of the system in the absence of random potential, is in turn. The substitution (15) corresponds to zero total quasimomentum. To get nonzero values of the total quasimomentum, one should use a slightly different ansatz, Ansatz (22) leads to the eigenvalue equation of the form (16) but with different coefficients a n and b n . In particular, considering the thermodynamic limit, equation (17) changes to 2(E k + g)(n + m)c n + g n(n + m)c n−1 + g (n + 1)(n + m + 1)c n+1 = Ec n .
We note, in passing, that if the spectra associated with different single-particle excitation energy and different total quasimomentum are superimposed, one typically finds a multiple degeneracy of the levels at U = 0 (see figure 8 below). It is interesting to compare the discussed Bogoliubov spectrum of an infinite system with the low-energy spectrum of a finite system. For this reason we calculate numerically the spectrum of N = 25 atoms in a lattice with L = 5 sites. 8 A few first levels of this system are depicted in the left panel of figure 7 where, to facilitate a comparison, we subtract the energy E 0 = −JN + UN(N − 1)/2L and rescale energy axis on the basis of the frequency ω 1 = ω 1 (U). The Bogoliubov spectrum, calculated by using equations (23) and (24), is depicted in the right panel of figure 7. Both similarities and differences are evident. The first two levels are seen to coincide in the whole interval 0 U 1. On the other hand, multiple degeneracy of the levels around U = 0.13 in the right panel is removed in the left panel. This is, in fact, not surprising. Indeed, let us consider the lowest group of levels, showing the degeneracy. These levels are associated with the quasimomentum Fock states |0, 2, 21, 2, 0 , |0, 2, 22, 0, 1 , |1, 0, 22, 2, 0 , and |1, 0, 23, 0, 1 . (Here we use a different notation for the Fock states, corresponding to the Brillouin zone −π < κ π.) These states belong to different spectra, labelled by k and m in equation (24), and are decoupled within the Bogolubov approach. However, for the considered finite system these states are coupled by interaction, where the coupling matrix elements are of the order of g/N.

Persistent current
In this subsection, we critically review the result of section 2 about the persistent current, carefully checking validity of the perturbative approach. The left panel in figure 8 shows the low-energy levels for N = 7 and L = 9 (the whole spectrum is shown, i.e., no symmetry selection according to the total quasimomentum). Remarkably, even for such a small number of atoms one still has a qualitative agreement with the Bogoliubov spectrum. Our states of interest in figure 8 are the supercurrent and the intermediate states |κ(m) , which originate from the point marked by an asterisk. It is seen, by comparing with figure 2(a), that (i) the splitting between these levels matches well equation (6) and (ii) the coupling of these states to the other states of the system is negligible, which is indicated by the absence of the avoided crossings. In addition to the case = 0, the right panel in figure 8 shows the spectrum of the atoms in the presence of a weak scattering potential (should be compared with figure 2(b)). Again, no avoided crossings with the other levels are seen. Hence, the approach of section 2 is well justified.
The above visual analysis of the spectrum can be made quantitative by considering the overlap of the supercurrent state |κ with the exact eigenstates, For κ = 0 the quantity (25) is obviously maximized by the ground state | 0 (U) , which is expected to coincide with the Bogoliubov state ground. The solid line in figure 9 shows the overlap of the state |κ = 0 with the ground Bogoliubov state. A monotonic decrease of Q = Q(U), seen in the figure, is due to population of the single-particle quasimomentum states with κ = 0, and is often referred to as the quantum or Bogoliubov depletion of the BEC [18,22]. Additionally, the dashed and dash-dotted lines in figure 9 depict the overlap of the states |κ = 2π/L and |κ = 4π/L with the Bogoliubov states, originating from these supercurrent states, which we calculate by substituting the single-particle excitation energy E k = J(1 − cos(2πk/L)) in equation (17) by E k = 0.5J[cos(κ + 2πk/L) + cos(κ − 2πk/L) − 2 cos κ] ∼ (2πk/L) 2 cos κ.
Finally, the series of dots correspond to the quantity (25). It is seen that for κ = 0 dots perfectly follow the solid continuous line. Thus the ground state of the system is indeed well approximated by the Bogoliubov state. With exception of two narrow avoided crossings this is also the case for the state of our interest κ = 2π/L. However, for higher initial quasimomentum κ = 4π/L, the Bogoliubov state is seen to be completely destroyed by the large number of avoided crossings.

Conclusions
Within the formalism of the Bose-Hubbard model, we have considered time evolution of the atomic supercurrent in a ring optical lattice with weak on-site disorder. For vanishing atomatom interactions, weak disorder induces Rabi oscillations of the atomic current, where the atoms periodically change their velocity to the opposite one. For non-vanishing atom-atom interactions, the supercurrent dynamics depend crucially on the initial quasimomentum κ (i.e., the initial velocity of the atoms). Namely, for a high quasimomentum κ ∼ π/2 the supercurrent exponentially decays as the interaction constant U exceeds some critical value, while for a low quasimomentum κ π/2 the oscillatory behaviour of the supercurrent changes to a persistent current.
The explanation for these effects is found in the structure of low-and high-energy spectra of the Bose-Hubbard model. It is shown that the low-energy spectrum of the system is regular, and the positions of the energy levels can be found by using a Bogoliubov approach. In contrast, the high-energy spectrum shows a transition from regular to a chaotic one if U exceeds its critical value. Using the results of the random matrix theory, we show that this transition is reflected in the exponential decay of the supercurrent with the decay constant proportional to U 2 .