Many-Body Localization in a Quasiperiodic System
Recent theoretical and numerical evidence suggests that localization can survive in disordered many-body systems with very high energy density, provided that interactions are sufficiently weak. Stronger interactions can destroy localization, leading to a so-called many-body localization transition. This dynamical phase transition is relevant to questions of thermalization in extended quantum systems far from the zero-temperature limit. It separates a many-body localized phase, in which localization prevents transport and thermalization, from a conducting (“ergodic”) phase in which the usual assumptions of quantum statistical mechanics hold. Here, we present numerical evidence that many-body localization also occurs in models without disorder but rather a quasiperiodic potential. In one dimension, these systems already have a single-particle localization transition, and we show that this transition becomes a many-body localization transition upon the introduction of interactions. We also comment on possible relevance of our results to experimental studies of many-body dynamics of cold atoms and non-linear light in quasiperiodic potentials.
In one-dimensional systems of non-interacting particles, an arbitrarily weak disordered potential generically localizes all quantum eigenstates Anderson (1958); Abrahams et al. (1979). Such a system is always an insulator, with a vanishing conductivity in the thermodynamic limit. The question of how this picture is modified by interactions remained unclear in the decades following Anderson’s original work on localization Fleishman and Anderson (1980); Thouless and Kirkpatrick (1981). Relatively recently, Basko, Aleiner, and Altshuler have argued that an interacting many-body system can undergo a so-called many-body localization (MBL) transition in the presence of quenched disorder. At low energy density and/or strong disorder, interactions are insufficient to thermalize the system, so the system remains a “perfect” insulator (i.e. with zero DC conductivity despited being excited); at higher energy density and/or weaker disorder, the conductivity can become nonzero and the system thermalizes, leading to a conducting phase Basko et al. (2006, 2007a).
The MBL transition is rather unique for several reasons. First, in contrast to more conventional quantum phase transitionsSachdev (2007), this is not a transition in the ground state. Instead, the MBL transition involves the localization of highly excited states of a many-body system, with finite energy density. This means that the transition differs from most metal-insulator transitions, which are sharp only at zero temperatureDobrosavljevic (2012). Furthermore, this MBL transition is of fundamental interest in the context of statistical mechanics. Local subsystems of interacting, many-body systems are generically expected to equilibrate with their surroundings, with statistical properties of these subsystems reaching thermal values after sufficient time. Studies of how this occurs in quantum systems have led to the so-called eigenstate thermalization hypothesis (ETH), which states that individual eigenstates of the interacting quantum system already encode thermal distributions of local quantities Deutsch (1991); Srednicki (1994). However, the many-body localized phase provides an example of a situation in which the ETH is false, and the ergodic hypothesis of quantum statistical mechanics is violated Oganesyan and Huse (2007); Pal and Huse (2010). Since the work of Basko et al., these intriguing aspects of MBL have motivated many studies aimed at locating and understanding this transition in disordered systems Oganesyan and Huse (2007); Žnidarič et al. (2008); Karahalios et al. (2009); Oganesyan et al. (2009); Pal and Huse (2010); Monthus and Garel (2010); Berkelbach and Reichman (2010); Canovi et al. (2011); Bardarson et al. (2012); Vosk and Altman (2012); De Luca and Scardicchio (2013); Khatami et al. (2012); Biroli et al. (2012).
On the other hand, it is important to note that single-particle localization does not require disorder. In 1980, Aubry and André studied a 1D single-particle tight-binding model that omits disorder in favor of a potential that is periodic, but with a period that is incommensurate with the underlying latticeAubry and André (1980). Harper had studied a similar model much earlier, but he had focused on a special ratio of hopping to potential strengthHarper (1955). Aubry and André showed that this point actually lies at a localization transition. It separates a weak potential phase, where all single-particle eigenstates are extended, from a strong potential phase, where all eigenstates are localized. In the 1980s and 1990s, physicists continued to study this quasiperiodic localization transition for its own peculiarities and because it mimics the situation in disordered systems in , where there is also a single-particle localization transition Thouless and Niu (1983); Chaves and Satija (1997); Evangelou and Katsanos (1997); Abanov et al. (1998); Eilmes et al. (1999); Siebesma and Pietronero (1987); Hashimoto et al. (1992); Vidal et al. (1999). The AA model was also actively investigated in the mathematical physics literature, because it involves a Schrödinger operator (i.e. the “almost Mathieu” operator) with particularly rich spectral properties. The contributions of mathematical physicists put the initial work on Aubry and André on more rigorous footing Simon (1982); Bellissard et al. (1983); Jitomirskaya and Simon (1994); Jitomirskaya (1999). More recently, the AA model has been directly experimentally realized in cold atom experimentsFallani et al. (2007); Lucioni et al. (2011) and also in photonic waveguidesLahini et al. (2009). The possibility of engineering quasiperiodic systems in the laboratory has inspired new theoretical and numerical work aimed at understanding the localization properties of such systems and how they differ from those with true disorder Albert and Leboeuf (2010); Cestari et al. (2011); Nessi and Iucci (2011); Tezuka and García-García (2012); He et al. (2012); Ribeiro et al. (2012); Gramsch and Rigol (2012).
Statement of the problem and summary of the results
In this paper, we ask whether there can be a MBL transition in an interacting extension of the AA model. More concretely, suppose we begin with a half-filled, one-dimensional system of fermions or hardcore bosons in a particular randomly chosen many-body Fock state, with some sites occupied and others empty. Such a configuration of particles is typically far from the ground state of the system. Instead, by sampling the initial configuration uniformly at random (i.e. without regard to its energy content), we are actually working in the so-called infinite temperature limit. If the particles are allowed to hop and interact for a sufficiently long time, the standard expectation is that the system should thermalize: that is, all microscopic states that are consistent with conservation laws should become equally likely and local observables should thereby assume some thermal distributionRigol et al. (2008). Can this expectation be violated in the presence of a quasiperiodic potential? In other words, can the system fail to serve as a good heat bath for itself? If so, can this be traced to the persistence of localization even in the presence of interactions?
The answer to both of these questions appears to be “yes.” We use numerical simulations of unitary evolution
of a many-body quasiperiodic system to measure three kinds of observables in the limit of very late times: the correlation between the initial and time-evolved particle density profiles, the many-body participation ratio, and the Rényi entropy.
Our observations are consistent with the existence of two phases in the parameter space of our model that differ qualitatively in ergodicity.
At finite interparticle interaction strength and large hopping ,
there exists a phase in which the usual assumptions of statistical mechanics appear to hold. The initial state evolves into a superposition of a finite fraction of the
total number of possible configurations, and consequently, local observables approximately assume their thermal values. This is the many-body ergodic phase.
However, at small hopping , there is a phase in which particle transport away from the initial configuration is not strongly enhanced by interactions. The system
explores only an exponentially small fraction of configuration space, and local observables do not even approximately thermalize. This is the many-body localized phase.
Figure 1 presents a schematic illustration of the proposed phase diagram. Although interactions induce an expansion of the ergodic regime, the localized phase
survives at finite , and consequently, there is evidence for a quasiperiodic MBL transition
There has certainly been substantial previous work on localization in many-body quasiperiodic systems. For instance, Vidal et al.Vidal et al. (1999) adapted the approach of Giamarchi and SchulzGiamarchi and Schulz (1988) to study the effects of a perturbative quasiperiodic potential on the low-energy physics of interacting fermions in one dimension. Very recently, He et al.He et al. (2012) studied the ground state Bose glass to superfluid transition for hardcore bosons in a 1D quasiperiodic lattice. Our work differs fundamentally from these and many other studies precisely because it focuses on non-equilibrium behavior in the high-energy (infinite temperature) limit and argues that a localization transition can even occur in this regime.
Organization of the paper
We begin our study in Section II by introducing our interacting extension of the standard AA model. Since the MBL transition is a non-equilibrium phase transition, our goal is to follow the real-time dynamics. To simplify this task, we describe a method of modifying the dynamics of our model, such that numerical integration of the new dynamics is somewhat easier than the original problem. In Section III, we introduce the quantities that we measure in our simulations and present the numerical results. Then, in Section IV, we argue that our data points to the existence of many-body localized and many-body ergodic phases by proposing model late-time states for each of these regimes and comparing to the numerical results from Section III. Next, in Section V, we extract estimates for the phase boundary from our data, motivating the phase diagram in Figure 1. Finally, we conclude in Section VI by summarizing our results, drawing connections to theory and experiment, and suggesting avenues for future extensions of our work.
We relegate two exact diagonalization studies to Appendix A. First, we examine the impact of our modified dynamics upon the single-particle and many-body problems. Second, we study the many-body level statistics of the interacting model. We find evidence for a crossover between Poisson and Wigner-Dyson statistics, consistent with the usual expectation in the presence of a localization transition Shklovskii et al. (1993).
Ii Model and Methodology
In this section, we motivate and introduce our model and our numerical methodology for studying real-time dynamics.
ii.1 The “Parent” Model
We would like to consider one-dimensional lattice models of the following general form:
Here, is a fermion annihilation operator, and is the corresponding fermion number operator. The three terms in the Hamiltonian (1) then correspond to an on-site potential, nearest-neighbor hopping, and nearest-neighbor interaction respectively. For now, we leave the boundary conditions unspecified. In 1D, the Hamiltonians (1) for hardcore bosons and fermions differ only in the matrix elements describing hopping over the boundary. With open boundary conditions, the Hamiltonians (and consequently all properties of the spectra) are identical.
If we set in the Hamiltonian (1) and take to be genuinely disordered, we recover the non-interacting Anderson Hamiltonian. If we then turn on a finite , we obtain a model that is related to the spin models that have been studied in the context of MBLPal and Huse (2010); Bardarson et al. (2012). Alternatively, suppose we set again and take:
With a generic irrational wavenumber and an arbitrary offset , we obtain the non-interacting AA modelAubry and André (1980). For our purposes, we would like to use an incommensurate potential of the form (2), with and and left as tuning parameters to explore different phases of the model (1).
Before proceeding, we should briefly review what is known about the single-particle AA model. With periodic boundary conditions and , this model is self-dual Aubry and André (1980); Albert and Leboeuf (2010). The self-duality can be realized by switching to Fourier space () and then performing a rearrangement of the wavenumbers such that the real-space potential term looks like a nearest-neighbor hopping in Fourier space and vice versa. On a finite lattice of length with periodic boundary conditions, such a rearrangement is possible whenever the wavenumber of the potential such that and are mutually prime. The duality construction reveals that, if the AA model has a transition, it must occur at . In the thermodynamic limit, there is indeed a transition at this value for nearly all irrational wavenumbers Thouless and Niu (1983). When , all single-particle eigenstates are spatially extended, and by duality, localized in momentum space; when , all single-particle eigenstates are spatially localized, and by duality, extended in momentum space. Exactly at , the eigenstates are multifractal Siebesma and Pietronero (1987); Hashimoto et al. (1992). The spatially extended phase of the AA model is characterized by ballistic, not diffusive, transport Aubry and André (1980). Recently, Albert and Leboeuf have argued that localization in the AA model is a fundamentally more classical phenomenon than disorder-induced Anderson localization, and that the AA transition at is most simply viewed as the classical trapping that occurs when the maximum eigenvalue of the kinetic (or hopping) term crosses the amplitude of the incommensurate potentialAlbert and Leboeuf (2010).
ii.2 Numerical Methodology and Modification of the Quantum Dynamics
Probing the MBL transition necessarily involves studying highly excited states of the system, and this precludes the application of much of the extensive machinery that has been developed for investigating low-energy physics. Consequently, several studies of MBL have resorted to exact diagonalization or other methods involving similar numerical cost Oganesyan and Huse (2007); Pal and Huse (2010); Monthus and Garel (2010). We too use a numerical methodology that scales exponentially in the size of the system. However, in order to access longer evolution times in larger lattices, we introduce a modification of the quantum dynamics. This modification is inspired by a scheme used previously by two of us in a study of classical spin chainsOganesyan et al. (2009). There, at any given time, either the even spins in the chain were allowed to evolve under the influence of the odd spins or vice versa. This provided access to late times that would have been too difficult to access by direct integration of the standard classical equations of motion.
By analogy, we propose allowing hopping on each bond in turn. At any given time, the instantaneous Hamiltonian looks like:
We will specify the value of in Section II.C below, where we discuss our choice of boundary conditions. The state of the system is allowed to evolve under this Hamiltonian for a time , and this evolution can be implemented by applying the unitary operator:
One full time-step of duration consists of cycling through all the bonds:
Note that, in (3), the hopping is enhanced by because the hopping on any given bond is activated only once per cycle, while the potential and interaction terms always act. Therefore, the factor of ensures that the average Hamiltonian over a time has the form (1). The advantage of employing the modified dynamics is that the only couple pairs of configurations, so preparing the reduces to exponentiating order two-by-two matrices, where is the size of the Hilbert space. This is generally a simpler task than exponentiating the original Hamiltonian (1). Our scheme only constitutes a polynomial speedup over exact diagonalization, but that speedup can increase the range of accessible lattice sizes by a few sites.
The modified dynamics raise several important issues that should be discussedMaricq (1982). The periodic time-dependence of the Hamiltonian induces so-called “multi-photon” (or “energy umklapp”) transitions between states of the “parent” model (1) that differ in energy by , reducing energy conservation to quasienergy conservation modulo . We need to question whether this destroys the physics of interest: does the single-particle Aubry-André transition survive, or do the multi-photon processes destroy the localized phase?
We take up this question in Appendix A, where we present a Floquet analysis of the single-particle and many-body problems. We find that, for sufficiently small , the universal behavior of the single-particle AA model is preserved. At larger , multi-photon processes can strongly mix eigenstates of the single-particle parent model, increasing the single-particle density-of-states and destroying the AA transition. In the spirit of the earlier referenced work on classical spin chainsOganesyan et al. (2009), our perspective in this paper is to identify whether MBL can occur in a model qualitatively similar to our parent model (1). Therefore, to explore dynamics on long time scales, we avoid destroying the single-particle transition, but still choose to be quite large within that constraint.
In Appendix A, we also examine the consequences of our choice of for the quasienergy spectrum of the many-body model. Our results suggest that multi-photon processes do not, in fact, strongly modify the parent model’s spectrum for much of the parameter range that we explore in this paper
Finally, we note in passing that several recent studies have focused on the localization properties of time-dependent modelsKitagawa et al. (2012); Martinez and Molina (2006); DÕAlessio and Polkovnikov (2013), including one on the quasiperiodic Harper modelKolovsky and Mantica (2012), but that the intricate details of this topic are somewhat peripheral to our main focus.
ii.3 Details of the Numerical Calculations
In studies of the 1D AA model, it is conventional to approach the thermodynamic limit by choosing lattice sizes according to the Fibonacci series (, , , , ) and wavenumbers for the potential (2) as ratios of successive terms in the seriesThouless and Niu (1983). These values of respect periodic boundary conditions while converging to the inverse of the golden ratio . For any finite lattice, the potential is only commensurate with the entire lattice (since successive terms in the Fibonacci series are mutually prime), and the duality mapping of the AA model is always exactly preserved. For our purposes however, this approach offers too few accessible system sizes and complicates matters by generating odd values of .
Instead, we found empirically that finite-size effects are least problematic if we use exclusively even , always keep the wavenumber of the potential fixed at , and set:
in equation (3), thereby forbidding hopping over the boundary
Using the approach described above, we have simulated systems up to size at half-filling. Our simulations always begin with a randomly chosen configuration (or Fock) state, so that the initial state has no entanglement across any spatial bond in the lattice (i.e. each site is occupied or empty with probability 1). Except in the exact diagonalization studies of Appendix A, we always set . We integrate out to and ultimately average the evolution of measurable quantities over several samples, where a sample is specified by the choice of the initial configuration and offset phase to the potential (2). The sample counts used in the numerics are provided in Table 1.
Iii Numerical Measurements
We now introduce the quantities that we measure to characterize the different regimes of our model. We also present the numerical data along with some qualitative remarks about the observed behavior. However, we largely defer quantitative phenomenology and modeling of the data to Section IV.
iii.1 Temporal Autocorrelation Function
One signature of localization is the system’s retention of memory of its initial state. Since we simulate the reversible evolution of a closed system, the quantum state of the entire system retains full memory of its past. Nevertheless, we may still ask if the information needed to deduce the initial state is preserved locally or if it propagates to distant parts of the system. A diagnostic measure with which to pose this “local memory” question is the temporal autocorrelator of site :
Here, the angular brackets refer to an expectation value in the quantum state. This single-site autocorrelator may be averaged over sites and then over samples (as defined in Section II.C) to obtain:
The sample average is indicated here with the large square brackets. Typically, to reduce the effects of noise, we also average over a few time steps within each sample (i.e. perform time binning) before taking the sample average.
We can discriminate three qualitatively different behaviors of vs. in our interacting model. Figure 2 shows examples of each of these behaviors at interaction strength .
Panel (a) is characteristic of the low regime, where the autocorrelator stays invariant over several orders-of-magnitude of time, and there is no statistically significant difference between time series for
different . At higher , as in panel (b), the time series show approximately power-law decay culminating in saturation to a late-time asymptote. For the largest systems, the power law is roughly consistent
with the diffusive expectation of decay. The late-time asymptote decays with (as expected from energy conservation
The difference between these two regimes is brought out more clearly in Figure 3. We focus on a late time and probe as a function of for different lattice sizes. Panels (a)-(c) show data for , , and respectively. All the panels show a “splaying” point of the vs. curves, separating a high regime in which decays with from a low regime in which it does not. The value of at this feature decreases monotonically with . Most importantly, in each case, this value is robust to changing ; if we halve from the value that appears in Figure 3, the feature appears at approximately the same value of . This property of the data is very fortunate: in Section IV.C below, we will use the splaying feature in these plots to put a numerical lower bound on the transition value of for different interaction strengths. Since time scales get very long near the transition, it is difficult to simulate out to convergence in this regime. Nevertheless, the fact that the value of at the splaying feature remains fixed in time implies that we can deduce the phase structure from our finite-time observations.
iii.2 Normalized participation ratio
One of the commonly used diagnostics for studying single-particle localization is the inverse participation ratio (IPR). This quantity is intended to probe whether quantum states explore the entire volume of the system and is often defined as the sum over sites of the amplitude of the state to the fourth power: . Typically, the IPR is inversely proportional to the localization volume in a single-particle localized phase and decays to zero as the inverse of the system volume in an extended phase.
We now describe how this quantity can be fruitfully exploited in the many-body context. Let denote some specific configuration of particles in sites. Then, we can write the state of the system in the configuration basis as:
The configuration-basis IPR is simply:
where the square brackets, as usual, denote a sample average. Interpreting as the inverse of the number of configurations on which has support, we now define the normalized participation ratio (NPR):
The quantity then represents the fraction of configuration space that the system explores. We expect to be independent of at late times in the ergodic phase. In the many-body localized phase, we expect to decay exponentially with .
In Figure 4, we plot vs. for , , and . The figure reveals an important difference between the non-interacting and interacting models. At low , both with and without interactions, decays exponentially with :
with . More surprisingly, also decays with at large in the non-interacting case; all that happens is that becomes essentially independent of . With even small interactions however, becomes system-size independent in the large regime, following our ansatz for an ergodic phase. We bring out this point more clearly in Figure 5, in which we extract estimates for the decay coefficient for various values of the interaction strength. Thus, the extended phase of the non-interacting AA model appears to be a special, non-ergodic limit.
Before proceeding, we should caution that, in panels (b) and (c) of Figure 4, the collapse at high looks very appealing because of the use of a semilog plot and would not be so striking on a normal scale. The axes have been chosen to highlight the exponential scaling at low , which would not be as apparent if we simply plotted vs. . However, regarding the absence of perfect collapse at high , note that the raw data for the IPR differ by several orders-of-magnitude for different values of the lattice size . Given this, the coincidence of the order-of-magnitude of for different values of is already a good indication of the proposed scaling, and some corrections to this scaling should be expected given the modest accessible system sizes.
iii.3 Rényi Entanglement Entropy
Unlike the normalized participation ratio, which provides a global characterization of the time-evolved state, bipartite entanglement is arguably a better proxy for whether a part of the system can act as a good heat bath for the rest. In the many-body ergodic phase, we expect the bipartite entanglement entropy to be a faithful reflection of the thermodynamic entropy. This implies an extensive entropy, pinned to its thermal infinite temperature value throughout the phase
Let subsystem A refer to lattice sites , and let subsystem B refer to the remaining sites in the chain. We can compute the reduced density matrix of subsystem A by beginning with the full density matrix and tracing out the degrees of freedom associated with subsystem B:
The sample-averaged order-2 Rényi entropy of subsystem A is then given by:
Both and the standard von Neumann entropy are expected to attain the same values in the ergodic phase; we choose to focus on the former to save on the computational cost of diagonalizing the reduced density matrix (13).
Our first task is to examine whether the putative localized phase of our model exhibits the same behavior that was observed with tDMRG Žnidarič et al. (2008); Bardarson et al. (2012). In panel (a) of Figure 6, we focus on a low value of and plot vs. for lattices. At very early times, the time series all tend to coincide, reflecting the formation of short-range entanglement at the cut between the subsystems. Afterwards, the non-interacting time series saturates for several orders-of-magnitude of time, while the interacting time series show behavior that is consistent with logarithmic growth. In order to clearly establish the saturation that follows the slow growth, we have had to focus on small lattices. Panel (b) of Figure 6 shows data for large . Here, the most striking difference between the non-interacting and interacting models lies in the saturation value of the entropy: the interacting model is substantially more entangled, but the saturation value does not appear to depend on the value of . We will see below that this is another indication that thermalization only occurs in the interacting, large regime.
Figure 7 shows late-time values of the Rényi entropy density plotted against the tuning parameter . We first focus on the high regime. In panel (a), , and for large . However, the entropy density is less than , which is the thermal result when the system has ergodic access to all configurations consistent with particle number conservation. The situation is dramatically different in panels (b) and (c), where and respectively. At high , the entropy actually looks superextensive. This is just a finite-size effect, because the entropy is well fit to a linear growth of the form:
where is a constant deficit, typically around . In Figure 8, we show that the slope at large in the interacting problem. This implies that the entropy is thermal in the limit, where the deficit is negligible.
Now, we turn to the low regime. Without interactions, the off-diagonal elements in the reduced density matrix (13) typically contain only a few frequencies originating from localized single particle orbitals immediately adjacent to the cut. The number of relevant orbitals is finite in . As a result, the off-diagonal elements cannot fully vanish, and the reduced density matrix never thermalizes. The resulting entanglement entropy is independent of as shown in the inset of panel (a). In the interacting problem, while the orbitals immediately adjacent to the cut still have roughly the same frequencies, the “spectral drift” (i.e. the spread of these lines due to sensitivity to the configuration of distant particles) allows for a much larger number of distinct and mutually incoherent contributions to offdiagonal elements of the reduced density matrix. These off-diagonal elements can dephase more efficiently, leading to a partial thermalization. This is the mechanism that likely underlies the extensive but subthermal entropy observed by Bardarson et al.Bardarson et al. (2012). For small , our numerical results in the low regime agree well with this expectation. For larger , the slow dynamics of the entropy formation makes it difficult to observe saturation, both in our work and in the tDMRG study of Bardarson et al.
If the entropy eventually becomes extensive for all , then the “crossing” feature that is present in panels (b) and (c) of Figure 7 would become a “splaying” feature, with the entropy density independent of at small . In any case, an interesting property of the data is that the values of at the crossing features of the vs. plots are consistent with the locations of the splaying features in the corresponding vs. plots of Figure 3. This seems to be the case for all . Thus, these features may be useful in locating the transition.
Iv Modeling the Many-Body Ergodic and Localized Phases
Above, we presented numerical evidence that our interacting AA model contains two regimes that show qualitatively distinct behavior of the autocorrelator, normalized participation ratio, and Rényi entropy. Next, we will propose and characterize model quantum states that qualitatively (and sometimes quantitatively) reproduce the numerically observed late-time behavior in the two regimes. These model states expose more clearly why the two regimes of our model are appropriately identified as many-body ergodic and localized phases.
iv.1 The Many-Body Ergodic Phase
To model the behavior of the putative ergodic phase, we begin by writing down a generic model state in the configuration basis:
Here, the refer to configurations of the full chain, whereas the and refer to configurations of the subsystems A and B, as defined in Section III.C above. The superscripts on the configurations and expansion coefficients refer to the number of particles in subsystem A. Writing the state in terms of the subsystem configurations will be useful shortly, but for now we focus on the statistical properties of the amplitude . We assume that this amplitude is distributed as a complex Gaussian random variable:
Within this distribution, and . From these average values, it is possible to deduce that:
for normalization and that the IPR is . This, in turn, implies:
This result is reproduced quantitatively in the numerics in Figure 4.
Next, suppose we compute the reduced density matrix of subsystem A in the state :
To find the Rényi entropy, we need to compute the trace of the square of this operator:
When we average over our distribution of amplitudes (17), only the coherent terms survive
The final term accounts for the double counting of terms where and simultaneously. We now introduce the notation:
and evaluate the expectation values in equation (21) to obtain:
Finally, using a Stirling approximation to the combination function and a saddle-point approximation for the sum, we find the entropy:
This is the same form observed in the numerics (15), and the deficit lies in the observed range. Asymptotically in , the entropy (25) is maximal, and this is exactly the expected behavior when the particle number thermalizes.
There is an important caveat to note here: we have argued above that, if multi-photon processes do not completely destroy energy conservation, then this can lead to relic autocorrelations at late times. This implies that the assumption of independent random amplitudes cannot be exactly correct on a finite lattice. However, the numerically-observed relic autocorrelations decay with , suggesting that our assumptions get better as the system size grows. Therefore, in the thermodynamic limit, this phase is truly thermal.
iv.2 The Many-Body Localized Phase
Our model for the time-evolved state in the localized regime is founded upon the following intuition: there exists a length scale , which is analogous to the single-particle localization length and beyond which particles are unlikely to stray from their positions in the initial state. Then, if we partition our lattice of length into blocks of size , exchange of particles between blocks is less important than rearrangements of the particles within each block. Consequently, the total number of configurations accessed by the state of the full system is approximately the product of the number of configurations accessed within each block. This multiplicative assumption should be very safe in a localized phase. We additionally assume that, within each block, the dynamics completely scramble the particle configuration. If a certain block of length contains particles in the initial state, then the time-evolved state contains equal amplitude for each of the possible ways of arranging particles in those sites. In keeping with our numerical protocol, we randomly select the initial state from the space of all possible Fock states of a certain global particle number. Then, a block of sites contains particles with probability:
We will consider the limit , where we can approximate the probability by the first term. The assumptions proposed above motivate writing down a state of the form:
where the tilde on the sum indicates that it should only run over configurations that are consistent with the initial distribution of particles among the blocks. The factors are complex phases which depend upon the configuration, and is a normalization which is equal to the total number of configurations represented in the state .
Before beginning our analysis of the state , we should note that, in contrast to our calculations in the ergodic phase, our goal in the localized regime will be to qualitatively tie the numerically observed large behavior to the existence of the length scale . Unfortunately, we cannot achieve the quantitative accuracy of the ergodic model state with the localized toy-model described above.
We begin by estimating the autocorrelator between the initial state and the model time-evolved state . A non-zero autocorrelator emerges, because each block is only at half-filling on average. Fluctuations away from half-filling (in either direction) yield a positive typical value of the autocorrelator within a block. Indicating an average over the distribution (26) with an overline, we find the block value . This is also the average value for the whole system when :
Next, to estimate the IPR, we need to compute the normalization factor . We begin by estimating the number of explored configurations in each block. The average of the logarithm of the number of explored configurations within a block is:
Then, using , we can estimate itself as:
Using this normalization, we can estimate the NPR :
This qualitatively agrees with the numerically observed behavior (12) up to subleading corrections, and in the large- limit:
Note that equations (28) and (32) imply a relationship between the scaling behaviors of and in the localized regime. This relationship is not reflected in our numerical data, in part because we cannot truly attain the limit . The numerically computed value of , for example, can contain finite-size corrections of order or . Also, we must keep in mind that the state is just a toy model that does not capture fine details of the time-evolved states in this regime. Thus, we must be content with reproducing the qualitative behavior of each measurable quantity individually, without expecting the relationships between these quantities in to be exactly reproduced in the data.
We now turn to the Rényi entropy, the quantity which most strikingly distinguishes between the non-interacting and interacting localized phases. To examine this quantity, we revert to partitioning the system in half, instead of into blocks of size . As long as , the assumptions that we made above about the blocks of size hold even better for the subsystems A and B. For example, we can assume that there are “explored sets” of configurations in subsystem A and configurations in subsystem B respectively, with . We consider computing the reduced density matrix , exactly as in equation (20) above. If the off-diagonal elements of this density matrix remain perfectly phase-coherent, it can easily be shown that . In reality, there will be a local contribution to the entropy from particles straying over the cut between subsystems A and B. This mimics the situation in non-interacting localized phases. Alternatively, suppose that dephasing is sufficiently strong that we can proceed by analogy with the ergodic phase, beginning with equation (21) and keeping only coherent terms as in equation (22). Thereafter, the calculation for the model localized state differs from the calculation for . We need to consider the statistics of the configuration probabilities . For , we need the configurations on both subsystems to lie within their respective explored sets; this occurs in subsystem A, for example, with probability . This reasoning leads to the “dephased” entropy:
where we have additionally made the approximation that typically . With only partial loss of coherence, the entropy would lie between these two limiting cases: . Thus, dephasing alone, without additional particle transport, can induce an extensive entropy.
Indeed, our numerics support the view that the main difference between the non-interacting and many-body localized phases is the amount of dephasing. There does not seem to be a qualitative difference in particle transport. The particle configuration stays trapped near its initial state, even with interactions, and the system does not thermalize.
V Tracing the Phase Boundary
in this section, we use the data from Section III to extract estimates of the phase boundary between the localized and ergodic phases. Estimating the location of the MBL transition is extremely challenging. Given the numerically accessible lattice sizes, satisfying finite-size scaling analyses are difficult to perform. Nevertheless, rough estimates have been made in the disordered problem Oganesyan and Huse (2007); Monthus and Garel (2010); Pal and Huse (2010); De Luca and Scardicchio (2013), so we will now attempt to extract an approximate phase boundary for our model.
We first consider the autocorrelator. Above, we noted the “splaying” feature in the late-time plots of the autocorrelator vs. . The value of at this feature can be taken as a lower bound for the transition. For slightly greater than this value, it is possible that only decays with because for accessible lattice sizes. Considering two lattice sizes ( and ) and finding when their values of deviate, we find the values reported in the first column of Table 2.
Next, we consider the fitting parameter in equation (12). In Figure 5, we see that there is a region where for finite interaction strength. Since , finite-size effects are clearly dominating the estimate in this region. We can use the value of where is minimal to track how this region moves as is varied. This yields the second column of the table.
Finally, a similar approach can be applied to extract estimates of from the fits (15). There exists a region where , but this is mathematically inconsistent in the thermodynamic limit. Therefore, if we find the value of that maximizes , we can again estimate the location of the region dominated by finite-size effects, yielding the final column of Table 2.
The estimates of the transition value in Table 2 were obtained using data for the latest time that we simulated (the time bin for and and for ). However, we have also estimated for data obtained at a half and a quarter of this integration time, finding consistent results. Thus, the general phase structure of the model is invariant to changing the observation time, even though not all measurable quantities have converged to their asymptotic values. Consolidating the information from the estimates in Table 2, we propose that the phase diagram qualitatively resembles Figure 1.
Before proceeding, it is worth noting that our rough estimates of the phase boundary do not make assumptions regarding the character of the MBL transition (i.e. whether it is continuous or first order). In fact, some of our plots (e.g. panel (c) of Figure 7) hint at the possibility of a discontinuous change in as a function of in the thermodynamic limit. We are not aware of any results that rule out a first-order MBL transition, so we must keep this possibility in mind.
Recently, evidence has accumulated that Anderson localization can survive the introduction of sufficiently weak interparticle interactions, giving rise to a many-body localization transition in disordered systems Basko et al. (2006, 2007a); Oganesyan and Huse (2007); Pal and Huse (2010); De Luca and Scardicchio (2013). The MBL transition appears to be a thermalization transition: in the proposed many-body localized phase, the fundamental assumption of statistical mechanics breaks down, and the system fails to serve as its own heat bath Oganesyan and Huse (2007); Pal and Huse (2010). We have presented numerical evidence that this type of transition can also occur in systems lacking true disorder if they instead exhibit “pseudodisorder” in the form of a quasiperiodic potential.
From one perspective, this may be an unsurprising claim. For the localized single-particle eigenstates of the quasiperiodic Aubry-André model have the same qualitative structure as those of the Anderson model, so the effects of introducing interactions ought to be similar. By this reasoning, perhaps it is even possible to guess the phase structure of an interacting AA model using knowledge of an interacting Anderson model: we simply match lines of the two phase diagrams that correspond to the same non-interacting, single-particle localization length.
However, this perspective misses important effects in all regions of the phase diagram. Most obviously, the AA model has a transition at , and it is interesting to see how this transition gets modified as it presumably evolves into the MBL transition at finite . It is also important to remember that quasiperiodic potentials are completely spatially correlated. This means that the AA model lacks rare-regions (Griffiths) effects, and this may have subtle consequences for the dynamics. Finally, the AA model contains a phase that is absent in the one-dimensional Anderson model, the extended phase, and we have seen above that interactions have a profound effect upon this regime.
Understanding MBL in the quasiperiodic context is especially pertinent given the current experimental situation. Some experiments that probe localization physics in cold atom systems use quasiperiodic potentials, constructed from the superposition of incommensurate optical lattices, in place of genuine disorder. The group of Inguscio, in particular, has recently explored particle transport for interacting bosons within this setup Fallani et al. (2007); Lucioni et al. (2011). Meanwhile, the AA model has also been realized in photonic waveguides, and experimentalists have studied the effects of weak interactions on light propagation through these systems. They have also investigated “quantum walks” of two interacting photons in disordered waveguides Lahini et al. (2009, 2010). This protocol resembles the one we have implemented numerically, so similar physics may arise. Finally, we note that Basko et al. have predicted experimental manifestations of MBL in solid-state materials. In such systems, there is always coupling to a phononic bath, so the MBL transition is expected to become a crossover that nevertheless retains interesting manifestations of the MBL phenomenaBasko et al. (2007b). Whether there exist quasiperiodic solid-state systems to which the predictions of Basko et al. apply remains to be understood.
Given the current experimental relevance of localization phenomena in quasiperiodic systems, we hope that our study will motivate further attempts to understand these issues. Unfortunately, our ability to definitively identify and analyze the MBL transition is limited by the modest lattice sizes and evolution times that we can simulate. Vosk and Altman recently developed a strong-disorder renormalization group for dynamics in the disordered problemVosk and Altman (2012), but the reliability of such an approach in the quasiperiodic context is unclear. A time-dependent density matrix renormalization (tDMRG) group study of this problem would be a valuable next step. Tezuka and García-García have published tDMRG results on localization in an interacting AA model, but their focus was not on the thermalization questions of many-body localizationTezuka and García-García (2012). It would be worthwhile to pose these questions using a methodology that allows access to much larger lattices. However, even tDMRG may have difficulty capturing the highly-entangled ergodic phase Žnidarič et al. (2008); Bardarson et al. (2012), so an effective numerical approach for definitively characterizing the transition remains elusive.
Acknowledgements.We thank E. Altman, M. Babadi, E. Berg, S.-B. Chung, K. Damle, D. Fisher, M. Haque, Y. Lahini, A. Lazarides, M. Moeckel, J. Moore, A. Pal, S. Parameswaran, D. Pekker, S. Raghu, A. Rey and J. Simon for helpful discussions. This research was supported, in part, by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grants CNS-0855217 and CNS-0958379. S.I. thanks the organizers of the 2010 Boulder School for Condensed Matter and Materials Physics. S.I. and V.O. thank the organizers of the Cargesè School on Disordered Systems. S.I. and G.R. acknowledge the hospitality of the Free University of Berlin. V.O. and D.A.H are grateful to KITP (Santa Barbara), where this research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. V.O. thanks NSF for support through award DMR-0955714, and also CNRS and Institute Henri Poincaré (Paris, France) for hospitality. D.A.H. thanks NSF for support through award DMR-0819860.
Appendix A Exact Diagonalization Results for the Single-Particle and Many-Body Problems
This appendix collects exact diagonalization results that supplement the real-time dynamics study in the main body of the paper.
a.1 Floquet Analysis of the Modified Dynamics
The goal of the first part of this appendix is to examine the consequences of the modifications to the quantum dynamics described in Section II.B above. We first verify that the AA transition survives by diagonalizing the single-particle AA Hamiltonian (i.e. the Hamiltonian (1) with ) and the single-particle unitary evolution operators (5) for various choices of the time step . Subsequently, we employ the same approach to examine how varying impacts the quasienergy spectrum of the interacting, many-body model.
Robustness of the Single-Particle Aubry-André Transition
To study the single-particle transition, we focus on the inverse participation ratio:
Here, denotes the amplitude of the wave function at site of an site lattice. We enclose the sum in equation (34) in parentheses to indicate important differences in the averaging procedure with respect to the many-body inverse participation ratio (10). In the many-body case, we computed the IPR as a sum over configurations in the quantum state at a particular time in the real-time evolution. Then, we averaged over samples, where a sample was specified by a choice of the offset phase to the potential (2) and an initial configuration. Throughout this appendix, we instead specify a “sample” solely by the offset phase , and we average over eigenstates within each sample before averaging over samples.
As noted previously, the usual AA model has a transition that must occur, by duality, at . Near the transition, the localization length is known to diverge with exponent Thouless and Niu (1983). Our exact diagonalization results indicate that, at the transition, . Hence, we can make the following scaling hypothesis for the IPR:
In panel (a) of Figure 9, we show that we can use this scaling hypothesis to collapse data for the standard AA model. We show data for to , with potential wavenumber and open boundary conditions. For all lattice sizes, we average over samples.
To establish the stability of the AA transition to the modified dynamics, we must ask: can the IPR obtained from diagonalizing the unitary evolution operators (5) be described using the scaling hypothesis (35)? Panel (b) of Figure 9 shows that this is indeed the case for . The only parameter that needs to be changed is , which decreases slightly as is raised. This implies that there is a transition in the Floquet spectrum of the system that can be tuned by varying . It would be a worthwhile exercise to map out the phase diagram of this single-particle problem in the plane. We leave this for future work.
Properties of the Many-Body Quasienergy Spectrum
We now turn our attention back to the effects of the modified dynamics upon the full, many-body model. In Section II.B above, we emphasized that our time-dependent model lacks energy conservation, with multi-photon processes inducing transitions between states of the parent model (1) that differ in energy by . In this part of the appendix, we will examine how varying impacts the quasienergy spectrum of the time-dependent model, using the approach that we applied to the single-particle case above: we diagonalize the time-independent Hamiltonian as well as the unitary evolution operator for one time step of the time-dependent model.
In Figure 10, we plot the density-of-states in quasienergy space of the parent model and time-dependent models for different values of . We focus on systems at half-filling with fermions (or, since we continue to use the boundary conditions described in Section II.C, hardcore bosons). We fix the interaction strength to and tune to explore different regimes of the model. In panels (a)-(c), we plot data for , , and . According to Table 2, these values of put the system in the localized phase, near the transition, and in the ergodic phase respectively.
We first consider the consequences of varying while holding the other parameters fixed. For sufficiently small , the quasienergy spectrum faithfully reproduces all the structure of the energy spectrum of the parent model. This is unsurprising, because if is greater than the bandwidth of the parent model’s spectrum, direct multi-photon processes will not take place. If we now tune so that it is less than this bandwidth, the quasienergy spectrum begins to deviate from the parent model’s spectrum at its edges. This effect can be seen, for instance, by examining the trace for in panels (a) or (b). For even higher values of (i.e. lower values of ), multi-photon processes strongly mix the states of the parent model, resulting in a flat quasienergy spectrum.
The effect of multi-photon processes can also be enhanced by broadening the parent model’s spectrum, which can be achieved by raising or . In panel (c) of Figure 10 for instance, multi-photon processes have significantly flattened the spectrum for , and deviations from the parent model are even visible for . Since we always use in our real-time dynamics simulations, it is perhaps fortunate that is well within the proposed ergodic phase for and that, near the critical point (i.e. in panel (b)), the quasienergy spectrum for still retains much of the structure of the parent model’s spectrum.
However, there is one more caveat to keep in mind: the energy content of the system also grows with . At fixed , , and , the properties of the parent and time-dependent models deviate from one another as the system size grows. If we truly want to faithfully reproduce the dynamics of the parent model with the modified dynamics, it may be necessary to scale down as we raise . However, recall that our goal is simply to find MBL in a model qualitatively similar to the parent model (1). Even with this more modest goal in mind, there is still the danger that, on sufficiently large lattices, multi-photon processes might couple a very large number of localized states and thereby destroy the many-body localized phase of the parent model. Our numerical observations indicate that this does not happen for the system sizes that we can simulate. We can keep fixed at unity for without issues, accepting the possibility that the sequence of models that we would in principle simulate on still larger lattices may require progressively smaller values of .
a.2 Level Statistics of the Many-Body Parent Model
Localization transitions are often characterized by transitions in the level statistics of the energy spectrumShklovskii et al. (1993). Two of us previously looked at the level statistics of the disordered problem and identified a crossover from Poisson statistics in the many-body localized phase to Wigner-Dyson statistics in the many-body ergodic phaseOganesyan and Huse (2007). The intuition that underlies this crossover is the following: in a localized phase, particle configurations that have similar potential energy are too far apart in configuration space to be efficiently mixed by the kinetic energy term in the Hamiltonian. Therefore, level repulsion is strongly suppressed, and Poisson statistics hold. Conversely, in an ergodic phase, there is strong level repulsion which lifts degeneracies, leading to Wigner-Dyson (i.e. random matrix) statistics.
Along the lines of the aforementioned study of the disordered problem, we focus on the gaps between successive eigenstates of the spectrum of the many-body parent model (1):
and a dimensionless parameter that captures the correlations between successive gaps in the spectrum:
For a Poisson spectrum, the are distributed as with mean ; meanwhile, when random matrix statistics hold, the mean value of has been numerically determined to be approximately Oganesyan and Huse (2007).
In Figure 11, we present exact diagonalization results for lattices at half-filling with potential wavenumber and the boundary conditions described in Section II.C above. We show data for the same parameter range examined in the body of this paper and average over samples for each value of and . For the largest value of , the mean value of interpolates between the expected values as is raised, consistent with the existence of a localization transition. We have also checked that the distributions of have the expected forms in the small and large limits in this regime. For smaller values of , we can speculate that grows with at large and approaches the expected value for very large . To argue for a MBL transition on the basis of exact diagonalization, we would need to study this sharpening of the crossover as is raised. This would indeed be an interesting avenue for future work. For our present purposes however, we only want to check consistency with our real-time dynamics data, as we have done in Figure 11.
- Both the many-body ergodic and localized phases differ qualitatively from their counterparts in the non-interacting AA model. The non-interacting extended phase is not ergodic, indicating that interactions are necessary for thermalization. Meanwhile, the many-body localized phase is expected to exhibit logarithmic growth of the bipartite entanglement entropy to an extensive value, albeit with subthermal entropy density. Such behavior is in fact consistent with the recent observations in the disordered problemŽnidarič et al. (2008); Bardarson et al. (2012). This growth is absent in the AA localized phase without interactions. Despite this difference, the interacting and non-interacting localized phases are similar in their inability to thermalize the particle density.
- There is an exception to this statement: multi-photon processes do seem to play an important role deep in the ergodic phase, where the energy content of the system is especially high. See Appendix A and the discussion of the time-dependence of the autocorrelator in Section III.A for more details.
- To appropriately realize open boundary conditions, we should also turn off interactions over the boundary. When exploring different options for the boundary conditions, we varied over the boundary and neglected to vary . This is unfortunate in that it makes the model somewhat stranger. However, our boundary conditions are chosen for convenience anyway, and the numerics suggest that the choice of boundary conditions does not impact the essential physics discussed in this paper.
- The statistical fluctuation of the total energy of the randomly chosen initial configuration is of order . Suppose the total energy is conserved by the dynamics. We can write . Here, the subscripts and refer to the initial and late-time states, and are bounded random numbers capturing the expectation value of interactions (and hopping at late-times), is the non-random amplitude of the quasiperiodic potential, and and are positive bounded amplitudes of the Fourier components at the wavevector of the quasiperiodic potential. This ansatz implies a finite correlation between the random phases and . Therefore, one of the Fourier modes of remains correlated as , and we expect in the ergodic phase. Note that this argument truly applies only to the energy-conserving parent model. In fact, in our numerics, there is only partial energy conservation, and energy non-conserving events become more prevalent as , , or is raised. This means that will generically decay faster than at large in the ergodic phase.
- This statement should be interpreted with some care. Quantum entanglement entropy measures, such as the Rényi entropy that we define in equation (14), carry information about the off-diagonal elements in the reduced density matrix. These terms have no classical analogue and would not be considered in a thermodynamic calculation. This difference can result in discrepancies in the subleading behavior. For instance, consider our calculation of the bipartite Rényi entropy of the model state in Section IV.A: the quantum Rényi entropy is one bit lower than the Rényi entropy calculated by classical counting of configurations. A more precise analogue of the classical entropy would thus be a “diagonal” entropy in which all off-diagonal elements of the reduced density matrix were neglected.
- Only the first term on the right-hand side of equation (22) would appear in a “classical counting” derivation of the thermodynamic entropy. The other two terms account for off-diagonal elements in the reduced density matrix (20). Please see footnote 53 for more details.
- P. Anderson, Physical Review 109, 1492 (1958).
- E. Abrahams, P. Anderson, D. Licciardello, and T. Ramakrishnan, Physical Review Letters 42, 673 (1979).
- L. Fleishman and P. Anderson, Physical Review B 21, 2366 (1980).
- D. Thouless and S. Kirkpatrick, Journal of Physics C: Solid State Physics 14, 235 (1981).
- D. Basko, I. Aleiner, and B. Altshuler, Annals of physics 321, 1126 (2006).
- D. Basko, I. Aleiner, and B. Altshuler, Problems of Condensed Matter Physics 1, 50 (2007a).
- S. Sachdev, Quantum phase transitions (Wiley Online Library, 2007).
- V. Dobrosavljevic, Conductor Insulator Quantum Phase Transitions p. 1 (2012).
- J. Deutsch, Physical Review A 43, 2046 (1991).
- M. Srednicki, Physical Review E 50, 888 (1994).
- V. Oganesyan and D. Huse, Physical Review B 75, 155111 (2007).
- A. Pal and D. Huse, Physical Review B 82, 174411 (2010).
- M. Žnidarič, T. c. v. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
- A. Karahalios, A. Metavitsiadis, X. Zotos, A. Gorczyca, and P. Prelovšek, Phys. Rev. B 79, 024425 (2009).
- V. Oganesyan, A. Pal, and D. A. Huse, Physical Review B 80, 115104 (2009).
- C. Monthus and T. Garel, Physical Review B 81, 134202 (2010).
- T. Berkelbach and D. Reichman, Physical Review B 81, 224429 (2010).
- E. Canovi, D. Rossini, R. Fazio, G. Santoro, and A. Silva, Physical Review B 83, 094431 (2011).
- J. H. Bardarson, F. Pollmann, and J. E. Moore, Physical Review Letters 109, 17202 (2012).
- R. Vosk and E. Altman, Arxiv preprint arXiv:1205.0026 (2012).
- A. De Luca and A. Scardicchio, EPL 101, 37003 (2013).
- E. Khatami, M. Rigol, A. Relaño, and A. Garcia-Garcia, Physical Review E 85, 050102 (2012).
- G. Biroli, A. Ribeiro-Teixeira, and M. Tarzia, ArXiv preprint arXiv:1211.7334 (2012).
- S. Aubry and G. André, Ann. Israel Phys. Soc 3, 1 (1980).
- P. Harper, Proceedings of the Physical Society. Section A 68, 874 (1955).
- D. Thouless and Q. Niu, Journal of Physics A: Mathematical and General 16, 1911 (1983).
- J. Chaves and I. Satija, Physical Review B 55, 14076 (1997).
- S. N. Evangelou and D. E. Katsanos, Physical Review B 56, 12797 (1997).
- A. G. Abanov, J. C. Talstra, and P. B. Wiegmann, Physical Review Letters 81, 2112 (1998).
- A. Eilmes, U. Grimm, R. Roemer, and M. Schreiber, The European Physical Journal B-Condensed Matter and Complex Systems 8, 547 (1999).
- A. Siebesma and L. Pietronero, EPL (Europhysics Letters) 4, 597 (1987).
- Y. Hashimoto, K. Niizeki, and Y. Okabe, Journal of Physics A: Mathematical and General 25, 5211 (1992).
- J. Vidal, D. Mouhanna, and T. Giamarchi, Physical Review Letters 83, 3908 (1999).
- B. Simon, Advances in Applied Mathematics 3, 463 (1982).
- J. Bellissard, R. Lima, and D. Testard, Communications in Mathematical Physics 88, 207 (1983).
- S. Jitomirskaya and B. Simon, Communications in Mathematical Physics 165, 201 (1994).
- S. Jitomirskaya, Annals of Mathematics-Second Series 150, 1159 (1999).
- L. Fallani, J. Lye, V. Guarrera, C. Fort, and M. Inguscio, Physical Review Letters 98, 130404 (2007).
- E. Lucioni, B. Deissler, L. Tanzi, G. Roati, M. Zaccanti, M. Modugno, M. Larcher, F. Dalfovo, M. Inguscio, and G. Modugno, Physical Review Letters 106, 230403 (2011).
- Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, N. Davidson, and Y. Silberberg, Physical Review Letters 103, 13901 (2009).
- M. Albert and P. Leboeuf, Physical Review A 81, 013614 (2010).
- J. Cestari, A. Foerster, M. Gusmão, and M. Continentino, Physical Review A 84, 055601 (2011).
- N. Nessi and A. Iucci, Physical Review A 84, 063614 (2011).
- M. Tezuka and A. M. García-García, Physical Review A 85, 031602 (2012).
- K. He, I. I. Satija, C. W. Clark, A. M. Rey, and M. Rigol, Physical Review A 85, 013617 (2012).
- P. Ribeiro, M. Haque, and A. Lazarides, ArXiv preprint arXiv:1211.6012 (2012).
- C. Gramsch and M. Rigol, Physical Review A 86, 053615 (2012).
- M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
- T. Giamarchi and H. Schulz, Physical Review B. 37, 325 (1988).
- B. Shklovskii, B. Shapiro, B. Sears, P. Lambrianides, and H. Shore, Physical Review B 47, 11487 (1993).
- M. Maricq, Physical Review B 25, 6622 (1982).
- T. Kitagawa, T. Oka, and E. Demler, Annals of Physics (2012).
- D. Martinez and R. Molina, The European Physical Journal B-Condensed Matter and Complex Systems 52, 281 (2006).
- L. DÕAlessio and A. Polkovnikov, Annals of Physics (2013).
- A. R. Kolovsky and G. Mantica, Physical Review B 86, 054306 (2012).
- Y. Lahini, Y. Bromberg, D. Christodoulides, and Y. Silberberg, Physical Review Letters 105, 163905 (2010).
- D. Basko, I. Aleiner, and B. Altshuler, Physical Review B 76, 052203 (2007b).