Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Stability of mobility edges in disordered interacting systems

Stability of mobility edges in disordered interacting systems 1 2 1 Pietro Brighi , Dmitry A. Abanin , and Maksym Serbyn IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria and Department of Theoretical Physics, University of Geneva, 24 quai Ernest-Ansermet, 1211 Geneva, Switzerland (Dated: August 28, 2020) Many-body localization provides a mechanism to avoid thermalization in isolated interacting quantum systems. The breakdown of thermalization may be complete, when all eigenstates in the many-body spectrum become localized, or partial, when the so-called many-body mobility edge separates localized and delocalized parts of the spectrum. Previously, De Roeck et al. [Phys. Rev. B 93, 014203 (2016)] suggested a possible instability of the many-body mobility edge in energy den- sity. The local ergodic regions | so called \bubbles" | resonantly spread throughout the system, leading to delocalization. In order to study such instability mechanism, in this work we design a model featuring many-body mobility edge in particle density : the states at small particle density are localized, while increasing the density of particles leads to delocalization. Using numerical sim- ulations with matrix product states we demonstrate the stability of MBL with respect to small bubbles in large dilute systems for experimentally relevant timescales. In addition, we demonstrate that processes where the bubble spreads are favored over processes that lead to resonant tunneling, suggesting a possible mechanism behind the observed stability of many-body mobility edge. We conclude by proposing experiments to probe particle density mobility edge in Bose-Hubbard model. Introduction.|Many-body localization (MBL) pro- gested in Ref. [18], which equally applies to MBME in vides a mechanism to avoid thermalization in isolated any extensive conserved quantity. First, using numer- quantum interacting systems [1, 2]. Despite intensive ical simulation with matrix product states (MPS), we theoretical [3, 4] and experimental [5{11] studies, only demonstrate that uniform dilute states remain localized fully-MBL phase in one spatial dimension is relatively even at system sizes of L = 40 sites up to 250 tunneling well understood. The fate of MBL in higher dimen- times (i.e. more than two orders of magnitude larger than sions [13{17] and the possibility of the coexistence of the inverse local hopping). Next, we use a region with localized and delocalized eigenstates in the same many- large particle density to reproduce the bubble described body spectrum [18] remain debated. in [18] and track its in uence on the dilute remainder of the system in a quantum quench. We do not nd any Similarly to the case of Anderson localization [19], the evidence of resonant tunneling of the bubble, at least on MBL and delocalized eigenstates cannot coexist at the experimentally relevant timescales. same energy suggesting the existence of many-body mo- In summary, the study of the particle density MBME bility edge (MBME) | a certain energy in the spectrum facilitates the state preparation and analysis and allows separating localized and delocalized eigenstates [2]. In us to access the dynamics of much larger systems using contrast to the non-interacting case, the energy of MBME time evolution with MPS. We report the stability of the scales extensively with system size. In the absence of a particle density mobility edge on long timescales and sug- coupling to a bath, this leads to an exactly vanishing gest that similar physics may be experimentally probed conductivity (in contrast to an exponentially small but using Bose-Hubbard model. nite value in Anderson insulator) until a certain critical Correlated hopping model.|We consider hard-core temperature [2]. bosons on an open chain of size L, with the following Recently De Roeck et al. [18] suggested a possible Hamiltonian, mechanism that may destroy MBME in large systems: a nite region with local energy density above the mobility L1 L X X edge | a \bubble" | may resonantly spread through- H = t (c c + h:c:) +  n ^ 1 i i i i+1 out the system thereby destroying localization every- i=1 i=1 where. However, recent experiments [11] and MPS sim- L1 ulations [12] gave evidence of MBME, at least on inter- + t (c n ^ c + h:c:): (1) 2 i i+1 i1 mediate timescales. In addition, a number of numerical i=2 studies observed a mobility edge [20{24] using exact di- agonalization (ED). Unfortunately, the ED is limited to The rst line corresponds to the non-interacting Ander- relatively small system sizes; experiments with MBME son's model [25], where random on-site potential has a in energy density are also challenging since they require uniform distribution,  2 [W; W ]. The facilitated hop- energy resolution. ping in the second line enables motion of a pair of par- In order to overcome the above challenges, we pro- ticles with amplitude t , $, making the model pose to study MBME in particle density. This allows interacting. The Hamiltonian (1) has two channels for us to directly probe the mechanism of instability sug- dynamics: the single particle hopping prevails in dilute arXiv:2005.02999v2 [cond-mat.dis-nn] 26 Aug 2020 2 states, while the pair hopping is dominant at larger den- sities. We note that a similar model was discussed in Ref. [18] in two dimensions, although only with two particles. The enhancement of localization length in the case of two interacting particles also received signi cant atten- tion [26, 27]. In a di erent direction, the fate of the single particle mobility edge in the presence of interactions was studied [10, 28]. In contrast, we study model (1) that does not have a single particle mobility edge and con- sider the nite particle density regime. We x the value of the hopping parameters t = 0:5 Figure 1. Scaling of level spacing ratio demonstrates that and t = 2 so that the localization length of a single at density  = 1=5 (solid lines, L = 15; 20; 25 with 3; 4; 5 particle  . 1 and at the same time a single pair has SP particles) the system enters MBL phase for W  6:3. In a localization length  & 2:5 for 2:5 . W . 6 [29]. contrast, at half- lling  = 1=2 (dashed curves, L = 10; : : : 18) For such a choice, our model does not su er from nite the critical disorder strengths is much larger and in the entire size e ects [30] and we establish MBME using eigenstates range of disorder r approaches thermal value with increasing av probes. system size. Data is generated from ED/SI simulations with Eigenstate probes of localization.|We use exact di- at least 10 disorder realizations using approximately 2% of agonalization and shift-invert (SI) numerical techniques eigenstates in the center of the spectrum. to provide evidence for MBME in Hamiltonian (1). We analyze the average ratio of level spacings,  = E E , in the middle of the spectrum, r = the dense and dilute cases. While in the dilute case the i+1 i av hmin( ;  )= max( ;  )i. This is a commonly used system retains memory of the initial state, see Fig. 2(a), i i+1 i i+1 probe of the MBL transition [21, 31] that attains the at  = 1=2 quantum dynamics leads to a progressively value r ' 0:39 for the Poisson level statistics, charac- more uniform density pro le with increasing system size, teristic of the MBL phase, and r ' 0:53 for the case Fig. 2(b). In order to quantify the di erence in the form GOE of random Gaussian orthogonal ensemble (GOE), typical of the density pro le at late times, in Fig. 2(c) we plot for chaotic Hamiltonians with time-reversal symmetry. the average deviation of the density from the equilibrium Figure 1 displays that at half- lling,  = N=L = 1=2, thermal value , n = (1=L) jhn ^ (T )i j. The i max i=1 where N is the total number of particles and L is the deviation of late-time density from the thermal value, chain length, the level statistics approaches GOE with n, in the dense regime decays exponentially with the L= increasing system size, which is consistent with the de- system size as n  e , where  ' 6:27. In con- localized phase. In contrast, at  = 1=5 lling r ows trast, for the dilute case n shows no dependence on the av towards r at strong disorder. In what follows we x system size, as is apparent in the density pro les. The the disorder strength to be W = 6:5, since at this value characteristic length  extracted in the dense case gives the dilute limit is localized while the dense limit clearly the minimum size for genuine egodic bubbles that can ows towards delocalization. The scaling of entangle- destroy the MBME according to Ref. [18]. ment entropy also con rms the coexistence of localized Having con rmed the coexistence of localized and de- and delocalized phases for disorder W = 6:5 [29]. localized states at di erent values of particle density Quench dynamics.|Having provided numerical evi- for the same disorder strength, we proceed with a more dence for the coexistence of localized and delocalized detailed study of the e ect of a bubble, whose behavior is phases in small systems, we turn to quantum quench dy- central to the mechanism proposed in [18]. Figure 2(d) namics that distinguishes MBL from ergodic phase [5, illustrates the evolution of a non-uniform initial state, 32]. We consider quenches where the system is initially where a dense region represents the bubble. The bubble prepared in a product state and then evolved with the region consists of 8 sites with two pairs of particles and Hamiltonian (1). Starting with a density wave of pe- has a local density of  = 1=2. The bubble is followed riod 1=, a con guration that contains no pairs, we cal- by a period-5 density wave that occupies L 10 sites culate the density pro le at late times. For the dilute and two additional empty sites at the end of the chain. case,  = 1=5, we use the time-evolved block decimation Although having  = 8=30 > 1=5, this state is still in (TEBD) with MPS [33, 34] (see [29] for additional details a localized sector, as shown in [29]. The bubble leaks and benchmarks). This allows to monitor dynamics of only weakly into the dilute region even at late times, see systems as large as L = 40 sites up to times T = 500. Fig. 2(d), with particles far away from the bubble not be- max In the dense case ( = 1=2) we use ED and Krylov sub- ing a ected. In contrast, in the dense case, Fig. 2(g), the space time evolution method. While ED allows to access bubble with average density of  = 2=3 successfully melts the in nite-time density pro le, with the Krylov method, the period-3 density wave state throughout the system. we simulate quantum dynamics up to T = 1000. max Next, in panels Fig. 2(e) and (h) we further illustrate The density pro les at late times look very di erent in the di erences between the density dynamics in the dense 3 •◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦••◦••◦◦◦•◦◦•◦◦• ••◦◦••◦◦◦◦•◦◦◦◦•◦◦ (a) (d) 1 . 0 (g) L =40 L =25 L =15 1. 00 1. 00 L =30 L =20 L =30 L =20 L =10 0 . 8 (a) 0. 75 0. 75 0 . 6 i i i n n 0. 50 0. 50 n 0 . 4 0. 25 0. 25 0 . 2 0. 00 0. 00 0 . 0 5 10 15 20 25 30 35 40 1 4 7 10 13 16 19 22 25 28 31 1 3 5 7 9 11 13 15 17 i i ••◦••◦••◦◦◦•◦◦•◦◦• •◦•◦•◦•◦•◦•◦•◦•◦•◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦ ◦◦ (b) (e) (h) 0. 6 L =10 L =14 L =18 1. 00 0. 7 L =12 L =16 0. 5 0. 75 0. 6 0. 4 n ˜ 0. 50 0. 5 0. 3 0. 25 0. 4 0. 2 0. 00 −1 0 1 2 −1 0 1 2 3 1 3 5 7 9 11 13 15 17 10 10 10 10 10 10 10 10 10 ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦••◦••◦◦◦•◦◦•◦◦• ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦ ◦◦ (c) (f) (j) 0.4 2. 0 c =4−5 c =3 c =8−9 c =6 1. 5 0.2 c =12−13 c =9 1. 0 c =18−19 c =12 ν = 1 / 5 0.1 0. 5 c =24−25 0.08 c =15 ν = 1 / 2 0.06 0. 0 0 −1 0 1 2 −1 0 1 2 3 10 15 20 25 30 35 40 10 10 10 10 10 10 10 10 10 L t Figure 2. (a)-(c) The quantum quench from the uniform density wave with period 1= reveals memory of the initial state at  = 1=5 in (a), whereas in the dense case  = 1=2 (b) the charge pattern relaxes to zero exponentially in the system size as is shown in (c). (d)-(f) Stability of the dilute system against the bubble consisting of a half- lled region with 4 particles is illustrated in panel (d) by the density pro le at T . (e) The time dynamics of density in the coarse grained regions (see max the legend at the top) shows the absence of signi cant relaxation in regions away from the bubble. The density in the region at the boundary with the bubble increases logarithmically in time. (f ) The onset of logarithmic entanglement dynamics after a transient is visible for all cuts (see the legend at the top) away from the bubble. (g)-(j) In contrast, the bubble delocalizes the system when the overall density  = 1=2. The residual density pro le at T = 1000 in panel (g) has only weak memory max of the initial state. In addition, densities coarse grained over 3-site regions in (f ) all tend to the equilibrium value of 1/2 and entanglement entropy in (j) displays faster than logarithmic growth for all cuts. The data are generated using TEBD and 4 4 3 Krylov (ED) dynamics using between 5 10 {100 (dilute) and 3 10 {10 (dense) disorder realizations. and dilute cases in presence of a bubble. In both cases change of density [see Fig. 2(e)], and after saturation of we plot the density of particles within subregions of small density dynamics, we expect an onset of the logarithmic i+k1 growth of entanglement. In contrast, the entanglement size k,  ~ = (1=k) hn i, that are shown at the top i j j=i dynamics in Fig. 2(j) is always faster than logarithmic. of the plot. In the dilute case, Fig. 2(e), we observe that In [29] we provide more details on the contribution of ~ remains far from its thermal value even at late times, particle transport to entanglement [38, 39], demonstrat- in contrast with [18], where an ergodic region larger than ing that it is responsible for logarithmic entanglement is expected to delocalize the system. The densities of increase, in agreement with [40], whereas the con gura- regions in the bubble and adjacent to the bubble seem tional entanglement grows faster than logarithmic, and to saturate, while the regions far away from the bubble total entropy shows power-law increase. show very slow dynamics. In contrast, the dense case, Bubble tunneling vs. decay processes.|The quench dy- Fig. 2(h), shows that all expectation values evolve to- namics discussed above suggests that a bubble is not able wards equilibrium, although the regions far away from to spread through the entire localized chain and remains the center of the chain display slow, logarithmic in time, in the vicinity of its initial position. At the same time, growth of density. most of our quench simulations are restricted to nite, Finally, we study the dynamics of the bipartite en- albeit long, times. In order to give a complementary evi- tanglement entropy, S , see Fig. 2(f ) and (j). The en- vN dence for the bubble localization, we return to eigenstate tanglement is de ned as S = tr  ln , where  is vN properties that e ectively probe the in nite time limit. the density matrix of the left subregion calculated from We start with an initial product state in the half- lled iHt j (t)i = e j i. Di erent entanglement cuts shown case illustrated for L = 12, at the top of Fig. 2(f ) and (j) are encoded by their color. Consistent with MBL, the increase of entanglement in j i =  ; (2) the region close to the bubble is logarithmic in time in Fig. 2(f ) [35{38]. The entanglement across the cuts fur- that contains a bubble of k = L=2 sites with  = 2=3 ll- ther away from the bubble begins to grow at signi cantly ing (boxed region), followed by a sparser region with the later times. For these more distant cuts, the initial up- same number of sites and density  = 1=3. To quantify rise in entanglement corresponds to a slow logarithmic the relation between the probability of the bubble tunnel- ( ) S(t) 4 ing to the opposite end of the system and the probability where the bubble spreads throughout the system, calling of the bubble spreading throughout the system, we use into question the applicability of the resonance argument a spatial re ection of j i and uniform density wave as of [18]. a representative of the state with bubble tunneling and Experimental realization.| Finally, we discuss a pos- spreading, respectively: sible way to observe the physics related to MBME in ex- periments with ultracold atoms. Within the disordered =   ; (3) Aubry-Andr e bosonic Hamiltonian, j i = ; (4) h i 2 X H = t(a a + h.c.) +  n + Un (n 1) ; (6) i+1 i i; i i illustrated for L = 12 and  = 1=2 lling. For dilute con gurations at  = 1=5 we de ne the bubble as a region that is actively used to study MBL physics [38, 42], the of size 2(L1) with density  = 1=2, joined with a dilute remainder. For L = 20 such a state is: bubbles can be represented by regions with ha a i =  > j i =  . 1 bosons per site. A particle within such region has a It is straightforward to show that the in nite time av- hopping matrix element enhanced by the Bose-factor of erage probability of nding the system with the wave hi, thus playing the role of hopping t in model (1). iHt In the regime of densities and disorder strengths such function e j i in the product state j i is given by 1 2 that the enhanced hopping hit corresponds to local- X ization lengths signi cantly larger than lattice spacing, P (j i ;j i) = jhE j ihE j ij ; (5) 1 2 1 2  > a, whereas a single boson localization length is dense =1 . a, this model will implement similar physics to our toy model. Note that at the same time it is important where jE i are the complete set of eigenstates of H . to keep interaction U low enough, U  t, to avoid the Eq. (5) quanti es the similarity in the expansion of two formation of minibands related to long-lived doublons. di erent states j i over the basis jE i and reduces to 1;2 By initializing the system in a product state with a the conventional participation ratio when j i = j i. 1 2 dense region of bosons in the center of the trap along In order to reveal the relation between bubble decay with low density of bosons away from such a region, the and tunneling processes, we calculate the ratio of proba- dynamics under Hamiltonian (6) will probe the ability of s s bilities of bubble decaying, P = P (j i ;j i), with j i s 1 2 2 the bubble to melt the imbalance [5] away from its origi- from Eq. (4), to bubble tunneling, P = P (j i ;j i) t 1 nal position. From our simulations we expect the absence with j i from Eq. (3). In the dense case, these two of imbalance relaxation far away from the bubble. In a probabilities are of the same order and moreover tend to di erent direction, doublons [43, 44] or second species of identity with increasing system size as expected in the particles not subject to disorder [9] are also promising delocalized phase, see Fig. 3. In the dilute case, the ra- candidates that can play a role of the bubble. tio P =P is rapidly increasing with both disorder, and s t Discussion.|We presented a model with MBME in system size. This proves that the bubble tunneling pro- particle density and investigated its properties numeri- cesses are strongly suppressed compared to the processes cally using ED and time evolution with MPS. We nd strong evidence of the persistence of localization at in- nite times for small systems and also observe memory of initial con guration until times of T = 500 for sys- 10 max L =12 L =10 tems with up to L = 40 sites. These times are at least L =14 L =15 two orders of magnitude larger compared to the inverse L =16 L =20 10 local hopping, ~=t , and are achievable with cold atoms experiments. While we cannot rule out a residual very slow delocalization at much later times, the constructed model allows us to bound the timescale up to which the localization remains stable in very large systems that are − 1 0 1 2 3 4 5 6 7 8 beyond the reach for ED. The model with MBME in particle density presented in this work allows for direct tests of the arguments about the instability of MBME [18]. In order for bubble to move Figure 3. The rapid increase of the ratio of P =P with s t throughout the system it is important that the bubble system size and disorder strength reveals that in the dilute does not disappear by spreading and that con gurations case,  = 1=5 the probability for bubble spreading is strongly with bubbles situated at di erent locations are e ectively enhanced compared to the probability of bubble tunneling coupled to each other. Our simulations reveal that dilute to the opposite end of the system. For the dense case these systems have no trace of bubble reemerging at a di erent two probabilities are of the same order and approach each location within the system. Moreover, even the expecta- other with increasing system size in a broad range of disorders. tion value of the pair density hn n i (pairs are building Averaging is done over at least 2:5 10 disorder realizations. i i+1 blocks of the bubble) is exponentially suppressed away P /P t s 5 from the original location of the dense bubble [29]. In an supported by the European Unions Horizon 2020 research alternative approach, we directly test the probability of and innovation programme under the Marie Sklodowska- the bubble to emerge at the opposite end of the system Curie Grant Agreement No. 665385. D.A. was supported at in nite time and nd it to be strongly suppressed. by the Swiss National Science Foundation. M.S. was To conclude, we expect that the proposed model will supported by European Research Council (ERC) un- enable further investigations of particle density MBME. der the European Union's Horizon 2020 research and Studies of the structure of matrix elements, extension of innovation programme (grant agreement No. 850899) the theory of LIOMs [45, 46] to systems with MBME in This work bene tted from visits to KITP, supported particle density [22], and studies of the e ect of a small by the National Science Foundation under Grant No. bath on a localized system [47{50] using our model rep- NSF PHY-1748958 and from the program \Thermal- resent promising avenues for future work. ization, Many body localization and Hydrodynamics" Acknowledgments.|We acknowledge useful discus- at International Centre for Theoretical Sciences (Code: sions with W. De Roeck and A. Michailidis. P.B. was ICTS/hydrodynamics2019/11). [1] L. Fleishman and P. W. Anderson, \Interactions and the uan Song, Hekang Li, Zhen Wang, Wenhui Ren, anderson transition," Phys. Rev. B 21, 2366{2377 (1980). Hang Dong, Dongning Zheng, Yu-Ran Zhang, Rubem [2] D.M. Basko, I.L. Aleiner, and B.L. Altshuler, \Metal{ Mondaini, Heng Fan, and H. Wang, \Observation of en- insulator transition in a weakly interacting many-electron ergy resolved many-body localization," arXiv e-prints , system with localized single-particle states," Annals of arXiv:1912.02818 (2019), arXiv:1912.02818 [quant-ph]. Physics 321, 1126 { 1205 (2006). [12] Titas Chanda, Piotr Sierant, and Jakub Zakrzewski, [3] Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and \Many-body localization transition in large quantum Maksym Serbyn, \Colloquium: Many-body localization, spin chains: The mobility edge," arXiv e-prints thermalization, and entanglement," Rev. Mod. Phys. 91, , arXiv:2006.02860 (2020), arXiv:2006.02860 [cond- 021001 (2019). mat.dis-nn]. [4] Rahul Nandkishore and David A. Huse, \Many-body lo- [13] Wojciech De Roeck and Francois Huveneers, \Stability calization and thermalization in quantum statistical me- and instability towards delocalization in many-body lo- chanics," Annual Review of Condensed Matter Physics calization systems," Phys. Rev. B 95, 155129 (2017). 6, 15{38 (2015). [14] Thorsten B. Wahl, Arijeet Pal, and Steven H. Simon, [5] Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, \Signatures of the many-body localized regime in two Henrik P. Lusc  hen, Mark H. Fischer, Ronen Vosk, Ehud dimensions," Nature Physics 15, 164{169 (2019). Altman, Ulrich Schneider, and Immanuel Bloch, \Obser- [15] Ionut-Dragos Potirniche, Sumilan Banerjee, and Ehud vation of many-body localization of interacting fermions Altman, \Exploration of the stability of many-body lo- in a quasirandom optical lattice," Science 349, 842{845 calization in d > 1," Phys. Rev. B 99, 205149 (2019). (2015). [16] Augustine Kshetrimayum, Marcel Goihl, and Jens [6] Pranjal Bordia, Henrik P. Lusc  hen, Sean S. Hodgman, Eisert, \Time evolution of many-body localized sys- Michael Schreiber, Immanuel Bloch, and Ulrich Schnei- tems in two spatial dimensions," arXiv e-prints, der, \Coupling identical one-dimensional many-body lo- arXiv:1910.11359 (2019), arXiv:1910.11359 [cond- calized systems," Phys. Rev. Lett. 116, 140401 (2016). mat.str-el]. [7] Jae-yoon Choi, Sebastian Hild, Johannes Zeiher, Peter [17] Elmer V. H. Doggen, Igor V. Gornyi, Alexander D. Schau, Antonio Rubio-Abadal, Tarik Yefsah, Vedika Mirlin, and Dmitry G. Polyakov, \Slow many- Khemani, David A. Huse, Immanuel Bloch, and body delocalization beyond one dimension," arXiv e- Christian Gross, \Exploring the many-body localization prints, arXiv:2002.07635 (2020), arXiv:2002.07635 [cond- transition in two dimensions," Science 352, 1547{1552 mat.dis-nn]. (2016). [18] Wojciech De Roeck, Francois Huveneers, Markus Muller,  [8] B. Chiaro et al., \Growth and preservation of entan- and Mauro Schiulaz, \Absence of many-body mo- glement in a many-body localized system," arXiv e- bility edges," Physical Review B 93, 1{14 (2016), prints, arXiv:1910.06024 (2019), arXiv:1910.06024 [cond- arXiv:1506.01505. mat.dis-nn]. [19] N.F. Mott, \Electrons in disordered structures," Ad- [9] Antonio Rubio-Abadal, Jae Yoon Choi, Johannes Zeiher, vances in Physics 16, 49{144 (1967). Simon Hollerith, Jun Rui, Immanuel Bloch, and Chris- [20] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Cri- tian Gross, \Many-Body Delocalization in the Presence terion for many-body localization-delocalization phase of a Quantum Bath," Physical Review X 9, 41014 (2019), transition," Phys. Rev. X 5, 041047 (2015). arXiv:1805.00056. [21] David J. Luitz, Nicolas La orencie, and Fabien Alet, [10] Thomas Kohlert, Sebastian Scherg, Xiao Li, Henrik P. \Many-body localization edge in the random- eld heisen- Lusc  hen, Sankar Das Sarma, Immanuel Bloch, and berg chain," Phys. Rev. B 91, 081103 (2015). Monika Aidelsburger, \Observation of many-body local- [22] Scott D. Geraedts, R. N. Bhatt, and Rahul Nandkishore, ization in a one-dimensional system with a single-particle \Emergent local integrals of motion without a complete mobility edge," Phys. Rev. Lett. 122, 170403 (2019). set of localized eigenstates," Phys. Rev. B 95, 064204 [11] Qiujiang Guo, Chen Cheng, Zheng-Hang Sun, Zix- (2017). 6 [23] Marcel Goihl, Jens Eisert, and Christian Krumnow, \Ex- localization transition," Nature 573, 385{389 (2019). ploration of the stability of many-body localized systems [43] Ulrich Krause, Th eo Pellegrin, Piet W. Brouwer, in the presence of a small bath," Phys. Rev. B 99, 195145 Dmitry A. Abanin, and Michele Filippone \Nucle- (2019). ation of ergodicity by a single doublon in supercooled [24] Ruixiao Yao, and Jakub Zakrzewski, \Many-body lo- insulators," arXiv e-prints , arXiv:1911.11711 (2019), calization in Bose-Hubbard model: evidence for the mo- arXiv:1911.11711 [cond-mat.dis-nn]. bility edge ," arXiv e-prints , arXiv:2002.00381 (2020), [44] Micheal Schecter, Thomas Iadecola, and Sankar Das arXiv:2002.00381 [condmat.dis-nn]. Sarma, \Con guration-controlled many-body localiza- [25] P. W. Anderson, \Absence of di usion in certain random tion and the mobility emulsion," Phys. Rev. B 98, lattices," Phys. Rev. 109, 1492{1505 (1958). 174201 (2018). [26] D. L. Shepelyansky, \Coherent propagation of two inter- [45] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Lo- acting particles in a random potential," Phys. Rev. Lett. cal conservation laws and the structure of the many-body 73, 2607{2610 (1994). localized states," Phys. Rev. Lett. 111, 127201 (2013). [27] Y Imry, \Coherent propagation of two interacting parti- [46] David A. Huse, Rahul Nandkishore, and Vadim cles in a random potential," Europhysics Letters (EPL) Oganesyan, \Phenomenology of fully many-body- 30, 405{408 (1995). localized systems," Phys. Rev. B 90, 174202 (2014). [28] Xiaopeng Li, Sriram Ganeshan, J. H. Pixley, and [47] Rahul Nandkishore, \Many-body localization proximity S. Das Sarma, \Many-body localization and quantum e ect," Phys. Rev. B 92, 245141 (2015). nonergodicity in a model with a single-particle mobility [48] David J. Luitz, Francois Huveneers, and Wojciech edge," Phys. Rev. Lett. 115, 186601 (2015). De Roeck, \How a small quantum bath can thermal- [29] \Supplemental Online material,". ize long localized chains," Phys. Rev. Lett. 119, 150602 [30] Z. Papi c, E. Miles Stoudenmire, and Dmitry A. Abanin, (2017). \Many-body localization in disorder-free systems: The [49] Katharine Hyatt, James R. Garrison, Andrew C. Potter, importance of nite-size constraints," Annals of Physics and Bela Bauer, \Many-body localization in the presence 362, 714 { 725 (2015). of a small bath," Phys. Rev. B 95, 035132 (2017). [31] Arijeet Pal and David A. Huse, \Many-body localization [50] Philip J. D. Crowley and Anushya Chandran, \Avalanche phase transition," Phys. Rev. B 82, 174411 (2010). induced co-existing localised and thermal regions in [32] Maksym Serbyn, Z. Papi c, and D. A. Abanin, \Quantum disordered chains," arXiv e-prints , arXiv:1910.10812 quenches in the many-body localized phase," Phys. Rev. (2019), arXiv:1910.10812 [cond-mat.dis-nn]. B 90, 174302 (2014). [51] Elmer V. H. Doggen, Frank Schindler, Konstantin S. [33] Guifr e Vidal, \Ecient classical simulation of slightly en- Tikhonov, Alexander D. Mirlin, Titus Neupert, tangled quantum computations," Phys. Rev. Lett. 91, Dmitry G. Polyakov, and Igor V. Gornyi, \Many-body 147902 (2003). localization and delocalization in large quantum chains," [34] ITensor Library (version 3.1.1) http://itensor.org. Phys. Rev. B 98, 174202 (2018). [35] M. Znidaric, T. Prosen, and P. Prelovsek, \Many-body localization in the heisenberg xxz magnet in a random eld," Phys. Rev. B 77, 064426 (2008). [36] Jens H. Bardarson, Frank Pollmann, and Joel E. Moore, \Unbounded growth of entanglement in models of many- body localization," Phys. Rev. Lett. 109, 017202 (2012). [37] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Universal slow growth of entanglement in interacting strongly disordered systems," Phys. Rev. Lett. 110, 260601 (2013). [38] Alexander Lukin, Matthew Rispoli, Robert Schittko, M. Eric Tai, Adam M. Kaufman, Soonwon Choi, Vedika Khemani, Julian L eonard, and Markus Greiner, \Prob- ing entanglement in a many-body{localized system," Sci- ence 364, 256{260 (2019). [39] Rajibul Islam, Ruichao Ma, Philipp M Preiss, M Eric Tai, Alexander Lukin, Matthew Rispoli, and Markus Greiner, \Measuring entanglement entropy in a quantum many-body system," Nature 528, 77 (2015). [40] Maximilian Kiefer-Emmanouilidis, Razmik Unanyan, Michael Fleischauer, and Jesko Sirker, \Evidence for unbounded growth of the number entropy in Many-Body Localized phases," Phys. Rev. Lett. 124, 243601 (2020). [41] David J. Luitz, and Yevgeny Bar Lev \Is there slow particle transport in the MBL phase?," arXiv e-prints , arXiv:2007.13767 (2020), arXiv:2007.13767 [cond-mat.dis-nn]. [42] Matthew Rispoli, Alexander Lukin, Robert Schittko, Sooshin Kim, M. Eric Tai, Julian L eonard, and Markus Greiner, \Quantum critical behaviour at the many-body S1 Supplementary material for \Particle density mobility edge" In this supplementary material we present additional data and details of the methods used in the main text. First, we show data for additional probes of ETH breakdown such as entanglement entropy. After this we present benchmarks of our time-evolving block decimation simulation of dynamics. Finally, we explore the behavior of the mutual IPR de ned in the main text. I. LOCALIZATION LENGTH AND ized, while at lower densities   1=5, when the typical PARAMETER CHOICE distance between pairs is large, we expect MBL phase. This motivates the choice t = 2, since at this value of t 2 2 (W ) approaches 2 at disorder strength around W  6. The dynamics generated by the constrained Hamilto- We note that we avoided further increase of t to keep the nian, Eq. (1), strongly depends on the choice of the hop- model away from the constrained limit: in the case when ping parameters t . In order to choose the most suitable 1;2 t dominates over t , the model would approximately re- 2 1 parameters for the study of MBME in particle density, we duce to a kinematically constrained model that has many explore localization lengths for a single particle  and SP disconnected sectors in the Hilbert space. for one pair of particles  . These localization lengths In order to rule out the presence of strong nite size are evaluated using ED. We calculate the in nite-time e ects, we studied the density of state in individual dis- average of the occupation number at each site for an ini- order realizations. In the regime when t  t the strong 2 1 tial state where either a single particle or a single pair nite size e ects would give rise to the presence of the are initialized at the rst site of the chain. We extract mini-bands and the DOS would become non-monotonous the localization lengths  ( ) from an exponential t SP P with numerous peaks corresponding to mini-band struc- of the density curve hn i . ture [30]. Figure S2 con rms that at our choice of pa- Resulting values of  for xed t = 0:5 and dif- P,SP 1 rameters even individual disorder realizations have a rel- ferent disorder values and di erent values of hopping t atively smooth density of states with Gaussian envelope, are shown in Fig. S1. The single particle hopping local- thus ruling out the presence of strong nite size e ects. ization length (dashed line in Fig. S1) does not depend on t , and becomes smaller than one lattice spacing for W & 4. The pair localization length is monotonously II. ED PROBES OF LOCALIZATION increasing with t at xed value of disorder strength, W . Our aim is to have  in the range between 2 and 5. In While in the main text we focused on the two values of this regime, the half- lling case is expected to be delocal- lling,  = 1=2 and 1=5, here we demonstrate the density dependence of critical disorder. For this purpose we cal- culate the average ratio of level spacings r for a single av system size L = 18 at varying values of density. Fig- ξ ξ , t =2 SP p 2 ξ , t =1 ξ , t =2.5 p 2 p 2 ξ , t =1.5 ξ , t =3 p 2 p 2 0 . 08 0 . 06 0 . 04 0 . 02 1 2 3 4 5 6 0 . 00 −20 −15 −10 −5 0 5 10 15 20 Figure S1. The localization length decreases, as expected, E with the disorder strength for all the values of t . For every constrained hopping amplitude t it is possible to locate the region of disorder where we expect to see a MBME in parti- Figure S2. The DOS from single disorder realizations show a cle density as the area among the two dashed lines. As the relatively smooth behavior and a Gaussian shape, thus con- curve crosses the rst dashed line, systems with typical par- rming the absence of strong nite size e ects. DOS refers ticle spacing 5 will be localized. Nevertheless denser states to a chain with L = 20 and  = 1=4. Disorder strength is will still be delocalized, having smaller distance among par- W = 5:0. Green, blue and orange curves correspond to dif- ticles. Data were obtained on a lattice of length L = 50 and ferent disorder realizations, while the black dashed line shows averaged over 5000 disorder realizations. disorder-averaged DOS. S2 0.53 L = 10 L = 12 0.30 0.8 0.51 L = 14 0.7 0.25 0.49 L = 16 0.6 0.20 0.47 ν 0.5 av 0.15 0.45 0.4 0.43 0.10 0.3 0.41 L = 10 0.2 0.05 L = 15 L = 20 0.39 0.00 2 4 6 8 1 2 3 4 5 6 7 8 Figure S3. The sharp di erence of r obtained for di erent av Figure S4. The behavior of half-chain entanglement entropy at the same disorder W clearly shows the MBME in our shows very distinct behavior for dilute (blue-shaded curves) model. Interestingly, the mobility edge curve W () is not and dense (red-shaded curves) states. The crossing in the symmetric, but is peaked around  = 2=3, implying that the dilute states implies that they entered the MBL phase, and states with the maximum number of pairs for xed size are thus have area-law entanglement entropy. On the other hand, the hardest to localize. The data is obtained for a system of dense states do not show a similar crossing in this range of size L = 18, using shift-invert method with 10 10 states disorder, suggesting that they are still in the ergodic phase. 4 3 from the middle of the spectrum and 5  10 10 disorder The data are obtained with shift-invert method for 10 10 realizations. eigenstates in the middle of the spectrum and averaged over 4 3 5 10 5 10 disorder realizations. ure S3 allows to estimate the dependence of the critical disorder on the lling, . At low densities ( <  (W )) sector of the Hamiltonian Eq. (1). Given the U (1) sym- states have r approaching value characteristic for Pois- av metry, we would expect that a xed lling sector presents son distribution of level spacings. In contrast, for dense MBME in energy density, similarly to the case of the ran- con gurations ( >  (W )) the level spacing ratio is close dom eld XXZ spin chain. In the middle of the spectrum to GOE prediction. the density of states is large and eigenstates may remain Figure S3 reveals that the most delocalized lling is delocalized, while at the same disorder strength the states = 2=3, which corresponds to the case when the best at the edges of the many-body spectrum are localized. packing of pairs in the chain,  , can be achieved. To explore the eventual presence of MBME in energy At this lling the Poisson values of r would be achieved av density in the half lling sector, we studied the energy beyond the upper limit of the considered disorder range. resolved level spacing ratio r (; W ). The results, dis- Av Decreasing particle density away from this value causes earlier onset of localization. For instance, xing disorder value W = 6:5 we observe that  = 1=2 and  = 1=5 are 1.0 situated well in delocalized and localized regions. 0.52 L =10 In addition to the level statistics indicator presented in 0.8 0.50 L =12 the main text, we studied other commonly used probes of L =14 0.48 0.6 ergodicity. In particular, Fig. S4 illustrates the behavior L =16 0.46 of bipartite entanglement entropy for di erent disorder 0.4 0.44 strengths and di erent llings. On the one hand, the - nite size scaling of entanglement entropy of eigenstates in 0.42 0.2 the middle of the spectrum shows that for  = 1=5 and 0.40 disorder W > W  6 the entanglement is consistent 0.0 3 4 5 6 with area-law. On the other hand, the entanglement of dense systems,  = 1=2, does not show a similar behav- ior. The nite size scaling, indeed, shows no crossing Figure S5. The plot, presenting the energy resolved level at these disorder values, thus suggesting volume-law of spacing ratio for L = 16 in the half lling sector, shows clear evidence of MBME. In the center of the band r approaches entanglement entropy for  = 1=2. av the GOE value, while at large and small energy density it is close to the Poisson value. Furthermore, the MBME curves obtained for smaller systems seem to converge at increasing A. Energy Density MBME system size, thus suggesting the persistence of the MBME in the thermodynamic limit. The plot was obtained averaging 4 3 2 2 over n = 10 ; 2 10 ; 5 10 ; 10 disorder realizations for In this section we brie y discuss the presence of many- increasing system size from L = 10 to L = 16. body mobility edge in energy density in a single density S/L Av S3 (a) (c) 10 − 12 δt = 0.1 ε = 10 − 9 0.005 δt = 0.1 ε = 10 − 3 0.000 − 6 − 12 δt = 0.05 ε = 10 − 9 − 9 − 0.005 10 δt = 0.05 ε = 10 t =1 t =200 − 9 δt = 0.02 ε = 10 t =100 t =300 − 12 − 0.010 0.02 0.03 0.05 0.08 0.1 1 4 7 10 13 16 19 δt i (b) (d) 0.02 − 12 δt =0.02 δt =0.1 δt = 0.1 ε = 10 − 2 − 9 0.01 δt =0.05 10 δt = 0.1 ε = 10 δ − 5 0.00 − 12 − 8 − 0.01 δt = 0.05 ε = 10 − 9 δt = 0.05 ε = 10 − 11 − 0.02 10 − 9 δt = 0.02 ε = 10 0 1 2 10 10 10 1 4 7 10 13 16 19 t i Figure S6. (a-b) Deviation of ground state delity from 1 in Trotter time evolution, F (t; t), shows power-law behavior both in time and in time-step, as expected. The data is obtained at density  = 1=5, system size L = 20 and disorder W = 6:5 for a particular disorder realization. The plots for other disorder realizations are qualitatively similar. (c-d) The comparison between ED and TEBD time evolution reveals that the most e ective way to increase accuracy of TEBD is to decrease the time-step t. 9 12 Indeed, the change in the truncation between " = 10 and " = 10 does not have much e ect on the the di erence between density pro les of exact diagonalization and TEBD. At the same time, the decrease of time step brings the local density pro le closer to ED results. The density pro les are calculated by propagating uniform density wave (c) and uniform pair-density wave (d) initial states to time t = 500 for a particular disorder realization with L = 20,  = 1=5, W = 6:5. played in gure S5, show evidence of many-body mobility related to the nite size of the time step in the p-th order edge; the level spacing ratio, as a function of the energy Trotter expansion grows as t . The other source of error EE min density  = , where E and E are the is the nite cuto , ", that governs the truncation of sin- min max E E max min ground state and the most excited state respectively, in- gular values in the singular value decomposition (SVD). creases from the Poisson to the GOE value as  goes While the instantaneous errors related to the trunca- from the lower edge to the center of the spectrum and tion and nite time step are known, understanding the decreases again from the center to the upper edge. This propagation of these errors with time and their possi- variation is such that at a xed disorder strength the ble interference is challenging. First we tested TEBD low and high energy states are localized, while the cen- algorithm by evolving the ground state of the same ter of the band is delocalized, thus de ning a many-body model. Provided that the time evolution is numerically mobility edge. The scaling of the MBME curves for dif- exact, the overlap between the TEBD-evolved ground TEBD ferent system size shows signs of convergence, suggesting state, j (t)i = U (t)jGSi and the exact time evo- {E t stability of the MBME in the thermodynamic limit. lution of the ground state, jGS(t)i = e jGSi, is supposed to give the identity h (t)jGS(t)i = 1 at all times. For the fourth-order Trotterization the behavior III. MPS SIMULATIONS OF QUENCH of F = 1jh (t)jGS(t)ij is known to be proportional DYNAMICS to (t) . The numerical results plotted in Fig. S6(a), con rm these expectations. In our MPS simulation, we time evolve dilute states in Next, we performed a benchmarking of TEBD algo- large systems L  30 up to time T = 500. For this rithm against ED time evolution for several disorder re- max we use the time-evolving block decimation (TEBD) algo- alizations and simulation parameters. An illustration of rithm with a fourth-order Trotter evolution based on the such benchmarking is shown in Fig. S6(c) and (d). In ITensor library [34]. The main parameter involved in the particular, we observed that time step t = 0:05 and cut- time evolution algorithm is the time step t used to split o " = 10 result in a good agreement between ED and the unitary evolution into a sequence of gates. The error TEBD dynamics. Smaller values of t; " would improve F (t, t ) F (t, t ) ED TN ED TN hn i h n i i i hn i h n i i i S4 IV. ADDITIONAL RESULTS FROM QUENCH • • ◦ ◦ • • ◦ ◦ ◦ ◦ • ◦ ◦ . . . 0.0005 DYNAMICS • ◦ ◦ ◦ ◦ • ◦ ◦ ◦ ◦ • ◦ ◦ . . . 0.0004 In the main text we discussed dynamics in quenches that begin from a uniform state or a bubble joined to a 0.0003 more dilute remainder. Below we present details for the 0.0002 quenches in presence of bubble. In addition, we discuss the dynamics resulting from the initial state containing 0.0001 density wave of particle pairs. 0.0000 0 100 200 300 400 500 A. Pair density and entanglement dynamics in presence of a bubble Figure S7. The normalized absolute value of the energy di er- Since particle pairs are the most mobile objects, we ence from the initial energy E(t) = jhH (t)i E(0)j=jE(0)j consider the pair density in quenches that are initialized remains very small for both the density wave, blue curve, with the bubble (see Fig. 2(d),(g) in the main text). The and the non-uniform, red curve, con gurations, con rming pair density is of special interest in these quenches as in the good accuracy of our numerical simulations beyond the Ref. [18] suggested that the instability of the system is ED benchmark. The larger deviation displayed by the non- ascribed to the ability of the bubble to move. In our uniform con guration is understood as a result of the pres- ence of the bubble in the lattice, that increases entanglement model the bubble consists of several pairs, thus motion growth. The results here shown are obtained averaging over of the bubble throughout the system would imply the 100 disorder realizations, for the system sizes, L = 30, and spreading of pairs. initial states described in the main text. The pair density de ned as hn n i measured at late i i+1 or in nite times is shown in Fig. S9. In the dense case the late time pair density pro le supports delocalization: the agreement but would result in a dramatic slowdown at late times the density of pairs becomes homogeneous of the evolution time. Therefore, we decided to use these throughout the formerly more dilute region of the system. parameters in the simulations presented in the main text. We note, that the pair density is not a conserved quantity, When larger system sizes are involved, as it is the and it can increase in the process of unitary dynamics. case for the simulations actually used in the main text, In contrast, for the dilute case the pair density pro le comparison with exact results is not available. There- has a pronounced exponential tail away from the initial fore, other indicators for the accuracy must be stud- ergodic region. This shows that pairs spreading away ied. Among these, energy conservation through the time from the initial bubble do not delocalize when encoun- evolution is a straightforward probe. The energy de- tering additional particles on their way. Indeed, while the viation E(t) = jE(0) E(t)j=E(0), where E(t) = late time pair density pro le has small peaks around the initial position of particles, these peaks are not very pro- h (t)j Hj (t)i, allows to control the propagation of the nounced. In addition, the study of the pair density pro le error during the Trotter time evolution. In Fig. S7 we show the results for E(t) in the two quenches presented in the main text, for L = 30. The two plots highlight that the average energy deviation is very small in both 17.5 con gurations. In spite of that, a clear di erence can be 15.0 observed among the two quenches, noticing that the non- 12.5 uniform state has larger error. This is probably due to the enhanced entanglement caused by the presence of the 10.0 bubble in the lattice. Nevertheless, E(t) remains very 7.5 small even at long times, thus con rming the reliability 5.0 of our long-time numerical simulations. In all the simulations performed using ITensor [34], we 2.5 used the U (1) symmetry implementation. In particular, 0.0 200 400 500 1000 1500 2000 to obtain the numerical results presented in Fig. 2(a) and Bond Dimension Bond Dimension (d) we set the maximum bond dimension to be 500 and 3000 respectively. As the histograms in Fig. S8 show, all the disorder realizations remained well below the maxi- Figure S8. The histogram representing the maximum bond mum threshold. This fact ensures that we have a control dimension of di erent disorder realizations show that the on the error encountered in the evolution, in contrast to threshold values of 500 and 3000 for the uniform density wave time evolution with TDVP with xed bond dimension, (left) and bubble states (right) were never saturated in our where error estimation is more challenging [51]. simulations. E(t) Counts Counts S5 • • ◦ • • ◦ • • ◦ ◦ ◦ • ◦ ◦ • ◦ ◦• L = 10 L = 12 vN L = 14 − 1 L = 16 L = 18 − 2 L = 20 L = 30 − 3 1 4 7 10 13 16 19 22 25 28 − 1 0 1 2 3 10 10 10 10 10 Figure S9. The nite size scaling of the pair density hn n i i i+1 shows opposite trend for the dense and dilute cases. The red- Figure S10. The di erent contributions to the entanglement shaded curves represent  = 1=2 con gurations: increasing entropy of the bubble show that the overall behavior of the the system size (from yellow to dark red) the pair density be- von Neumann entropy is faster than logarithmic. Neverthe- comes more uniform and approaches the thermal value, hence less this behavior can be ascribed to the sole con gurational in the thermodynamic limit the probability of nding a pair entropy S , while the particle transport contributes to the far from the bubble is almost the same as nding it in the bub- purely logarithmic growth. The curves are obtained through ble. On the contrary, blue curves ( = 1=5) show exponential Krylov evolution up to T = 1000 averaged over 1000 dis- max vanishing of the pair density and, furhtermore, increasing sys- order realizations for L = 18 and W = 6:5. tem size (from light blue to dark blue) the density decreases, suggesting that at the thermodynamic limit there will be no pair outside the thermal region. Data were obtained with ED, con gurations with the same particle number. Interest- Krylov (T = 1000) and TEBD (T = 500) algorithms max max ingly, Fig. S10 shows that while the overall entanglement averaging over 100 disorder samples for the largest MPS sim- entropy grows faster than logarithmic, this is due only to 4 4 3 3 ulations (L = 20; 30), 3  10 ; 10 ; 5  10 and 10 for ED 3 the con guration part (yellow curve) and the entangle- (from L = 10 to L = 16) and over 10 for Krylov algorithm ment due to particle transport has logarithmic growth. (L = 18). The logarithmic growth of S is consistent with the log- arithmic particle transport presented in Fig. 2(h) in the main text and with other transport measures presented in the uniform density wave at  = 1=5 reveals an almost in the next section. We identify this behavior as a hall- constant behavior, centered around hn n i  10 , i i+1 mark of MBME, and note that it happens on long, yet which corresponds to the values reached at the end of the experimentally accessible timescales t  50(~=t ). exponential tail in the system with L = 30 in Fig. S9. Recent work [40] demonstrated that the logarithmic Next, we focus on understanding di erent contribu- tions to entanglement growth. Exploiting the U (1) sym- metry of our model and following Refs. [38, 39], we split the von Neumann entanglement entropy into a con g- vN uration and a particle transport contributions. Indeed, due to conservation of the total number of bosons the full reduced density matrix  must have a block-diagonal (n) form. Individual blocks within  can be written as p  , where p gives the probability to have n particles in the (n) (n) 0.5 susbsystem A and  is normalized as tr  = 1. Us- ing such representation of the reduced density matrix we 0.2 can split the full entropy into S = S + S as: vN C n 0.1 − 1 0 1 2 3 (n) (n) 10 10 10 10 10 S = tr  log  = p tr  log p vN n n X X (n) (n) (S1) = p log p p tr  log n n n Figure S11. The log-log plot of the entropies discussed in n n Eq. (S1) con rms that only the total entropy is growing as a = S + S : power-law, while both S and S grow slower. In particular, n C n c S grows as a rst degree polynomial in log(t) (dashed green In this way the entanglement growth is split into two line) and S as a second degree polynomial in log(t) (orange dashed line). The sum of these two behaviors (red dashed contributions: one coming from the particle transport, line) agrees with the power-law behavior (dashed blue line). and another originating from dephasing between di erent hn n i i i+1 S(t) S(t) S6 L/2 L growth of the number entropy is expected in the thermal phase, provided there is particle transport over distances that increase as a power-law in time, l / t . The au- (a) thors also predict a power-law scaling of the con guration entropy, which we do not observe, as shown in gure S11. 0.015 We attribute the slower than power-law growth of con g- 0.010 uration entanglement to the localized nature of the right half of the chain, which in turn reduces the number of 0.005 possible con gurations until the particle transport from the left half leads to delocalization. 0.000 0 1 Further analysis of the entanglement dynamics shows 10 10 that S grows in a power-law fashion S  at , cor- vN vN responding to the dashed blue t in Fig. S11, over a rel- (b) evant time interval. On the other hand, both S and 0.2 S behave as polynomials in log(t), of rst and second degree respectively, and their sum (dashed red line) rea- sonably approximates the power-law behavior of S , as vN 0.1 b 2 one can expand t  1 + b log(t) + log (t). Finally, to support our interpretation of logarithmic in- 0.0 crease of S as due to transport, we study the dynamics − 1 0 1 2 3 10 10 10 10 10 of density correlation functions and uctuations. Fig- ure S12 presents the dynamics of connected correlation (c) functions C (i; t) = hn ^ n ^ ihn ^ihn ^ i with respect to i i L=2 L=2 the central site of the chain, the local density uctua- 1.0 2 2 tions n = hn ^ i hn ^ i and the density uctuations in i i 2 2 the dilute part of the chain n = hn ^ ihn ^ i , where R R n ^ = n ^ . The logarithmic dynamics of these R i 0.5 i=L=2 W =5.0 W =6.0 quantities is consistent with the behavior of number en- W =5.5 W =6.5 tropy, thus proving the further support for the existence of slow transport in the dilute part of the chain. 0.0 − 1 0 1 2 3 10 10 10 10 10 Figure S12. Time dynamics of correlation functions (a), lo- cal density uctuations (b) and density uctuations in the B. Quench dynamics from a pair density wave state dilute half of the chain (c) all show, after an initial power-law growth, logarithmic increase with time. In particular, corre- lations far away from the central site (blue curves, as encoded Below we consider quench from a pair-density wave of in the legend above) show signs of a logarithmic light-cone. period 2=. These con gurations accommodate the max- Similarly, local density uctuations deep in the localized re- imal possible number of pairs in the uniform state. Fig- gion present slower dynamics. Finally, panel (c) shows how ure S13 con rms that such state is localized at  = 1=5 increasing disorder slows the growth of the global density uc- and is relaxing in the dense case. Dense systems display tuation of the dilute half. These results were obtained with the Krylov method on system size L = 18 for W = 6:5 (a), strong dependence on the system size and increased ten- with ED on system size L = 16 and W = 6:5 (b) and for dency towards relaxation at larger system sizes, L. In di erent disorder values (c), averaging over 200 disorder real- contrast, at  = 1=5 the late time density pro le has al- izations. most no dependence on the size of the system. In particu- lar, even at very large lengths the curves do not approach the average density represented by the dashed black line. V. P AS A FUNCTION OF BUBBLE DISTANCE In the text we introduced a quantity P (j i ;j i) that 1 2 measures how similar the expansion of over eigen- 1;2 states is. We bire y mentioned that the case where We note, that although the authors of Ref. [40] report unbounded j i = j i corresponds to the usual participation ratio, 1 2 growth of the number entanglement in the MBL phase, the suc- hence we refer to the inverse of P as mutual inverse par- cessive work of Luitz and Bar Lev [41] shows that this is due to ticipation ratio (mIPR). The mutual IPR assumes very rare particle uctuations around the boundary between the two subsystems and the growth disappears at large enough system di erent values depending on the nature of the two states: size. values of mIPR O(N ) correspond to two vectors that n (t) C(i, t) R n (t) i S7 ••◦◦◦◦◦◦◦◦••◦◦◦◦◦◦ ◦◦ (a) we illustrate the behavior of mIPR between pair of states L =40 L =20 L =10 1 . 00 which correspond to a smaller bubble displacement. L =25 L =15 In our analysis we measure the mIPR, P = 0 . 75 P ( ; ), between the following states in the dense L d 0 . 50 limit (half- lling, L = 12), 0 . 25 j i =  ; (S2) 0 . 00 j i =    : (S3) 5 10 15 20 25 30 35 40 |{z} i d ••◦◦••◦◦••◦◦•• ◦◦ (b) 1 . 0 Here we use the bubble that contains all particles to max- imize the range of achievable displacements. For the di- L = 10 L = 14 L = 18 0 . 8 lute case, L = 15, we use similar pair of states with 0 . 6 bubble containing 3 particles ( = 1=5). In the thermal phase, eigenstates are approximately 0 . 4 given by random vectors in the Hilbert space and 0 . 2 L = 12 L = 16 their average overlap with other normalized vectors ap- 0 . 0 proaches the value predicted by random matrix theory, 1 3 5 7 9 11 13 15 17 irrespective of the state or the eigenstate. In the weak disorder limit, we then expect P to be independent (c) on the distance between the two bubbles and to have n = 1 / 5 the same behavior as the conventional IPR: P  N . n = 1 / 2 This expectation is con rmed by the results presented in g. S14(a) and (b) for W = 0:5. On the other hand, in the MBL phase eigenstates are −1 not similar to random vectors, but instead are character- ized by a set of local integrals of motion that have a nite 10 15 20 25 30 35 40 overlap with the local particle density. Thus, two prod- uct states with globally di erent arrangement of particles are expected to have drastically di erent expansion over eigenstates. Therefore, we expect P / exp [d=]. As Figure S13. The late time density pro les of the pair density presented in gure S14(a) and (b), at strong disorder our waves at dense,  = 1=2, and dilute,  = 1=5, llings show results support this hypothesis for both dilute (a) and very di erent behavior. (a) Dilute con gurations are essen- dense (b) states. tially frozen, and do not approach the thermal density repre- At intermediate disorder strength, we observe a quali- sented by the black dashed line. (b) In contrast, at  = 1=2 tative di erence between dense and dilute cases. Dilute relaxation is enhanced at larger L. (c) The deviation of late time density from the thermal value decay exponentially with con gurations, Fig. S14(a), show exponential behavior pair L= pair T already at W = 4:5, whereas dense states in Fig. S14(b) system size as e with   8:1. In contrast for = 1=5, the residual density remains nearly constant with need much stronger disorder to clearly present the same system size and na ve t to the exponential gives an order trend. This result con rms the presence of mobility dilute of magnitude larger scale,   84. Data at  = 1=5 is T edge and is consistent with the observed absence of pair obtained via ED (L = 10, L = 15, and L = 20 with 5 10 , spreading reported in Figure S9 and also with the nite 4 3 10 , and 2  10 disorder realizations), Krylov time evolu- size scaling of mIPRs shown in the main text, Fig. 3. 3 3 tion (L = 25, T = 10 and 10 disorder realizations) and max TEBD (L = 40, T = 300 and 100 disorder realizations). max For  = 1=2 we used ED (L = 10, 12, 14, and 16 with 3 10 , 4 3 3 10 , 5 10 , and 10 disorder realizations) and Krylov time VI. DYNAMICAL PROBE OF THE ABSENCE 3 3 evolution (L = 18, T = 10 , and 10 disorder realizations). max OF RESONANCES The discussion on mutual IPR showed how tunneling processes are strongly suppressed in the dilute case of our have similar expansion over eigenstates, while very large model. In addition to eigenstates analysis, we also stud- values of mIPR imply that the expansion is very di er- ied long time dynamics of states with a thermal bubble. ent. In the main text we analyzed the mIPR between In this way, it was possible to verify whether a bubble two product states where the bubble is located at the initialized at a certain position can dynamically give rise left and right end of the system respectively, see Eqs (2)- to a dense region somewhere else in the chain. In order to (3). Such pair of states corresponds to the maximum study this process we de ned a projector onto the subset possible displacement of the bubble in the chain. Below of Hilbert space that has large density in a certain region. Δ S8 (a) (b) 15 15 10 10 13 W =0.5 W =6.5 13 W =0.5 W =6.5 W =2.5 W =8.5 W =2.5 W =8.5 11 11 W =4.5 W =10.5 W =4.5 W =10.5 9 9 − − 7 7 2 4 6 8 10 12 1 2 3 4 5 6 d d Figure S14. The mutual IPR, I , that quanti es the inverse probability of bubble tunneling d sites, increases exponentially with d at strong disorder. At weak disorder the mIPR approaches the Hilbert space dimension, N , shown by a dashed line. In the dilute system in (b), W = 4:5 marks the onset of the exponential growth, suggesting that the thermal bubble is frozen at its initial position. On the other hand, for  = 1=2, in (a), the clear exponential behavior emerges only at larger disorder. I was calculated for system sizes L = 15 and L = 12 in dilute and dense case respectively and averaged over 10 disorder realizations. More speci cally, we de ne size of the region L , we use the lengthscale extracted from the decay of n. Fit in Fig. 2(c) in the main text P (L ; i) = j ih j ; (S4) yields L ' 6 7, while t in Fig. S13(c) gives a some- c 0 what larger scale. We de ne an initial state j i that has j i2C an entangled dense region of size approximately L (de- where states ji are all possible product states that sat- scribed by a linear superposition of product states j i) isfy the condition    in the region [i; i + L ]. This c 0 followed by a product state: projector selects all con gurations where the system is lo- ^ N cally above the mobility edge. We notice that P (L ; i) X j i = p j i ji : (S5) takes into account all possible con gurations, thus con- 0 i i=1 sidering also the entropic factor. In order to understand what is the minimal required Below, we x the overall density to  = 1=4 and W = 6:5, which still corresponds to a localized system. The dense region is obtained as a superposition of di erent con gurations with N 1 particles in L = 2(N 1) L =12 sites. The remaining particle is initialized in the middle L =16 of the last segment of the chain. For instance, for L = 16 L =20 − 1 this results into following initial state: j i =   + − 2 (S6) 10   + + : : : ; 0 2 4 6 8 10 12 where the boxed area contains a dense entangled bubble and the remainder is in the dilute state. The initial state j i is then evolved through the Hamiltonian (1) in a quench protocol. After time evo- Figure S15. The late time evolution of hP (L ; d)i shows ex- n 0 lution up to a maximum time T = 1000, we measure max ponential decay for all the system sizes studied (L = 12; 16; 20 hP (L ; d)i = h (t)j P (L ; d)j (t)i, which quanti es 0  0 c c at density  = 1=4). Furthermore, we notice that increasing the probability of encountering a bubble shifted by d sites the system size the exponential vanishing becomes more se- from the initial position of the bubble. vere, suggesting that in the thermodynamic limit there would be no motion of the bubble at all. These results were obtained Finally, averaging over all di erent product states in 4 3 3 using 10 ; 5 10 and 10 disorder realizations for the system the dilute part of the chain and over disorder we ob- sizes from smaller to larger. tain the data in Fig. S15. This plot reveals that the hP i d S9 probability of having a dense ( >  ) region decays onantly but rather tunnels throughout the system. More- exponentially with the distance d from its initial loca- over, the nite size scaling analysis shows that increasing tion. This is in agreement with our long-time TEBD the system size the decay of hP (L ; d)i with distance d dynamics, Fig. S9, that reveals localization of individual is enhanced. Therefore in the dilute regime of our model pairs. Thus, we conclude that bubble does not spread res- the bubble remains localized around its initial position. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Quantum Physics arXiv (Cornell University)

Stability of mobility edges in disordered interacting systems

Brighi, Pietro; Abanin, Dmitry; Serbyn, Maksym
Quantum Physics , Volume 2020 (2005) – May 6, 2020
Download PDF
15 pages

Loading next page...
 
/lp/arxiv-cornell-university/stability-of-mobility-edges-in-disordered-interacting-systems-nguwYUb0XH

References (54)

ISSN
2469-9950
eISSN
ARCH-3342
DOI
10.1103/PhysRevB.102.060202
Publisher site
See Article on Publisher Site

Abstract

1 2 1 Pietro Brighi , Dmitry A. Abanin , and Maksym Serbyn IST Austria, Am Campus 1, 3400 Klosterneuburg, Austria and Department of Theoretical Physics, University of Geneva, 24 quai Ernest-Ansermet, 1211 Geneva, Switzerland (Dated: August 28, 2020) Many-body localization provides a mechanism to avoid thermalization in isolated interacting quantum systems. The breakdown of thermalization may be complete, when all eigenstates in the many-body spectrum become localized, or partial, when the so-called many-body mobility edge separates localized and delocalized parts of the spectrum. Previously, De Roeck et al. [Phys. Rev. B 93, 014203 (2016)] suggested a possible instability of the many-body mobility edge in energy den- sity. The local ergodic regions | so called \bubbles" | resonantly spread throughout the system, leading to delocalization. In order to study such instability mechanism, in this work we design a model featuring many-body mobility edge in particle density : the states at small particle density are localized, while increasing the density of particles leads to delocalization. Using numerical sim- ulations with matrix product states we demonstrate the stability of MBL with respect to small bubbles in large dilute systems for experimentally relevant timescales. In addition, we demonstrate that processes where the bubble spreads are favored over processes that lead to resonant tunneling, suggesting a possible mechanism behind the observed stability of many-body mobility edge. We conclude by proposing experiments to probe particle density mobility edge in Bose-Hubbard model. Introduction.|Many-body localization (MBL) pro- gested in Ref. [18], which equally applies to MBME in vides a mechanism to avoid thermalization in isolated any extensive conserved quantity. First, using numer- quantum interacting systems [1, 2]. Despite intensive ical simulation with matrix product states (MPS), we theoretical [3, 4] and experimental [5{11] studies, only demonstrate that uniform dilute states remain localized fully-MBL phase in one spatial dimension is relatively even at system sizes of L = 40 sites up to 250 tunneling well understood. The fate of MBL in higher dimen- times (i.e. more than two orders of magnitude larger than sions [13{17] and the possibility of the coexistence of the inverse local hopping). Next, we use a region with localized and delocalized eigenstates in the same many- large particle density to reproduce the bubble described body spectrum [18] remain debated. in [18] and track its in uence on the dilute remainder of the system in a quantum quench. We do not nd any Similarly to the case of Anderson localization [19], the evidence of resonant tunneling of the bubble, at least on MBL and delocalized eigenstates cannot coexist at the experimentally relevant timescales. same energy suggesting the existence of many-body mo- In summary, the study of the particle density MBME bility edge (MBME) | a certain energy in the spectrum facilitates the state preparation and analysis and allows separating localized and delocalized eigenstates [2]. In us to access the dynamics of much larger systems using contrast to the non-interacting case, the energy of MBME time evolution with MPS. We report the stability of the scales extensively with system size. In the absence of a particle density mobility edge on long timescales and sug- coupling to a bath, this leads to an exactly vanishing gest that similar physics may be experimentally probed conductivity (in contrast to an exponentially small but using Bose-Hubbard model. nite value in Anderson insulator) until a certain critical Correlated hopping model.|We consider hard-core temperature [2]. bosons on an open chain of size L, with the following Recently De Roeck et al. [18] suggested a possible Hamiltonian, mechanism that may destroy MBME in large systems: a nite region with local energy density above the mobility L1 L X X edge | a \bubble" | may resonantly spread through- H = t (c c + h:c:) +  n ^ 1 i i i i+1 out the system thereby destroying localization every- i=1 i=1 where. However, recent experiments [11] and MPS sim- L1 ulations [12] gave evidence of MBME, at least on inter- + t (c n ^ c + h:c:): (1) 2 i i+1 i1 mediate timescales. In addition, a number of numerical i=2 studies observed a mobility edge [20{24] using exact di- agonalization (ED). Unfortunately, the ED is limited to The rst line corresponds to the non-interacting Ander- relatively small system sizes; experiments with MBME son's model [25], where random on-site potential has a in energy density are also challenging since they require uniform distribution,  2 [W; W ]. The facilitated hop- energy resolution. ping in the second line enables motion of a pair of par- In order to overcome the above challenges, we pro- ticles with amplitude t , $, making the model pose to study MBME in particle density. This allows interacting. The Hamiltonian (1) has two channels for us to directly probe the mechanism of instability sug- dynamics: the single particle hopping prevails in dilute arXiv:2005.02999v2 [cond-mat.dis-nn] 26 Aug 2020 2 states, while the pair hopping is dominant at larger den- sities. We note that a similar model was discussed in Ref. [18] in two dimensions, although only with two particles. The enhancement of localization length in the case of two interacting particles also received signi cant atten- tion [26, 27]. In a di erent direction, the fate of the single particle mobility edge in the presence of interactions was studied [10, 28]. In contrast, we study model (1) that does not have a single particle mobility edge and con- sider the nite particle density regime. We x the value of the hopping parameters t = 0:5 Figure 1. Scaling of level spacing ratio demonstrates that and t = 2 so that the localization length of a single at density  = 1=5 (solid lines, L = 15; 20; 25 with 3; 4; 5 particle  . 1 and at the same time a single pair has SP particles) the system enters MBL phase for W  6:3. In a localization length  & 2:5 for 2:5 . W . 6 [29]. contrast, at half- lling  = 1=2 (dashed curves, L = 10; : : : 18) For such a choice, our model does not su er from nite the critical disorder strengths is much larger and in the entire size e ects [30] and we establish MBME using eigenstates range of disorder r approaches thermal value with increasing av probes. system size. Data is generated from ED/SI simulations with Eigenstate probes of localization.|We use exact di- at least 10 disorder realizations using approximately 2% of agonalization and shift-invert (SI) numerical techniques eigenstates in the center of the spectrum. to provide evidence for MBME in Hamiltonian (1). We analyze the average ratio of level spacings,  = E E , in the middle of the spectrum, r = the dense and dilute cases. While in the dilute case the i+1 i av hmin( ;  )= max( ;  )i. This is a commonly used system retains memory of the initial state, see Fig. 2(a), i i+1 i i+1 probe of the MBL transition [21, 31] that attains the at  = 1=2 quantum dynamics leads to a progressively value r ' 0:39 for the Poisson level statistics, charac- more uniform density pro le with increasing system size, teristic of the MBL phase, and r ' 0:53 for the case Fig. 2(b). In order to quantify the di erence in the form GOE of random Gaussian orthogonal ensemble (GOE), typical of the density pro le at late times, in Fig. 2(c) we plot for chaotic Hamiltonians with time-reversal symmetry. the average deviation of the density from the equilibrium Figure 1 displays that at half- lling,  = N=L = 1=2, thermal value , n = (1=L) jhn ^ (T )i j. The i max i=1 where N is the total number of particles and L is the deviation of late-time density from the thermal value, chain length, the level statistics approaches GOE with n, in the dense regime decays exponentially with the L= increasing system size, which is consistent with the de- system size as n  e , where  ' 6:27. In con- localized phase. In contrast, at  = 1=5 lling r ows trast, for the dilute case n shows no dependence on the av towards r at strong disorder. In what follows we x system size, as is apparent in the density pro les. The the disorder strength to be W = 6:5, since at this value characteristic length  extracted in the dense case gives the dilute limit is localized while the dense limit clearly the minimum size for genuine egodic bubbles that can ows towards delocalization. The scaling of entangle- destroy the MBME according to Ref. [18]. ment entropy also con rms the coexistence of localized Having con rmed the coexistence of localized and de- and delocalized phases for disorder W = 6:5 [29]. localized states at di erent values of particle density Quench dynamics.|Having provided numerical evi- for the same disorder strength, we proceed with a more dence for the coexistence of localized and delocalized detailed study of the e ect of a bubble, whose behavior is phases in small systems, we turn to quantum quench dy- central to the mechanism proposed in [18]. Figure 2(d) namics that distinguishes MBL from ergodic phase [5, illustrates the evolution of a non-uniform initial state, 32]. We consider quenches where the system is initially where a dense region represents the bubble. The bubble prepared in a product state and then evolved with the region consists of 8 sites with two pairs of particles and Hamiltonian (1). Starting with a density wave of pe- has a local density of  = 1=2. The bubble is followed riod 1=, a con guration that contains no pairs, we cal- by a period-5 density wave that occupies L 10 sites culate the density pro le at late times. For the dilute and two additional empty sites at the end of the chain. case,  = 1=5, we use the time-evolved block decimation Although having  = 8=30 > 1=5, this state is still in (TEBD) with MPS [33, 34] (see [29] for additional details a localized sector, as shown in [29]. The bubble leaks and benchmarks). This allows to monitor dynamics of only weakly into the dilute region even at late times, see systems as large as L = 40 sites up to times T = 500. Fig. 2(d), with particles far away from the bubble not be- max In the dense case ( = 1=2) we use ED and Krylov sub- ing a ected. In contrast, in the dense case, Fig. 2(g), the space time evolution method. While ED allows to access bubble with average density of  = 2=3 successfully melts the in nite-time density pro le, with the Krylov method, the period-3 density wave state throughout the system. we simulate quantum dynamics up to T = 1000. max Next, in panels Fig. 2(e) and (h) we further illustrate The density pro les at late times look very di erent in the di erences between the density dynamics in the dense 3 •◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦••◦••◦◦◦•◦◦•◦◦• ••◦◦••◦◦◦◦•◦◦◦◦•◦◦ (a) (d) 1 . 0 (g) L =40 L =25 L =15 1. 00 1. 00 L =30 L =20 L =30 L =20 L =10 0 . 8 (a) 0. 75 0. 75 0 . 6 i i i n n 0. 50 0. 50 n 0 . 4 0. 25 0. 25 0 . 2 0. 00 0. 00 0 . 0 5 10 15 20 25 30 35 40 1 4 7 10 13 16 19 22 25 28 31 1 3 5 7 9 11 13 15 17 i i ••◦••◦••◦◦◦•◦◦•◦◦• •◦•◦•◦•◦•◦•◦•◦•◦•◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦ ◦◦ (b) (e) (h) 0. 6 L =10 L =14 L =18 1. 00 0. 7 L =12 L =16 0. 5 0. 75 0. 6 0. 4 n ˜ 0. 50 0. 5 0. 3 0. 25 0. 4 0. 2 0. 00 −1 0 1 2 −1 0 1 2 3 1 3 5 7 9 11 13 15 17 10 10 10 10 10 10 10 10 10 ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦ ◦ ••◦••◦••◦◦◦•◦◦•◦◦• ••◦◦••◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦◦◦•◦◦ ◦◦ (c) (f) (j) 0.4 2. 0 c =4−5 c =3 c =8−9 c =6 1. 5 0.2 c =12−13 c =9 1. 0 c =18−19 c =12 ν = 1 / 5 0.1 0. 5 c =24−25 0.08 c =15 ν = 1 / 2 0.06 0. 0 0 −1 0 1 2 −1 0 1 2 3 10 15 20 25 30 35 40 10 10 10 10 10 10 10 10 10 L t Figure 2. (a)-(c) The quantum quench from the uniform density wave with period 1= reveals memory of the initial state at  = 1=5 in (a), whereas in the dense case  = 1=2 (b) the charge pattern relaxes to zero exponentially in the system size as is shown in (c). (d)-(f) Stability of the dilute system against the bubble consisting of a half- lled region with 4 particles is illustrated in panel (d) by the density pro le at T . (e) The time dynamics of density in the coarse grained regions (see max the legend at the top) shows the absence of signi cant relaxation in regions away from the bubble. The density in the region at the boundary with the bubble increases logarithmically in time. (f ) The onset of logarithmic entanglement dynamics after a transient is visible for all cuts (see the legend at the top) away from the bubble. (g)-(j) In contrast, the bubble delocalizes the system when the overall density  = 1=2. The residual density pro le at T = 1000 in panel (g) has only weak memory max of the initial state. In addition, densities coarse grained over 3-site regions in (f ) all tend to the equilibrium value of 1/2 and entanglement entropy in (j) displays faster than logarithmic growth for all cuts. The data are generated using TEBD and 4 4 3 Krylov (ED) dynamics using between 5 10 {100 (dilute) and 3 10 {10 (dense) disorder realizations. and dilute cases in presence of a bubble. In both cases change of density [see Fig. 2(e)], and after saturation of we plot the density of particles within subregions of small density dynamics, we expect an onset of the logarithmic i+k1 growth of entanglement. In contrast, the entanglement size k,  ~ = (1=k) hn i, that are shown at the top i j j=i dynamics in Fig. 2(j) is always faster than logarithmic. of the plot. In the dilute case, Fig. 2(e), we observe that In [29] we provide more details on the contribution of ~ remains far from its thermal value even at late times, particle transport to entanglement [38, 39], demonstrat- in contrast with [18], where an ergodic region larger than ing that it is responsible for logarithmic entanglement is expected to delocalize the system. The densities of increase, in agreement with [40], whereas the con gura- regions in the bubble and adjacent to the bubble seem tional entanglement grows faster than logarithmic, and to saturate, while the regions far away from the bubble total entropy shows power-law increase. show very slow dynamics. In contrast, the dense case, Bubble tunneling vs. decay processes.|The quench dy- Fig. 2(h), shows that all expectation values evolve to- namics discussed above suggests that a bubble is not able wards equilibrium, although the regions far away from to spread through the entire localized chain and remains the center of the chain display slow, logarithmic in time, in the vicinity of its initial position. At the same time, growth of density. most of our quench simulations are restricted to nite, Finally, we study the dynamics of the bipartite en- albeit long, times. In order to give a complementary evi- tanglement entropy, S , see Fig. 2(f ) and (j). The en- vN dence for the bubble localization, we return to eigenstate tanglement is de ned as S = tr  ln , where  is vN properties that e ectively probe the in nite time limit. the density matrix of the left subregion calculated from We start with an initial product state in the half- lled iHt j (t)i = e j i. Di erent entanglement cuts shown case illustrated for L = 12, at the top of Fig. 2(f ) and (j) are encoded by their color. Consistent with MBL, the increase of entanglement in j i =  ; (2) the region close to the bubble is logarithmic in time in Fig. 2(f ) [35{38]. The entanglement across the cuts fur- that contains a bubble of k = L=2 sites with  = 2=3 ll- ther away from the bubble begins to grow at signi cantly ing (boxed region), followed by a sparser region with the later times. For these more distant cuts, the initial up- same number of sites and density  = 1=3. To quantify rise in entanglement corresponds to a slow logarithmic the relation between the probability of the bubble tunnel- ( ) S(t) 4 ing to the opposite end of the system and the probability where the bubble spreads throughout the system, calling of the bubble spreading throughout the system, we use into question the applicability of the resonance argument a spatial re ection of j i and uniform density wave as of [18]. a representative of the state with bubble tunneling and Experimental realization.| Finally, we discuss a pos- spreading, respectively: sible way to observe the physics related to MBME in ex- periments with ultracold atoms. Within the disordered =   ; (3) Aubry-Andr e bosonic Hamiltonian, j i = ; (4) h i 2 X H = t(a a + h.c.) +  n + Un (n 1) ; (6) i+1 i i; i i illustrated for L = 12 and  = 1=2 lling. For dilute con gurations at  = 1=5 we de ne the bubble as a region that is actively used to study MBL physics [38, 42], the of size 2(L1) with density  = 1=2, joined with a dilute remainder. For L = 20 such a state is: bubbles can be represented by regions with ha a i =  > j i =  . 1 bosons per site. A particle within such region has a It is straightforward to show that the in nite time av- hopping matrix element enhanced by the Bose-factor of erage probability of nding the system with the wave hi, thus playing the role of hopping t in model (1). iHt In the regime of densities and disorder strengths such function e j i in the product state j i is given by 1 2 that the enhanced hopping hit corresponds to local- X ization lengths signi cantly larger than lattice spacing, P (j i ;j i) = jhE j ihE j ij ; (5) 1 2 1 2  > a, whereas a single boson localization length is dense =1 . a, this model will implement similar physics to our toy model. Note that at the same time it is important where jE i are the complete set of eigenstates of H . to keep interaction U low enough, U  t, to avoid the Eq. (5) quanti es the similarity in the expansion of two formation of minibands related to long-lived doublons. di erent states j i over the basis jE i and reduces to 1;2 By initializing the system in a product state with a the conventional participation ratio when j i = j i. 1 2 dense region of bosons in the center of the trap along In order to reveal the relation between bubble decay with low density of bosons away from such a region, the and tunneling processes, we calculate the ratio of proba- dynamics under Hamiltonian (6) will probe the ability of s s bilities of bubble decaying, P = P (j i ;j i), with j i s 1 2 2 the bubble to melt the imbalance [5] away from its origi- from Eq. (4), to bubble tunneling, P = P (j i ;j i) t 1 nal position. From our simulations we expect the absence with j i from Eq. (3). In the dense case, these two of imbalance relaxation far away from the bubble. In a probabilities are of the same order and moreover tend to di erent direction, doublons [43, 44] or second species of identity with increasing system size as expected in the particles not subject to disorder [9] are also promising delocalized phase, see Fig. 3. In the dilute case, the ra- candidates that can play a role of the bubble. tio P =P is rapidly increasing with both disorder, and s t Discussion.|We presented a model with MBME in system size. This proves that the bubble tunneling pro- particle density and investigated its properties numeri- cesses are strongly suppressed compared to the processes cally using ED and time evolution with MPS. We nd strong evidence of the persistence of localization at in- nite times for small systems and also observe memory of initial con guration until times of T = 500 for sys- 10 max L =12 L =10 tems with up to L = 40 sites. These times are at least L =14 L =15 two orders of magnitude larger compared to the inverse L =16 L =20 10 local hopping, ~=t , and are achievable with cold atoms experiments. While we cannot rule out a residual very slow delocalization at much later times, the constructed model allows us to bound the timescale up to which the localization remains stable in very large systems that are − 1 0 1 2 3 4 5 6 7 8 beyond the reach for ED. The model with MBME in particle density presented in this work allows for direct tests of the arguments about the instability of MBME [18]. In order for bubble to move Figure 3. The rapid increase of the ratio of P =P with s t throughout the system it is important that the bubble system size and disorder strength reveals that in the dilute does not disappear by spreading and that con gurations case,  = 1=5 the probability for bubble spreading is strongly with bubbles situated at di erent locations are e ectively enhanced compared to the probability of bubble tunneling coupled to each other. Our simulations reveal that dilute to the opposite end of the system. For the dense case these systems have no trace of bubble reemerging at a di erent two probabilities are of the same order and approach each location within the system. Moreover, even the expecta- other with increasing system size in a broad range of disorders. tion value of the pair density hn n i (pairs are building Averaging is done over at least 2:5 10 disorder realizations. i i+1 blocks of the bubble) is exponentially suppressed away P /P t s 5 from the original location of the dense bubble [29]. In an supported by the European Unions Horizon 2020 research alternative approach, we directly test the probability of and innovation programme under the Marie Sklodowska- the bubble to emerge at the opposite end of the system Curie Grant Agreement No. 665385. D.A. was supported at in nite time and nd it to be strongly suppressed. by the Swiss National Science Foundation. M.S. was To conclude, we expect that the proposed model will supported by European Research Council (ERC) un- enable further investigations of particle density MBME. der the European Union's Horizon 2020 research and Studies of the structure of matrix elements, extension of innovation programme (grant agreement No. 850899) the theory of LIOMs [45, 46] to systems with MBME in This work bene tted from visits to KITP, supported particle density [22], and studies of the e ect of a small by the National Science Foundation under Grant No. bath on a localized system [47{50] using our model rep- NSF PHY-1748958 and from the program \Thermal- resent promising avenues for future work. ization, Many body localization and Hydrodynamics" Acknowledgments.|We acknowledge useful discus- at International Centre for Theoretical Sciences (Code: sions with W. De Roeck and A. Michailidis. P.B. was ICTS/hydrodynamics2019/11). [1] L. Fleishman and P. W. Anderson, \Interactions and the uan Song, Hekang Li, Zhen Wang, Wenhui Ren, anderson transition," Phys. Rev. B 21, 2366{2377 (1980). Hang Dong, Dongning Zheng, Yu-Ran Zhang, Rubem [2] D.M. Basko, I.L. Aleiner, and B.L. Altshuler, \Metal{ Mondaini, Heng Fan, and H. Wang, \Observation of en- insulator transition in a weakly interacting many-electron ergy resolved many-body localization," arXiv e-prints , system with localized single-particle states," Annals of arXiv:1912.02818 (2019), arXiv:1912.02818 [quant-ph]. Physics 321, 1126 { 1205 (2006). [12] Titas Chanda, Piotr Sierant, and Jakub Zakrzewski, [3] Dmitry A. Abanin, Ehud Altman, Immanuel Bloch, and \Many-body localization transition in large quantum Maksym Serbyn, \Colloquium: Many-body localization, spin chains: The mobility edge," arXiv e-prints thermalization, and entanglement," Rev. Mod. Phys. 91, , arXiv:2006.02860 (2020), arXiv:2006.02860 [cond- 021001 (2019). mat.dis-nn]. [4] Rahul Nandkishore and David A. Huse, \Many-body lo- [13] Wojciech De Roeck and Francois Huveneers, \Stability calization and thermalization in quantum statistical me- and instability towards delocalization in many-body lo- chanics," Annual Review of Condensed Matter Physics calization systems," Phys. Rev. B 95, 155129 (2017). 6, 15{38 (2015). [14] Thorsten B. Wahl, Arijeet Pal, and Steven H. Simon, [5] Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, \Signatures of the many-body localized regime in two Henrik P. Lusc  hen, Mark H. Fischer, Ronen Vosk, Ehud dimensions," Nature Physics 15, 164{169 (2019). Altman, Ulrich Schneider, and Immanuel Bloch, \Obser- [15] Ionut-Dragos Potirniche, Sumilan Banerjee, and Ehud vation of many-body localization of interacting fermions Altman, \Exploration of the stability of many-body lo- in a quasirandom optical lattice," Science 349, 842{845 calization in d > 1," Phys. Rev. B 99, 205149 (2019). (2015). [16] Augustine Kshetrimayum, Marcel Goihl, and Jens [6] Pranjal Bordia, Henrik P. Lusc  hen, Sean S. Hodgman, Eisert, \Time evolution of many-body localized sys- Michael Schreiber, Immanuel Bloch, and Ulrich Schnei- tems in two spatial dimensions," arXiv e-prints, der, \Coupling identical one-dimensional many-body lo- arXiv:1910.11359 (2019), arXiv:1910.11359 [cond- calized systems," Phys. Rev. Lett. 116, 140401 (2016). mat.str-el]. [7] Jae-yoon Choi, Sebastian Hild, Johannes Zeiher, Peter [17] Elmer V. H. Doggen, Igor V. Gornyi, Alexander D. Schau, Antonio Rubio-Abadal, Tarik Yefsah, Vedika Mirlin, and Dmitry G. Polyakov, \Slow many- Khemani, David A. Huse, Immanuel Bloch, and body delocalization beyond one dimension," arXiv e- Christian Gross, \Exploring the many-body localization prints, arXiv:2002.07635 (2020), arXiv:2002.07635 [cond- transition in two dimensions," Science 352, 1547{1552 mat.dis-nn]. (2016). [18] Wojciech De Roeck, Francois Huveneers, Markus Muller,  [8] B. Chiaro et al., \Growth and preservation of entan- and Mauro Schiulaz, \Absence of many-body mo- glement in a many-body localized system," arXiv e- bility edges," Physical Review B 93, 1{14 (2016), prints, arXiv:1910.06024 (2019), arXiv:1910.06024 [cond- arXiv:1506.01505. mat.dis-nn]. [19] N.F. Mott, \Electrons in disordered structures," Ad- [9] Antonio Rubio-Abadal, Jae Yoon Choi, Johannes Zeiher, vances in Physics 16, 49{144 (1967). Simon Hollerith, Jun Rui, Immanuel Bloch, and Chris- [20] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Cri- tian Gross, \Many-Body Delocalization in the Presence terion for many-body localization-delocalization phase of a Quantum Bath," Physical Review X 9, 41014 (2019), transition," Phys. Rev. X 5, 041047 (2015). arXiv:1805.00056. [21] David J. Luitz, Nicolas La orencie, and Fabien Alet, [10] Thomas Kohlert, Sebastian Scherg, Xiao Li, Henrik P. \Many-body localization edge in the random- eld heisen- Lusc  hen, Sankar Das Sarma, Immanuel Bloch, and berg chain," Phys. Rev. B 91, 081103 (2015). Monika Aidelsburger, \Observation of many-body local- [22] Scott D. Geraedts, R. N. Bhatt, and Rahul Nandkishore, ization in a one-dimensional system with a single-particle \Emergent local integrals of motion without a complete mobility edge," Phys. Rev. Lett. 122, 170403 (2019). set of localized eigenstates," Phys. Rev. B 95, 064204 [11] Qiujiang Guo, Chen Cheng, Zheng-Hang Sun, Zix- (2017). 6 [23] Marcel Goihl, Jens Eisert, and Christian Krumnow, \Ex- localization transition," Nature 573, 385{389 (2019). ploration of the stability of many-body localized systems [43] Ulrich Krause, Th eo Pellegrin, Piet W. Brouwer, in the presence of a small bath," Phys. Rev. B 99, 195145 Dmitry A. Abanin, and Michele Filippone \Nucle- (2019). ation of ergodicity by a single doublon in supercooled [24] Ruixiao Yao, and Jakub Zakrzewski, \Many-body lo- insulators," arXiv e-prints , arXiv:1911.11711 (2019), calization in Bose-Hubbard model: evidence for the mo- arXiv:1911.11711 [cond-mat.dis-nn]. bility edge ," arXiv e-prints , arXiv:2002.00381 (2020), [44] Micheal Schecter, Thomas Iadecola, and Sankar Das arXiv:2002.00381 [condmat.dis-nn]. Sarma, \Con guration-controlled many-body localiza- [25] P. W. Anderson, \Absence of di usion in certain random tion and the mobility emulsion," Phys. Rev. B 98, lattices," Phys. Rev. 109, 1492{1505 (1958). 174201 (2018). [26] D. L. Shepelyansky, \Coherent propagation of two inter- [45] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Lo- acting particles in a random potential," Phys. Rev. Lett. cal conservation laws and the structure of the many-body 73, 2607{2610 (1994). localized states," Phys. Rev. Lett. 111, 127201 (2013). [27] Y Imry, \Coherent propagation of two interacting parti- [46] David A. Huse, Rahul Nandkishore, and Vadim cles in a random potential," Europhysics Letters (EPL) Oganesyan, \Phenomenology of fully many-body- 30, 405{408 (1995). localized systems," Phys. Rev. B 90, 174202 (2014). [28] Xiaopeng Li, Sriram Ganeshan, J. H. Pixley, and [47] Rahul Nandkishore, \Many-body localization proximity S. Das Sarma, \Many-body localization and quantum e ect," Phys. Rev. B 92, 245141 (2015). nonergodicity in a model with a single-particle mobility [48] David J. Luitz, Francois Huveneers, and Wojciech edge," Phys. Rev. Lett. 115, 186601 (2015). De Roeck, \How a small quantum bath can thermal- [29] \Supplemental Online material,". ize long localized chains," Phys. Rev. Lett. 119, 150602 [30] Z. Papi c, E. Miles Stoudenmire, and Dmitry A. Abanin, (2017). \Many-body localization in disorder-free systems: The [49] Katharine Hyatt, James R. Garrison, Andrew C. Potter, importance of nite-size constraints," Annals of Physics and Bela Bauer, \Many-body localization in the presence 362, 714 { 725 (2015). of a small bath," Phys. Rev. B 95, 035132 (2017). [31] Arijeet Pal and David A. Huse, \Many-body localization [50] Philip J. D. Crowley and Anushya Chandran, \Avalanche phase transition," Phys. Rev. B 82, 174411 (2010). induced co-existing localised and thermal regions in [32] Maksym Serbyn, Z. Papi c, and D. A. Abanin, \Quantum disordered chains," arXiv e-prints , arXiv:1910.10812 quenches in the many-body localized phase," Phys. Rev. (2019), arXiv:1910.10812 [cond-mat.dis-nn]. B 90, 174302 (2014). [51] Elmer V. H. Doggen, Frank Schindler, Konstantin S. [33] Guifr e Vidal, \Ecient classical simulation of slightly en- Tikhonov, Alexander D. Mirlin, Titus Neupert, tangled quantum computations," Phys. Rev. Lett. 91, Dmitry G. Polyakov, and Igor V. Gornyi, \Many-body 147902 (2003). localization and delocalization in large quantum chains," [34] ITensor Library (version 3.1.1) http://itensor.org. Phys. Rev. B 98, 174202 (2018). [35] M. Znidaric, T. Prosen, and P. Prelovsek, \Many-body localization in the heisenberg xxz magnet in a random eld," Phys. Rev. B 77, 064426 (2008). [36] Jens H. Bardarson, Frank Pollmann, and Joel E. Moore, \Unbounded growth of entanglement in models of many- body localization," Phys. Rev. Lett. 109, 017202 (2012). [37] Maksym Serbyn, Z. Papi c, and Dmitry A. Abanin, \Universal slow growth of entanglement in interacting strongly disordered systems," Phys. Rev. Lett. 110, 260601 (2013). [38] Alexander Lukin, Matthew Rispoli, Robert Schittko, M. Eric Tai, Adam M. Kaufman, Soonwon Choi, Vedika Khemani, Julian L eonard, and Markus Greiner, \Prob- ing entanglement in a many-body{localized system," Sci- ence 364, 256{260 (2019). [39] Rajibul Islam, Ruichao Ma, Philipp M Preiss, M Eric Tai, Alexander Lukin, Matthew Rispoli, and Markus Greiner, \Measuring entanglement entropy in a quantum many-body system," Nature 528, 77 (2015). [40] Maximilian Kiefer-Emmanouilidis, Razmik Unanyan, Michael Fleischauer, and Jesko Sirker, \Evidence for unbounded growth of the number entropy in Many-Body Localized phases," Phys. Rev. Lett. 124, 243601 (2020). [41] David J. Luitz, and Yevgeny Bar Lev \Is there slow particle transport in the MBL phase?," arXiv e-prints , arXiv:2007.13767 (2020), arXiv:2007.13767 [cond-mat.dis-nn]. [42] Matthew Rispoli, Alexander Lukin, Robert Schittko, Sooshin Kim, M. Eric Tai, Julian L eonard, and Markus Greiner, \Quantum critical behaviour at the many-body S1 Supplementary material for \Particle density mobility edge" In this supplementary material we present additional data and details of the methods used in the main text. First, we show data for additional probes of ETH breakdown such as entanglement entropy. After this we present benchmarks of our time-evolving block decimation simulation of dynamics. Finally, we explore the behavior of the mutual IPR de ned in the main text. I. LOCALIZATION LENGTH AND ized, while at lower densities   1=5, when the typical PARAMETER CHOICE distance between pairs is large, we expect MBL phase. This motivates the choice t = 2, since at this value of t 2 2 (W ) approaches 2 at disorder strength around W  6. The dynamics generated by the constrained Hamilto- We note that we avoided further increase of t to keep the nian, Eq. (1), strongly depends on the choice of the hop- model away from the constrained limit: in the case when ping parameters t . In order to choose the most suitable 1;2 t dominates over t , the model would approximately re- 2 1 parameters for the study of MBME in particle density, we duce to a kinematically constrained model that has many explore localization lengths for a single particle  and SP disconnected sectors in the Hilbert space. for one pair of particles  . These localization lengths In order to rule out the presence of strong nite size are evaluated using ED. We calculate the in nite-time e ects, we studied the density of state in individual dis- average of the occupation number at each site for an ini- order realizations. In the regime when t  t the strong 2 1 tial state where either a single particle or a single pair nite size e ects would give rise to the presence of the are initialized at the rst site of the chain. We extract mini-bands and the DOS would become non-monotonous the localization lengths  ( ) from an exponential t SP P with numerous peaks corresponding to mini-band struc- of the density curve hn i . ture [30]. Figure S2 con rms that at our choice of pa- Resulting values of  for xed t = 0:5 and dif- P,SP 1 rameters even individual disorder realizations have a rel- ferent disorder values and di erent values of hopping t atively smooth density of states with Gaussian envelope, are shown in Fig. S1. The single particle hopping local- thus ruling out the presence of strong nite size e ects. ization length (dashed line in Fig. S1) does not depend on t , and becomes smaller than one lattice spacing for W & 4. The pair localization length is monotonously II. ED PROBES OF LOCALIZATION increasing with t at xed value of disorder strength, W . Our aim is to have  in the range between 2 and 5. In While in the main text we focused on the two values of this regime, the half- lling case is expected to be delocal- lling,  = 1=2 and 1=5, here we demonstrate the density dependence of critical disorder. For this purpose we cal- culate the average ratio of level spacings r for a single av system size L = 18 at varying values of density. Fig- ξ ξ , t =2 SP p 2 ξ , t =1 ξ , t =2.5 p 2 p 2 ξ , t =1.5 ξ , t =3 p 2 p 2 0 . 08 0 . 06 0 . 04 0 . 02 1 2 3 4 5 6 0 . 00 −20 −15 −10 −5 0 5 10 15 20 Figure S1. The localization length decreases, as expected, E with the disorder strength for all the values of t . For every constrained hopping amplitude t it is possible to locate the region of disorder where we expect to see a MBME in parti- Figure S2. The DOS from single disorder realizations show a cle density as the area among the two dashed lines. As the relatively smooth behavior and a Gaussian shape, thus con- curve crosses the rst dashed line, systems with typical par- rming the absence of strong nite size e ects. DOS refers ticle spacing 5 will be localized. Nevertheless denser states to a chain with L = 20 and  = 1=4. Disorder strength is will still be delocalized, having smaller distance among par- W = 5:0. Green, blue and orange curves correspond to dif- ticles. Data were obtained on a lattice of length L = 50 and ferent disorder realizations, while the black dashed line shows averaged over 5000 disorder realizations. disorder-averaged DOS. S2 0.53 L = 10 L = 12 0.30 0.8 0.51 L = 14 0.7 0.25 0.49 L = 16 0.6 0.20 0.47 ν 0.5 av 0.15 0.45 0.4 0.43 0.10 0.3 0.41 L = 10 0.2 0.05 L = 15 L = 20 0.39 0.00 2 4 6 8 1 2 3 4 5 6 7 8 Figure S3. The sharp di erence of r obtained for di erent av Figure S4. The behavior of half-chain entanglement entropy at the same disorder W clearly shows the MBME in our shows very distinct behavior for dilute (blue-shaded curves) model. Interestingly, the mobility edge curve W () is not and dense (red-shaded curves) states. The crossing in the symmetric, but is peaked around  = 2=3, implying that the dilute states implies that they entered the MBL phase, and states with the maximum number of pairs for xed size are thus have area-law entanglement entropy. On the other hand, the hardest to localize. The data is obtained for a system of dense states do not show a similar crossing in this range of size L = 18, using shift-invert method with 10 10 states disorder, suggesting that they are still in the ergodic phase. 4 3 from the middle of the spectrum and 5  10 10 disorder The data are obtained with shift-invert method for 10 10 realizations. eigenstates in the middle of the spectrum and averaged over 4 3 5 10 5 10 disorder realizations. ure S3 allows to estimate the dependence of the critical disorder on the lling, . At low densities ( <  (W )) sector of the Hamiltonian Eq. (1). Given the U (1) sym- states have r approaching value characteristic for Pois- av metry, we would expect that a xed lling sector presents son distribution of level spacings. In contrast, for dense MBME in energy density, similarly to the case of the ran- con gurations ( >  (W )) the level spacing ratio is close dom eld XXZ spin chain. In the middle of the spectrum to GOE prediction. the density of states is large and eigenstates may remain Figure S3 reveals that the most delocalized lling is delocalized, while at the same disorder strength the states = 2=3, which corresponds to the case when the best at the edges of the many-body spectrum are localized. packing of pairs in the chain,  , can be achieved. To explore the eventual presence of MBME in energy At this lling the Poisson values of r would be achieved av density in the half lling sector, we studied the energy beyond the upper limit of the considered disorder range. resolved level spacing ratio r (; W ). The results, dis- Av Decreasing particle density away from this value causes earlier onset of localization. For instance, xing disorder value W = 6:5 we observe that  = 1=2 and  = 1=5 are 1.0 situated well in delocalized and localized regions. 0.52 L =10 In addition to the level statistics indicator presented in 0.8 0.50 L =12 the main text, we studied other commonly used probes of L =14 0.48 0.6 ergodicity. In particular, Fig. S4 illustrates the behavior L =16 0.46 of bipartite entanglement entropy for di erent disorder 0.4 0.44 strengths and di erent llings. On the one hand, the - nite size scaling of entanglement entropy of eigenstates in 0.42 0.2 the middle of the spectrum shows that for  = 1=5 and 0.40 disorder W > W  6 the entanglement is consistent 0.0 3 4 5 6 with area-law. On the other hand, the entanglement of dense systems,  = 1=2, does not show a similar behav- ior. The nite size scaling, indeed, shows no crossing Figure S5. The plot, presenting the energy resolved level at these disorder values, thus suggesting volume-law of spacing ratio for L = 16 in the half lling sector, shows clear evidence of MBME. In the center of the band r approaches entanglement entropy for  = 1=2. av the GOE value, while at large and small energy density it is close to the Poisson value. Furthermore, the MBME curves obtained for smaller systems seem to converge at increasing A. Energy Density MBME system size, thus suggesting the persistence of the MBME in the thermodynamic limit. The plot was obtained averaging 4 3 2 2 over n = 10 ; 2 10 ; 5 10 ; 10 disorder realizations for In this section we brie y discuss the presence of many- increasing system size from L = 10 to L = 16. body mobility edge in energy density in a single density S/L Av S3 (a) (c) 10 − 12 δt = 0.1 ε = 10 − 9 0.005 δt = 0.1 ε = 10 − 3 0.000 − 6 − 12 δt = 0.05 ε = 10 − 9 − 9 − 0.005 10 δt = 0.05 ε = 10 t =1 t =200 − 9 δt = 0.02 ε = 10 t =100 t =300 − 12 − 0.010 0.02 0.03 0.05 0.08 0.1 1 4 7 10 13 16 19 δt i (b) (d) 0.02 − 12 δt =0.02 δt =0.1 δt = 0.1 ε = 10 − 2 − 9 0.01 δt =0.05 10 δt = 0.1 ε = 10 δ − 5 0.00 − 12 − 8 − 0.01 δt = 0.05 ε = 10 − 9 δt = 0.05 ε = 10 − 11 − 0.02 10 − 9 δt = 0.02 ε = 10 0 1 2 10 10 10 1 4 7 10 13 16 19 t i Figure S6. (a-b) Deviation of ground state delity from 1 in Trotter time evolution, F (t; t), shows power-law behavior both in time and in time-step, as expected. The data is obtained at density  = 1=5, system size L = 20 and disorder W = 6:5 for a particular disorder realization. The plots for other disorder realizations are qualitatively similar. (c-d) The comparison between ED and TEBD time evolution reveals that the most e ective way to increase accuracy of TEBD is to decrease the time-step t. 9 12 Indeed, the change in the truncation between " = 10 and " = 10 does not have much e ect on the the di erence between density pro les of exact diagonalization and TEBD. At the same time, the decrease of time step brings the local density pro le closer to ED results. The density pro les are calculated by propagating uniform density wave (c) and uniform pair-density wave (d) initial states to time t = 500 for a particular disorder realization with L = 20,  = 1=5, W = 6:5. played in gure S5, show evidence of many-body mobility related to the nite size of the time step in the p-th order edge; the level spacing ratio, as a function of the energy Trotter expansion grows as t . The other source of error EE min density  = , where E and E are the is the nite cuto , ", that governs the truncation of sin- min max E E max min ground state and the most excited state respectively, in- gular values in the singular value decomposition (SVD). creases from the Poisson to the GOE value as  goes While the instantaneous errors related to the trunca- from the lower edge to the center of the spectrum and tion and nite time step are known, understanding the decreases again from the center to the upper edge. This propagation of these errors with time and their possi- variation is such that at a xed disorder strength the ble interference is challenging. First we tested TEBD low and high energy states are localized, while the cen- algorithm by evolving the ground state of the same ter of the band is delocalized, thus de ning a many-body model. Provided that the time evolution is numerically mobility edge. The scaling of the MBME curves for dif- exact, the overlap between the TEBD-evolved ground TEBD ferent system size shows signs of convergence, suggesting state, j (t)i = U (t)jGSi and the exact time evo- {E t stability of the MBME in the thermodynamic limit. lution of the ground state, jGS(t)i = e jGSi, is supposed to give the identity h (t)jGS(t)i = 1 at all times. For the fourth-order Trotterization the behavior III. MPS SIMULATIONS OF QUENCH of F = 1jh (t)jGS(t)ij is known to be proportional DYNAMICS to (t) . The numerical results plotted in Fig. S6(a), con rm these expectations. In our MPS simulation, we time evolve dilute states in Next, we performed a benchmarking of TEBD algo- large systems L  30 up to time T = 500. For this rithm against ED time evolution for several disorder re- max we use the time-evolving block decimation (TEBD) algo- alizations and simulation parameters. An illustration of rithm with a fourth-order Trotter evolution based on the such benchmarking is shown in Fig. S6(c) and (d). In ITensor library [34]. The main parameter involved in the particular, we observed that time step t = 0:05 and cut- time evolution algorithm is the time step t used to split o " = 10 result in a good agreement between ED and the unitary evolution into a sequence of gates. The error TEBD dynamics. Smaller values of t; " would improve F (t, t ) F (t, t ) ED TN ED TN hn i h n i i i hn i h n i i i S4 IV. ADDITIONAL RESULTS FROM QUENCH • • ◦ ◦ • • ◦ ◦ ◦ ◦ • ◦ ◦ . . . 0.0005 DYNAMICS • ◦ ◦ ◦ ◦ • ◦ ◦ ◦ ◦ • ◦ ◦ . . . 0.0004 In the main text we discussed dynamics in quenches that begin from a uniform state or a bubble joined to a 0.0003 more dilute remainder. Below we present details for the 0.0002 quenches in presence of bubble. In addition, we discuss the dynamics resulting from the initial state containing 0.0001 density wave of particle pairs. 0.0000 0 100 200 300 400 500 A. Pair density and entanglement dynamics in presence of a bubble Figure S7. The normalized absolute value of the energy di er- Since particle pairs are the most mobile objects, we ence from the initial energy E(t) = jhH (t)i E(0)j=jE(0)j consider the pair density in quenches that are initialized remains very small for both the density wave, blue curve, with the bubble (see Fig. 2(d),(g) in the main text). The and the non-uniform, red curve, con gurations, con rming pair density is of special interest in these quenches as in the good accuracy of our numerical simulations beyond the Ref. [18] suggested that the instability of the system is ED benchmark. The larger deviation displayed by the non- ascribed to the ability of the bubble to move. In our uniform con guration is understood as a result of the pres- ence of the bubble in the lattice, that increases entanglement model the bubble consists of several pairs, thus motion growth. The results here shown are obtained averaging over of the bubble throughout the system would imply the 100 disorder realizations, for the system sizes, L = 30, and spreading of pairs. initial states described in the main text. The pair density de ned as hn n i measured at late i i+1 or in nite times is shown in Fig. S9. In the dense case the late time pair density pro le supports delocalization: the agreement but would result in a dramatic slowdown at late times the density of pairs becomes homogeneous of the evolution time. Therefore, we decided to use these throughout the formerly more dilute region of the system. parameters in the simulations presented in the main text. We note, that the pair density is not a conserved quantity, When larger system sizes are involved, as it is the and it can increase in the process of unitary dynamics. case for the simulations actually used in the main text, In contrast, for the dilute case the pair density pro le comparison with exact results is not available. There- has a pronounced exponential tail away from the initial fore, other indicators for the accuracy must be stud- ergodic region. This shows that pairs spreading away ied. Among these, energy conservation through the time from the initial bubble do not delocalize when encoun- evolution is a straightforward probe. The energy de- tering additional particles on their way. Indeed, while the viation E(t) = jE(0) E(t)j=E(0), where E(t) = late time pair density pro le has small peaks around the initial position of particles, these peaks are not very pro- h (t)j Hj (t)i, allows to control the propagation of the nounced. In addition, the study of the pair density pro le error during the Trotter time evolution. In Fig. S7 we show the results for E(t) in the two quenches presented in the main text, for L = 30. The two plots highlight that the average energy deviation is very small in both 17.5 con gurations. In spite of that, a clear di erence can be 15.0 observed among the two quenches, noticing that the non- 12.5 uniform state has larger error. This is probably due to the enhanced entanglement caused by the presence of the 10.0 bubble in the lattice. Nevertheless, E(t) remains very 7.5 small even at long times, thus con rming the reliability 5.0 of our long-time numerical simulations. In all the simulations performed using ITensor [34], we 2.5 used the U (1) symmetry implementation. In particular, 0.0 200 400 500 1000 1500 2000 to obtain the numerical results presented in Fig. 2(a) and Bond Dimension Bond Dimension (d) we set the maximum bond dimension to be 500 and 3000 respectively. As the histograms in Fig. S8 show, all the disorder realizations remained well below the maxi- Figure S8. The histogram representing the maximum bond mum threshold. This fact ensures that we have a control dimension of di erent disorder realizations show that the on the error encountered in the evolution, in contrast to threshold values of 500 and 3000 for the uniform density wave time evolution with TDVP with xed bond dimension, (left) and bubble states (right) were never saturated in our where error estimation is more challenging [51]. simulations. E(t) Counts Counts S5 • • ◦ • • ◦ • • ◦ ◦ ◦ • ◦ ◦ • ◦ ◦• L = 10 L = 12 vN L = 14 − 1 L = 16 L = 18 − 2 L = 20 L = 30 − 3 1 4 7 10 13 16 19 22 25 28 − 1 0 1 2 3 10 10 10 10 10 Figure S9. The nite size scaling of the pair density hn n i i i+1 shows opposite trend for the dense and dilute cases. The red- Figure S10. The di erent contributions to the entanglement shaded curves represent  = 1=2 con gurations: increasing entropy of the bubble show that the overall behavior of the the system size (from yellow to dark red) the pair density be- von Neumann entropy is faster than logarithmic. Neverthe- comes more uniform and approaches the thermal value, hence less this behavior can be ascribed to the sole con gurational in the thermodynamic limit the probability of nding a pair entropy S , while the particle transport contributes to the far from the bubble is almost the same as nding it in the bub- purely logarithmic growth. The curves are obtained through ble. On the contrary, blue curves ( = 1=5) show exponential Krylov evolution up to T = 1000 averaged over 1000 dis- max vanishing of the pair density and, furhtermore, increasing sys- order realizations for L = 18 and W = 6:5. tem size (from light blue to dark blue) the density decreases, suggesting that at the thermodynamic limit there will be no pair outside the thermal region. Data were obtained with ED, con gurations with the same particle number. Interest- Krylov (T = 1000) and TEBD (T = 500) algorithms max max ingly, Fig. S10 shows that while the overall entanglement averaging over 100 disorder samples for the largest MPS sim- entropy grows faster than logarithmic, this is due only to 4 4 3 3 ulations (L = 20; 30), 3  10 ; 10 ; 5  10 and 10 for ED 3 the con guration part (yellow curve) and the entangle- (from L = 10 to L = 16) and over 10 for Krylov algorithm ment due to particle transport has logarithmic growth. (L = 18). The logarithmic growth of S is consistent with the log- arithmic particle transport presented in Fig. 2(h) in the main text and with other transport measures presented in the uniform density wave at  = 1=5 reveals an almost in the next section. We identify this behavior as a hall- constant behavior, centered around hn n i  10 , i i+1 mark of MBME, and note that it happens on long, yet which corresponds to the values reached at the end of the experimentally accessible timescales t  50(~=t ). exponential tail in the system with L = 30 in Fig. S9. Recent work [40] demonstrated that the logarithmic Next, we focus on understanding di erent contribu- tions to entanglement growth. Exploiting the U (1) sym- metry of our model and following Refs. [38, 39], we split the von Neumann entanglement entropy into a con g- vN uration and a particle transport contributions. Indeed, due to conservation of the total number of bosons the full reduced density matrix  must have a block-diagonal (n) form. Individual blocks within  can be written as p  , where p gives the probability to have n particles in the (n) (n) 0.5 susbsystem A and  is normalized as tr  = 1. Us- ing such representation of the reduced density matrix we 0.2 can split the full entropy into S = S + S as: vN C n 0.1 − 1 0 1 2 3 (n) (n) 10 10 10 10 10 S = tr  log  = p tr  log p vN n n X X (n) (n) (S1) = p log p p tr  log n n n Figure S11. The log-log plot of the entropies discussed in n n Eq. (S1) con rms that only the total entropy is growing as a = S + S : power-law, while both S and S grow slower. In particular, n C n c S grows as a rst degree polynomial in log(t) (dashed green In this way the entanglement growth is split into two line) and S as a second degree polynomial in log(t) (orange dashed line). The sum of these two behaviors (red dashed contributions: one coming from the particle transport, line) agrees with the power-law behavior (dashed blue line). and another originating from dephasing between di erent hn n i i i+1 S(t) S(t) S6 L/2 L growth of the number entropy is expected in the thermal phase, provided there is particle transport over distances that increase as a power-law in time, l / t . The au- (a) thors also predict a power-law scaling of the con guration entropy, which we do not observe, as shown in gure S11. 0.015 We attribute the slower than power-law growth of con g- 0.010 uration entanglement to the localized nature of the right half of the chain, which in turn reduces the number of 0.005 possible con gurations until the particle transport from the left half leads to delocalization. 0.000 0 1 Further analysis of the entanglement dynamics shows 10 10 that S grows in a power-law fashion S  at , cor- vN vN responding to the dashed blue t in Fig. S11, over a rel- (b) evant time interval. On the other hand, both S and 0.2 S behave as polynomials in log(t), of rst and second degree respectively, and their sum (dashed red line) rea- sonably approximates the power-law behavior of S , as vN 0.1 b 2 one can expand t  1 + b log(t) + log (t). Finally, to support our interpretation of logarithmic in- 0.0 crease of S as due to transport, we study the dynamics − 1 0 1 2 3 10 10 10 10 10 of density correlation functions and uctuations. Fig- ure S12 presents the dynamics of connected correlation (c) functions C (i; t) = hn ^ n ^ ihn ^ihn ^ i with respect to i i L=2 L=2 the central site of the chain, the local density uctua- 1.0 2 2 tions n = hn ^ i hn ^ i and the density uctuations in i i 2 2 the dilute part of the chain n = hn ^ ihn ^ i , where R R n ^ = n ^ . The logarithmic dynamics of these R i 0.5 i=L=2 W =5.0 W =6.0 quantities is consistent with the behavior of number en- W =5.5 W =6.5 tropy, thus proving the further support for the existence of slow transport in the dilute part of the chain. 0.0 − 1 0 1 2 3 10 10 10 10 10 Figure S12. Time dynamics of correlation functions (a), lo- cal density uctuations (b) and density uctuations in the B. Quench dynamics from a pair density wave state dilute half of the chain (c) all show, after an initial power-law growth, logarithmic increase with time. In particular, corre- lations far away from the central site (blue curves, as encoded Below we consider quench from a pair-density wave of in the legend above) show signs of a logarithmic light-cone. period 2=. These con gurations accommodate the max- Similarly, local density uctuations deep in the localized re- imal possible number of pairs in the uniform state. Fig- gion present slower dynamics. Finally, panel (c) shows how ure S13 con rms that such state is localized at  = 1=5 increasing disorder slows the growth of the global density uc- and is relaxing in the dense case. Dense systems display tuation of the dilute half. These results were obtained with the Krylov method on system size L = 18 for W = 6:5 (a), strong dependence on the system size and increased ten- with ED on system size L = 16 and W = 6:5 (b) and for dency towards relaxation at larger system sizes, L. In di erent disorder values (c), averaging over 200 disorder real- contrast, at  = 1=5 the late time density pro le has al- izations. most no dependence on the size of the system. In particu- lar, even at very large lengths the curves do not approach the average density represented by the dashed black line. V. P AS A FUNCTION OF BUBBLE DISTANCE In the text we introduced a quantity P (j i ;j i) that 1 2 measures how similar the expansion of over eigen- 1;2 states is. We bire y mentioned that the case where We note, that although the authors of Ref. [40] report unbounded j i = j i corresponds to the usual participation ratio, 1 2 growth of the number entanglement in the MBL phase, the suc- hence we refer to the inverse of P as mutual inverse par- cessive work of Luitz and Bar Lev [41] shows that this is due to ticipation ratio (mIPR). The mutual IPR assumes very rare particle uctuations around the boundary between the two subsystems and the growth disappears at large enough system di erent values depending on the nature of the two states: size. values of mIPR O(N ) correspond to two vectors that n (t) C(i, t) R n (t) i S7 ••◦◦◦◦◦◦◦◦••◦◦◦◦◦◦ ◦◦ (a) we illustrate the behavior of mIPR between pair of states L =40 L =20 L =10 1 . 00 which correspond to a smaller bubble displacement. L =25 L =15 In our analysis we measure the mIPR, P = 0 . 75 P ( ; ), between the following states in the dense L d 0 . 50 limit (half- lling, L = 12), 0 . 25 j i =  ; (S2) 0 . 00 j i =    : (S3) 5 10 15 20 25 30 35 40 |{z} i d ••◦◦••◦◦••◦◦•• ◦◦ (b) 1 . 0 Here we use the bubble that contains all particles to max- imize the range of achievable displacements. For the di- L = 10 L = 14 L = 18 0 . 8 lute case, L = 15, we use similar pair of states with 0 . 6 bubble containing 3 particles ( = 1=5). In the thermal phase, eigenstates are approximately 0 . 4 given by random vectors in the Hilbert space and 0 . 2 L = 12 L = 16 their average overlap with other normalized vectors ap- 0 . 0 proaches the value predicted by random matrix theory, 1 3 5 7 9 11 13 15 17 irrespective of the state or the eigenstate. In the weak disorder limit, we then expect P to be independent (c) on the distance between the two bubbles and to have n = 1 / 5 the same behavior as the conventional IPR: P  N . n = 1 / 2 This expectation is con rmed by the results presented in g. S14(a) and (b) for W = 0:5. On the other hand, in the MBL phase eigenstates are −1 not similar to random vectors, but instead are character- ized by a set of local integrals of motion that have a nite 10 15 20 25 30 35 40 overlap with the local particle density. Thus, two prod- uct states with globally di erent arrangement of particles are expected to have drastically di erent expansion over eigenstates. Therefore, we expect P / exp [d=]. As Figure S13. The late time density pro les of the pair density presented in gure S14(a) and (b), at strong disorder our waves at dense,  = 1=2, and dilute,  = 1=5, llings show results support this hypothesis for both dilute (a) and very di erent behavior. (a) Dilute con gurations are essen- dense (b) states. tially frozen, and do not approach the thermal density repre- At intermediate disorder strength, we observe a quali- sented by the black dashed line. (b) In contrast, at  = 1=2 tative di erence between dense and dilute cases. Dilute relaxation is enhanced at larger L. (c) The deviation of late time density from the thermal value decay exponentially with con gurations, Fig. S14(a), show exponential behavior pair L= pair T already at W = 4:5, whereas dense states in Fig. S14(b) system size as e with   8:1. In contrast for = 1=5, the residual density remains nearly constant with need much stronger disorder to clearly present the same system size and na ve t to the exponential gives an order trend. This result con rms the presence of mobility dilute of magnitude larger scale,   84. Data at  = 1=5 is T edge and is consistent with the observed absence of pair obtained via ED (L = 10, L = 15, and L = 20 with 5 10 , spreading reported in Figure S9 and also with the nite 4 3 10 , and 2  10 disorder realizations), Krylov time evolu- size scaling of mIPRs shown in the main text, Fig. 3. 3 3 tion (L = 25, T = 10 and 10 disorder realizations) and max TEBD (L = 40, T = 300 and 100 disorder realizations). max For  = 1=2 we used ED (L = 10, 12, 14, and 16 with 3 10 , 4 3 3 10 , 5 10 , and 10 disorder realizations) and Krylov time VI. DYNAMICAL PROBE OF THE ABSENCE 3 3 evolution (L = 18, T = 10 , and 10 disorder realizations). max OF RESONANCES The discussion on mutual IPR showed how tunneling processes are strongly suppressed in the dilute case of our have similar expansion over eigenstates, while very large model. In addition to eigenstates analysis, we also stud- values of mIPR imply that the expansion is very di er- ied long time dynamics of states with a thermal bubble. ent. In the main text we analyzed the mIPR between In this way, it was possible to verify whether a bubble two product states where the bubble is located at the initialized at a certain position can dynamically give rise left and right end of the system respectively, see Eqs (2)- to a dense region somewhere else in the chain. In order to (3). Such pair of states corresponds to the maximum study this process we de ned a projector onto the subset possible displacement of the bubble in the chain. Below of Hilbert space that has large density in a certain region. Δ S8 (a) (b) 15 15 10 10 13 W =0.5 W =6.5 13 W =0.5 W =6.5 W =2.5 W =8.5 W =2.5 W =8.5 11 11 W =4.5 W =10.5 W =4.5 W =10.5 9 9 − − 7 7 2 4 6 8 10 12 1 2 3 4 5 6 d d Figure S14. The mutual IPR, I , that quanti es the inverse probability of bubble tunneling d sites, increases exponentially with d at strong disorder. At weak disorder the mIPR approaches the Hilbert space dimension, N , shown by a dashed line. In the dilute system in (b), W = 4:5 marks the onset of the exponential growth, suggesting that the thermal bubble is frozen at its initial position. On the other hand, for  = 1=2, in (a), the clear exponential behavior emerges only at larger disorder. I was calculated for system sizes L = 15 and L = 12 in dilute and dense case respectively and averaged over 10 disorder realizations. More speci cally, we de ne size of the region L , we use the lengthscale extracted from the decay of n. Fit in Fig. 2(c) in the main text P (L ; i) = j ih j ; (S4) yields L ' 6 7, while t in Fig. S13(c) gives a some- c 0 what larger scale. We de ne an initial state j i that has j i2C an entangled dense region of size approximately L (de- where states ji are all possible product states that sat- scribed by a linear superposition of product states j i) isfy the condition    in the region [i; i + L ]. This c 0 followed by a product state: projector selects all con gurations where the system is lo- ^ N cally above the mobility edge. We notice that P (L ; i) X j i = p j i ji : (S5) takes into account all possible con gurations, thus con- 0 i i=1 sidering also the entropic factor. In order to understand what is the minimal required Below, we x the overall density to  = 1=4 and W = 6:5, which still corresponds to a localized system. The dense region is obtained as a superposition of di erent con gurations with N 1 particles in L = 2(N 1) L =12 sites. The remaining particle is initialized in the middle L =16 of the last segment of the chain. For instance, for L = 16 L =20 − 1 this results into following initial state: j i =   + − 2 (S6) 10   + + : : : ; 0 2 4 6 8 10 12 where the boxed area contains a dense entangled bubble and the remainder is in the dilute state. The initial state j i is then evolved through the Hamiltonian (1) in a quench protocol. After time evo- Figure S15. The late time evolution of hP (L ; d)i shows ex- n 0 lution up to a maximum time T = 1000, we measure max ponential decay for all the system sizes studied (L = 12; 16; 20 hP (L ; d)i = h (t)j P (L ; d)j (t)i, which quanti es 0  0 c c at density  = 1=4). Furthermore, we notice that increasing the probability of encountering a bubble shifted by d sites the system size the exponential vanishing becomes more se- from the initial position of the bubble. vere, suggesting that in the thermodynamic limit there would be no motion of the bubble at all. These results were obtained Finally, averaging over all di erent product states in 4 3 3 using 10 ; 5 10 and 10 disorder realizations for the system the dilute part of the chain and over disorder we ob- sizes from smaller to larger. tain the data in Fig. S15. This plot reveals that the hP i d S9 probability of having a dense ( >  ) region decays onantly but rather tunnels throughout the system. More- exponentially with the distance d from its initial loca- over, the nite size scaling analysis shows that increasing tion. This is in agreement with our long-time TEBD the system size the decay of hP (L ; d)i with distance d dynamics, Fig. S9, that reveals localization of individual is enhanced. Therefore in the dilute regime of our model pairs. Thus, we conclude that bubble does not spread res- the bubble remains localized around its initial position.

Journal

Quantum Physics – arXiv (Cornell University)

Published: May 6, 2020

There are no references for this article.