Localization of weakly interacting Bose gas in quasiperiodic potential

We study the localization properties of weakly interacting Bose gas in a quasiperiodic potential commonly known as Aubry-Andr\'e model. Effect of interaction on localization is investigated by computing the `superfluid fraction' and `inverse participation ratio'. For interacting Bosons the inverse participation ratio increases very slowly after the localization transition due to `multisite localization' of the wave function. We also study the localization in Aubry-Andr\'e model using an alternative approach of classical dynamical map, where the localization is manifested by chaotic classical dynamics. For weakly interacting Bose gas, Bogoliubov quasiparticle spectrum and condensate fraction are calculated in order to study the loss of coherence with increasing disorder strength. Finally we discuss the effect of trapping potential on localization of matter wave.


I. INTRODUCTION
In recent years 'Anderson localization' of particles and waves has regained interest in quantum many particle physics. 'Anderson localization' is a remarkable quantum phenomenon for which the propagating wave becomes exponentially localized in the presence of disorder [1,2] and is a well studied subject in the context of electronic systems in presence of disorder. Like particles, waves can also localize in disordered medium. Localization of 'matter wave' has recently been observed in experiments on ultracold Bose gases in presence of speckle potential [3] and in bichromatic optical lattices [4]. For certain parameters, Bosons in bichromatic lattice can be mapped on to Aubry-André (AA) model [5] with 'quasiperiodic' potential. Also localization of light in AA model has been observed experimentally [6]. Unlike Anderson model with random disorder, the AA model exhibits localization transition in one dimension at a critical strength of the potential. Apart from that the quasiperiodicity of AA model gives rise to various interesting spectral properties [7]. In recent experiment [8] effect of interaction on AA-model has been studied. Repulsive interaction plays an important role in the formation of correlated phases like 'Bose-glass' phase [9][10][11][12]. Also the dynamics and diffusion of interacting Bose gas in presence of disorder have been investigated both experimentally [13,14] and theoretically [15,16]. In recent years 'many body localization' [17] and localization at finite temperature [18][19][20] have generated an impetus to study disordered Bose gas.
In this work we investigate localization of weakly interacting Bose gas in AA potential. Apart from calculating the ground state properties, we also study localization of wavefunction by using a classical dynamical map. The paper is organized as follows: In section II, we discuss non-interacting Bose gas in quasiperiodic potential. In subsection A, we review the AA model and discuss the single particle localization properties. Localization in non-interacting system using dynamical map approach is presented in subsection B. In section III, various physical quantities like inverse participation ratio and superfluid fraction are computed to investigate localization of weakly interacting Bose gas within meanfield theory. Effect of non-linearity due to interactions on the dynamical map is studied. Further, we compute the Bogoliubov quasiparticle energies and amplitudes to investigate quantum fluctuations. We study the localization of both Bogoliubov amplitudes and non-condensate densities.
Effect of trapping potential on localization and effect of disorder on the center of mass motion are presented in section IV. Finally we summarize our results in section V.

II. NON-INTERACTING BOSONS IN QUASIPERIODIC POTENTIAL
In the original experiment [4] localization of ultracold Bosons has been studied using bichromatic optical lattice which can be mapped on to Aubry-André model within tight binding approximation and for certain parameter regime [21]. The Aubry-André model is defined by the Hamiltonian, where |n is a Wannier state at lattice site n, J is nearest neighbour hopping strength, λ is the strength of onsite potential and period of potential is determined by β. In the rest of the paper, we would be working in the units in which J = 1.
For λ = 2, the Hamiltonian (1) is equivalent to the Harper model [22] describing the motion of an electron in a square lattice in the presence of a perpendicular magnetic field, where the flux Φ through each plaquette in units of flux quantum Φ 0 = h/e is given by the parameter β = Φ/Φ 0 . Energy spectrum of this problem gives rise to the well known 'Hofstadter butterfly' [23]. The Hamiltonian given in Eq.(1), poses very interesting properties for irrational values of β. When β is chosen to be a 'Diophantine number', AA model undergoes a localization transition at a critical value of the potential strength λ = 2 [5,24]. On the contrary, localization transition is absent in one dimensional Anderson model with random disorder. The quasiperiodic potential generates a correlated disorder in the AA model.

A. Single particle localization
In order to study localization transition in AA model, we choose β = ( √ 5 − 1)/2 which is the inverse of 'golden mean' and a Diophantine number . This is particularly helpful since rational approximation of β can be done by Fibonacci series. Although incommensurability of the potential with the underlying lattice plays an crucial role in AA model, for numerical studies one can approximate β = F n−1 /F n , where F n is nth Fibonacci number for sufficiently large value n. This rational approximation of β fixes the lattice size N s = F n to impose periodic boundary condition [25]. Duality in AA model can be shown by introducing new basis states in momentum space, The dual model is obtained by substituting Eq.2 in Eq.1, (3) It is important to note that AA model becomes self-dual at a critical coupling λ = 2, where the localization transition occurs. To obtain the eigenvalues and eigenfunctions of the Hamiltonian we expand the state vector in terms of Wannier states, |ψ = n ψ n |n , where ψ n is the wavefunction at nth lattice site. The eigenvalue equation of the Hamiltonian Eq.1 is reduced to a discrete Schrödinger equation, where is the energy eigenvalue. The degree of localization of a normalized state |ψ can be quantified by 'inverse participation ratio' (IPR) I, The wavefunction of an extremely localized particle at a site n 0 is given by ψ n = δ n,n0 , for which the IPR becomes unity. On the other hand, for the completely delocalized wavefunction ψ n = 1/ √ N s , IPR is 1/N s which vanishes in the thermodynamic limit. For AA model all energy eigenfunctions in real space are exponentially localized and IPR sharply increases to unity above the critical coupling λ = 2. Due to the duality of AA model, the localization of wavefunction in real space and in dual momentum space shows opposite behavior. In real space the wavefunctions are localized above λ = 2, wheras localization in dual momentum space occurs for λ ≤ 2. The IPR of the ground state wavefunction in real space and in dual monentum space as a function of λ are shown in Fig.1.
It is interesting to note that IPR in real and momentum space intersects at the self dual point λ = 2.
The spectrum of AA model also shows various interesting properties. For λ = 0 it has usual single band energy spectrum of periodic lattice. Addition of the potential generates quasiperiodic structure and destroys the band like dispersion. Successive rational approximation of β generates more periodicity over the underlying lattice, which in turn breaks the original energy band into many subbands and it leads to opening of band gaps. The energy spectrum of AA model is obtained by numerical diagonalization and is depicted in Fig.2a for increasing values of λ. The variation of energy gaps with the coupling strength can be noticed in this figure. The spectral statistics and distribution of energy gaps are analyzed in [7]. Mathematically it can be shown that the energy spectrum of AA model forms a Cantor set [26]. Self-similarity in energy levels can be understood from the integrated level density, where i is the i'th eigenvalue and θ(x) is heaviside step function. The integrated density of states N ( ) shows 'Devil's staircase' like fractal structure, which is evident from Fig.2b where we plotted the normalized integrated density of states N ( ) as a function of scaled energy within the interval of zero to one. For different values of λ, the Devil's staircase structures of N ( ) do not overlap, which indicates that the fractal dimension changes with the strength of the potential λ. Transport properties also show significant changes in localization transition. Due to spatial localization of single particle states transport coefficients vanish in the thermodynamic limit. In electronic systems the conductivity vanishes exponentially with the system size due to 'Anderson localization'. For neutral superfluid corresponding physical quantity is 'superfluid fraction'(SFF), which is measured by generating a superflow by applying a phase twist [27]. In presence of the phase twist the original Hamiltonian in Eq.1 becomes, where Θ is arbitrarily small phase difference across the boundary. For one dimensional system with periodic boundary condition this is exactly similar to a quantum ring in presence of a flux which generates supercurrent through the ring. The superfluid fraction f s is defined as [28], where E 0 (Θ) is the ground state energy of the Hamiltonian H Θ with arbitrary small value of Θ, and E 0 (0) is the ground state energy of the original Hamiltonian given in Eq.1. To obtain the ground state energy upto Θ 2 order, the Hamiltonian H Θ is expanded as, where we define a current operatorĴ = −i n (|n n + 1| − h.c) and the usual kinetic energyT = − n (|n n + 1| + h.c). Using second order perturbation, the SFF in Eqn.8 can be written as, where i , |ψ i are eigenvalues and eigenfunctions of Hamiltonian given in Eq.1 and |ψ 0 is ground state wavefunction. Variation of SFF with the coupling strength of the potential is shown in Fig.3. For increasing values of λ the SFF decreases from unity and vanishes at the critical point λ = 2.

B. Hamiltonian map
Localization phenomena can also be studied by an alternative method of classical Hamiltonian map(CHM), which is very useful to obtain analytical estimate of localization length [29,30]. The discrete Schrödinger equation given in Eq.4 can be written as following dynamical map, where the classical dynamical variables are given by, x n = ψ n and p n = ψ n − ψ n−1 . Here the wavefunction ψ n plays the role of position in CHM and the number of lattice site becomes the number of iteration or time axis of the dynamics. We can choose initial real wavefunctions ψ 0 and ψ 1 (or equivalently x 1 and p 1 ) and evolve the dynamical variables by the transfer matrix, By iteration of this map we obtain the asymptotic behavior of the dynamical variables for a fixed energy eigenvalue . Mathematically it can be shown that the dynamical instability of the above Hamiltonian map corresponds to the localization phenomena. Exponentially localized wavefunctions asymptotically fall off as ψ n ∼ e −n/ξ , where ξ is the localization length. But for arbitrary choice of the initial wavefunction, the solution of Eq.4 also has an exponentially growing component ∼ e n/ξ . Hence the localization is manifested by the exponential growth of the dynamical variables of the Hamiltonian map which gives rise to chaos in the classical phase space. Whereas the extended states represents periodic motion in the phase space. To understand the localization transition in AA model, we choose an eigenvalue close to the band center and calculate the phase space trajectories of the Hamiltonian map for an ensemble of random initial values of x and p. The phase portraits of Eq.12 are shown in Fig.3 for increasing values of λ. In the delocalized regime with small disorder strength λ, the phase portrait consists of closed elliptic orbitals as seen from Fig.3a. Increasing the strength of quasiperiodic potential λ leads to diffusive behavior at the outer part of the phase space and regular portion of phase space with periodic orbits is reduced (see Fig.3b and Fig.3c). Finally for λ close to the critical value, an instability sets in the CHM and phase space trajectories in Fig.3d show chaotic behavior which indicates localization transition. In dynamical systems 'Lyapunov exponent' (LE) is a measure to quantify chaos. Since the number of lattice site plays the role of time in the classical map, the LE corresponds to the inverse localization length of the wavefunction. For periodic motion LE vanishes and in the chaotic regime the non-vanishing LE gives finite localization length ξ. To calculate the LE we construct the matrix U Ns = T Ns T Ns−1 .....T 2 T 1 by multiplying the transfer matrices T n sequentially for N s iteration. The LE l can be obtained by using the formula, where λ + is the largest eigenvalue of the matrix U Ns and LE is obtained from ξ = 1/l. The localization length can also be calculated by using Thouless formula [31]. Using the duality of AA model the analytical expression of LE is given by ξ = log(λ/2) [5], which is independent of energy. Within CHM approach we numerically compute the localization length ξ from Eq.14 using QR decomposition method [32]. Numerically obtained localization length as a function of λ for different energy eigenvalues are compared with the analytical result in Fig.5. For the localization transition in AA model with quasiperiodic potential, the choice of β as irrational number (particularly a Diophantine number) plays a crucial role [24]. To understand this mathematical condition, we analyze the phase space of CHM for successive rational approximation of β, which is shown in Fig.6. Using Fibonacci series, the rational approximation of β can be written as β n = F n−1 /F n and the potential has periodicity F n for successive integer n. For sufficiently large value of n, β n approaches to the inverse of 'golden mean' and the potential becomes quasiperiodic. In the localized regime, for λ = 2.2 the phase space trajectories of CHM with increasing order of rational approximation of β are presented in Fig.6 (a) to (c) for fixed initial condition. Even in the localized regime phase space contains periodic orbits for β 5 . As shown in Fig.6b to Fig.6d, the periodic orbits break as β approaches to the inverse of 'golden mean' by successive rational approximation and finally dynamics become chaotic which is consistent with the localization phenomena.

III. LOCALIZATION OF WEAKLY INTERACTING BOSE GAS
Interacting Bosons in presence of quasiperiodic potential can be described by Bose-Hubbard model [9], where a † i (a i ) are creation (annihilation) operators for Bosons at site i,n i = a † i a i , and U is the strength of onsite repulsive interaction. Above quantum manybody Hamiltonian can capture various correlated phases of strongly interacting Bosons which undergoes quantum phase transition. It is known that Bose-Hubbard model in the presence of random disorder can give rise to 'Bose glass' phase [9]. For sufficiently weak interaction strength U and for high average density of Bosons a 'quasi-condensate' may form in one dimensional systems [33]. In this regime, one can replace the bosonic operators a i by a classical field ψ i which represents the macroscopic wavefunction of the 'quasi-condensate'. Minimization of the classical energy corresponding to Eq.15 leads to the 'discrete nonlinear Schrödinger equation' (DNLS), where µ is the chemical potential, and the normalization of the wavefunction of N b number of Bosons gives We numerically find out the ground state wavefunction of Eq.16 by imaginary time propagation method. For weak disorder the wavefunction is extended and numerical convergence is fast, whereas close to the localization transition more care is needed to obtain the ground state since many metastable states appear in this region.
To obtain the degree of localization of the ground state wavefunction in presence of repulsion we calculate the IPR in real space using Eq.5. The variation of IPR of condensate wavefunction with increasing strength of disorder λ for different values of effective interaction U N b is shown in Fig.7(a). The IPR becomes nonzero for λ ≥ 2 and increases with a slower rate compared to the noninteracting system, indicating spreading of the wavefunction due to the repulsive interaction. The degree of localization decreases with increasing strength of repulsive interaction U N b . In Fig.7(b) the spatial variation of the wavefunction with disorder strength λ is represented by color scale plot. It is evident from this figure that single site localization is not favorable energetically and the wavefunction is localized at almost degenerate but spa-tially separated sites. This particular feature of the fragmented condensate has also been studied in [34]. Due to the multisite localization the IPR is much less than unity and increases very slowly with λ. A change in the slope of IPR with increasing λ occurs when the number of localized sites decreases and we have checked that finally at very large value of λ the wavefunction becomes localized at a single site.
To study the interplay between disorder and interaction in the transport properties of dilute Bose gas, we calculate the superfluid fraction f s . To generate a superflow we introduce a small amount of phase twist (Θ ∼ 0.1) in the hopping term of DNLS (in Eq.16) similar to Eq.7 and then the superfluid fraction is computed using the formula, where E cl (Θ) is the classical energy corresponding to the Hamiltonian in Eq.15, The wavefunction φ i minimize E cl (Θ) and is normalized to unity. The superfluid fraction f s as a function of λ for different values of repulsive interactions U N b is shown in Fig.8. The SFF of weakly interacting Bose gas obtained from 'density matrix renormalization group' also shows similar behavior [35]. Due to the repulsive interaction, f s vanishes at larger strength of quasiperiodic potential λ > 2, however the IPR rises from zero at λ ≈ 2. This behavior is different from the noninteracting AA model. We also investigated the localization properties of DNLS by Hamiltonian map approach. Eq.16 can be writ-ten in the form of nonlinear classical map, (20) where, nonlinearity gives rise to an unstable classical potential ∼ − U N b 4 x 4 i due to which the classical trajectories become unstable at large values of phase space variables. To avoid this problem we choose relatively small region of phase space within which the potential is metastable, and study the effect of disorder in the phase space dynamics. In Fig.9, for a fixed value of nonlinearity U N b we show the phase portrait with an ensemble of initial configurations for increasing values of disorder strength λ. Similar qualitative features like non interacting system are also seen in the phase space dynamics of DNLS, however the classical periodic orbits are modified due to the nonlinearity.

A. Fluctuations within Bogoliubov approximation
So far we studied the weakly interacting Bosons within mean-field approximation using macroscopic wavefunction for the condensate. It is also important to analyze the quantum fluctuations induced by the quasiperiodic disorder potential. Within Bogoliubov approximation, the quantum field operators can be approximated by, where ψ i is the macroscopic wavefunction of the condensate satisfies Eq.16, u ν i ,v ν i are amplitudes corresponding to νth eigenmode with bosonic operators b ν , b † ν . Bogoliubov quasiparticle energies ω ν can be obtained from, and normalization condition gives For a given interaction strength U and average density of Bosons n 0 = N b /N s , we first numerically obtain the ground state macroscopic wavefunction ψ i , then diagonalize Eq.23 to calculate Bogoliubov quasiparticle energies and amplitudes u ν i , v ν i . The Bogoliubov energy spectrum for increasing strength of disorder λ is depicted in Fig.10(a) for a fixed value of interaction U N b = 1 and N s = 144. The normalized integrated density of states N (ω) also shows Devil's staircase like structure as shown in Fig.10(b). Due to the quantum fluctuations, depletion of the condensate occurs and the non-condensate density ρ nc at zero temperature can be obtained from the Bogoliubov theory, The condensate fraction is given by N c /N b = 1 − i,ν |v ν i | 2 /N b . In one dimensional system the noncondensate fraction diverges as log(N s ), which prohibits the formation of condensate in the thermodynamic limit. However for sufficiently weak interaction a quasi condensate can form in quasi one dimensional and finite system [33]. To investigate the disorder induced quantum fluctuation in a quasi-condensate in a finite lattice with N s = 144, we calculate the condensate fraction with increasing disorder strength λ, which is shown in Fig.11(a). In absence of disorder the quantum depletion is small in weakly interacting gas of Bosons and increases with the interaction strength. With the increase in disorder strength λ the condensate fraction remains close to unity in the delocalized regime and then rapidly decreases around the critical point λ ≈ 2. This qualitative feature (as shown in Fig.11(a)) indicates that enhanced quantum fluctuations near the localization transition can destroy the quasi condensate above λ ≈ 2 and strongly correlated phases can appear. The variation of non-condensate density ρ nc with disorder strength λ is depicted in Fig.11(c) by color scale plot. It is interesting to note that the distinct feature of multisite localization for λ > 2 is observed even for non-condensate density. Also the IPR of normalized non-condensate density shows similar behavior as that of the condensate wavefunction and increases from λ ≈ 2 (as seen from Fig.11(c)). Localization of Bogoliubov quasiparticles and enhancement of quantum fluctuations in presence of quasiperiodic potential particularly for λ ≥ 2 are clearly evident from this analysis.

IV. LOCALIZATION IN THE PRESENCE OF A TRAP
Although the localization transition in AA model occurs in thermodynamic limit, in real experimental setup a weak trapping potential is always present in order to confine the ultracold atoms. Main features of the localization can also be observed in trapped system provided the length scale of the trapping potential is larger compared to the localization length. Additionally some interesting effects due to the trap can also be seen. A harmonic trap is introduced by adding a potential V i = 1 2 ω 2 HO (i − i c ) 2 in the Hamiltonian (Eq.15), where ω HO is related to the trapping frequency, i is the site index, and i c is the center of the trap. Although in a trap the wavefunction is width of the wavepacket is sharply reduced due to disorder induced localization. The density profile of the condensate in presence of a harmonic trap for different disorder strength λ is shown in Fig.12b. We calculate the IPR of the ground state wavefunction with increasing disorder strength λ and an increase of IPR around λ ≈ 2 is observed as expected. However the IPR does not vanish in the delocalized regime λ < 2 and takes a small value due to the tapping potential. We have noticed earlier that in the absence of a trap the IPR increases very slowly after the localization transition due to the repulsive interaction and the wavefunction is localized at spatially separated sites with quasi-degenerate energies. These quasi-degeneracy of onsite energies can be lifted by introducing a harmonic trap which leads to enhancement of the degree of localization which is elucidated in Fig.12a where a rapid increase of IPR to unity is shown by increasing the trap frequency by a small amount.
Collective oscillation of a Bose-Einstein condensate in a trap can also show its superfluid properties. A small displacement of the condensate from the center of the trap generates center of mass (COM) oscillation which is a well studied collective mode of Bose-Einstein condensate. Here we numerically study the motion of COM of trapped condensate in presence of quasi periodic disorder. Since the superfluidity is affected by the disorder, we calculate the coherence factor Ψ of the oscillating condensate which is given by [36] We can notice that the above expression contains information of relative phase difference between neighboring sites and |Ψ| 2 gives a quantitative measure of overall phase coherence of the oscillating condensate. In Fig.13 (a) and (b), we have shown the COM motion of the trapped condensate and coherence factor of corresponding macroscopic wavefunction for increasing disorder strength λ. It is clear that the coherence of the time dependent wavefunction decreases with increasing disorder. It is also interesting to study the variation of condensate fraction with disorder. In one dimensional trapped condensate the divergence of non-condensate fraction is less severe due to the finiteness of the trap. Like the homogeneous system, we expect enhancement of quantum fluctuation due to disorder which may destroy the condensate around a critical disorder strength λ ≈ 2.

V. CONCLUSIONS
In conclusion, we have investigated localization of both non-interacting Bosons as well as weakly interacting quasi-condensate in presence of a quasiperiodic potential. Apart from calculating various physical properties, understanding localization transition through the approach of classical Hamiltonian map is one of the main result of this work.
In the non-interacting Aubry-André model the localization transition can be identified by the vanishing 'superfluid fraction' and the rise of IPR at the critical strength of quasiperiodic potential λ = 2. A classical Hamiltonian map is constructed from the Schrödinger equation. The phase space trajectories of the corresponding classical map shows periodic orbits for small disorder strength. With the increasing disorder strength λ chaotic behavior is observed at the outer region of the phase space and finally near the critical value λ = 2 all periodic orbits are destroyed due to the onset of chaos in the localized regime. The parameter β = ( √ 5 − 1)/2 is chosen to be an irrational Diophantine number(inverse of golden mean) which is an essential requirement for localization transition. This has been elucidated by means of Hamiltonian map for λ > 2 where the destruction of periodic orbits by successive rational approximation of β indicates onset of localization.
We have investigated the localization of weakly interacting Bose gas in quasiperiodic potential by mean-field approach. Unlike non-interacting case, the vanishing of SFF and rise of IPR does not take place at same strength of disorder λ. Due to the repulsive interactions, the macroscopic wavefunction is localized at many sites with quasi-degenerate energies for λ > 2. The multisite localization of the interacting system is manifested by a much slower increase of IPR starting from λ = 2. With increasing the strength of quasiperiodic potential, the number of sites over which the condensate is localized decreases and finally the wavefunction becomes localized at a single site for very large value of λ and IPR approaches to unity. The SFF decreases with increasing disorder strength λ and vanishes at λ > 2 due to the repulsive interaction. The repulsive interaction gives rise to an unstable nonlinear potential in the CHM approach, due to which the stable region of phase space decreases. In the phase portrait the stable region containing the periodic orbits decreases with increasing λ and finally the onset of chaos at λ ≈ 2 signifies localization of the wavefunction. Further, the Bogoliubov quasiparticle spectrum has been calculated numerically. We notice that disorder enhances the quantum fluctuations due to which the condensate fraction of a finite system decreases rapidly around λ ≈ 2. This indicates the possible formation of glassy phase and multisite localized insulators.
Finally we also considered the effect of trapping potential on the localization transition. In the localized regime, the number of sites over which the wavefunction is localized reduces due to the presence of a trap which has been shown from the rapid increase in the IPR by tuning the trap frequency. The center of mass motion of the condensate in a harmonic trap also shows the signature of localization. The center of mass oscillations become incoherent with increasing disorder. To summarize, the present study provides a clear picture of localization of non-interacting and weakly interacting Bose gas in the presence of a quasiperiodic disorder and it reveals various interesting features which are interesting for both academic point of view, as well for future experiments.