# Early breakdown of area-law entanglement at the many-body delocalization transition

###### Abstract

We introduce the numerical linked cluster (NLC) expansion as a controlled numerical tool for the study of the many-body localization (MBL) transition in a disordered system with continuous non-perturbative disorder. Our approach works directly in the thermodynamic limit, in any spatial dimension, and does not rely on any finite size scaling procedure. We study the onset of many-body delocalization through the breakdown of area-law entanglement in a generic many-body eigenstate. By looking for initial signs of an instability of the localized phase, we obtain a value for the critical disorder, which we believe should be a lower bound for the true value, that is higher than current best estimates from finite size studies. This implies that most current methods tend to overestimate the extent of the localized phase due to finite size effects making the localized phase appear stable at small length scales. We also study the mobility edge in these systems as a function of energy density, and find that our conclusion is the same at all examined energies.

Introduction— The eigenstate thermalization hypothesis (ETH) is a powerful statement relating observables of the high energy eigenstates of a quantum many-body system with their thermal expectation values eth (); rigol-eth1 (). However, this principle can be violated in certain systems with strong enough disorder, where even the high energy eigenstates possess only local entanglement mbl-review (); rigol-eth2 (). Anderson localization is a one-body example of this. An area of key interest is how far this localization persists in a many-body state in presence of interactions mbl (). At what point are interactions strong enough that the localization is destroyed and the system obeys ETH? This is the problem of many-body localization (MBL) transition, which is a topic of active research both theoretically and experimentally vosk-review (); muller-review (); laumann-sem (); vidal-tn (); abanin-iom (); muller-dynamics (); vishwanath-loc (); abanin-periodic (); nayak-qc (); rigol-per (); demler-probe (); scardicchio-edge (); mucciolo-stats (); abanin (); huse-quasiper (); eisert-mbl (); altman-random (); sondhi-mbl-top (); grover-liq (); sondhi-order (); liangsheng (); demler (); gross-heis (); demarco-hub (); bloch (); bloch2 ().

The surge of interest in many-body localized systems has motivated many numerical studies. Most studies have focused on exact diagonalization or Lanczos methods which are able to address both sides of the transition in small systems abanin-mbl (); alet-mbl (); huse-rfheisenberg (); scardicchio-mbl (); pollmann-isingmbl (); bardarson-mbl (); reichman-mbl (); bhatt-mbl (); abanin-per (); santos-dyn (); moore-mbl (); moore-mi (); rigol-fluc (); moore-growth (). However, since much about this phase transition is still not well understood, extension of finite size results to the thermodynamic limit can prove difficult. We would like to examine this phase transition using expansion methods, which provide an alternate way of addressing the thermodynamic limit. While standard perturbative series expansions are very powerfuloitmaa (); series (); imbrie (), they suffer from small energy denominators in models with continuous non-perturbative disorder. Thus, we turn to the numerical linked cluster (NLC) expansion nlc (); nlc-ent (); rigol-nlce (), which does not suffer from this problem of small energy denominators.

In this paper, we provide evidence that for a prototypical model of MBL, approaching the critical disorder from the localized side, the localized phase actually becomes unstable only at increasingly long length scales inaccessible to most numerical techniques. This implies that finite size numerical studies on the MBL transition tend to overestimate the extent of the localized phase.

Model— The system we study explicitly is the spin-1/2 Heisenberg Hamiltonian with random fields along the direction,

(1) |

on the 1d chain, where the sum is over adjacent pairs. The random field is picked from a uniform distribution . This is one of the simplest models to study many-body localization on, and has been studied numerically in great detail huse-rfheisenberg (); alet-mbl (); abanin-mbl (); nayak (); scardicchio-mbl (); huse-entspread (). At low , the system is in a thermalizing phase obeying the ETH, while at high , the system is in the localized (MBL) phase.

To identify these different phases, we focus on the entanglement properties of the eigenstates. The typical measure for entanglement in a pure state bipartitioned into two parts and is the von Neumann entropy, defined for some state as

(2) |

where is the reduced density matrix, obtained by tracing over all external degrees of freedom from the density matrix.

A typical eigenstate in an ETH obeying system will exhibit thermal volume law entanglement. The entanglement entropy will approach the classical thermal entropy (required by ETH) and scale with the volume of the regions and . In the localized phase, the eigenstates will instead obey an area-law, scaling with the area between and . This can be understood by regarding them as simultaneous eigenstates of many local operators mbl-review (): only due to mixing contained in operators near the boundary will one get contributions to the entanglement, which therefore grows with the area of the bipartitioning.

Numerical Linked Cluster Expansions— NLC is similar to perturbative series expansions in that interactions within clusters of increasing sizes must be considered, but rather than perturbatively treating each interaction within a cluster, we solve them numerically, typically by exact diagonalization. Our treatment of disorder in NLC is different from usual spin-glass (); rigol-nlce (), allowing us to deal with continuous non-perturbative disorder. The procedure is outlined briefly here.

Let be the order to which we wish to do the calculation. The order of the calculation defined as the number of spins in the largest cluster considered. We identify a finite size region of the infinite system to work with. For the chain, this is simply a length chain, with a bipartitioning cut in the middle. Each of the sites is assigned a field , which is held fixed until the calculation is complete. The reason for this choice of system size is so that the results remain correct for the infinite system, to the desired order, as explained later.

We define a cluster to be a set of sites. A Hamiltonian can be obtained considering only the spins in , and diagonalized numerically to obtain the eigenstates with eigenvalues labeled by . Our quantity of interest, the eigenstate averaged entanglement entropy, can then be calculated as

(3) |

where is the normalization factor and is the inverse temperature.

The entropy for the infinite lattice can be expressed as a sum over the weight of all clusters that can be embedded in the lattice: . The weight of a cluster is then defined recursively by the principle of inclusion and exclusion:nlc ()

(4) |

One can show that only connected clusters which cross the boundary can have a nonzero weight. First, if a cluster does not cross the boundary there can obviously be no entanglement in it or its subclusters, so the weight is trivially zero. Second, proving that only connected clusters can contribute simply amounts to proving that obeys the linked cluster property, that is, for a cluster with two disconnected components and , . This follows from the fact that . Thus, we must simply consider connected clusters of up to size that have sites on both sides of the partition. Our finite size representation was chosen to contain all the necessary clusters of the infinite system up to order . The count for the clusters crossing the boundary scales with the area thus guaranteeing an area-law as long as NLC converges.

Using this, we can obtain a series whose sum gives the total eigenstate averaged entanglement entropy per unit area ,

(5) |

where is simply for the chain. Finally, the entire calculation must be repeated for different realizations of to obtain a disorder averaged value for . footnote1 ()

The NLC scheme used here is slightly different from what has typically been used for random systems in the past, where different embeddings of the same graph are treated identically and one does not need a consistent finite system spin-glass (); rigol1 (); rigol2 (). The more standard scheme has been applied to study MBL systems with discrete disorder rigol-nlce (), which allows one to perform disorder averaging before subgraph subtraction. When there is continuous disorder, partial disorder averaging over a finite number of realizations means that the linked cluster property is only approximate, and thus large errors will build up at high orders. In our approach, the linked cluster property is guaranteed and one is free to average over many realizations of the system. This is a key aspect of our calculation which allows us to treat disordered systems with continuous non-perturbative randomness.

Does it converge?— We first examine . If entanglement satisfies a thermal volume-law, interpreting as a proxy length scale nlc-ent (), we expect to eventually saturate to the volume-law constant at high enough footnote2 (). We should note that our model (Eq 1) does not possess a strongly thermalizing regime, due to the integrability at , and therefore we do not yet see this saturation to the thermal value within our range of footnote3 (). In the localized phase, the additional entanglement due to the addition of one site far away from the cut should become exponentially small with distance, so we expect to decay exponentially to 0 once is larger than some localization length . We define the MBL phase in our study to be one in which the sum of converges exponentially.

Fig 1 shows for a range of values. To estimate convergence or divergence, a linear extrapolation to can be performed. If the extrapolation predicts , we argue that this corresponds to a breakdown of area-law. Although we expect to eventually go to 0 exponentially in the localized phase, this would only happen when our cluster sizes are much greater than the localization length scale. This can be more clearly shown by examining the ratio of terms, which we will discuss next. Note that the case of an area-law with logarithmic corrections would correspond to one where heads linearly to , which we consider in this analysis the boundary between convergence and divergence.

Let us define the ratio of the (disorder averaged) series terms . In the MBL phase, we can say more about the overall trend of . Again interpreting as a proxy length scale, for large , we expect to decrease exponentially with potentially power-law prefactors. The leading contributions at should be of the form , where is some positive number, is the localization length, and is some arbitrary constant mbl-review (). Therefore, in the large limit, discarding terms smaller than , we expect

(6) |

Therefore, plotting versus , should approach from below, with a slope of . However, near the transition can become very large and we do not actually see this behavior within our range of attainable . We can, however, predict whether this kind of exponential convergence is possible given the behavior of the series at finite .

Fig 2 shows the behavior of for our range of and . The trend seems to be for to increase steadily with (although eventually must approach 1 in the delocalized phase). If the system is in the MBL phase and the series is to converge exponentially, we expect that at some , will begin decreasing exponentially and will begin heading towards with slope . Barring bizarre behavior such as increasing to some high value and then suddenly decreasing before finally increasing again towards , this places a restriction on what can be when begins its exponential decay. Because the slope of the approach is negative or 0, . But also , which means this decay can only begin occurring when

(7) |

Therefore, once has increased above , cannot converge exponentially. This is not a rigorous claim, but should be valid as long as behaves in a regular manner. This clearly shows (in Fig 2) that the series for cannot converge exponentially, and thus is not in the MBL phase. However, the series for are still within this region, and thus may diverge or converge. Hence, our result should serve as a lower bound, with our best estimate being at . Going to higher order in NLC can further refine this value.

Discussion— A way of viewing our result huse-thanks () is that we are seeing an instability to thermalization of an almost-localized regime nandkishore-nearly-localized (). That is (in Fig 1), initially acts quite localized in that it is much smaller than the thermal value and is getting smaller as is increased, but may start increasing at higher , signaling the “onset” of thermalization. This onset of thermalization moves to higher as is increased, but goes beyond our range of accessible after . However, by looking in a sensitive way for initial signs of an instability, we are able to place a lower bound for at .

To understand how our analysis is more sensitive to this transition than other methods, let us focus on the entanglement per unit volume . decreasing with system size is often associated with area-law entanglement and therefore localization alet-mbl (). In our study, for a system of size would correspond to the quantity . So even if had already “turned up” and was increasing (clearly thermalizing), would not begin increasing until had increased above the mean of all the previous terms in the series. Our analysis predicts this upturn, which itself would precede estimates from finite size systems using .

Other methods, which are more focused on seeing the full onset of thermalization, do not observe the transition near our bound huse-rfheisenberg (); alet-mbl (); abanin-mbl (); scardicchio-mbl (). Results from Lanczos on systems of up to sites with finite size scaling show evidence for a transition at alet-mbl (). However, the scaling exponent obtained from finite size scaling strongly violates the Harris-Chayes bound in 1d chayesharris (); anushya-fss (), evidence that perhaps they are still far from the true critical point. This implies that finite size effects are significant, even in the systems accessible to Lanczos, and some corrections to the finite size scaling is needed. These effects cause an overestimate of the stability of the MBL phase, which actually becomes unstable to thermalization earlier only at much longer length scales. Note that an interesting alternate possibility is the existence of an intermediate phase between the ETH and MBL phase, with neither thermal nor area-law entanglement grover (); scardicchio-nonerg-extended (); xiaopeng (), which we do not pursue further.

This onset of thermalization at high order is what one would expect from a long lengthscale delocalization mechanism. In studies of this transition by a renormalization group approach, one also finds that the transition is driven by rare metallic inclusions potter (); vosk (). Near the transition, these are rare enough that small systems look localized, but actually become thermalizing at long length scales.

The mobility edge— Finally, we can observe the transition at different energy windows by varying in Eq. 3 of our NLC calculation, thus probing states at a given energy defined by the thermal ensemble. Following the same arguments as at , we can obtain estimates for at a given temperature or energy density. Fig 3 shows our estimates for various values, with the scaled energy on the vertical axis: , where and are the lowest and highest energies in the energy spectrum. The shape of our estimates are very similar to those obtained in previous numerical calculations alet-mbl (); abanin-mbl (), along with the slight asymmetry expected around laumann-edge (). As with the case at , we find our estimates are consistently higher than previous numerical calculations.

There is much debate on whether a mobility edge exists in the thermodynamic limit. Numerical results suggest the presence of such an edge alet-mbl (); abanin-mbl (); laumann-edge (), but there are also arguments against it schuilaz-bubbles (). While our phase diagram shows a similar shape as previous numerical studies, our analysis gives a lower bound for , and hence does not negate the claims of the absence of a mobility edge. If the transition actually occurs at a single for all energies, the fact that our estimates are lower away from the center of the spectrum would indicate that much larger length scales would be needed to observe delocalization and that finite size effects would be much stronger in those regions.

Conclusions— In conclusion, we have studied the MBL transition in the random field Heisenberg model using NLC expansions. We focus on the breakdown of the area-law of entanglement in the eigenstates of the Hamiltonian. Our approach works directly in the thermodynamic limit and does not rely on any finite size scaling. By looking for signs of instability in the MBL phase, we are able to estimate a lower bound for the critical disorder in the thermodynamic limit. At all energies examined, our estimates are consistently higher than found by finite size studies. This implies that numerical methods which look for the full onset of thermalization tend to overestimate the extent of the MBL phase, which actually becomes unstable earlier but only at much longer length scales. Near the transition, finite size effects are significant and hence caution must be taken when relating to the infinite system.

Our result can be readily verified by cold atoms experiments, which are able to present very well characterized systems gross-heis (); demarco-hub (); bloch (); bloch2 (). If the 1d random field Heisenberg model is experimentally realized, measurements of the critical disorder for large systems should lie above our estimate. We may also extend our result for other similar models of the disorder driven MBL transition, where delocalization also occurs over a long lengthscale potter (); vosk (), and suggest that the true critical point would be higher than finite size scaling estimates from small systems. Also of interest are quantum chaotic Wannier-Stark systems wannier1 (); wannier2 (), which are experimentally accessible and possess a localization-delocalization transition in the absence of disorder.

###### Acknowledgements.

We would like to thank David Huse for many valuable discussions. This work is supported in part by NSF grant number DMR-1306048.## References

- (1) J. M. Deutsch, PRA 43, 2146 (1991); M. Srednicki, PRE 50, 888 (1994).
- (2) Marcos Rigol, Vanja Dunjko and Maxim Olshanii, Nature 452, 854–858 (17 April 2008)
- (3) R. Nandkishore and D. A. Huse, Annual Review of Condensed Matter Physics, Vol. 6: 15–38 (2015).
- (4) Marcos Rigol, Phys. Rev. Lett. 103, 100403 (2009)
- (5) L. Fleishman and P. W. Anderson, Phys. Rev. B 21, 2366 (1980); B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Phys. Rev. Lett. 78, 2803 (1997); I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, PRL 95, 206603 (2005); D. M. Basko, I. L. Aleiner and B. L. Altshuler, Annals of Physics 321, 1126 (2006); V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
- (6) Ehud Altman and Ronen Vosk Annual Review of Condensed Matter Physics Vol. 6: 383–409 (2015)
- (7) Xiaoquan Yu and Markus MÃ¼ller Annals of Physics 337 55 (2013)
- (8) Anushya Chandran and C. R. Laumann, Phys. Rev. B 92, 024301 (2015)
- (9) A. Chandran, J. Carrasquilla, I. H. Kim, D. A. Abanin, and G. Vidal, Phys. Rev. B 92, 024201 (2015)
- (10) Anushya Chandran, Isaac H. Kim, Guifre Vidal, and Dmitry A. Abanin, Phys. Rev. B 91, 085425 (2015)
- (11) Mauro Schiulaz, Alessandro Silva, and Markus MÃ¼ller, Phys. Rev. B 91, 184202 (2015)
- (12) Yasaman Bahri, Ronen Vosk, Ehud Altman and Ashvin Vishwanath, Nature Communications 6, Article number: 7341 (2015)
- (13) Pedro Ponte, Anushya Chandran, Z. PapiÄ, Dmitry A. Abanin Annals of Physics Volume 353, February 2015, Pages 196â204
- (14) Bela Bauer and Chetan Nayak, Phys. Rev. X 4, 041021 (2014)
- (15) Luca DâAlessio and Marcos Rigol, Phys. Rev. X 4, 041048 (2014)
- (16) M. Serbyn, M. Knap, S. Gopalakrishnan, Z. PapiÄ, N.âY. Yao, C.âR. Laumann, D.âA. Abanin, M.âD. Lukin, and E.âA. Demler, Phys. Rev. Lett. 113, 147204 (2014)
- (17) C.âR. Laumann, A. Pal, and A. Scardicchio, Phys. Rev. Lett. 113 200405 (2014)
- (18) Claudio Chamon, Alioscia Hamma, and Eduardo R. Mucciolo, Phys. Rev. Lett. 112, 240501 (2014)
- (19) Maksym Serbyn, Z. PapiÄ, and D. A. Abanin, Phys. Rev. B 90, 174302 (2014); ibid. Phys. Rev. Lett. 111, 127201 (2013); ibid. Phys. Rev. Lett. 110, 260601 (2013)
- (20) Shankar Iyer, Vadim Oganesyan, Gil Refael, and David A. Huse, Phys. Rev. B 87, 134202 (2013)
- (21) Christian Gogolin, Markus P. MÃ¼ller, and Jens Eisert, Phys. Rev. Lett. 106, 040401 (2011)
- (22) Ronen Vosk and Ehud Altman, Phys. Rev. Lett. 112, 217204 (2014)
- (23) Anushya Chandran, Vedika Khemani, C. R. Laumann, and S. L. Sondhi, Phys. Rev. B 89, 144201 (2014)
- (24) Tarun Grover and Matthew P A Fisher J. Stat. Mech. (2014) P10010
- (25) David A. Huse, Rahul Nandkishore, Vadim Oganesyan, Arijeet Pal, and S. L. Sondhi, Phys. Rev. B 88, 014206 (2013)
- (26) Liangsheng Zhang, Hyungwon Kim, David A. Huse, Phys. Rev. E 91, 062128 (2015).
- (27) Kartiek Agarwal, Sarang Gopalakrishnan, Michael Knap, Markus Muller, and Eugene Demler, PRL 114, 160401 (2015)
- (28) Sebastian Hild, Takeshi Fukuhara, Peter SchauÃ, Johannes Zeiher, Michael Knap, Eugene Demler, Immanuel Bloch, and Christian Gross, Phys. Rev. Lett. 113, 147205 (2014)
- (29) S.âS. Kondov, W.âR. McGehee, W. Xu, and B. DeMarco, Phys. Rev. Lett. 114, 083002 (2015)
- (30) Michael Schreiber, Sean S. Hodgman1, Pranjal Bordia, Henrik P. LÃ¼schen, Mark H. Fischer, Ronen Vosk, Ehud Altman, Ulrich Schneider, Immanuel Bloch, Science 349 (6250): 842–845
- (31) Pranjal Bordia, Henrik P. LÃ¼schen, Sean S. Hodgman, Michael Schreiber, Immanuel Bloch, Ulrich Schneider, arXiv:1509.00478.
- (32) M. Serbyn, Z. Papic, D. A. Abanin, arXiv:1507.01635.
- (33) D. J. Luitz, N. Laflorencie, F. Alet, Phys. Rev. B 91, 081103 (2015)
- (34) A. Pal, D. A. Huse, Phys. Rev. B 82, 174411 (2010).
- (35) A. De Luca, A. Scardicchio, EPL 101 (2013) 37003
- (36) J. Kjall, J. Bardarson, F. Pollman, Phys. Rev. Lett. 113, 107204
- (37) Soumya Bera, Henning Schomerus, Fabian Heidrich-Meisner, and Jens H. Bardarson, Phys. Rev. Lett. 115 046603 (2015)
- (38) Yevgeny Bar Lev, Guy Cohen, and David R. Reichman, Phys. Rev. Lett. 114, 100601 (2015)
- (39) Sonika Johri, Rahul Nandkishore, and R.âN. Bhatt, Phys. Rev. Lett. 114, 117401 (2015)
- (40) Pedro Ponte, Z. PapiÄ, FranÃ§ois Huveneers, and Dmitry A. Abanin, Phys. Rev. Lett. 114, 140401 (2015)
- (41) E. J. Torres-Herrera and Lea F. Santos, Phys. Rev. B 92, 014208 (2015)
- (42) R. Vasseur, S. A. Parameswaran, and J. E. Moore, Phys. Rev. B 91, 140202(R) (2015)
- (43) Yichen Huang and Joel E. Moore, Phys. Rev. B 90, 220202(R) (2015)
- (44) Lea F Santos and Marcos Rigol, 2012 Phys. Scr. 2012 014033
- (45) Jens H. Bardarson, Frank Pollmann, and Joel E. Moore, Phys. Rev. Lett. 109, 017202 (2012)
- (46) J. Oitmaa, C. Hamer and W. Zheng, Series Expansion Methods for strongly interacting lattice models (Cambridge University Press, 2006).
- (47) R. R. P. Singh, R. G. Melko and J. Oitmaa, Phys. Rev. B 86, 075106 (2012). T. Devakul and R. R. P. Singh, Phys. Rev. B 90, 125136 (2014); ibid. Phys. Rev. B 90, 054415 (2014).
- (48) J. Z. Imbrie, arXiv:1403.7837.
- (49) M. Rigol, T. Bryant and R.R.P. Singh Phys. Rev. Lett. 97, 187202 (2006).
- (50) A. B. Kallin, K. Hyatt, R. G. Melko and R. R. P. Singh, Phys. Rev. Lett. 110, 135702 (2013).
- (51) B. Tang, D. Iyer, M. Rigol, Phys Rev B 91, 161109 (2014)
- (52) Bela Bauer and Chetan Nayak, J. Stat. Mech. (2013) P09005
- (53) Arun Nanduri, Hyungwon Kim, and David A. Huse, PRB 90, 064201 (2014)
- (54) R. R. P. Singh and S. Chakravarty, Phys. Rev. Lett. 57, 245 (1986).
- (55) M. Rigol, Phys. Rev. Lett. 112, 170601 (2014)
- (56) Deepak Iyer, Mark Srednicki, and Marcos Rigol, Phys. Rev. E 91, 062142 (2015)
- (57) One could consider the median entropy rather than the disorder averaged entropy. Median is not a linear operation, however, but one could still define a series as the difference between the median value of the partial sums up to each . Doing so gives similar results as the mean.
- (58) This was confirmed by performing the expansion for random pure states.
- (59) In a model with a strongly thermalizing phase, such as the Floquet model studied in liangsheng (), we do see saturating at the thermal value.
- (60) We thank David Huse for this insight.
- (61) Sarang Gopalakrishnan and Rahul Nandkishore, Phys. Rev. B 90, 224203 (2014)
- (62) J. T. Chayes, L. Chayes, Daniel S. Fisher, and T. Spencer, Phys. Rev. Lett. 57, 2999 (1986)
- (63) A. Chandran, C. R. Laumann, V. Oganesyan, arXiv:1509.04285
- (64) Tarun Grover, arXiv:1405.1471.
- (65) Andrea De Luca, B. L. Atshuler, V. E. Kravtsov, A. Scardicchio, Phys Rev Lett 113, 046806
- (66) Xiaopeng Li, Sriram Ganeshan, J. H. Pixley, S. Das Sarma, arXiv:1504.00016
- (67) W. de Roeck, F. Huveneers, M. Muller, M. Schiulaz, arXiv:1506.01505
- (68) Andrew C. Potter, Romain Vasseur, S.A. Parameswaran, arXiv:1501.03501
- (69) Ronen Vosk, David A. Huse, Ehud Altman, arXiv:1412.3117
- (70) Ian Mondragon-Shem, Arijeet Pal, Taylor L. Hughes, and Chris R. Laumann, arxiv:1501.03824
- (71) Carlos A. Parra-Murillo, Javier Madronero, and Sandro Wimberger, Phys. Rev. A 88, 032119 (2013); ibid. Phys. Rev. A 89, 053610 (2014).
- (72) Andrea Tomadin, Riccardo Mannella, and Sandro Wimberger Phys. Rev. Lett. 98, 130402 (2007).