Many-body localization in Bose-Hubbard model: evidence for the mobility edge
MMany-body localization in Bose-Hubbard model: evidence for the mobility edge
Ruixiao Yao and Jakub Zakrzewski
2, 3 School of Physics, Peking University, Beijing 100871, China Institute of Theoretical Physics, Jagiellonian University in Krakow, (cid:32)Lojasiewicza 11, 30-348 Krak´ow, Poland Mark Kac Complex Systems Research Center, Jagiellonian University in Krakow, Krak´ow, Poland. ∗ (Dated: May 15, 2020)Motivated by recent experiments on interacting bosons in quasi-one-dimensional optical lattice[Nature , 385 (2019)] we analyse theoretically properties of the system in the crossover betweendelocalized and localized regimes. Comparison of time dynamics for uniform and density wave likeinitial states enables demonstration of the existence of the mobility edge. To this end we definea new observable, the mean speed of transport at long times. It gives us an efficient estimate ofthe critical disorder for the crossover. We also show that the mean velocity growth of occupationfluctuations close to the edges of the system carries the similar information. Using the quantumquench procedure we show that it is possible to probe the mobility edge for different energies. Many-body localization (MBL) despite numerous ef-forts of the last 15 years (for some reviews see [1–4]) isstill a phenomenon not fully understood. Recently evenits very existence in the thermodynamic limit has beenquestioned [5] which created a vivid debate [6–8]. Simu-lations of large systems dynamics are also not fully con-clusive [9, 10]. It seems that, as suggested by [8], presentday computer resources prevent us from drawing a defi-nite conclusion on this point. The issue may be addressedin experiments via quantum simulator approach [11] al-though the required precision may be also prohibitive.The experimental studies of MBL are much less nu-merous. Early work [12] provided indications of MBL ina large fermionic system in an optical lattice. Subsequentstudies showing the lack of thermalization and a mem-ory of the initial configuration in time evolution were re-ported for interacting fermions in quasiperiodic disorderin optical quasi-one-dimensional (1D) lattice [13]. Ex-periments also considered rather large systems with ei-ther fermions [14] or bosons [15] (the latter in two di-mensions). Only recently quite small systems attractedexperimental attention first for interacting photons [16]where even the attempt at level spacing statistics mea-surement was made as well as for bosonic atoms wherelogarithmic spreading of entanglement was observed [17]as well as long-range correlations analysed [18].On the theory side bosonic system were not the primechoice for the analysis for a simple reason. While forspin-1/2 (or spinless fermion) the local Hilbert space di-mension is two, twice that for spinful fermions, for bosonsthe strict limit is set by a total number of particles in thesystem. This severely limits the number of sites that canbe included in any simulation performed within exactdiagonalization-type studies. In effect only few papersaddressed MBL with bosons. By far the most notableis a courageous attempt to treat bosonic system in twodimensions by approximate tensor network approach [19](still restricting the model to maximally double occupa-tions of sites). Earlier treatment [20] considered bothsmall and large systems in 1D revealing the existence of
FIG. 1. Mean gap ratio r as a function of the energy for U = 1 (left) and U = 2 .
87 right, for a system of 8 bosons on 8lattice sites. Observe that for the higher U case, high energystates are localized at lower disorder amplitudes revealing theinverted mobility edge. Red and yellow lines represent theenergy of the uniform initial state (with unit occupation ofall sites) and the density wave with even states being dou-bly occupied while odd states being empty. Observe that forlarger U the transition to localized situation as revealed by¯ r statistics depends strongly on the initial state energy, thedensity wave initial state should localize more easily. the so called reverse mobility edge in the energy spectra.With assumed 3/2 filling of the system it has been shownthat higher lying in energy states are localized for loweramplitude of the disorder for sufficiently strong interac-tions. Many body localization with superconducting cir-cuits were discussed in [21] while the role of doublons forthermalization in [22]. Very recently the same system wasdiscussed also at half-filling [23] confirming the existenceof the mobility edge. Let us mention for completenessalso works on MBL of bosons with random interactions[24, 25] instead of a random on-site potential.In these works the tight-binding model describing thephysics is given by 1D Bose-Hubbard Hamiltonian: H = − J L − (cid:88) i (cid:16) ˆ b † i ˆ b i +1 + h.c. (cid:17) + U L (cid:88) i ˆ n i (ˆ n i − L (cid:88) i µ i ˆ n i , (1)where ˆ b i (ˆ b † i ) are bosonic annihilation (creation) opera-tors on site i with [ˆ b i , ˆ b † j ] = δ ij and ˆ n i = ˆ b † i ˆ b i . We assumeuniform interactions U across the lattice while the on-site a r X i v : . [ c ond - m a t . d i s - nn ] M a y energies (chemical potential) µ i is site dependent. While[20, 23] considered a random uniform on-site disorder,following the experiments [13, 17, 18] from now on weassume the quasi-periodic disorder µ i = W cos(2 πβi + φ )where, for a given realization, φ is fixed and disorderaveraging corresponds to an average over uniformly dis-tributed φ ∈ [0 , π ). It is known that the localizationproperties of the system strongly depend on the value of β [26, 27] , from now on we take β = 1 / .
618 and a unitfilling following [17, 18] and work with open boundaryconditions. We consider sufficiently deep optical lattices,so (1) may be used, for shallow quasiperiodic potentialcase see [28].With the model defined we study first its spectral prop-erties. The localization may be probed looking at themean gap ratio, r [29] defined as an average of r n -the minimum of the ratio of consequtive level spacings s n = E n +1 − E n : r n = min { s n +1 s n , s n s n +1 } . (2) r is a simple, dimensionless probe of level statistics [29]with r ≈ .
38 for Poisson statistics (PS) (correspondingto localized, integrable case) and r ≈ .
53 for Gaussianorthogonal ensemble (GOE) of random matrices well de-scribing statistically the ergodic case [30, 31]. Diagonaliz-ing the model for 8 bosons on 8 sites and calculating ¯ r asa function of the disorder amplitude W and the relativeenergy (cid:15) n = E n − E min E max − E min we obtain the plot representedin Fig. 1. While for U = 1 the energy dependence of ¯ r seems only weakly dependent on energy, for U = 2 .
87 aclear inverted mobility edge emerges: at the same disor-der values while states of higher energies seem localized,the low lying states have ¯ r corresponding to extendedstates. This is in agreement with earlier studies [20, 23]at different densities.A more qualitative analysis might be carried out withfinite size scaling (FSS) [32]. Recent works have indi-cated [33–35] , however, that MBL transition may be ofKosterlitz-Thouless type rendering symmetric FSS ques-tionable. In view of that we show r as a function ofdisorder for different system sizes with crossing of thecurves being an indicator of the crossover disorder value W c - compare Fig. 2. We may observe that the transitionfor U = 2 .
87 clearly depends on the initial state energysuggesting the existence of the mobility edge. The cross-ing between L = 8 and L = 9 data occurs for W c = 6 . W c = 4 .
0) for (cid:15) corresponding to uniform - UN (densitywave - DW) state. On the other hand for U = 1 W c ≈ . FIG. 2. Left: Mean gap ratio r reveals that transition toMBL, as indicated by crossing of curves for different systemsizes (indicated in the figure) occurs for energies correspond-ing to a uniform state [(a) panel], (cid:15) ≈ .
15] for much strongerdisorder than for (cid:15) = 0 . x (see text for the definition) as obtained from time dynamicsstarting from those two states (right column) which revealssubdiffussive growth on the localized side (straight line be-haviour in log-log plots with the slope well below 0.5). Thecurves correspond to W = 1 , , , , , , W = 1in both cases saturates at intermediate times due to a finitesystem size - the time dependence is shown for L = 8. used for different DW type states [20] we shall considerthe transport distance (for other studies of the relatedquantity for spin systems see [37, 38]) defined as [18]:∆ x = 2 | (cid:80) d d × (cid:104) G (2) c ( i, i + d ) (cid:105) i | - i.e. the modulus ofthe disorder averaged (as denoted by the overbar) andsite-averaged (as denoted by a subscript i ) second ordercorrelation function of density G (2) c ( i, i + d ) = (cid:104) ˆ n i ˆ n i + d (cid:105) −(cid:104) ˆ n i (cid:105)(cid:104) ˆ n i + d (cid:105) . We concentrate mostly on U = 2 .
87 casewhich reveals mobility edge in above spectral analysisand confront them with data obtained for U = 1 whenno mobility edge evidence is expected for UN or DWstates.The time dependence of the transport distance ob-tained for L = 8 system is shown in Fig. 2 in rightpanels for different values of the disorder. Data are ob-tained using exact diagonalization approach. On the de-localized side one observes a fast almost ballistic growth(for W = 1) followed by a saturation when the finitesize of the system dominates the dynamics. Upon ap-proaching the transition the power-law growth becomesapparent for intermediate times, the corresponding power β = d ln ∆ x/d ln t decreasing smoothly within the inter-val of [0,0.35] indicating a subdiffusive dynamics for bothinitial states considered (similar behavior is observed for U = 1 [36]). For sufficiently large W the motion freezeswith very small β . The data are averaged over 200 dis- FIG. 3. The mean transport velocity V x (averaged in t ∈ [100 , U = 1(left) and U = 2 .
87 (right) for UN state (top) and DW (bot-tom). The derivative seems to indicate the transition to MBLquite clearly, again manifesting different critical disorder value(when it vanishes on the localized side) for UN and for DWstate at U = 2 .
87. Statistical errors of V x are smaller thanthe symbols. order realizations, the residual fluctuations apparent inthe figures form an estimate of the statistical error.Let us inspect in more detail the system properties insubdiffusive region. Apart from L = 8 case, we con-sider also L = 10 ,
12 to observe the effects due to thesystem size. For L = 10 ,
12 time evolution is carriedour using the standard Chebyshev technique [39, 40].Due to a moderate number of 200 disorder realizationsto minimize local fluctuations one may average ∆ x ( t )obtained over some time interval, for sufficiently largetimes. On the localized side one expects such an aver-age, ∆ x L to be small, signalling the transition. Whilesuch an approach signals the mobility edge for U = 2 . x ( t ) over some large interval.Such velocity of transport V x can be readily obtained as V x = [∆ x ( t ) − ∆ x (( t )] / ( t − t ). It is shown for t = 100and t = 250 in Fig. 3. The larger the system size themore pronounced are the maxima of the derivative. Onthe delocalized side V x is small as ∆ x saturates due tofinite sizes considered. On the localized side V x practi-cally vanishes as transport for long times is prohibited.Thus the pronounced maximum of the derivative in thecritical transition regime is a robust and aposteriori ex-pected phenomenon. The critical disorder is read out asthe point where V x practically vanishes (being say 1% ofthe maximal value). While it is important to analyse V x at large times, after the initial growth, the results arerobust against varying the choice of the time interval aswell as the method of estimating the mean velocity [36].For U = 1 (left column) the uniform (a) and the DW (c)initial state lead to similar V x dependence as a function of W . V x approaches 0 and becomes almost independent ofthe system size for W ≈ ± . FIG. 4. Same as Fig. 3 but for the edge fluctuation veloc-ity V F . For U = 2 .
87 (right) a clear difference between UNand DW initial states also appears, the transition to MBLoccurring for larger W for UN state, in agreement with thegap ratio analysis. Statistical errors in V F determination arecomparable to symbols’ sizes. critical disorder values obtained from gap ratio analysis.For U = 2 .
87 case, the position of the peaks of V x for UN (b) and DW (d) states in Fig. 3 strongly differ.Using the same criterion of vanishing V x in the localizedregime, we get the critical disorder W c ≈ W c ≈ . x re-quires measurements of all second order correlations. Asit turns out this may not be necessary. Consider on-sitefluctuations F i = (cid:104) n i (cid:105)−(cid:104) n i (cid:105) . While [41] considered fluc-tuations in the middle of the chain, it is advantageous toconcentrate on F = ( F + F L ) / F , at long times shows indications of thecrossover to localized phase [36], we present in Fig. 4 itsderivative V F = dF/dx averaged over “long times” in amanner completely analogous to V x . Remarkably, it alsoreveals the mobility edge when varying W in a similarmanner to V x plotted just above.Our analysis up till now have been restricted to initialFock-like states and the corresponding energies. It wassuggested, however, that different energy regimes may beaddressed by a “quantum quench spectroscopy” [41] toreveal the mobility edge energy dependence. The pro-posed scheme assumes a preparation of the ground stateof the system at some value of a parameter, then the fastquench of that parameter to the final investigated valuetransfers the system into an excited wavepacket with ex-cess energy dependent on the change of this parameter.By changing the initial parameter value one may, hope-fully, scan the final energy. This method was tested in[41] on spinless fermions case with quenching the disordervalue and observing the time dependence of the entan- FIG. 5. (a): Quantum quench from initial disorder ampli-tude W i to final values W (as indicated in the figures) allowsto prepare initial states of different energy (cid:15) . (b): The meantransport velocity at long times V s obtained for such pre-pared states reveals a W -dependent crossover to MBL as afunction of (cid:15) visualizing the mobility edge seen in Fig. 1(b).The crossover epsilon values agree with gap ratio statistics forcorresponding W values shown in panel (c). Edge-fluctuationvelocity V F [panel (d)] is also sensitive to the crossover in amanner similar to V s . glement entropy as well as on site density fluctuations.Let us apply the same method to the current prob-lem. We assume that the ground state is prepared atdifferent disorder amplitudes W i and the rapid quenchbrings the disorder amplitude to a desired final value W .Changing W i we may change the energy ( (cid:15) ) of the pre-pared wavepacket as shown in Fig. 5(a). Note that forreaching significant excitation (cid:15) we start with “negative”amplitude W i . In this way we can reach different final W values. The obtained energies are characterized by, un-fortunately, quite significant error. For a given preparedinitial W i different realizations of the disorder (differentphase φ ) lead to different excitation energies. As notedin [41] the resulting error diminishes with the increasingsystem size, here we consider a small system L = 8 mostly(see [36] for larger systems) as then gap ratio statisticsallows us to verify the results obtained. The error barsshown in Fig. 5 are due to this effect and limit the accu-racy of the final energy obtained by quenching.The initial state is then evolved in time determining∆ x ( t ) and its mean derivative V x = d ∆ x/dt exactly asbefore. At different final disorder values a transition fromdelocalized to localized regimes is observed when chang-ing “smoothly” (cid:15) – compare Fig. 5(b). Again the tailof V x , when its value becomes close to zero is an indi-cator of MBL regime. For comparison, Fig. 5(c) showsthe corresponding gap ratio statistics at these disordervalues. The agreement is spectacular showing that thequantum quench spectroscopy combined with the trans-port distance measurements allows to continuously moni-tor the mobility edge in our system. The drawback of the method is the inherent limitation of the energy epsilon resolution discussed above.To complete the picture, Fig. 5(d) presents the edge-fluctuation derivative V F = dF/dx averaged over “longtimes”. As for UN and DW states in Fig. 4 it also revealsthe mobility edge when varying final disorder strength W . This is quite promising for possible applications ofquench spectroscopy to systems of larger size where mea-surements leading to transport distance (and V x ) maybecome costly while measurements of edge fluctuationsrequire site resolution at the edges only.To conclude, by considering transport properties in thetransition between extended and localized states in Bose-Hubbard Hamiltonian describing bosons in optical latticewith diagonal quasiperiodic disorder we have shown thatobservables directly accessible to the experiment [18] re-veal the existence of the mobility edge in the system.The mobility edge has been convincingly shown to exitsfor the finite size Heisenberg chain [32] studying differ-ent spectra measures (most notably the gap ratio statis-tics) as well as using quantum quench spectroscopy andtime dependence of the entanglement and number fluc-tuations [41]. It has also been postulated to exist forBose-Hubbard system with a larger density by consider-ing different density wave states [20]. Here we show thatwithin the realm of state of the art, current experiment[18] the mobility edge may be verified, for the first time,experimentally, via readily accessible observables. Weconsidered different initial states, either prepared as Fockstates with different density patterns or states obtainedvia quantum quench of the disorder amplitude. The lat-ter approach allows us (modulo inherent uncertainties)to follow the mobility edge in energy. Both the globaltrends and values of the critical disorder obtained by usquantitatively agree across different observables taken foranalysis and agree with the spectral gap ratio analysis.When this work was completed we have learnt that en-ergy resolved MBL was also considered in spin quantumsimulator [43].We are grateful to Sooshin Kim and Markus Greinerfor discussing the details of simulations accompanyingthe experiment [18] and to Piotr Sierant for remarks onthe manuscript. JZ thanks Dominique Delande and TitasChanda for discussions. The numerical computationshave been possible thanks to High-Performance Comput-ing Platform of Peking University. Support of PL-Gridinfrastructure is also acknowledged. This research hasbeen supported by National Science Centre (Poland) un-der project 2016/21/B/ST2/01086 (J.Z.). ∗ [email protected][1] D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys.Rev. B , 174202 (2014). [2] R. Nandkishore and D. A. Huse, Ann. Rev. Cond. Mat.Phys. , 15 (2015).[3] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018), quantum simulation / Simulation quan-tique.[4] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev.Mod. Phys. , 021001 (2019).[5] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, arXive-prints , arXiv:1905.06345 (2019), arXiv:1905.06345[cond-mat.str-el].[6] D. A. Abanin, J. H. Bardarson, G. D. Tomasi,S. Gopalakrishnan, V. Khemani, S. A. Parameswaran,F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur,(2019), arXiv:1911.04501 [cond-mat.str-el].[7] P. Sierant, D. Delande, and J. Zakrzewski, (2019),arXiv:1911.06221 [cond-mat.dis-nn].[8] R. K. Panda, A. Scardicchio, M. Schulz, S. R. Tay-lor, and M. ˇZnidariˇc, (2019), arXiv:1911.07882 [cond-mat.dis-nn].[9] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D.Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi,Phys. Rev. B , 174202 (2018).[10] T. Chanda, P. Sierant, and J. Zakrzewski, Phys. Rev. B , 035148 (2020).[11] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Dem-ler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriks-son, K.-M. C. Fu, M. Greiner, K. R. A. Hazzard, R. G.Hulet, A. J. Kollar, B. L. Lev, M. D. Lukin, R. Ma,X. Mi, S. Misra, C. Monroe, K. Murch, Z. Nazario, K.-K.Ni, A. C. Potter, P. Roushan, M. Saffman, M. Schleier-Smith, I. Siddiqi, R. Simmonds, M. Singh, I. B. Spielman,K. Temme, D. S. Weiss, J. Vuckovic, V. Vuletic, J. Ye,and M. Zwierlein, arXiv e-prints , arXiv:1912.06938(2019), arXiv:1912.06938 [quant-ph].[12] S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco,Phys. Rev. Lett. , 083002 (2015).[13] M. Schreiber, S. S. Hodgman, P. Bordia, H. P. L¨uschen,M. H. Fischer, R. Vosk, E. Altman, U. Schneider, andI. Bloch, Science , 842 (2015).[14] H. P. L¨uschen, P. Bordia, S. Scherg, F. Alet, E. Altman,U. Schneider, and I. Bloch, Phys. Rev. Lett. , 260401(2017).[15] J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch,and C. Gross, Science , 1547 (2016).[16] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Basti-das, A. Megrant, R. Barends, Y. Chen, Z. Chen,B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen,M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus,M. Neeley, C. Quintana, D. Sank, A. Vainsencher,J. Wenner, T. White, H. Neven, D. G. Ange-lakis, and J. Martinis, Science , 1175 (2017),https://science.sciencemag.org/content/358/6367/1175.full.pdf.[17] A. Lukin, M. Rispoli, R. Schittko, M. E. Tai,A. M. Kaufman, S. Choi, V. Khemani, J. L´eonard,and M. Greiner, Science , 256 (2019),https://science.sciencemag.org/content/364/6437/256.full.pdf.[18] M. Rispoli, A. Lukin, R. Schittko, S. Kim, M. E. Tai,J. L´eonard, and M. Greiner, Nature , 385 (2019).[19] T. B. Wahl, A. Pal, and S. H. Simon, Nature Physics , 164 (2019).[20] P. Sierant and J. Zakrzewski, New Journal of Physics ,043032 (2018).[21] T. Orell, A. A. Michailidis, M. Serbyn, and M. Silveri, Phys. Rev. B , 134504 (2019).[22] U. Krause, T. Pellegrin, P. W. Brouwer, D. A. Abanin,and M. Filippone, arXiv e-prints , arXiv:1911.11711(2019), arXiv:1911.11711 [cond-mat.dis-nn].[23] M. Hopjan and F. Heidrich-Meisner, arXiv e-prints, arXiv:1912.09443 (2019), arXiv:1912.09443 [cond-mat.str-el].[24] P. Sierant, D. Delande, and J. Zakrzewski, Phys. Rev.A , 021601 (2017).[25] P. Sierant, D. Delande, and J. Zakrzewski, Acta Phys.Polon. A , 1707 (2017).[26] V. Guarrera, L. Fallani, J. E. Lye, C. Fort, and M. In-guscio, New J. Phys. , 107 (2007).[27] E. V. H. Doggen and A. D. Mirlin, Phys. Rev. B ,104203 (2019).[28] H. Yao, H. Khoudli, L. Bresque, and L. Sanchez-Palencia, Phys. Rev. Lett. , 070405 (2019).[29] V. Oganesyan and D. A. Huse, Phys. Rev. B , 155111(2007).[30] M. L. Mehta, Random Matrices (Elsevier, Amsterdam,1990).[31] F. Haake, Quantum Signatures of Chaos (Springer,Berlin, 2010).[32] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B , 081103 (2015).[33] N. Mac´e, N. Laflorencie, and F. Alet, SciPost Phys. ,50 (2019).[34] N. Laflorencie, G. Lemari´e, and N. Mac´e, arXive-prints , arXiv:2004.02861 (2020), arXiv:2004.02861[cond-mat.dis-nn].[35] J. ˇSuntajs, J. Bonˇca, T. Prosen, and L. Vidmar, arXive-prints , arXiv:2004.01719 (2020), arXiv:2004.01719[cond-mat.dis-nn].[36] See Supplemental Material at [URL will be inserted bypublisher] for additional details.[37] S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Phys.Rev. Lett. , 196801 (2017).[38] F. Weiner, F. Evers, and S. Bera, Phys. Rev. B ,104204 (2019).[39] H. Tal-Ezer and R. Kosloff, The Jour-nal of Chemical Physics , 3967 (1984),https://doi.org/10.1063/1.448136.[40] H. Fehske and R. Schneider,Computational many-particle physics (Springer, Ger-many, 2008).[41] P. Naldesi, E. Ercolessi, and T. Roscilde, SciPost Phys. , 010 (2016).[42] V. Khemani, D. N. Sheng, and D. A. Huse, Phys. Rev.Lett. , 075702 (2017).[43] Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang,W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, R. Mondaini,H. Fan, and H. Wang, arXiv e-prints , arXiv:1912.02818(2019), arXiv:1912.02818 [quant-ph]. SUPPLEMENTARY MATERIALComments on gap ratio analysis and finite sizescaling
In the letter we did not show the results of finite sizescaling analysis. While several authors used a simple,single parameter scaling of the form W → ( W − W c ) L /ν to extract the critical disorder value W c , this procedurehas been put recently under critique. Already [32] no-ticed that the obtained values of W c and ν dangerouslydepend on system sizes taken for the analysis. More-over, the exponent ν appears to be close to unity whichviolates the so called Harris bound (for details see [32]and references therein). Having at our disposal systemswith L = 7 , , W c for U = 2.87 for (cid:15) values correspondingto UN and DW states is significant - giving an argu-ment towards the existence of the mobility edge. On theother hand variations of ν are large (20%) while there isno reason to believe that the character of the transitionchanges. Moreover, ν values, as for spin systems, violatethe Harris bound. FIG. S.1. 1 (left) and U = 2.87 (right) in the scaled en-ergy regions corresponding to two interesting states. Top rowcorresponds to (cid:15) = 0 . / .
15 (left/right) energy region appro-priate for uniform initial state for L = 8 (compare to Fig. 1in main text). Bottom row shows scaling for (cid:15) = 0 . / . Recent works on spin systems tend to describe MBLtransition as being of Kosterlitz-Thouless (KT) type[34, 35] and perform FSS with appropriate KT corre-lation function. We refrain from doing so as for bosonsmuch smaller system sizes are available and such a pro-cedure would also be doubtful. For that reason we usea simple crossing of r curves for different system sizes asa reasonable estimate of the characteristic disorder whenMBL sets in. One should also keep in mind that forsuch small sizes (as also used in the experiment [18]) thecrossover region must take a finite range of disorder am-plitudes, W , so determining a precise number for W c in the thermodynamic limit is not needed really. Details on accuracy issues and numericalsimulations
Both spectral data (gap ratio) and time dynamics areobtained with “exact” methods such as an exact diag-onalization and/or Chebyshev propagation scheme (forsizes L = 10 , φ from a randomuniform distribution from [0 , π ] interval. Those errorsmay be minimized by simply increasing the number ofdisorder realizations. Secondly, we define “average” or“mean” quantities at long times such as the mean trans-port distance, the mean derivative of the transport dis-tance etc. which are the quantities averaged over sometime interval. These definitions inherently induce addi-tional errors that we describe below. FIG. S.2. Left: The transport distance ∆ x for U = 1 casefor a uniform density (a) and the density wave (c) states for W = 3 , , , , , W . The dashed lines show alinear fit (in the log-log plot) indicating a power law time-dependence for intermediate times. Panel (c) enlarges (in lin-ear scale) the region of panel (a) used to determine “the meantransport distance ∆ x L (see text). Different curves are shiftedvertically for a better comparison. The fluctuations observedare due to a finite 200 disorder realizations (discrete pointscorrespond to the finite step of the output from the numericalcode. Panel (d) visualizes errors in the mean velocity of thetransport V x - obtained as a mean slope in [100,250] interval.Dashed line corresponds to a linear fit, dashed-dotted lineconnects simply the t = 100 and t = 250 points for W = 3.As discussed further in the text both can be used to estimate V x . Let us recall from the main text the definition of thetransport distance ∆ x = 2 | (cid:80) d d × (cid:104) G (2) c ( i, i + d ) (cid:105) i | [18]. FIG. S.3. The transport distance ∆ x L obtained around t =250 as a function of the disorder for all cases studied (the orderof panels is as in the previous figure). Data are for differentsystem sizes. Observe that for U = 1 (left column) the datafor uniform and DW initial states are similar. For U = 2 . W for the uniform state in agreement withgap ratio analysis. Here the second order correlation function of density G (2) c ( i, i + d ) = (cid:104) ˆ n i ˆ n i + d (cid:105) − (cid:104) ˆ n i (cid:105)(cid:104) ˆ n i + d (cid:105) is averaged overdisorder realizations and different sites i . As shown inFig. 2 of the letter, compare also Fig. S.2 the transportdistance first rapidly grows (on the time scale of few tun-nelings) redistributing particles. On a longer time scaletwo main features are observed. For small disorder (de-localized regime) ∆ x saturates due to a finite small sizeof the system. With increasing disorder this saturationshifts to longer times, see W = 3 curves which bends indi-cating a saturation at times t > x with time follows to a good ac-curacy a power law behaviour (as shown by dashed linesin left panels of Fig. S.2) with the power β ∈ (0 , . W values - showing the same effectof disorder. On a localized side both the growth and thevalues reached are quite small allowing for identificationof the localized regime.One could attempt to use these log-log fits to estimatethe transition to MBL (e.g. when β decays to zero). Sucha procedure has many drawbacks, as regions of approx-imate power law growth change with W , the smaller β the errors become more important. Instead, also in viewof short time-scale fluctuations due to a finite number ofdisorder realizations, it may be useful to consider a meantransport distance, ∆ x L at some experimentally reach-able time of few hundreds of tunneling times. We findthe mean value of the transport distance in the interval t ∈ [220 , x practically vanishes (reaches, say, 1/1000of the maximal value) - then the error of ∆ x L is verysmall. For small disorder values, e.g. W = 3 the errorcomes not from the disorder induced fluctuations (whichare the main reason for averaging over a finite time in-terval) but from the remaining growth. Still this errorfor W = 3 which may be estimated from Fig. S.2(b) tobe about 0.05 is less than one percent of the mean valuein this interval ∆ x L ≈ . W may be compared and fluctuations visu-alized.The main observable we use to identify the critical dis-order strength for a given energy is, however, not themean distance but rather the mean velocity of trans-port V x for intermediate and large times. For two times t and t > t it can be defined as V x = [∆ x ( t ) − ∆ x (( t )] / ( t − t ) and it is nothing else as d∆ x ( t ) / d t averaged over [ t , t ] interval. Alternatively one may fita linear slope to ∆ x ( t ) in the same interval. Both pro-cedures are illustrated in Fig. S.2(d) for W = 3 curveand lead, practically, for a sufficiently large interval, toalmost identical slopes. The error of of the slope fittingprocedure in [100 , FIG. S.4. The mean derivative of the transport distance V x (averaged in t ∈ [120 , The choice of the interval [ t , t ] is not important, aslong as we consider the interval for t sufficiently large,say t >
80. The mean velocity obtained for [120 , x ( t ) with equal weights we may calcu-late the weighted linear fit in which squared deviationsare weighted by a gaussian centered at the center of theinterval considered with a standard deviation σ = 25.Transport velocities obtained by such a procedure for thecenter at t = 175 in the same interval yield essentiallythe same data as in Fig. S.4 i.e. within the size of thesymbols. This indicates that long time mean transportvelocity is a robust measure. Quantum quench –additional results
Here we supplement the results presented in Fig. 5 ofthe main text for the quantum quench scenario.
FIG. S.5. Transport distance ∆ x (left: panels (a-b))and the edge-fluctuation F (right (panels (c-d)) as a func-tion of time for W = 1 (top) and W = 3 bottom. Dif-ferent curves correspond to different initial starting disor-der amplitude W i . The curves are naturally ordered by W i = − , − , − , − , − , −