The Aubry-Andr\'e model as the hobbyhorse for understanding localization phenomenon

TL;DR

This paper provides a comprehensive pedagogical analysis of localization in a one-dimensional quasiperiodic lattice using the Aubry-André model, exploring stationary and dynamical properties to understand the phenomenon.

Contribution

It offers a detailed analysis of localization phenomena in the Aubry-André model, combining stationary and dynamical perspectives to enhance understanding.

Findings

Hofstadter butterfly pattern observed in energy spectrum

Localization characterized by IPR and NPR metrics

Dynamical spreading distinguishes localized from delocalized states

Abstract

We present a thorough pedagogical analysis of the single particle localization phenomenon in a quasiperiodic lattice in one dimension. Description of disorder in the lattice is represented by the Aubry-Andr\'e model. Characterization of localization is performed through the analysis of both, stationary and dynamical properties. The stationary properties investigated are the inverse participation ratio (IPR), the normalized participation ratio (NPR) and the energy spectrum as a function of the disorder strength. As expected, the distinctive Hofstadter pattern is found. Two dynamical quantities allow discerning the localization phenomenon, being the spreading of an initially localized state and the evolution of population imbalance in even and odd sites across the lattice.

Figures

Equations

$$\hat{H}=\int dx{\hat{\psi }}^{†}(x)\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+V(x)\right)\hat{\psi }(x)$$

$$+\frac{1}{2}\frac{4\pi a_s{\hslash }^2}{m}\int dx{\hat{\psi }}^{†}(x){\hat{\psi }}^{†}(x)\hat{\psi }(x)\hat{\psi }(x),$$

$$V_{opt}(x)$$

$$\text{eqalign}\hat{\psi }(x)=\sum _i{\hat{b}}_{i}w(x-x_i),\text{cr}{\hat{\psi }}^{†}(x)=\sum _i{\hat{b}}_i^{†}{w}^{\ast }(x-x_i).$$

$$\hat{H}=-\sum _{i,j}{J}_{ij}{\hat{b}}_i^{†}{\hat{b}}_j+\sum _{ij}{ϵ}_{ij}{\hat{b}}_i^{†}{\hat{b}}_j+\sum _{ij}{\Delta }_{ij}{\hat{b}}_i^{†}{\hat{b}}_j+\sum _{i,j,l,v}{U}_{i,j,l,v}{\hat{b}}_i^{†}{\hat{b}}_j^{†}{\hat{b}}_l{\hat{b}}_v,$$

$$\text{eqalign}{J}_{ij}=-\int dx{w}^{\ast }(x-x_i)\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+s_1{E}_{R_1}{\sin }^2(k_{1}x)\right)w(x-x_j)\text{cr}{ϵ}_{ij}=\int dx{w}^{\ast }(x-x_i)V_T(x)w(x-x_j)\text{cr}{\Delta }_{ij}=s_2{E}_{R_1}{\beta }^2\int dx{w}^{\ast }(x-x_i){\sin }^2(\beta k_{1}x+\phi )w(x-x_j)\text{cr}{U}_{ijlv}=\frac{4\pi a_s{\hslash }^2}{2m}\int dx{w}^{\ast }(x-x_i)w^{\ast }(x-x_j)w(x-x_l)w(x-x_v).$$

$${\Delta }_{ij}=-\frac{s_2{E}_{R_1}{\beta }^2}{2}\int dx{w}^{\ast }(x-x_i)\cos (2\beta k_{1}x+{\phi }^{\prime })w(x-x_j).$$

$${\Delta }_{ii}=-\frac{s_2{E}_{R_1}{\beta }^2}{2}\cos (2\pi \beta i+{\phi }^{\prime })\int dy\cos (2\beta k_{1}y)|w(y){|}^2.$$

$$\displaystyle\eqalign{\cos(2\beta k_{1}y+2\beta k_{1}x_{i}+\varphi^{\prime})=&\cos(2\pi\beta i+\varphi^{\prime})\cos(2\beta k_{1}y)-\cr&\sin(2\pi\beta i+\varphi^{\prime})\sin(2\beta k_{1}y),}$$

$${\Delta }_{ij}=\Delta \cos (2\pi \beta i+\varphi ){\delta }_{ij},$$

$$\Delta =\frac{s_2{E}_{R_1}{\beta }^2}{2}\int dy\cos (2\beta k_{1}y)|w(y){|}^2.$$

$$U=U_{i,i,i,i}=\frac{4\pi a_s{\hslash }^2}{2m}\int |w(x){|}^{4}dx.$$

$$\hat{H}=-J\sum _{⟨i,j⟩}{\hat{b}}_i^{†}{\hat{b}}_j+\Delta \sum _i\cos (2\pi \beta i+\varphi ){\hat{n}}_i+U\sum _i{\hat{n}}_i({\hat{n}}_i-1),$$

$$\hat{H}=-J\sum _{⟨i,j⟩}{\hat{b}}_i^{†}{\hat{b}}_j+\Delta \sum _i\cos (2\pi \beta i+\varphi ){\hat{n}}_i.$$

$$\hat{H}=-J\sum _j(|w_j⟩⟨w_{j+1}|+|w_{j+1}⟩⟨w_j|)+\Delta \sum _j\cos (2\pi \beta j+\varphi )|w_j⟩⟨w_j|.$$

$$|k_s⟩=\sum _j{e}^{2\pi i\beta k_{s}j}|w_j⟩$$

$$H=-\frac{\Delta }{2J}J\sum _s(|k_s⟩⟨k_{s+1}|+|k_{s+1}⟩⟨k_s|)+\frac{2J}{\Delta }\Delta \sum _s\cos (2\pi \beta s)|k_s⟩⟨k_s|,$$

$$l(E{)}^{-1}={\int }_{-\infty }^{\infty }\ln |E-ϵ|d\rho (ϵ).$$

$$G(y)=\frac{{\Delta }^2}{2\pi }\cos (2\pi \beta y).$$

$$\hat{H}=-J\sum _i(|x_i⟩⟨x_{i+1}|+|x_{i+1}⟩⟨x_i|)+\Delta \sum _i\cos (2\pi \beta i+\varphi )|x_i⟩⟨x_i|$$

$$-({\psi }_{j+1}+{\psi }_{j-1})+\frac{\Delta }{J}\cos (2\pi \beta j+\varphi ){\psi }_j=\frac{E}{J}{\psi }_j,$$

$$\displaystyle\mathbf{H}=\left(\begin{array}[]{ccccc}\Delta/J\cos(2\pi\beta+\varphi)&-1&\cdots&-1\\ -1&\Delta/J\cos(4\pi\beta+\varphi)&\cdots&0\\ \vdots&\ddots&\ddots&-1\\ -1&0&-1&\Delta/J\cos(2\Omega\pi\beta+\varphi)\end{array}\right),$$

$$\mathrm{IPR}(|\psi ⟩)=\sum _i|{\psi }_i{|}^4.$$

$$\mathrm{NPR}(|\psi ⟩)=\frac{1}{\Omega }\frac{1}{{\sum }_i|{\psi }_i{|}^4}.$$

$$i\hslash \frac{\partial }{\partial t}|\psi ⟩=H|\psi ⟩.$$

$$\frac{d}{dt}{\psi }_i=-\frac{i}{\hslash }\sum _j{H}_{ij}{\psi }_j,$$

$$\text{eqalign}\overrightarrow{k_1}=\mathbf{H}({\tau }_j)\vec{\psi }({\tau }_j)\text{cr}\overrightarrow{k_2}=\mathbf{H}({\tau }_j+\Delta \tau /2)(\vec{\psi }({\tau }_j)+{\vec{k}}_1/2)\text{cr}\overrightarrow{k_3}=\mathbf{H}({\tau }_j+\Delta \tau /2)(\vec{\psi }({\tau }_j)+{\vec{k}}_2/2)\text{cr}\overrightarrow{k_4}=\mathbf{H}({\tau }_j+\Delta \tau )(\vec{\psi }({\tau }_j)+{\vec{k}}_3)\text{cr}\psi ({\tau }_{j+1})=\psi ({\tau }_j)+\frac{\Delta \tau }{6}(k_1+2k_2+2k_3+k_4).$$

$$\sigma (\tau )={\left[\sum _i{i}^2|{\psi }_i(\tau ){|}^2\right]}^{1/2}.$$

$$\mathrm{I}(\tau )=\frac{n_e(\tau )-n_o(\tau )}{n_e(\tau )+n_o(\tau )},$$

$$\mathrm{I}=\underset{\tau \to \infty }{\lim }\frac{1}{\tau }\int d{\tau }^{\prime }\frac{n_e({\tau }^{\prime })-n_o({\tau }^{\prime })}{n_e({\tau }^{\prime })+n_o({\tau }^{\prime })}$$

Full text

• December 2018

1 Introduction

The understanding of electronic mobility in quasiperiodic or disordered media is one of the fundamental issues in the condensed matter domain. Disorder, together with the interparticle interactions intrinsically present in every macroscopic sample, are the responsible ones of the physical behavior and response properties of solids. A material becomes insulator as a result of either electron-electron or electron-ion interactions. While the study of electron-electron interactions demands the use of a many body theory and leads to identify a Mott insulating phase [1], electron-ion interactions are addressed within the single electron theory and allow to discern among several types of insulators. Among them, band insulators arising from periodicity in the lattice [2], Peierls insulators associated with lattice distortions, and Anderson insulators resulting from lattice imperfections or impurities, also known as lattice disorder. The purpose of this manuscript is to show how this loss of long range order in the lattice leads to localize single electrons and thus cause the absence of its diffusion across the lattice. Although localization phenomenon has been extensively investigated, is up to now a central topic which still have a manyfold of open questions to be addressed. For instance, the interplay between disorder and dimensionality, interactions and inhomogeneity created by external fields among others.

Characterization of localization phenomenon requires as a first step the election of an effective model representing the disorder, and then the use of the standard quantum mechanics techniques to analyze their effects. There are two general schemes from which localization has been envisaged, the Anderson model [3], in which disorder is represented by a random amplitude of the on-site energies, and the Aubry-André model [4, 5], where disorder arises from the superposition of two lattice potentials with incommensurate wavelengths. These models capture the metal-insulator transition in disordered lattices and allow to characterize such transition by tracking different properties as we describe below. Here, we shall use the Aubry-André model as the hobbyhorse for studying and characterizing the influence that disorder has in producing localized states. This model introduced in 1980, has shown to be very successful in describing such transition not only in the single electron case but when interparticle interactions are considered [6, 7, 8, 9, 10]. Localization in a lattice can be recognized through several signatures, either of stationary or dynamical character. What it is important to stress is that localization can result from both, destructive interference associated with the multiple scattering of the wave function traveling along the disordered medium and the spectral properties of the Schrödinger equation [11]. Destruction of wave coherence or loss of mobility is quantified in terms of several properties that can be extracted from the wave function. The first distinctive signature of localization, identified in the seminal work of Anderson, was the localization length that measures the size of the exponentially localized single particle state as a function of disorder strength. Here we concentrate on analyzing the properties enunciated in the following lines. First, we analyze the properties of the Aubry-André model, and then we investigate the inverse participation ratio (IPR) and its opposite, the normalized participation ratio (NPR), that quantify the fraction of sites contributing to the state along the lattice. Next, we investigate the energy spectrum that also allows monitoring the transition to localization. In addition to these quantities characterizing localization of the stationary states, there are dynamical parameters that also allow to track the evolution of a given initial state in the presence of disorder. Among them, the spreading of the initial state and the imbalance between the density probability of even and odd sites in the lattice, as a function of time.

Current experiments with ultracold neutral atoms realized in the laboratory represent the ideal scenario where the spatial quasiperiodicity of the Aubry-André model can be recreated. Optical lattice potentials produced by standing waves resulting from interfering laser fields emulate such a non crystalline environment seen by electrons moving across the ion cores. Nowadays, such large ensembles of fermionic or bosonic atoms loaded in optical lattices offer advantages with respect to experiments performed in solids since ultracold atoms can be prepared to analyze isolated effects present in solids, without the influence of further outcomes [12]. As a matter of fact, 10 years ago the Aubry-André model was experimentally set in a laboratory for the first time [13].

The aim of the present manuscript is to present a pedagogical description of the Aubry-André model to understand and characterize the localization phenomenon. Advanced undergraduate and graduate students should be able to follow this article with no difficulty. We believe that this material should give the appropriate tools and techniques to face and approach forefront problems including the many-body localization phenomenon. The manuscript is organized in 6 sections. First in section 2 we derive the Aubry-André model demonstrating how quasiperiodicity in the potential gives rise to a cosine function incommensurate with the underlying periodic tight-binding 1D lattice. Then, in section 3 the properties characterizing the Aubry-André model are delineated. Sections 4 and 5 account for the time independent and time dependent analysis that characterizes the localization transition. Finally, in section 6 a summary of results is presented.

2 Model

Our starting point is the 1D Hamiltonian operator for ultracold bosonic atoms with mass $m$ confined in an external potential $V(x)$ and interacting via a contact potential written, as usual, in terms of the $s$- wave scattering length $a_{s}$,

[TABLE]

$\hat{\psi}^{\dagger}(x)$ and $\hat{\psi}(x)$ are the bosonic creation and annihilation field operators satisfying the standard commutations rules for bosons, the external potential $V(x)$ is given by the superposition $V(x)=V_{T}(x)+V_{opt}(x)$, in which $V_{T}(x)$ is a slowly varying magnetic harmonic trap and, $V_{opt}(x)$ an optical lattice potential. Figure 1(a) shows a sketch of the resulting potential $V(x)$. In presence of disorder, $V_{opt}(x)$ consists of two optical lattices [14], the main lattice $V_{1}(x)=s_{1}E_{R_{1}}\sin^{2}(k_{1}x)$ which is used to create a tight-binding environment for the atoms and a secondary one $V_{2}(x)=s_{2}E_{R_{2}}\sin^{2}(k_{2}x)$ which introduces an optical disorder [15]. Superimposing both lattices gives rise to the following bichromatic potential:

[TABLE]

where $k_{i}=2\pi/\lambda_{i}$ (i=1,2) are the wave vectors, with $\lambda_{i}$ the wavelength of the lasers fields, $s_{i}$ are the heights of the lattices in units of the recoil energy $E_{R_{i}}=h^{2}/(2m\lambda_{i}^{2})$, $\varphi$ is an arbitrary phase and $\beta=\lambda_{1}/\lambda_{2}$ the ratio between the wavelengths. When $s_{2}<<s_{1}$ and $\beta$ is an incommensurate number, the secondary lattice does not change considerably the positions of the potential minima generated by the main lattice [16]. Instead, as shown in figure 1(b) it has the effect of shifting the local site energy by an amount $\Delta_{i}$ only.

For single atoms and no net disorder $s_{2}=0$, the eigenstates of equation (1) are Bloch wave functions [17]. As it is well known, an appropriate linear combination of Bloch states yields a Wannier wave function $w_{\nu}(x-x_{i})$, characterized by large probability amplitude around lattice site $i$, that is, a localized wave function at each site $i$. Since the atoms under study are at ultracold temperatures, it is well justified the assumption that the energies involved in the system are smaller compared to the energy required to allow second and higher band populations. This consideration allows us to drop the band index $\nu$ in the Wannier functions and contemplate first band populations only. Having this assumption in mind, it is convenient to expand the field operators $\hat{\psi}(x)$ and $\hat{\psi}^{\dagger}(x)$ in the Wannier basis:

[TABLE]

Being $\hat{b}_{i}$ and $\hat{b}^{\dagger}_{i}$ the annihilation and creation operators for a particle in a Wannier state at the lattice site $i$ respectively. It is worthwhile to stress that the sums in equation (5) run over all lattice sites. As mentioned above, for weak disorder $s_{2}<<s_{1}$ the minima of the main lattice are not remarkably affected and we can safely substitute the latter expansion of the field operators in equation (1). After some straightforward algebra, one can obtain the following expression for the Hamiltonian (1):

[TABLE]

where we have defined the following constants:

[TABLE]

The first term in equation (6) describes the energy cost for an atom to hop from site $i$ to site $j$ with $i\neq j$, see figure 2.

Note that the hopping probability $J_{ij}$ is proportional to the overlap between the Wannier functions centered at different lattice sites. Within the so-called tight-binding approximation, this overlap is essential only for nearest neighbors [17], thus we can neglect the tunneling between next nearest neighbors and beyond. Also, since the main lattice potential is invariant under translations by one lattice period, $a=\lambda/2$, the hopping parameter $J_{ij}$ becomes a constant $J$ independent of the lattice site. For the case $i=j$ the $J_{ii}$ term represents an on-site energy shift which is equal for all sites and thus can be dropped. Also, the second term $\epsilon_{ij}$ represents an on-site shift of the energy. However, this time we have to assume that the frequency $\omega_{T}$ of the harmonic trap satisfies [18] $\sqrt{\hbar/(m\omega_{T}})<<\lambda_{1}/2$ which allow us to consider the variation of the energy from site to site as equal, thus to a first approximation, we can neglect the contribution of the harmonic trap, see figure 1(a). The third term in equation (6) is the responsible for the optical disorder in the lattice. In order to deal with it, we first use the trigonometric relation $\sin^{2}(\beta k_{1}x+\varphi)=(1-\cos(2\beta k_{1}x+\varphi^{\prime}))/2$ with $\varphi^{\prime}=2\varphi$. Inserting this relation into the third equality of equation (7) and dropping the constant term, we obtain:

[TABLE]

Again, for deep enough lattices, the leading contribution of equation (8) is the $i=j$ term, which corresponds to an on-site energy shift variation. Further, we can make the change of variable $y=x-x_{i}$, leading to:

[TABLE]

Where in the last equation we have identified $x_{i}\rightarrow i$, used the subsequent trigonometric identity:

[TABLE]

and symmetric properties to drop the sine integral. Following the above steps we finally get the usual disorder term [14]:

[TABLE]

where $\phi=\varphi^{\prime}+\pi$ and $\Delta$ is defined as the following constant parameter:

[TABLE]

The term associated with the interaction energy can also be simplified by taking into account the tight binding approximation. The dominant term of the overlap of four Wannier functions is due to the term $i=j=l=v$, which corresponds to an onsite interaction, where the atoms only “see each other” whenever they are in the same lattice site.

[TABLE]

Summarizing all the above approximations, we end with the following interacting Hamiltonian:

[TABLE]

where the notation $\langle i,j\rangle$ indicates that the sum runs over nearest neighbors only and $\hat{n}_{i}=\hat{b}_{i}^{\dagger}\hat{b}_{i}$ is the number operator at site $i$. One of the most outstanding advantages of the experiments with ultracold atomic gases is the possibility of tuning the strength of the pairwise interactions between atoms, via an external magnetic field. This procedure called Feshbach resonance [12] allows the experimentalist to prepare a gas of atoms with a zero scattering length $a_{s}$, and consequently $U=0$. Such a non-interacting system constitutes an experimental realization of the non-interacting Harper [5] or Aubry-André [4] model:

[TABLE]

This model and also the interacting version has been realized in experiments with ultracold atoms in bichromatic optical lattices potentials [13, 19]. Written in the Dirac notation, the above Hamiltonian takes the form:

[TABLE]

It is important to emphasize that recent investigations [11, 20] pointed out that the Hamiltonian in equation (18) is strictly valid in the extreme tight-binding limit of a very deep main lattice potential.

3 Properties of the Aubry-André model

Having set the Aubry-André Hamiltonian it is worthwhile to expose some basic properties of this model. Beginning with its duality in space and momentum representations, one can transform the Hamiltonian (18) written in the Wannier representation to one in the momentum representation, via the following transformation:

[TABLE]

After substitution and straightforward algebra, we find the dual Hamiltonian:

[TABLE]

which has the same structure as that of (18) except that the tunneling rate has changed from $J\rightarrow\Delta/2$ and the disorder strength from $\Delta\rightarrow 2J$. Also, for simplicity we set $\phi=0$. Now, according to the Heisenberg uncertainty principle the Hamiltonian (20) has localized states where the first Hamiltonian (18) has extended states, and viceversa. Obviously, the transition from extended to localized states or from localized to extended states must take place at the same set of parameters for both Hamiltonians. This reasoning leads us to impose a condition over the two parameters characterizing both Hamiltonians, the condition at the transition point is $\Delta/J=2$. Furthermore, the localization transition is sensitive to the value of the parameter $\beta$ [14]. For instance, an integer value $\beta$ would not display any localization transition. In order to observe such a transition, $\beta$ must have some degree of incommensurability. One way to achieve this requirement is to choose $\beta$ as the ratio of two adjacent Fibonacci numbers $F_{n-1}/F_{n}$ [21], with $F_{n-1}$ and $F_{n}$ two consecutive elements of the Fibonacci sequence. Such a procedure approaches the inverse Golden ratio $\varphi=(\sqrt{5}-1)/2$ for large enough Fibonacci numbers. Also, one can express $\beta$ as the ratio of two relative prime numbers $\beta=P/Q$ with $P$ and $Q$ larger than the system size of simulation [22]. However, in experimental realizations, $\beta$ is restricted to the available laser wavelengths. For instance, recent experiments [20] were performed with $\beta=532/738$. Summarizing, we can say that we must ensure that $\beta$ is such that the system remains aperiodic within the size of interest [11]. If these requirements are accomplished one would observe that for $\Delta/J=2$, the wave function develops peaks around certain lattice sites, such peaks become a single peak as the disorder amplitude $\Delta$ is increased. As shown by Thouless [23], there exists a link between the density of states $\rho$ and the inverse of the localization length $l$ of the wave function. In general, the localization length depends on the energy of the localized wave function. One can find such dependence by solving the following integral:

[TABLE]

For the Aubry-André model considered here, the latter equation gives an energy independent localization length of $l=d/\ln(\Delta/2J)$ [22] which implies that all eigenstates would display the same localization length.

Another remarkable property of the Aubry-André model is related to the type of disorder generated, which is called correlated disorder. This is because, the correlation function of the disorder potential $G(y)=\langle V(x)V(x+y)\rangle$ takes the value:

[TABLE]

This disorder correlation is in contrast with the uncorrelated Anderson disorder [3] where any random disorder localizes the system.

4 Time independent results

In this section, we describe some of the most distinctive single particle time independent results of the Aubry-André model. In the deep tight-binding approximation, the Wannier functions are highly localized and thus can be represented by the site basis [17], that is, $|w_{i}\rangle\rightarrow|x_{i}\rangle$. In this scenario, the equation (18) is reduced to the Schrödinger equation $\hat{H}|\psi\rangle=E|\psi\rangle$ with $|\psi\rangle$ being an eigenstate of the one particle Hamiltonian $\hat{H}$:

[TABLE]

By expanding the wave function $|\psi\rangle$ in terms of the site basis $|\psi\rangle=\sum_{i}\psi_{i}|x_{i}\rangle$ and calculating the product $\langle x_{j}|\hat{H}|\psi\rangle$, we obtain the following difference equation for the coefficients $\psi_{i}$:

[TABLE]

where we have divided by $J$ in order to have dimensionless equations and at the same time to set the energy scale. It is easy to see that equation (24) can be rewritten in a matrix form $\mathbf{H}\vec{\psi}=\frac{E}{J}\vec{\psi}$ being $\vec{\psi}=(\psi_{1},\psi_{2},\cdots,\psi_{\Omega})^{T}$ the state vector and $\mathbf{H}$ the Hamiltonian matrix:

[TABLE]

where $\Omega$ is the number of sites in the lattice. In the following calculations we consider $\Omega=987$ sites, although similar results are found for larger lattice sizes. At this point, it is instructive to state two aspects that were considered in our numerical calculations. The first one is related to the boundary conditions of the problem. As can be seen from the Hamiltonian matrix (29), we take periodic boundary conditions, which means that the 1D lattice closes itself. Also, for the subsequent calculations we consider $\beta=610/987\approx 0.618034$.

Now, we are left to find the eigenvalues and eigenvectors of the matrix $\mathbf{H}$. In the case of vanishing disorder $\Delta/J=0$, the equation (24) is easily solved with the ansatz $\psi_{j}=e^{ikaj}$, which displays the energy spectrum of a free particle in a 1D lattice $E_{k}=-2J\cos(ka)$ [2]. In figure 3(a) we plot the ground state density $|\psi_{i}|^{2}$ as a function of the lattice site $i$ for zero disorder $\Delta/J=0$ and nonzero disorder $\Delta/J=1.5$ in figure 3(b).

As one would expect, the ground state profile in the absence of disorder is a normalized constant at each site. This means that the particle is completely delocalized in the lattice. For the case $\Delta/J=1.5$, the density profile of the ground state displays multiple peaks which indicates the presence of different potential depths across the sites. However, the wave function is still extended over all the lattice.

With aim of sketching the ground state density for two different values of the disorder amplitude satisfying $\Delta/J\geq 2$, in figure 4 we exhibit two cases, the left one associated to $\Delta/J=2$ and the right one corresponding to $\Delta/J=4$. As can be seen from these density profiles, $\Delta/J=2$ exhibits the transition from an extended to a localized state, while the case $\Delta/J=4$ shows that when the disorder amplitude is increased, the localization becomes sharply.

An important quantity, that arises in describing the localization transition, and is widely used in the literature, is the inverse participation ratio (IPR). For a normalized state $|\psi\rangle=\sum_{i}\psi_{i}|x_{i}\rangle$, it is defined as:

[TABLE]

The IPR gives us the inverse of the number of sites occupied by the wave function. For instance, it approaches zero as $1/\Omega$, for an extended wave function, while it goes to $1$ for a localized state on a single lattice site.

The great advantage of using the IPR parameter is that instead of looking at the wave function in each realization, we just have to check a single parameter to confirm the nature of the wave function. In figure 5 we illustrate the IPR associated with the ground state as a function of the disorder strength $\Delta/J$. Each point in this figure corresponds to an average over 10 realizations of the phase $\phi$ which is set randomly in the interval $\phi\in[0,2\pi)$. As it is shown, the IPR parameter becomes different from zero for $\Delta/J=2$ and it approaches the unity as the disorder increases. This peculiar behavior makes the IPR a good parameter for studying the localization transition. Even though the above definition of the IPR is related to the localization in real space, one could extend the idea to momentum space or more exotic basis as the Floquet basis in periodic driven optical lattices [24, 25].

Another interesting quantity is the normalized participation ratio (NPR), which plays the opposite role of the IPR parameter [11]. That is, the NPR parameter remains finite for spatially extended states, while approaching zero for a localized one. For a given normalized state $|\psi\rangle=\sum_{i}\psi_{i}|x_{i}\rangle$, the NPR parameter is defined as follows:

[TABLE]

In order to illustrate the utility of the NPR parameter, first, one should notice that the above definitions for the IPR and NPR parameters are related to a single eigenstate $|\psi\rangle$. However, one can calculate these two quantities for the full eigenstate spectrum, and display the average of the IPR and NPR parameters. In figure 6 we plot such averages as a function of the disorder strength $\Delta/J$. As it is shown in this figure, for $\Delta/J<2$ we obtain, as expected, IPR $=0$ and NPR $\neq 0$, while for $\Delta/J>2$, IPR $\neq 0$ and NPR$=0$. This indicate that the spectrum is either, completely delocalized ($\Delta/J<2$) or completely localized ($\Delta/J>2$), but not a mixture of both, localized and extended states.

We should emphasize that this peculiar behavior of IPR and NPR curves, is a direct consequence of the tight-binding approximation used in deriving the Aubry-André model [11]. Nevertheless, when hopping to next nearest neighbors is considered, a noticeable overlap between both curves emerges [11, 26, 27, 28]. In such scenario extended and localized states take place in the same spectrum, the value of the energy that separates localized and delocalized eigenstates is called mobility edge energy. At the transition point $\Delta/J=2$ all the eigenstates exhibit a multifractal structure [29], a subject which is out of the scope of this work. Summarizing, we can organize the latter result in the diagram shown in figure 7 which displays an absence of mixture between localized and extended eigenstates in the full spectrum.

To conclude this section, in figure 8 we show the energy spectrum of the Aubry-André model as a function of the incommensurate parameter $\beta$ for two different disorder strengths $\Delta/J$. As described above, the values of $\beta$ must be chosen such that they have some degree of incommensurability. Here we have used $\beta\in[0,1)$ such that $\beta=i/\Omega$ with $i=1,2,3,...,987$. This spectrum has been studied in numerous works [30, 31, 32, 33, 34, 35] since it displays very rich structure in both the extended and localized regimes.

For values of disorder amplitude above and below of the transition point, the energy spectrum structure is completely different from that associated to $\Delta/J=2$, where it shows the Hofstadter butterfly spectrum [32]. The case $\Delta/J=2$ has been widely investigated, since it describes the quantum physics of an electron moving on a two dimensional square lattice in the presence of a transverse magnetic field [5].

5 Time dependent results

In this section, we discuss the time evolution of a given initial condition, in presence of disorder. This is one of the ways in which experimentalists measure how much a system apart from an initially delocalized or extended state [13, 20]. First, we briefly summarize some basic concepts related to the evolution in time of single particle problems. According to quantum mechanics, the time evolution of a ket $|\psi\rangle$ is given by the time dependent Schrödinger equation:

[TABLE]

Again, we expand the ket in the site basis $|\psi\rangle=\sum_{i}\psi_{i}|x_{i}\rangle$ and after substitution in the equation (32), we obtain the time dependent equation for the coefficients $\psi_{i}$

[TABLE]

where $H_{ij}$ are the matrix elements of $\mathbf{H}$ given in equality (29). Equation (33) represents a system of $\Omega$ coupled ordinary differential equations which can be easily solved by using the Runge-Kutta 4 (RK4) method:

[TABLE]

Where we have set $\tau=Jt/\hbar$ as our dimensionless unit of time. The advantage of using the RK4 method resides on one side in its accuracy $(\Delta\tau)^{4}$, and on the other side, in the relatively simple way in which the above equations can be implemented. In our calculations, we set $\Delta\tau=0.01$ which displays conservation of both, the norm and the energy within the whole numerical time evolution.

As it is well known, the time evolution of an eigenstate of $H$ would give trivial results. Nevertheless, the evolution of an arbitrary state can yield signatures of the presence of disorder. For this reason, we first study the evolution in time of an initially fully localized state in the middle of the lattice $|\psi(\tau=0)\rangle=|x_{0}\rangle$. This initial condition mimics ”designs ” prepared in current experiments performed with ultracold atoms [13]. In figure 9(a) and figure 9(b), we plot the initial density profile and the spreading of such initial state in the absence of disorder for a time of $\tau=100$ respectively.

In order to quantify the spreading of the initial wave function, we determine the root mean square of the displacement (RMSD) in each time step, this latter quantity is defined as

[TABLE]

In figure 10 we plot the RMSD, in logarithmic scale, as a function of the time $\tau$ for five different disorder strengths. The shaded area in each curve represents the standard error over $10$ realizations of random phase $\phi$. As one can observe, for zero disorder $\Delta/J=0$ the wave packet propagates ballistically, that is $\sigma\propto\tau$, showing that the RMSD grows linearly in time. Due to the system finite size, the RMSD reach a maximum value and oscillates around it. For this reason, we let the system evolve until the RMSD reaches its maximum value for zero disorder.

The time dependence of the RMSD can be fitted with a power law ansatz $\sigma\propto\tau^{\gamma}$. This fit must be done at intermediate times scales, where one neglects the contribution of the transient behavior at short times and the maximum spreading at later ones [25]. Figure 11 shows the exponent $\gamma$ of the above fit as a function of the disorder strength $\Delta/J$. This plot allows us to identify the ballistic regime $\gamma=1$, the superdiffusive $1/2<\gamma<1$, the subdiffusive $0<\gamma<1/2$, and the localized one $\gamma=0$, associated to the diffusion of an initially localized wave packet.

Another interesting observable that allows us to discern between a localized and extended phase, and can also be detected in current experiments [19, 20] is the imbalance $I(\tau)$. For a single particle problem, the imbalance parameter is defined as follows:

[TABLE]

where $n_{e}(\tau)=\sum_{i\in even}|\psi_{i}(\tau)|^{2}$ is the total probability density of the particle on even sites and $n_{o}(\tau)=\sum_{i\in odd}|\psi_{i}(\tau)|^{2}$ corresponds to the total probability density of the particle on odd sites of the lattice. In order to condense the whole time evolution of the imbalance, we defined $I$ as the asymptotic value of $I(\tau)$ [36]

[TABLE]

For calculation purposes, our numerical simulations were performed until $\tau=1000$, which corresponds to a real time of $t=1000\hbar/J$ that is much bigger than the hopping time in the lattice.

Besides being a measurable observable, the most significant advantage of using the imbalance parameter as an order parameter for probing localization is that, it can also provide signatures of many body localization [37], when interactions are present. To show the dynamical behavior of $I(\tau)$, we start by considering the evolution of a density-wave like pattern in which only even sites are initially occupied. To have meaningful calculations we have to impose the same number of odd and even sites in the full lattice. Here, we consider $\Omega=1000$ and the incommensurate parameter $\beta=987/1597$. Then, we calculate the value of the imbalance $I(\tau)$ at every time step until $\tau=1000$ at which we observe that the imbalance oscillates around its asymptotic value. For illustration purposes in figure 12, we plot the imbalance as a function of time $\tau$ for a disorder $\Delta/J=4.0$

In figure 13 we show the asymptotic value of the imbalance as a function of the disorder strength $\Delta/J$. Since for $\Delta/J<2$ all the states are extended, the particle can easily tunnel to nearby sites leading to a zero value of the imbalance in a few tunneling times. The vanishing of the imbalance must be interpreted as an indication of ergodicity since the system completely loses any previous information associated with the initial state. However, for $\Delta/J>2$ the imbalance reaches a finite value which is closer to the initial value as the disorder is increased. This suggests that the system is non-ergodic as it retains certain memory of the initial configuration.

6 Final Remarks

The main purpose of the manuscript is to introduce the undergraduate student to one of the most studied topics in condensed matter, the localization phenomenon. For this purpose, we have presented a comprehensive study of a single particle moving in a disordered lattice in one dimension. In particular, the disorder here analyzed corresponds to a quasiperiodic one. Considering as a starting point a quantum analog of such condensed matter system, namely a weakly interacting ultracold Bose gas confined in a 1D quasiperiodic lattice, a straightforward derivation of the Aubry-André Hamiltonian representing this type of disorder was presented. Then, we proceeded to review the essential results of the one-dimensional Aubry-André model, as well as the form in which localization transition can be recognized. First, we focused on the description of the Aubry-André model properties and then, to characterize the localization transition as a function of the disorder amplitude, we investigated stationary and dynamical properties. Particularly, we determine the inverse participation ratio (IPR) and its opposite, the normalized participation ratio (NPR) that provide information of how localized across the lattice a wave function is, as a function of the disorder amplitude. These parameters were calculated for both, the ground state and the full energy spectrum exhibiting that the spectrum is either, completely delocalized or completely localized, but not a mixture of both, localized and extended states. Another stationary property here studied was the energy spectrum that exhibits the classical Hofstadter pattern when the disorder amplitude reaches a critical value.

All the formalism and techniques used here are at the level of an advanced undergraduate student or equivalent. We believe that this kind of readings brings a student to become closer to the comprehension of current research on single particle and many particle localization phenomena. With the tools used in this paper, the interested reader can address the study of vanguard problems related to the central topic of this manuscript, for instance, the dependence of localization on dimensionality, the competition of disorder and interactions, the effects of next nearest neighbors in the localization phenomenon, and the response of the system in the novel driven disordered lattices among others. Indeed, the effects of disorder combined with strong interactions remains an outstanding challenge to the theory.

Finally, in this paragraph we want to briefly summarize some of the reported predictions of the Aubry-André model. These include spin-orbit coupling effects [42], closed expressions for the energy separating localized and non-localized states [43] and coexistence of localized and extended states in interacting quasiperiodic systems [44] among others. At the many body level, localization of the ground state established rigorously in the weakly interacting regime for both, repulsive and attractive interactions [45] and many-body localization vs. thermalization and onset of equilibrium [46], which can have implications for quantum devices and quantum computation. All of these phenomena can be analyzed with the same tools and properties here presented.

Acknowledgments

We acknowledge useful discussion with G. García Naumis. GADC acknowledges CONACYT scholarship. This work was partially funded by grant INI105217 DGAPA (UNAM) and 255573 CONACYT.