Quantum dynamics in the interacting Fibonacci chain
Cecilia Chiaracane, Francesca Pietracaprina, Archak Purkayastha, John Goold
QQuantum dynamics in the interacting Fibonacci chain
Cecilia Chiaracane, ∗ Francesca Pietracaprina, † Archak Purkayastha, ‡ and John Goold § School of Physics, Trinity College Dublin, Dublin 2, Ireland (Dated: February 2, 2021)Quantum dynamics on quasiperiodic geometries has recently gathered significant attention in ultra-cold atomexperiments where non trivial localised phases have been observed. One such quasiperiodic model is the socalled Fibonacci model. In this tight-binding model, non-interacting particles are subject to on-site energiesgenerated by a Fibonacci sequence. This is known to induce critical states, with a continuously varying dy-namical exponent, leading to anomalous transport. In this work, we investigate whether anomalous diffusionpresent in the non-interacting system survives in the presence of interactions and establish connections to a pos-sible transition towards a localized phase. We investigate the dynamics of the interacting Fibonacci model bystudying real-time spread of density-density correlations at infinite temperature using the dynamical typicalityapproach. We also corroborate our findings by calculating the participation entropy in configuration space andinvestigating the expectation value of local observables in the diagonal ensemble.
I. INTRODUCTION
Anderson localisation is the phenomenon where electronsundergo quantum coherent scattering with random impuritiesand give rise to a metal-insulator transition. In three spa-tial dimensions, a critical energy depending on the disorderstrength, the so called mobility edge, separates localized fromextended states associated with diffusion [1]. In one dimen-sion, instead, all the eigenstates get spatially localized wheneven an infinitesimal amount of uncorrelated disorder is in-troduced, and the wires switch from ballistic conductors toinsulators [2]. Yet if the random disorder is replaced bya quasiperiodic potential, incommensurate with the underly-ing periodicity of the lattice, a wider variety of behavioursarises even in one dimension. The paradigmatic example isAubry-Andr´e-Harper (AAH) model, where the cosine poten-tial, modified by an irrational factor in its argument, inducesa transition from a completely delocalized phase to a com-pletely localized phase as the strength of the quasiperiodicpotential is increased [3, 4].At the critical point of the non-interacting AAH model boththe spectrum and the eigenfunctions show a fractal structure,and transport becomes sub-diffusive [5–8]. Recently, therehave been several theoretical works investigating the effect ofmany-body interactions on the AAH model [9–18]. More-over, the control over the Hamiltonian and the initial condi-tions in ultra-cold atoms set-ups has given the platform to re-alize quasiperiodic models and probe their non-trivial trans-port properties [19] from the perspective of dynamics. Inthese experiments, by tuning the relative depths of the opticallattices trapping the atoms, it is possible to investigate boththe non-interacting limit of the models as well as, the effectof many-body interactions on the dynamics. For example, thesingle particle localization in the AAH model is known to giverise to many-body localization (MBL) [20–23] in presence ofinteractions. ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] Another well known example of a quasiperiodic systemis the Fibonacci model, characterized by a potential gener-ated by the Fibonacci substitution rule. Though this model istopologically related to the AAH model [24, 25], its trans-port properties are known to be extremely different. Forexample there is no delocalization-localization transition inthe Fibonacci model. Instead, the Fibonacci potential in-duces critical behavior of all eigenstates at every potentialstrength [6, 18, 26–29]. The transport exponents show asmooth crossover from ballistic to sub-diffusive with increasein the potential strength [27, 30, 31]. A natural question, then,is what happens to the transport behavior of the Fibonaccimodel in presence of interactions?Different answers to this question have been proposed inthe literature, ranging from a transition towards MBL [32], tometal-insulator transitions at low energies [33, 34], to persis-tence of the anomalous diffusion [35]. Motivated by the lackof experimental results at this moment, in this paper, we fo-cus on characterizing the transport of the interacting Fibonaccimodel with a further approach, exploiting dynamical quantumtypicality [36, 37]. We study the real-time broadening of theexpectation values of local number operators, starting from anon-equilibrium typical state, that can act as representative ofthe equilibrium ensemble. Via dynamical quantum typicality,the quantities are directly related to spread of density-densitycorrelations and thereby to classification of transport via theGreen-Kubo formula in the isolated system [38–40]. Thisapproach allows us to access larger systems sizes and muchlonger time scales than otherwise possible. We find strongevidence of sub-diffusive transport at large enough poten-tial strengths, which precedes a crossover to a possible MBLphase. We further corroborate the crossover to this phase bymeans of a study of participation entropy and calculations inthe diagonal ensemble.This article is structured as follows. In Section II, wepresent the non-interacting Fibonacci model and display itsnon-trivial transport properties by reproducing the known re-sults of the spreading of an initially localized wave-packet. InSection III, we add many-body interactions to the Hamilto-nian, and review the previous works in the literature about theinteracting Fibonacci model. Dynamical quantum typicalityis introduced in Section IV, where we also show explicitly a r X i v : . [ c ond - m a t . d i s - nn ] J a n how to exploit it in order to compute the infinite temperaturedensity correlations we use to classify transport and presentour main results. In Section V and VI, we complete theobservations from the study of the dynamics with a furtherinvestigation respectively on the participation entropy of thesystem and the expectation values of both the local occupa-tion and imbalance in diagonal ensemble. We conclude andsummarize in Section VII. II. NON-INTERACTING FIBONACCI MODEL
The Fibonacci model is an one-dimensional system of non-interacting fermions, described by the Hamiltonian ˆ H NI = N − (cid:88) l =1 t h (ˆ a † l ˆ a l +1 + h . c) + N (cid:88) l =1 u l ˆ a † l ˆ a l , (1)where ˆ a l is the annihilation operator of a fermion on site l , t h is the tunnelling constant and u l the on-site energy of site l .The on-site potential is binary u l = ( u A , u B ) , and the chainof values on the N sites is obtained by repeatedly applying theFibonacci substitution rule, given by u A → u A u B (2) u B → u A . (3)Alternatively, the sequence of N on-site energies can bebuilt connecting together two smaller sequences. Startingfrom two initial chains C = [ u B ] , C = [ u A ] , one gets C = C C = [ u A , u B ] . Increasingly longer chains are gen-erated by concatenation of the segments from the two pre-vious generations C n = [ C n − , C n − ] . Consequently, thelength of every chain C n belongs to the Fibonacci sequence F n ∈ { , , , , , , . . . } . Unlike a periodic sequence, gen-erated by smaller parts of it and presenting the same rateof u A to u B even in the indefinitely extended limit, the Fi-bonacci chain exhibits a u A /u B ratio equal to F n /F n − atthe n th generation, which goes to /τ for n → ∞ , with τ = (1 + √ / the golden ratio [41]. Experiments or sim-ulations can be strongly limited in system size and usuallyinvolve chains of generic length N that do not belong to theFibonacci chain. To treat the model with the Fibonacci po-tential within such small system sizes, we adopt the averagingprocedure used in Refs. [28, 31]. We consider an “infinite ”sequence with N ∞ (cid:29) N and cut out finite samples of length N . It is possible to prove that there exist N + 1 nonequiva-lent samples, among which one (two) is reflection symmetricaround the center of the chain for N even ( N odd) and eachof the remaining configurations has a symmetric partner, withsame eigenvalues and eigenfunctions [41]. Therefore, afterdiscarding the reflection symmetric examples and the sym-metric partners of the samples already considered, N/ (or ( N − / ) distinct samples are available to average over.This averaging procedure also restores effective translationalinvariance in the thermodynamic limit. The quasiperiodicityof the potential gives rise to a multi-fractal spectrum at every u A and u B [26, 28]. Therefore, we assume without loss of ∆ x ( t ) t u = 0 . .. u = 1 . .. u = 1 . .. u = 2 . .. u = 5 . .. (a) . . . α u (b) Figure 1. (a) The mean squared displacement ∆ x of a state ini-tially localized at the middle of a non-interacting Fibonacci chain of N = 1001 sites is computed in time for different example values ofthe potential strength u . The fits shown by the dashed lines in thelog-log plot follows Eq. (5). We notice the curves saturating at low u as the size of the system is reached, where the fast dynamics makesthe finite size of the system visible at shorter times. (b) The extractedexponent α varies continuously with u , indicating anomalous diffu-sion. generality to control a single parameter u A = − u B = u inunits of t h . It has been shown that the multi-fractality of thespectrum induces anomalous behaviour in the transport prop-erties of the non-interacting system [27].For completeness, we start here by reproducing knownresult on the transport properties of the non-interacting Fi-bonacci model. We initialize the isolated system in openboundary condition (OBC) with a fermion on site N/ . Thusthe chosen initial state at t = 0 is | Ψ(0) (cid:105) = (cid:80) l Ψ l (0) ˆ a † l | (cid:105) ,with Ψ l (0) = δ lN/ the components over the site basis and | (cid:105) the vacuum state. We trace the time-evolution of the com-ponents Ψ l ( t ) , governed by the Schr¨odinger equation, and themean squared displacement of the wavefunction ∆ x ( t ) = (cid:88) l (cid:2) ( l − N/ | Ψ l ( t ) | (cid:3) . (4)The asymptotic time-dependence of this quantity can be writ-ten as ∆ x ( t ) ∼ t α , (5)where α = 0 implies localisation, α = 2 denotes ballis-tic transport, α = 1 implies diffusion, < α < im-plies super-diffusive transport, and < α < implies sub-diffusive transport. In Fig. 1a we show the ∆ x ( t ) as func-tion of time for a chain of N = 1001 , averaged over the ( N − / nonequivalent configurations, for different valuesof the potential strength u . Because of the finite system size,at long enough times, ∆ x ( t ) saturates. This saturation hap-pens at time scales where the initially localized wave-packethas spread over the entire system. For small u , this time scalegoes as ∼ N/ . As u increases, the transport slows down,and correspondingly the saturation happens at longer times.The dynamical exponent α corresponds to the thermodynamiclimit behavior, and needs to be obtained from the long time be-havior of the system before this saturation happens. In Fig. 1bwe show α , extracted from the fits of the ∆ x ( t ) curves, as afunction of u . As expected, α tends to and transport be-comes ballistic for u → . It then decreases continuouslytowards for increasing u , with two regimes: super-diffusive( α > for u (cid:46) . ) and sub-diffusive ( α < for u (cid:38) . ).The above characterization of transport in the non-interactingFibonacci model is well known. The main results of this workconcerns the the Fibonacci model in presence of interactions,which we describe below. III. THE INTERACTING FIBONACCI MODEL
The interacting Fibonacci model is realized by adding toEq. (1) a nearest neighbour density-density term ˆ H I = ˆ H NI + 2∆ N − (cid:88) l =1 ˆ n l +1 ˆ n l , (6)where ˆ n l = ˆ a † l ˆ a l is the fermionic number operator on site l ,and ∆ the strength of the many-body interaction.The interacting version of the Fibonacci model has re-cently started to receive attention with respect to how themany-body term affects the transport properties of the sys-tem. In the recent work by Varma and ˇZnidariˇc, the dynam-ics of polarized domain walls and the boundary-driven Lind-blad equation steady states reveal diffusion at small interac-tion strengths [31]. The spectral analysis in Ref. [32] pro-vides, instead, evidence for a localization transition at finitepotential strength, that would constitute a genuine many-bodyeffect, since the non-interacting model does not exhibit a lo-calized phase. Finally, a non equilibrium Green’s functionsapproach in Ref. [35] suggests in the Fermi-Hubbard realiza-tion of the model a slow sub-diffusive behaviour at high poten-tial strength, determined by the non-trivial spectral propertiesof the model.Here we aim to investigate the survival of anomalous diffu-sion in the interacting model. In presence of interactions, the classification of transport based on spread of a localized wave-packet is no longer possible. Instead, we classify transport byspread of an inhomogeneity on the infinite temperature stateof the system. Via dynamical typicality, as we show below,this is exactly analogous to spread of a localized wave-packet,and reduces to the same in absence of interactions. IV. DYNAMICAL QUANTUM TYPICALITY
In many-body problems, the numerical simulation of theSchr¨odinger equation is technically challenging due the ex-ponential growth of the Hilbert space dimension D with thenumber of degrees of freedom of the system. Popular tech-niques, such as the time-dependent density matrix renormal-ization group [42], can push the simulations to large systemsizes N ∼ , but are limited to short times due to the in-creasing of entanglement, and cannot generally reach the timescales required to study equilibrium properties. However, itis possible to exploit the concept of dynamical quantum typi-cality (DQT) to circumvent part of these difficulties. The ap-proximation tells that it is possible to infer the dynamics ofthe system from a single pure state | ψ (cid:105) drawn at random onan arbitrary basis {| φ k (cid:105)} Dk =1 , that can be considered as “typi-cal” representative of the statistical ensemble [36, 37], as weexplain below.We write explicitly the typical state as | ψ (cid:105) = ˆ R D (cid:88) k =1 c k | φ k (cid:105) , c k = a k + ib k , (7)with ˆ R an arbitrary linear operator, and a k and b k mutually in-dependent random variables from Gaussian distributions withzero mean and variance 1/2. It can be shown from the proper-ties of the coefficients that the averaged expectation value ofan arbitrary Hermitian operator ˆ O in the typical state is equiv-alent to the expectation value taken with respect to a densitymatrix ˆ ρ , as follows O = (cid:104) ψ | ˆ O | ψ (cid:105) = T r [ˆ ρ ˆ O ] , (8)where the overline indicates the average over the probabilitydistribution, and the density matrix is defined as ˆ ρ = ˆ R ˆ R † . (9)Since the density matrix is positive semi-definite, it can al-ways be written in the above form. Thus, any mixed state ˆ ρ can be represented in terms of an ensemble of typical purestates. When we further look at the sample to sample fluctu-ations in taking the average over the distribution, assuming ˆ O is Hermitian, we get σ O = (cid:0) (cid:104) ψ | ˆ O | ψ (cid:105) (cid:1) − ( O ) = T r [(ˆ ρ ˆ O ) ] . (10)The expression can be bounded from above by σ O ≤ T r [ˆ ρ ˆ O ] ≤ (cid:107) ˆ O (cid:107) T r [ˆ ρ ] , (11)with (cid:107) ˆ O (cid:107) = (cid:0) max λ n (cid:1) considering the eigenvalues { λ on } ofthe operator, and T r [ˆ ρ ] the purity of the state ˆ ρ . The aboveresult can also be generalized to cases where ˆ O is not Hermi-tian, by breaking ˆ O into Hermitian and anti-Hermitian parts.For highly mixed state in a high dimensional Hilbert space T r [ˆ ρ ] (cid:28) . In such cases, the sample to sample fluctuationsin doing the ensemble average also becomes small, so that, fora large enough system size, a small number of realizations isenough to calculate expectation values of operators. In partic-ular, in the infinite temperature limit, the state is completelymixed, and essentially one typical state realization can be usedas representative of the whole ensemble: ˆ ρ = 11 D , | ψ (cid:105) = 1 √ D D (cid:88) k =1 c k | φ k (cid:105) , (12)where D = 2 N for the Hamiltonian ˆ H I .The formulation described above does not depend on anyspecific property of the operator ˆ O . It can also be a combi-nation of operators in the Heisenberg picture so that both thedynamics, as well as two-time correlations can be obtained.Dynamical typicality can be used to connect two-time den-sity correlations in the infinite temperature state to dynamicsfollowing an initially localized quench over the infinite tem-perature state, as we explain in the following section. A. Density correlations at infinite temperature and spread of alocalized quench
Two-time density correlations are intimately connectedwith transport properties of an isolated system in the thermo-dynamic limit [43, 44]. The infinite temperature correlationfunction reads C lq ( t ) = (cid:104) ˆ n l ( t )ˆ n q (cid:105) − (cid:104) ˆ n l (cid:105) (cid:104) ˆ n q (cid:105) = T r [ˆ n l ( t )ˆ n q ]2 N − T r [ˆ n l ]2 N T r [ˆ n q ]2 N . (13)Let us choose q = N/ , and define C l ( t ) = C lN/ ( t ) . We canthen further simplify by using ˆ n N/ = ˆ n N/ and T r [ˆ n p ] =2 N − , as shown in the following C l ( t ) = T r [ˆ n l ( t )ˆ n N/ ]2 N − T r [ˆ n l ]2 N T r [ˆ n N/ ]2 N = T r [ˆ n N/ ˆ n l ( t )ˆ n N/ ]2 N − (cid:104) ψ N/ | ˆ n l ( t ) | ψ N/ (cid:105) − (14)In the last equality, we have exploited typicality with | ψ N/ (cid:105) = ˆ n N/ | ψ (cid:105) , | ψ (cid:105) = 12 N/ D (cid:88) k =1 c k | φ k (cid:105) (15)where, as seen previously in Eq. (12), | ψ (cid:105) is the typical stateassociated to the thermal state at infinite temperature. Since the dimension of the Hilbert space grows exponentially, atlarge enough sizes N (cid:38) sample to sample fluctuationsare negligible and it can be considered only one typical state.When we normalize the typical state to | ˜ ψ N/ (cid:105) ≈ √ | ψ N/ (cid:105) (shown in Appendix A), we finally obtain C l ( t ) ≈ C typl ( t ) = 12 (cid:16) n l ( t ) − (cid:17) , (16)with n l ( t ) = (cid:104) ˜ ψ N/ | ˆ n l ( t ) | ˜ ψ N/ (cid:105) . Thus, the density-densitycorrelation C l ( t ) is given by the dynamics of the expectationvalue of ˆ n l after a quench induced by the normalized projec-tion of a typical state onto the subspace where the site N/ isoccupied [43, 44]. The subtraction of 1/2 within the parenthe-ses amounts to subtracting the background initial occupationof sites away from the middle of the chain l (cid:54) = N/ . In gen-eral, as computed explicitly in Appendix A, we have C l ( t ) ≈ C typl ( t ) = 12 (cid:16) n l ( t ) − n l (0) (cid:17) , l (cid:54) = N/ (17)where, again, the expectation values are evaluated with re-spect to | ˜ ψ N/ (cid:105) .In order to classify transport, we define the spatial variance σ ( t ) = 4 N (cid:88) l =1 (cid:18) l − N (cid:19) C typl ( t ) . (18)As shown in Appendix. B, this quantity is related to the time-dependent diffusion coefficient at infinite temperature limit[38–40], dσ dt = 8 D ( t ) , (19) D ( t ) = lim β → lim N →∞ N (cid:90) t dt (cid:48) (cid:104) ˆ I ( t (cid:48) ) ˆ I (cid:105) β,µ . (20)For a diffusive system, the diffusion coefficient is constant, D ( t ) = D and hence σ ( t ) = 8 D t . If D ( t ) decreases withtime, it points to sub-diffusive transport. If D ( t ) is zero, thenthere is no transport. If D ( t ) ∼ t , then the transport is ballis-tic, while if D ( t ) diverges with a power less than , then thetransport is super-diffusive. Consequently, if, in general, σ ( t ) ∼ t α , (21)then, α = 1 corresponds to diffusive transport, α = 2 corre-sponds to ballistic transport, < α < corresponds to super-diffusive transport, < α < corresponds to sub-diffusivetransport, and α = 0 corresponds to lack of transport. Thisis exactly same as the transport classification in terms of ∆ x ( t ) for the non-interacting system given in Eq. (5). In-deed, as we show in Appendix C, for a non-interacting system, σ ( t ) = ∆ x ( t ) .Another way to characterize transport is via decay of den-sity autocorrelation with time. The infinite temperature den-sity autocorrelation at site N/ is given by C N/ ( t ) = T r [ˆ n N/ ( t )ˆ n N/ ]2 N − T r [ˆ n l ]2 N ≈ C typN/ ( t ) = 12 (cid:16) n N/ ( t ) − (cid:17) (22) σ ( t ) t u = 0 . u = 0 . u = 1 . u = 2 . u = 4 . u = 8 . (a) α u N = 20 N = 22 N = 24 (b) σ ( t ) t N = 20 N = 22 N = 24 f ( t ) ∝ t . ∝ t . ∝ t . (c) σ ( t ) t N = 20 N = 22 N = 24 (d) Figure 2. (a) Log-log plot of σ ( t ) vs. time t at different Fibonacci potential strengths u , for a chain of N = 24 spins. The data computedfrom the Krylov time-evolution are shown with continuous lines. We show the corresponding fits of the form of Eq. (21) in dotted lines; theshort-time fits are shown in pink, while the long-time ones are in black. We notice σ ( t ) growing faster than t at low u , and slowing down as u increases. (b) Exponent α extracted from the σ ( t ) as a function of the potential strength u . The crosses correspond to the fast dynamics(shown in pink in panel a), while the dots are relative to the long-time dynamics (shown in black in panel a). The errors of each data pointare smaller than the dot size. (c) σ ( t ) is shown at u = 1 . for three different system-sizes in linear scale, along with their correspondingpower-law fits. (d) The same quantity is displayed for u = 8 . at three different system sizes N = 20 , , for a longer time in linear scale. Thus, via typicality, it corresponds to how the occupation atthe middle site approaches its thermal value following thequench. We assume a general power-law decay of autocor-relation, C N/ ( t ) ≈ C typN/ ( t ) ∼ t − ν . (23)For a ballistic system, ν = 1 , for a diffusive system ν = 1 / ,while for a localized system ν = 0 . Correspondingly, <ν < / points to super-diffusive transport, and / < ν < points to sub-diffusive transport. For diffusive and ballisticsystems, this exponent ν is related to the exponent α as α =2 ν . But, for anomalous transport, these two exponents maynot be directly related. In the following, we present resultsfor transport classification of the interacting Fibonacci chainbased on calculation of both the exponents. B. Results
We numerically study the Hamiltonian, recast throughJordan-Wigner transformations into that of a spin 1/2 XXZmodel with external magnetic field, ˆ H I = N − (cid:88) l =1 [ t h (ˆ s + l ˆ s − l +1 +h . c)+ 2∆ ˆ s zl ˆ s zl +1 )]+ N (cid:88) l =1 u l ˆ s zl (24)with ˆ s + l , ˆ s − l , ˆ s zl respectively the raising, lowering and z spinoperators at site l , and u l the on-site potential following theFibonacci sequence. We restrict our calculations to sectors ofthe total Hilbert space with fixed magnetization, choosing, inparticular, the largest one with N/ spins up. Details on howto modify the expressions in the DQT approach are reportedin Appendix A. We generate a single typical state by tak-ing a normalized state-vector | ˜ ψ s (cid:105) with random coefficients,and apply the operator ˆ n N/ . We time-evolve this state using Krylov subspace method [45] and calculate the density pro-file at each time point, from which both σ ( t ) and C N/ ( t ) .The dynamical typicality approach, together with the Krylovsubspace methods, allows us to do a long time simulation ofa maximum system size of N = 24 . We fix the hopping term t h = 1 , and the interaction strength ∆ = 0 . , and investigatethe nature of transport as a function of the strength of the Fi-bonacci potential u . All the results shown are averaged overthe collection of the non-equivalent samples, as described inSec. II.In Fig. 2a, we show σ ( t ) as a function of t at significantvalues of u , for N = 24 . As with the non-interacting sys-tem, for small u , u (cid:46) ∆ , due to the finite size of the system,saturation occurs at time ∼ N/ . Power-law fitting of thedata before the saturation yields a super-diffusive exponent, > α > (see plot for u = 0 . in Fig. 2a). On increasing u ,the transport slows down, and therefore it takes much longertime to hit saturation. For u (cid:38) ∆ , we see a clear sub-diffusiveexponent, > α > (plots for u = 1 . , . in Fig. 2a)and saturation is not reached within our simulation time scalesand system sizes. The crossover from super-diffusive to sub-diffusive behavior seem to occur at u ≈ ∆ , where from ourresults, there does not seem to be a clear power-law behaviorbefore the saturation happens. It is possible, at best, to fit twodifferent power-laws at two different time regimes, between (cid:46) t (cid:46) with α ∼ . , and between (cid:46) t (cid:46) with α ∼ . . The exponents obtained from the power-law fitsare given in Fig. 2b, which shows the crossover from super-diffusive to sub-diffusive transport. At much higher values of u , u (cid:29) ∆ , σ ( t ) again to quickly saturates to a finite, lowvalue: this points to a lack of spreading of the initially lo-calized quench, α = 0 , and thereby pointing at a many-bodylocalized (MBL) regime ( u = 4 . , . in Fig. 2a); this is remi-niscent of the results of Ref. [32] and will be explored in moredetail in the following section. As often in literature, it is hardto pinpoint exactly at which values of u the crossover to MBL n N / ( t ) − / t u = 0 . u = 0 . u = 1 . u = 2 . u = 4 . u = 8 . Figure 3. Time evolution of C N/ ( t ) = n N/ ( t ) − evaluatedon the typical state projected over the subspace where the site at thecenter of the chain N/ is initially occupied by one particle. At u = 0 . , . , . , . , C N/ ( t ) decays as t − ν , with respectively ν =0 . , . , . , . . The fits are shown in dashed lines. At high u , C N/ ( t ) does not seem to show any decay up to longest simulationtime. regime happens from dynamical results.To highlight the differences between the sub-diffusive andthe MBL regime and to discuss finite-size effects, in Figs. 2cand 2d, we show plots of σ ( t ) for u = 1 . and u = 8 . re-spectively for different system sizes. In Fig. 2c, the long timepower-law growth of σ ( t ) with a sub-diffusive exponent isclear for all three system sizes N = 20 , , . With increasein system size, the time extent of the power-law growth in-creases, as expected, and the power-law exponent also con-verges (to α = 0 . ). However, the different system sizes no-ticeably do not overlap at any time scale. This is due to theeffect of the finite system size coupled with the quasiperiodicpotential: results for quasiperiodic systems, even for largesystem sizes are dependent on particular choice of system-sizes [5–7], particularly, for the Fibonacci potential, on howdifferent the system-sizes are from Fibonacci numbers. Thissystem-size dependence may be reduced by averaging oversamples, but the small number of available samples (equal tothe system size N ) limits the kind of averaging that is possibleto perform in our system sizes. We note that, while this be-havior holds for all values of u , this does not affect our abilityto obtain the power-law exponent and that, nevertheless, allthe results for the three different system sizes are of the sameorder of magnitude. In Fig. 2d, this same size-dependent ef-fect is shown in the localized regime for u = 8 . . Here, wehighlight the presence of oscillations, while at the same timeshowing no signs of a power-law growth trend.Next, we look at C N/ ( t ) and characterize transport interms of the exponent ν (Eq. (23)). The plots of C N/ ( t ) = n N/ ( t ) − are shown in Fig. 3. C N/ ( t ) shows a oscil-lations on top of a very clear power-law decay for u < . .For u (cid:46) ∆ , the power-law exponent is consistent with super-diffusive transport, > ν > . , for u (cid:38) ∆ , the power-lawexponent is consistent with sub-diffusive transport . > ν > . For u (cid:29) ∆ , corresponding to u = 4 . , . in Fig 3, we do P E S P / l n D u N = 14 N = 16 N = 18 N = 20 N = 22 Figure 4. Partecipation entropy S P / ln D associated to the centralregion of the spectrum for different Fibonacci potential strength u .The curves are displayed for multiple chain sizes. not see any power-law decay up to longest time scales that wesimulated, thereby suggesting localization. This is consistentwith our results from time scaling of σ ( t ) .In Ref. [31], it was shown that at small u at large enoughsystem sizes, the behavior becomes diffusive. Since our re-sults are limited to much smaller system sizes, we cannotcompletely rule out that possibility. Nevertheless, in Ref. [31],at larger values of u one parameter was reported showing sub-diffusive transport. This is completely consistent with the sub-diffusive behavior we observe for u (cid:38) ∆ .Our results strongly suggest that anomalous transport sur-vives in the Fibonacci model in presence of interactions.Moreover, at u (cid:29) ∆ , our investigations on dynamics sug-gests a crossover to MBL. This is very interesting because, inabsence of interactions, there is no localized phase. MBL inthe Fibonacci model was previously reported in Ref.[32], at adifferent strength of interaction ∆ . For the remainder of thestudy, we investigate the existence of MBL for our choice ofparameters from spectral properties of the Hamiltonian, andcalculations in the diagonal ensemble. V. EIGENSTATE PROPERTIES ACROSS DIFFERENTTRANSPORT REGIMES
The spectral properties of the Hamiltonian (6) have alreadyshown evidence of a many-body localization transition at fi-nite critical potential strength at ∆ = 1 . [32]. This phasewould be introduced uniquely by the interplay of quasidis-order and many-body interactions, since localization is notpresent in the noninteracting limit of the model. We performhere an analysis similar to [32], by computing the R´enyi- participation entropy of the Fibonacci chain with ∆ = 0 . through exact diagonalization (ED) [46]. These quantitieshave been used to characterize localization both in single-particle and in many-body interacting systems [47–49].Let | ψ E (cid:105) represent a many-body energy eigenstate. Thiscan be expanded in an arbitrary basis, which we choose to bethe configuration-space basis, as | ψ (cid:105) = (cid:80) Dk =1 d k | χ k (cid:105) . Theprobability p k = | d k | indicates the “participation” of the el-ement | χ k (cid:105) from the arbitrary basis {| χ k (cid:105)} k in the state | ψ E (cid:105) .The second R´enyi participation entropy (PE) is given by S P = − ln (cid:16) D (cid:88) k =1 p k (cid:17) . (25)If the eigenstate is completely delocalized, S P / log( D ) → .On the other hand, if S P / log( D ) → D , then the eigenstateis fractal with a fractal dimension of D . For a system show-ing MBL, the mid-spectrum energy eigenstates, the regionof the Hilbert space that is sampled by the isolated systemat infinite temperature, are expected to be fractal with a lowfractal dimension. For systems which are neither completelydelocalized nor in MBL, S P / log( D ) for the mid-spectrumeigenstates may not converge to a constant. The study of S P ,thereby allows to capture crossover to MBL.In Fig. 4, we plot S P / log( D ) as a function of the potentialstrength for different system sizes, N = 14 , , , , .All the points are obtained from an average over ∼ mid-specturm eigenstates, with the exception of the data for N = 22 that are averaged over ∼ eigenstates. Finally, thePE are averaged over the nonequivalent realizations of the Fi-bonacci potential. The mid-spectrum eigenvalues and eigen-states are obtained through the shift-invert algorithm [50].At very low u , the PE S P / log( D ) is close to , but stillshows dependence on the system size. At larger values of u , S P / log( D ) decays rapidly with u , and eventually showsa collapse for the different system sizes. Thus, two regimescan be identified, corresponding to the transport (either su-perdiffusive or subdiffusive) and no transport regimes foundin Sec. IV B, the latter reminiscent of the many-body local-ized phase identified in Ref. [32]. More definitive statementsabout the transition in the thermodynamic limit would requirea systematic study of the finite size scaling which is beyondthe purpose of the present work. Instead, in the following, weexplore yet another way of characterizing the MBL transitionfrom finite system sizes. VI. DIAGONAL ENSEMBLE
In Sec. IV B, we obtained finite-time results for the dynam-ics of the system at different potential strengths u . In thissection, we instead focus on the asymptotic results by usingfull ED and the diagonal ensemble, or infinite time averagedstate, to investigate the time infinite limit of n N/ ( t ) of theisolated system initialized in a non-equilibrium state | ψ (cid:105) , inorder to understand if the system reaches absence of transportat high values of u or rather exhibits a region of slow dynam-ics. Given the use of full ED, we are limited to smaller systemsizes than in Sec. IV B.The time infinite limit of an arbitrary observable ˆ O reads as O ( t → ∞ ) = lim T →∞ T (cid:90) T (cid:104) ψ ( t ) | ˆ O | ψ ( t ) (cid:105) dt, (26) . . . n N / ( t →∞ ) - /N u = 0 . u = 0 . u = 1 . u = 1 . u = 2 . u = 4 . u = 8 . u = 10 (a) . . . I ( t →∞ ) /N u = 0 . u = 0 . u = 1 . u = 1 . u = 2 . u = 4 . u = 8 . u = 10 (b) Figure 5. (a) Expectation value of the occupation at half chain ˆ n N/ in the diagonal ensemble associated to the initial typical state, whichgives the infinite time limit of the operator. The dotted lines indicateextrapolation for N → ∞ for the first 3 values of u , described by /N γ with γ = 2 . , . , . for increasing u . (b) Expectationvalue of the imbalance ˆ I in the diagonal ensemble for the initial N´eelstate, with the same colour code of (a). Again, the dotted lines rep-resent the fits /N γ we use to extrapolate the value of I ( t → ∞ ) inthe thermodynamic limit, with γ = 8 . , . , . for increasing u . and, if the spectrum is not degenerate, it can be easily re-written as O ( t → ∞ ) = O diag = (cid:88) k (cid:104) φ k | ˆ O | φ k (cid:105) | (cid:104) φ k | ψ (cid:105) | , (27)where O diag indicates the expectation value of the operator ˆ O in the diagonal ensemble relative to the initial state | ψ (cid:105) .However, the computation of O diag requires full ED of theHamiltonian, so our results are limited to the system sizes N = 10 , , , , , up to a maximum of obtained onlyat u ≥ . .We focus on the occupation number at half chain ˆ n N/ ,considering the diagonal ensemble for the typical state. Theoccupation of the initial state is 1 by construction, and it willeventually reach the equilibrium value of ∼ . in case ofthermalization. The results are shown in Fig. 5a for differentpotential strengths as a function of the inverse of the systemsize. At low u (cid:46) , the value of the observable at infinite timedecreases with N and we are able to extrapolate the infinite-size limit result through a fit of the form ∼ /N γ ; the fits areshown with dotted lines in Fig. 5a and extrapolate to . , indi-cating that the system thermalizes in the thermodynamic limit.At larger u , we do not assume a form for the finite size scalingof n N/ and thus we do not extrapolate the infinite-size limit.We also consider the imbalance [19, 20], a density correla-tion function defined as I ( t ) = 4 N N (cid:88) j =1 (cid:104) ψ (0) | (ˆ n j ( T ) − / n (0) − / | ψ (0) (cid:105) , (28)which in the case of initial N´eel state | ψ N (cid:105) can be written asthe following operator ˆ I = [ˆ n e − ˆ n o ] /N, (29)where ˆ n e/o = (cid:80) l e/o ˆ n l is the number of particles at theeven ( e ) or odd ( o ) sites. We compute the initial imbalance I ( t = 0) = (cid:104) ψ N | ˆ I | ψ N (cid:105) , and derive its infinite time limit fromits expectation values in the diagonal ensemble associated tothe N´eel state. The initial value at t = 0 is and will eventu-ally reach the equilibrium value of if there is thermalization.We show the infinite-time limit of the imbalance in Fig. 5b asa function of /N . The results are similar to those from ˆ n N/ and the random typical state. At low u , it is possible to ex-trapolate the imbalance in the thermodynamic limit, giving .However, at larger potential strength, namely for u > , thedata shows a lack of decay with N , up to the system sizes wehave access to, and supports our results obtained in Sec. IV Bpointing to absence of transport in the system and localization. VII. CONCLUSIONS
In this work, we have studied the dynamics of density-density correlations at infinite temperature of the Fibonaccimodel in presence of nearest neighbour many-body interac-tions via direct numerical simulation using the DQT approach.The DQT approach, coupled with Krylov subspace methods[38–40], have allowed us to obtain the density correlations forlarger system sizes and much longer time scales than other-wise possible. This allowed us to extract the dynamical expo-nents corresponding to the transport properties of the model.We have further correlated our results with calculations of theparticipation entropy of the mid-spectrum states, and with ex-act diagonalization calculations in the diagonal ensemble cor-responding to non-equilibrium initial states, at smaller sys-tem sizes. We have focused on a fixed interaction strength, ∆ = 0 . t h , and have characterized the transport as functionof the strength of the Fibonacci potential u . The following pic-ture emerges from our investigation. For u (cid:46) ∆ , the transportis relatively fast. We find some evidence of super-diffusionin this regime, although the fast transport and the finite sys-tem sizes do not allow us to extract a long time transport ex- ponent. On increasing u , the transport slows down, allow-ing us to extract long time exponents. For u (cid:38) ∆ , we finda strong evidence of sub-diffusive transport. The crossoverfrom super-diffusive to sub-diffusive behavior seems to occurat u ∼ ∆ , where we are unable to extract a single dynami-cal exponent. On further increasing u , i.e, for u (cid:29) ∆ , we findstrong evidence that the system crosses over to an MBL phase.The MBL phase is then further corroborated with studies ofthe participation entropy and the diagonal ensemble, both ofwhich complement the results from study of the dynamics.The above picture that emerges from our study contributestowards filling a gap in our present understanding of in-teracting quasiperiodic systems. Most studies of interact-ing quasiperiodic systems have focused on the AAH poten-tial, both in theory and in experiment [9–14, 16, 17, 20–23].Though related with the AAH model, the non-interacting Fi-bonacci model is known to have very different transport prop-erties, which continuously crosses over from ballistic to sub-diffusive as a function of the strength of the potential [27, 30].There has been only few works exploring the Fibonacci modelin presence of interactions [31, 32, 35]. In Ref. [32], the spec-tral properties of the Fibonacci model were studied as a func-tion of the potential strength at a fixed interaction strength of ∆ = t h . At this interaction strength, in absence of the Fi-bonacci potential, i.e, in the ordered XXZ-chain, the transportis known to be super-diffusive [40, 51, 52]. A transition toMBL phase was predicted. This is very interesting, because,in absence of interactions, the Fibonacci potential shows nolocalization. The question, then, is whether this MBL canbe seen at lower interaction strengths. The infinite tempera-ture transport properties at small interaction were investigatedin Ref. [31]. This study gave strong evidence that presenceof a small interaction makes the transport diffusive at all po-tential strengths. This again, is very non-trivial, because, inabsence of interactions there is smooth crossover from ballis-tic to sub-diffusive. This crossover between diffusive to sub-diffusive transport also occur in some disordered interactingsystems [53]. The question, then, becomes whether trans-port can become anomalous again at intermediate interactionstrengths. One parameter point was shown in favour of thisin Ref. [31]. Our choice of interaction strength, ∆ = 0 . t h ,is intermediate. At this choice of interaction strength, in ab-sence of the Fibonacci potential, i.e, in ordered XXZ-chain,the transport is known to be ballistic [40, 51, 52]. Our findingsprovide strong evidence that crossover to MBL with increasein the strength of Fibonacci potential can happen even at thisinteraction strength, and is preceded by a regime a anomaloussub-diffusive transport. This answers both the above ques-tions. This is also very different from a third work, Ref. [35],which studied transport properties of a different model, thespinful Fibonacci model with Fermi-Hubbard interaction, andshowed that a localized phase cannot occur in that system, andthere will always be slow sub-diffusive transport at large in-teractions, and large potential strengths.More definitive results on the MBL phase would requirestudy of larger systems up to longer times, which is beyondcurrent state-of-the-art numerical techniques. Given the pe-culiar spectral properties of this class of models, a study ofthe energy dependence of the transport properties is a verypromising direction for a subsequent investigation, for exam-ple using open systems techniques [54, 55]. Moreover, allpresent studies of transport properties are limited to infinitetemperature and zero temperature [33, 34]. The finite temper-ature transport properties, as well as the study of thermoelec-tric behavior [56], are other interesting but challenging direc-tions for future work. ACKNOWLEDGEMENTS
We thank Mark. T. Mitchison, Marlon Brenes, N. Lo Gulloand N. Laflorencie for insightful discussions. The spin con-figurational basis, the Hamiltonian and the operators are gen-erated by the open source python package QuSpin [57]. Weacknowledge the provision of computational facilities by theDJEI/DES/SFI/HEA Irish Centre for High-End Computing(ICHEC). This project received funding from the EuropeanResearch Council (ERC) under the European Union’s Hori-zon 2020 research and innovation program (Grant AgreementNo. 758403). J.G. is supported by a SFI-Royal Society Uni-versity Research Fellowship. FP has received funding fromthe European Union’s Horizon 2020 research and innovationprogramme under the Marie Sklodowska-Curie grant agree-ment No 838773. AP is supported by the European Union’sHorizon 2020 research and innovation programme under theMarie Sklodowska-Curie grant agreement No. 890884.
APPENDIXAppendix A: Normalization and restriction to half-filled sector
In Eq. (16), we use the normalized typical state | ˜ ψ N/ (cid:105) = 1 √ C | ψ N/ (cid:105) , (A1)where C is the normalization constant, given by C = (cid:104) ψ N/ | ψ N/ (cid:105) = (cid:104) ψ | ˆ n N/ | ψ (cid:105)≈ (cid:104) ψ | ˆ n N/ | ψ (cid:105) = T r [ˆ n N/ ]2 N = 12 . (A2)Moreover, the subtraction of the factor 1/2 from the termwithin the parentheses amounts to subtracting the backgroundinitial occupation of sites away from the site N/ , where thetypical state is initially localized. In order to verify it, we no-tice that for r (cid:54) = N/ , n r (0) = (cid:104) ˜ ψ N/ | ˆ n r | ˜ ψ N/ (cid:105) ≈ (cid:104) ψ N/ | ˆ n r | ψ N/ (cid:105) = 2 (cid:104) ψ | ˆ n N/ ˆ n r ˆ n N/ | ψ (cid:105) = 2 (cid:104) ψ | ˆ n r ˆ n N/ | ψ (cid:105)≈ (cid:104) ψ | ˆ n r ˆ n N/ | ψ (cid:105) = 2 T r [ˆ n r ˆ n N/ ]2 N = 12 , (A3)where in the last equality we use that T r [ˆ n r ˆ n N/ ] = 2 N − .The above results and those described in Sec. IV A do notmake use of the fact that the system is number conserving. For a large enough number conserving system, the biggest contri-bution to Eq. (17) comes from the half-filled sector. It is plau-sible that in such case, one can completely restrict the calcula-tion to the half-filled sector, starting from a typical state in thesector, thus saving computational resources and pushing theforward the system size. In complete analogy to Sec. IV A,we define a typical state in the half-filled subsector | ψ s (cid:105) = 1 D s D s (cid:88) k =1 c k | φ sk (cid:105) , D s = N !( N/ N/ , (A4)where {| φ sk (cid:105)} D s k =1 is an orthonormal basis in the half-filled sec-tor. The new normalization constant C s in | ˜ ψ sN/ (cid:105) = 1 √ C s | ψ sN/ (cid:105) , | ψ sN/ (cid:105) = ˆ n N/ | ψ s (cid:105) , (A5)is given by C s = (cid:104) ψ sN/ | ψ sN/ (cid:105) = (cid:104) ψ s | ˆ n N/ | ψ s (cid:105) ≈ (cid:104) ψ s | ˆ n N/ | ψ s (cid:105) = T r [ˆ n N/ ] D s = ( N − N − ! N − ! N ! N ! N ! = 12 . (A6)As before, we have | ˜ ψ sN/ (cid:105) ≈ √ | ψ sN/ (cid:105) . (A7)However, the background occupation of sites q (cid:54) = N/ is nowless than 1/2, as it is possible to notice by reproducing theresult of Eq. (A3) in the half-filled sector: n sr (0) = (cid:104) ˜ ψ sN/ | ˆ n r | ˜ ψ sN/ (cid:105) ≈ (cid:104) ψ sN/ | ˆ n r | ψ sN/ (cid:105) = 2 (cid:104) ψ s | ˆ n r ˆ n N/ | ψ s (cid:105) ≈ (cid:104) ψ s | ˆ n r ˆ n N/ | ψ s (cid:105) = 2 T r [ˆ n r ˆ n N/ ] D s = ( N − N − ! N ! N ! N ! N != 12 (cid:16) − N − (cid:17) . (A8)Finally, in analogy with Eq. (16) we are able to define C sl ( t ) = 12 ( n sl ( t ) − n sl (0)) , l (cid:54) = N/ ≈ (cid:104) n sl ( t ) − (cid:16) − N − (cid:17)(cid:105) , (A9)where n sl ( t ) is the expectation value of the operator ˆ n l at time t , starting from the initial state given by | ˜ ψ sN/ (cid:105) . For largeenough system, we expect C sl ( t ) ≈ C l ( t ) . We have con-firmed this by directly comparing simulations performed forshort time interval on a chain of size N = 20 both in the totalHilbert space and in the largest sector at half-filling confirmsour conjecture. Appendix B: Relation with Green-Kubo conductivity
The Green-Kubo formula for particle conductivity at finitetemperature can be written as σ GK = β lim t →∞ lim N →∞ N Re (cid:18)(cid:90) t dt (cid:48) (cid:104) ˆ I ( t (cid:48) ) ˆ I (cid:105) β,µ (cid:19) , (B1)0where ˆ I is the particle current operator, and (cid:104) ... (cid:105) β,µ denotesaverage taken over the thermal state of the system with tem-perature β and chemical potential µ . In above the order oflimits cannot be interchanged. For one-dimensional systemswith open boundary condition, the particle current operator isgiven by ˆ I = d ˆ xdt , (B2)where ˆ x = N (cid:88) p =1 p ˆ n p (B3)is the position operator. This definition gives, (cid:104) ˆ I ( t ) ˆ I ( t ) (cid:105) = ddt ddt (cid:32) N (cid:88) p,q =1 pq (cid:104) ˆ n p ( t )ˆ n q ( t ) (cid:105) β,µ (cid:33) . (B4)Using time translational invariance of the thermal state andchanging variable to τ = t − t , we have (cid:104) ˆ I ( τ ) ˆ I (cid:105) = − d dτ (cid:32) N (cid:88) p,q =1 pqC pq ( β, µ, t ) (cid:33) , (B5)where C pq ( β, µ, τ ) = (cid:104) ˆ n p ( τ )ˆ n q (cid:105) β,µ . (B6)Now we use the relation, pq = p + q − ( p − q ) , along withthe assumption that the Hamiltonian is number conserving, sothat ddτ (cid:16)(cid:80) Np =1 ˆ n p ( τ ) (cid:17) = 0 , to obtain (cid:104) ˆ I ( τ ) ˆ I (cid:105) β,µ = 12 d dτ (cid:32) N (cid:88) p,q =1 ( p − q ) C pq ( β, µ, τ ) (cid:33) . (cid:90) τ dt (cid:48) (cid:104) ˆ I ( t (cid:48) ) ˆ I (cid:105) β,µ = 12 ddτ (cid:32) N (cid:88) p,q =1 ( p − q ) (cid:104) C pq ( β, µ, τ ) (cid:33) . (B7)Using above equation in Eq.(B1), we have, σ GK = β t →∞ ddt m nn ( t ) m nn ( t ) = lim N →∞ N Re (cid:32) N (cid:88) p,q =1 ( p − q ) C pq ( β, µ, t ) (cid:33) . (B8)Further simplification of m nn ( t ) is possible if the systemhas translational invariance in the thermodynamic limit. Inthat case, C pq ( t ) becomes almost independent of q for largeenough system sizes. So, we can fix q = N/ , to obtain, m nn ( t ) = lim N →∞ Re (cid:32) N (cid:88) p =1 (cid:18) p − N (cid:19) C p N ( β, µ, t ) (cid:33) , (B9) Writing the above expression in the high-temperature limit, β → , using the dynamical typicality and the fact that at β → , m nn ( t ) is real, we have, in presence of translationalinvariance in the thermodynamic limit, σ ( t ) = 4 lim β → m nn ( t ) . (B10)The scaling of the quantity σ ( t ) with time gives the nature ofhigh temperature transport. Let us define the time-dependentdiffusion coefficient at high temperature as D ( t ) = lim β → lim N →∞ N (cid:90) t dt (cid:48) (cid:104) ˆ I ( t (cid:48) ) ˆ I (cid:105) β,µ . (B11)Then, from above, dσ dt = 8 D ( t ) . (B12)This derivation relies on translational invariance of the systemin the thermodynamic limit. Though the Fibonacci model isnot translationally invariant in the thermodynamic limit, thetranslational invariance is effectively restored on averagingover the various realizations. Appendix C: Relation with the spread of wavepacket in thenon-interacting system
We go back again to the case of the non-interacting system,described by the Hamiltonian in Eq. (1), that we diagonalizeas Φ T H NI Φ = D , D = diag { (cid:15) ν } , (C1)with the single-particle eigenvectors given by the columns of Φ , and the eigenvalues by { (cid:15) ν } . In the diagonalized basis, theHamiltonian reads ˆ H NI = N (cid:88) ν =1 (cid:15) ν ˆ c † ν ˆ c ν , (C2)where ˆ c ν = (cid:80) Np =1 Φ pν ˆ a p are the fermionic annihilation op-erators in the eigenbasis. The two time density correlation atfinite temperature can be simplified as follows C pq ( β, t ) = (cid:104) ˆ n p ( t )ˆ n q (cid:105) − (cid:104) ˆ n p (cid:105) (cid:104) ˆ n q (cid:105) = N (cid:88) ν,α =1 Φ νp Φ νq Φ αp Φ αq e it/ (cid:126) ( (cid:15) ν − (cid:15) α ) [1 − n F ( (cid:15) α )] n F ( (cid:15) ν ) , (C3)where (cid:104)(cid:105) indicates the ensemble average, after applying theWick’s theorem (cid:104) ˆ a † p ( t p )ˆ a q ( t q )ˆ a † m ( t m )ˆ a n ( t n ) (cid:105) = (cid:104) ˆ a † p ( t p )ˆ a q ( t q ) (cid:105) (cid:104) ˆ a † m ( t m )ˆ a n ( t n ) (cid:105) + (cid:104) ˆ a † p ( t p )ˆ a n ( t n ) (cid:105) (cid:104) ˆ a q ( t q )ˆ a † m ( t m ) (cid:105) , (C4)1and the following relations (cid:104) ˆ a † p ( t p )ˆ a q ( t q ) (cid:105) = N (cid:88) ν =1 Φ νp Φ νq e i(cid:15) ν / (cid:126) ( t p − t q ) n F ( (cid:15) ν ) (cid:104) ˆ a p ( t p )ˆ a † q ( t q ) (cid:105) = N (cid:88) ν =1 Φ νp Φ νq e i(cid:15) ν / (cid:126) ( t p − t q ) (1 − n F ( (cid:15) ν )) (C5)with n F ( E ) = { β ( E − µ )] } − the Fermi-Dirac dis-tribution. Now, we take the infinite temperature limit β → of Eq. (C3), and shift the labels to consider the correlation be-tween the middle of the chain and the other sites q = N/ and p = l + N/ as in the previous section, C l ( t ) = 14 N (cid:88) ν,α =1 Φ νl Φ νN/ Φ αl Φ αN/ e it ( (cid:15) ν − (cid:15) α ) / (cid:126) (cid:16) N (cid:88) α =1 Φ αl Φ αN/ e − it(cid:15) α / (cid:126) (cid:17)(cid:16) N (cid:88) ν =1 Φ νl Φ νN/ e it(cid:15) ν / (cid:126) (cid:17) = 14 | Ψ l ( t ) | , (C6) where Ψ l ( t ) = N (cid:88) ν =1 Φ νl Φ νN/ e it(cid:15) ν / (cid:126) . (C7)The dynamics of each Ψ l ( t ) corresponds to evolution accord-ing to i d Ψ l ( t ) dt = N (cid:88) r =1 H lr Ψ l ( t ) , (C8)starting from the initial condition Ψ l ( t ) = δ lN/ . Thus, fromEq. (C6) we derive that the two time density correlation at infi-nite temperature we use to classify transport in the interactingsystem is directly proportional to | Ψ l ( t ) | in single-particlesystems: | Ψ l ( t ) | = 4 C l ( t ) ≈ C typl ( t ) . (C9)Physically, | Ψ l ( t ) | gives the probability of finding a particleat site l = p − N/ , after initializing the system with a singleparticle located at site N/ . From above, and Eqs.(4) and (18),we see that, for a non-interacting system, σ ( t ) = ∆ x ( t ) .But unlike ∆ x ( t ) , σ ( t ) is well-defined also in presence ofinteractions. [1] G. Semeghini, M. Landini, P. Castilho, S. Roy, G. Spagnolli,A. Trenkwalder, M. Fattori, M. Inguscio, and G. Modugno,Nat. Phys. , 554 EP (2015).[2] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V.Ramakrishnan, Phys. Rev. Lett. , 673 (1979).[3] S. Aubry and G. Andr´e, Ann. Israel Phys. Soc , 18 (1980).[4] H. Hiramoto and M. Kohmoto, Int. J. Mod. Phys. B , 281(1992).[5] A. Purkayastha, S. Sanyal, A. Dhar, and M. Kulkarni, Phys.Rev. B , 174206 (2018).[6] V. K. Varma, C. de Mulatier, and M. ˇZnidariˇc, Phys. Rev. E ,032130 (2017).[7] J. Sutradhar, S. Mukerjee, R. Pandit, and S. Banerjee, Phys.Rev. B , 224204 (2019).[8] A. Purkayastha, Journal of Statistical Mechanics: Theory andExperiment , 043101 (2019).[9] Y. Yoo, J. Lee, and B. Swingle, Phys. Rev. B , 195142(2020).[10] T. Cookmeyer, J. Motruk, and J. E. Moore, Phys. Rev. B ,174203 (2020).[11] M. ˇZnidariˇc and M. Ljubotina, Proceedings of the NationalAcademy of Sciences , 4595 (2018).[12] Y. B. Lev, D. M. Kennes, C. Kl¨ockner, D. R. Reichman, andC. Karrasch, EPL (Europhysics Letters) , 37003 (2017).[13] P. Naldesi, E. Ercolessi, and T. Roscilde, SciPost Phys. , 010(2016).[14] V. Mastropietro, Phys. Rev. Lett. , 180401 (2015).[15] J. Settino, N. Lo Gullo, A. Sindona, J. Goold, and F. Plastina,Phys. Rev. A , 033605 (2017).[16] S. Iyer, V. Oganesyan, G. Refael, and D. A. Huse, Phys. Rev.B , 134202 (2013). [17] M. Tezuka and A. M. Garc´ıa-Garc´ıa, Phys. Rev. A , 031602(2012).[18] J. X. Zhong and R. Mosseri, Journal of Physics: CondensedMatter , 8383 (1995).[19] H. P. L¨uschen, S. Scherg, T. Kohlert, M. Schreiber, P. Bordia,X. Li, S. Das Sarma, and I. Bloch, Phys. Rev. Lett. , 160404(2018).[20] T. Kohlert, S. Scherg, X. Li, H. P. L¨uschen, S. Das Sarma,I. Bloch, and M. Aidelsburger, Phys. Rev. Lett. , 170403(2019).[21] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen, M. H.Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Sci-ence , 842 (2015).[22] H. P. L¨uschen, P. Bordia, S. S. Hodgman, M. Schreiber,S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, andU. Schneider, Phys. Rev. X , 011034 (2017).[23] H. P. L¨uschen, P. Bordia, S. Scherg, F. Alet, E. Altman,U. Schneider, and I. Bloch, Phys. Rev. Lett. , 260401(2017).[24] Y. E. Kraus and O. Zilberberg, Phys. Rev. Lett. , 116404(2012).[25] V. Goblot, A. ˇStrkalj, N. Pernet, J. L. Lado, C. Dorow,A. Lemaˆıtre, L. Le Gratiet, A. Harouri, I. Sagnes, S. Ravets,A. Amo, J. Bloch, and O. Zilberberg, Nature Physics , 832(2020).[26] M. Kohmoto, B. Sutherland, and C. Tang, Phys. Rev. B ,1020 (1987).[27] H. Hiramoto and S. Abe, Journal of the Physical Society ofJapan , 230 (1988).[28] N. Mac´e, A. Jagannathan, and F. Pi´echon, Phys. Rev. B ,205153 (2016). [29] A. Jagannathan, (2020), arXiv:2012.14744 [cond-mat.stat-mech].[30] J. Zhong, R. B. Diener, D. A. Steck, W. H. Oskay, M. G. Raizen,E. W. Plummer, Z. Zhang, and Q. Niu, Phys. Rev. Lett. ,2485 (2001).[31] V. K. Varma and M. ˇZnidariˇc, Phys. Rev. B , 085105 (2019).[32] N. Mac´e, N. Laflorencie, and F. Alet, SciPost Phys. , 50(2019).[33] J. Vidal, D. Mouhanna, and T. Giamarchi, Phys. Rev. Lett. ,3908 (1999).[34] J. Vidal, D. Mouhanna, and T. Giamarchi, Phys. Rev. B ,014201 (2001).[35] J. Settino, N. W. Talarico, F. Cosco, F. Plastina, S. Maniscalco,and N. Lo Gullo, Phys. Rev. B , 144303 (2020).[36] C. Bartsch and J. Gemmer, Phys. Rev. Lett. , 110403 (2009).[37] P. Reimann, Phys. Rev. E , 062129 (2018).[38] R. Steinigeweg, H. Wichterich, and J. Gemmer, EPL (Euro-physics Letters) , 10004 (2009).[39] T. Heitmann, J. Richter, D. Schubert, and R. Steinigeweg,Zeitschrift f¨ur Naturforschung A , 421 (01 May. 2020).[40] R. Steinigeweg, F. Jin, D. Schmidtke, H. De Raedt,K. Michielsen, and J. Gemmer, Phys. Rev. B , 035155(2017).[41] G. R. Goodson, Chaotic dynamics: fractals, tilings and substi-tutions (Cambridge University Press, 2017).[42] S. Paeckel, T. K¨ohler, A. Swoboda, S. R. Manmana,U. Schollw¨ock, and C. Hubig, Annals of Physics , 167998(2019).[43] J. Richter, F. Jin, L. Knipschild, J. Herbrych, H. De Raedt,K. Michielsen, J. Gemmer, and R. Steinigeweg, Phys. Rev. B , 144422 (2019).[44] J. Richter, D. Schubert, and R. Steinigeweg, Phys. Rev. Re-search , 013130 (2020).[45] A. Nauts and R. E. Wyatt, Phys. Rev. Lett. , 2238 (1983).[46] N. Mac´e, F. Alet, and N. Laflorencie, Phys. Rev. Lett. ,180601 (2019).[47] B. Kramer and A. MacKinnon, Reports on Progress in Physics , 1469 (1993).[48] D. J. Luitz, N. Laflorencie, and F. Alet, Journal of StatisticalMechanics: Theory and Experiment , P08007 (2014).[49] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B ,081103 (2015).[50] F. Pietracaprina, N. Mac´e, D. J. Luitz, and F. Alet, SciPostPhys. , 45 (2018).[51] R. Steinigeweg, EPL (Europhysics Letters) , 67001 (2012).[52] M. ˇZnidariˇc, Phys. Rev. Lett. , 220601 (2011).[53] V. K. Varma, A. Lerose, F. Pietracaprina, J. Goold, andA. Scardicchio, Journal of Statistical Mechanics: Theory andExperiment , 053101 (2017).[54] J. J. Mendoza-Arenas, M. ˇZnidariˇc, V. K. Varma, J. Goold, S. R.Clark, and A. Scardicchio, Phys. Rev. B , 094435 (2019).[55] M. Brenes, J. J. Mendoza-Arenas, A. Purkayastha, M. T.Mitchison, S. R. Clark, and J. Goold, Phys. Rev. X , 031040(2020).[56] C. Chiaracane, M. T. Mitchison, A. Purkayastha, G. Haack, andJ. Goold, Phys. Rev. Research , 013093 (2020).[57] P. Weinberg and M. Bukov, SciPost Phys.2