# Many-body localization in infinite chains

## Abstract

We investigate the phase transition between an ergodic and a many-body localized phase in infinite anisotropic spin- Heisenberg chains with binary disorder. Starting from the Néel state, we analyze the decay of antiferromagnetic order and the growth of entanglement entropy during unitary time evolution. Near the phase transition we find that decays exponentially to its asymptotic value in the localized phase while the data are consistent with a power-law decay at long times in the ergodic phase. In the localized phase, shows an exponential sensitivity on disorder with a critical exponent . The entanglement entropy in the ergodic phase grows sub-ballistically, , , with varying continuously as a function of disorder. Exact diagonalizations for small systems, on the other hand, do not show a clear scaling with system size and attempts to determine the phase boundary from these data seem to overestimate the extent of the ergodic phase.

###### pacs:

75.10.Jm, 05.70.Ln, 72.15.Rn## I Introduction

It is by now well established that disorder can drive closed one-dimensional quantum many-body systems into a many-body localized (MBL) phase Imbrie2016 (); Review1 (); Review2 (). In such a phase the system fails to act as a bath for its own subsystems and thermalization does not occur. Instead, memory of the initial conditions is retained. The ‘drosophila’ to study properties of the MBL phase is the spin- Heisenberg chain

(1) |

with and a random variable drawn from a uniform box distribution with disorder strength . Here is the length of the system and is the component of the spin operator acting at site . Studies of this model have been based mainly on exact diagonalization (ED) for small systems PalHuse (); Luitz1 (); Luitz2 (); BarLev (); Agarwal (). These numerical results have then been used to determine a critical point between the ergodic and MBL phase by showing, for instance, that the level statistics changes from a Wigner-Dyson distribution at small but nonzero to a Poisson distribution at (MBL) with . Furthermore, deep in the MBL phase the entanglement entropy is shown to increase logarithmically during unitary time evolution BardarsonPollmann (), confirming results from an earlier density matrix renormalization group study ZnidaricProsen ().

ED studies of small systems are, however, ill-suited to address the properties of weakly disordered systems as well as the phase transition itself because in both cases the relevant length scale will be much larger than the achievable system sizes . This creates, in particular, a significant obstacle in understanding this novel type of dynamical phase transition where the entanglement entropy changes from volume law (ergodic) to area law (MBL), making it distinct from regular thermal transitions or ground state critical points. Two approaches have so far been used to tackle this problem: On the one hand, it has been tried to investigate the critical regime based on extrapolations from ED data to larger systems Luitz1 (); BarLev (); Agarwal (); Luitz2 (). Assuming that the transition is described by a single diverging length scale , the obtained results are mostly consistent with a critical exponent . This, however, would violate a Harris-type bound which demands in dimensions in order for the transition to be stable NandkishorePotter (); ChandranLaumann (). A second recent approach is based on a real-space renormalization group (RG) applied to effective minimal models assuming that only two energy scales exist VoskHuse (); PotterVasseur (). The length scale is then found to diverge with an exponent consistent with the Harris bound. However, it is important to stress that the RG approaches are not based on microscopic models and contradict the results from previous ED studies.

In this work we shed new light on this controversial point by studying a disordered interacting quantum chain directly in the thermodynamic limit (TDL). In this way we avoid the fundamental obstacle one faces in ED studies of the phase transition. In the following we focus on the anisotropic Heisenberg chain, Eq. (1), with binary disorder instead of the more commonly used box disorder. This naturally arises as an effective model for a bosonic system with a mobile and an immobile species in the limit of strong onsite Hubbard interactions and also exhibits a transition from an ergodic to an MBL phase.AndraschkoEnssSirker (); TangRigol () As in the noninteracting caseZhangSirker (), one expects that the chosen disorder distribution leads to quantitative changes while the qualitative features, in particular the properties of the transition, are universal. The goals of this work are to establish the phase diagram of the model (1) with binary disorder as a function of disorder strength and anisotropy (see Fig. 1) and to study the ergodic-MBL phase transition directly in the TDL. In order to obtain an exact disorder average in a single simulation, we introduce an ancilla spin-, at each site and replace . The state of then determines the local binary disorder ParedesVerstraete (); AndraschkoEnssSirker (). We consider the unitary time evolution starting from an initial product state in the Hilbert space of spins and ancillas, where represents a superposition of all possible disorder configurations. Following recent experiments SchneiderBloch1 (); SmithMonroe (); SchneiderBloch2 () we prepare the spins in the Néel state ( in the equivalent fermionic picture). We then study the exactly disorder averaged decay of the antiferromagnetic order

(2) |

where measures the staggered magnetization (imbalance) of the physical spins and the identity operator acts on the ancillas. The time evolved state is defined by where includes the coupling between spins and ancillas . In addition, we also study the growth of the disorder averaged entanglement entropy . Here denotes the reduced density matrix of half of the infinite chain consisting of spins and ancillas. Because the density matrix includes the ancillas, the entanglement entropy is quantitatively not the same as for a semi-infinite chain consisting of spins only. The ancillas are, however, completely static so that the entanglement entropies with and without the ancillas show the same scaling with time. We simulate the translationally invariant system of spins and ancillas using the light cone renormalization group (LCRG), a variant of the density matrix renormalization group which yields results directly in the TDL EnssSirker (); AndraschkoEnssSirker ().

We choose the LCRG bond dimension such that the truncation error always remains smaller than . By comparing with results obtained keeping the truncation error smaller than we make sure that our results are numerically exact for the times shown. This requires bond dimensions of up to states. The scales we are reaching in these simulations are unprecedented: at small disorder in the ergodic phase correlations spread approximately ballistically as where is the maximal velocity of excitations in the lattice. The maximal times in our simulations therefore test the system at length scales of at least . While we cannot exclude the possibility that the scaling of the quantities we study changes qualitatively at even larger scales, our data represent substantial progress compared to ED studies which are limited to scales of .

## Ii Decay of order parameter

In the clean free fermion case () the decay of the order parameter is given by with being the Bessel function of the first kind and time measured in units of . For interactions it has been shown that the asymptotic decay in the clean case is well described by the free fermion asymptotics multiplied by an exponential decay BarmettlerPunk (); BarmettlerPunk2 (). Turning on disorder introduces barriers between thermalizing clusters with equal Zeeman field, see inset of Fig. 1(a). In the ergodic phase, a finite thermalization time across such barriers must exist where is the number of jumps of the Zeeman field within the barrier and a function depending on disorder and anisotropy . The probability that a particular site is part of a barrier with jumps is given by . After time only clusters separated by barriers of size will not have thermalized, and the asymptotic decay in the ergodic phase follows

(3) |

up to logarithmic corrections. In the MBL phase, on the other hand, the staggered magnetization will not decay completely, . Combining the different limiting cases, we fit the LCRG data for anisotropies , disorder , and times by the functions

(4) |

with lifetime and exponent of a power-law decay. We perform fits using both fit functions and check for consistency, i.e., in the ergodic phase and in the MBL phase with . As shown in Fig. 2, this leads to excellent fits which allow to extract an estimate for in the MBL phase, in the ergodic phase, as well as the lifetime of the oscillations , see Fig. 3(a,b,c). In Appendix .1 we show that the fit parameters depend only weakly on the time window used for the fit, cf. Fig. 7. In particular, the asymptotic value of the order parameter is very robust in all fits, cf. Fig. 8. As in the clean case BarmettlerPunk (); BarmettlerPunk2 (), we cannot find any fitting function which describes the data for small disorder in the regime well. For , on the other hand, we find that the asymptotics is very well described by a pure non-oscillating exponential decay , see Fig. 2(c).

Based on the RG analysis of a minimal model, an exponential sensitivity of the residual imbalance in the MBL phase () has been predicted PotterVasseur ().

As shown in Fig. 3(a), we obtain an excellent data collapse for different with a critical exponent using and as fitting parameters. The critical values obtained from the data collapse lead to the phase boundary shown in Fig. 1(a). We note that violates the Harris bound, see below. For comparison, the power-law exponent is shown in Fig. 3(b): In a theory with a single length scale , one would expect that where is the dynamical critical exponent PotterVasseur (). However, the fits yield absolute statistical errors in the power-law exponent between making it impossible to extract a scaling close to . The values determined by nevertheless are consistent with, although slightly larger than, the values based on the data collapse for the magnetization. The relaxation time , on the other hand, can be extracted with statistical errors of less than and is shown in Fig. 3(c). For very small disorder we qualitatively find the same behavior as in the clean case BarmettlerPunk (); BarmettlerPunk2 (). The relaxation time decreases approximately as for and increases proportional to for . For and we find that the relaxation times immediately decrease when disorder is added; in a region around , however, the relaxation times remain stable at first before increasing at larger disorder strengths. When plotting the value of the smallest disorder where deviates substantially (by more than ) from the clean case for different anisotropies , we find that this change in relaxation time does occur when crossing from the ergodic to the MBL phase, see Fig. 1(a).

Similarly to the phase diagram for the XXZ chain with box disorder—obtained by ED in Ref. BarLev, —we observe reentrant behavior: for fixed and small in the MBL phase, increasing interactions can first drive the system into the ergodic phase before localization is again stabilized at large interactions.

## Iii Entanglement growth

To investigate the properties of the phase transition in more detail, we now turn to an analysis of the entanglement entropy . Using the same type of argument as for the decay of the order parameter, a power law is expected in the ergodic phase; a power law in the entanglement entropy is also found in the critical Harper model RooszIgloi (). If the RG theories of Refs. VoskHuse, ; PotterVasseur, do describe the transition correctly then holds. On the MBL side, on the other hand, we have shown previously in Ref. AndraschkoEnssSirker, that as is predicted on general grounds.VoskAltman (); SerbynPapicA (); SerbynPapicB (); NanduriKim (); HuseNandkishore () For we fit the LCRG data for to a power law and obtain excellent fits with statistical errors of less than . Furthermore, the exponent , shown in Fig. 4(a), is only weakly affected by a modification of the fit interval provided that the behavior for small times is excluded. For and intermediate disorder, on the other hand, we find two different regimes: a power law increase at intermediate times followed by a much slower increase for , see Fig. 4(b). Because of the limited time range available, it is not clear if the latter regime corresponds to a power-law increase with a smaller exponent or signals a crossover to logarithmic scaling.

Overall, we find a sub-ballistic spreading including, surprisingly, an extended region of disorder strengths for small where the entanglement spreads diffusively, and . Remarkably, for the parameters , where the entanglement spreads diffusively the system seems to be already deep in the MBL phase according to the phase diagram Fig. 1(a); the order parameter shown in Fig. 8 signals localization at least up to times . Physically, this intermediate diffusive regime might be explained by the existence of many relatively narrow barriers between thermalizing segments which lead to diffusion while rare wide barriers lead to an exponential enhancement of the entanglement time and finally prevent the system from fully thermalizing. Our findings might therefore possibly indicate that the transition is not described by a single length scale. In this case the scaling hypothesis is violated and a Harris criterion does not apply. The region with corresponds to classical diffusion implying, in particular, that the spin-spin autocorrelation function decays as with . Note that even in the clean case the presence or absence of diffusion at infinite temperatures in is an open and controversially discussed topic, with numerical results for short and intermediate times showing a power-law decay with an exponent depending on the fit interval FabriciusMcCoy (); SirkerDiff (). At small finite temperatures, on the other hand, has recently been established by field theoretical methods in the TDL and confirmed by numerical data KarraschPereiraSirker ().

## Iv Comparison with ED

While the LCRG data for infinite systems support a consistent interpretation of the MBL transition, it is instructive to compare to exact diagonalization results for finite systems. Two commonly used methods to establish the phase diagram of the disordered model (1) are calculating the level statistics and studying the time average of an order parameter.

To obtain the level statistics we define with the difference between adjacent energy eigenvalues. At the integrable point and also in the MBL phase where additional local conserved charges exist we expect Poisson statistics with an average value , while Wigner-Dyson statistics with is expected in the ergodic phase for oganesyan2007 (); PalHuse (), see Appendix .2 for details. In Fig. 5(a), results for model (1) with binary disorder, , and system sizes are shown where the disorder averages are exact for while inequivalent samples have been used for . Contrary to the box disorder case PalHuse () we do not find a point where appears to be close to stationary which has been interpreted as being indicative of the critical point in the thermodynamic limit. Note, however, that even in the box disorder case a ‘drifting’ of the crossing points of curves with different has been observed for increasing system size.

An alternative criterion to estimate the phase boundary is to fix as the disorder value where crosses the intermediate value between Wigner-Dyson and Poisson statistics. If there is a sharp transition between an ergodic and an MBL phase in the thermodynamic limit then will converge to the critical value in the limit . However, even using this alternative criterion the problem persists that no clear scaling with is obtained for the limited system sizes available. For , takes an intermediate value between Wigner-Dyson and Poisson statistics around disorder (see Fig. 9 in the appendix), which is an order of magnitude larger than the value for established above for the infinite chain.

With increasing anisotropy the system approaches the Ising limit where each local becomes approximately conserved. For small disorder it would then require very large systems to see level repulsion and Wigner-Dyson statistics. We therefore consider only anisotropies using ED, see Fig. 1(b).

A naive linear extrapolation in of the time averaged magnetizations also yields a critical for , see Fig. 5(b). LCRG, on the other hand, shows quite clearly that is already deep inside the MBL phase (see Fig. 2(b2)). Using both and to extract a phase boundary shows that the ED results can lead to a significantly larger extent of the ergodic phase for all , see Fig. 1. This might not be completely unexpected because any system with length much smaller than the localization length will look ergodic.

The difference between infinite and finite systems is exemplified clearly in the time evolution of the order parameter for , shown in Fig. 6. According to the finite-size scaling of the time averaged obtained by ED, this is far in the ergodic phase (cf. Figs. 5(b) and 12). However, we observe that ED for and LCRG for only agree up to . LCRG, instead, shows that saturates for at least up to , which corresponds to an effective system size . Since the LCRG data for the infinite chain test the dynamics at length scales which are a factor larger than the length scales reached in ED, the most plausible explanation for this discrepancy appears to be that the scaling of is non-monotonic. In order to check this tentative explanation, one would need to diagonalize much larger systems.

## V Conclusions

Using time-dependent density matrix renormalization group calculations we have established the phase diagram of the XXZ spin- chain with binary disorder in the TDL. For weak disorder in the ergodic phase we are able to test the dynamics on length scales of the order of lattice sites which is significantly larger than the lengths which can be studied in exact diagonalization. Our results generalize previous studies of the decay of Néel order (imbalance), , from clean to disordered systems which is highly relevant to interpret recent SchneiderBloch1 (); SmithMonroe (); SchneiderBloch2 () and future cold atomic gas experiments. We find that in the MBL phase shows an exponential sensitivity on disorder with a critical exponent near the ergodic-MBL phase transition of . For the entanglement entropy we find a power-law growth at intermediate times with an exponent which varies continuously as a function of disorder. For small we find, in particular, a diffusive growth of entanglement at intermediate times in the MBL phase near the transition while is expected at long times. This intermediate time behavior might indicate a second relevant length scale in the problem. In this case the scaling hypothesis is violated and a Harris bound does not apply.

###### Acknowledgements.

J.S. acknowledges support by the Natural Sciences and Engineering Research Council (NSERC, Canada) and by the Deutsche Forschungsgemeinschaft (DFG) via Research Unit FOR 2316. We are grateful for the computing resources and support provided by Compute Canada and Westgrid.## Appendix

In the Appendix we provide technical details regarding the fitting of the LCRG data and the exact diagonalizations. Furthermore, we present spectra of time averaged magnetizations for individual disorder realizations which show qualitative differences in the ergodic and deep in the MBL phase and might be a useful tool for experimental analysis.

### .1 Fits of the order parameter

Using the LCRG algorithm we have obtained data for the decay of directly in the thermodynamic limit. To analyze these data we have used the two fit functions given in Eq. (4). Here we want to show that these fits are quite stable with regard to the time window chosen provided that one excludes the initial fast decay which is not well described by the fit functions. We always start with the data for the smallest disorder, , using the values for the free fermion case without disorder , , , , and as initial guess. The fitting parameters obtained from the converged least square fit are then used as initial parameters for the fit of the data set with the next larger disorder.

As an example, we show here different fits and additional data for the particularly interesting case . In Fig. 7 we compare the parameters of three different fits.

For all three fits give parameters which are quite close to each other and which change smoothly as a function of disorder. For larger disorder values most parameters remain stable except for the phase shift . We want to stress again that the fit functions are based on the asymptotics for the clean free fermion case and are not expected to yield good fits for large disorder and/or interaction strengths. For we can obtain reasonable fits up to disorder .

For we extract the remaining magnetization at infinite times, , from the fits while for we simply take the average of in the specified time window. The results are shown in Fig. 8(b) and are almost independent of the time window.

We have found that there appears to be a diffusive entanglement spreading for and , see Fig. 4. From the magnetization data it appears, however, that for these disorder strengths we are already in the MBL phase, see Fig. 8(a). A possible interpretation is that the diffusive entanglement spreading only holds at intermediate times while a crossover to the expected logarithmic scaling will happen at larger times, inaccessible to our numerical calculations.

### .2 Exact Diagonalization (ED)

We use ED for small XXZ chains of length to study the level statistics of the disordered Hamiltonian ( values) as well as the time evolution of observables such as the staggered magnetization, .

For each disorder configuration the Hamiltonian (1) conserves the total spin quantum number . Here, we consider chains with no average magnetization, , for even . The energy spectrum determines the level statistics, while all eigenvectors are needed for expectation values such as the staggered magnetization. The computation is repeated for different disorder realizations and the results averaged. In particular, in the case of binary disorder, we explicitly average over all possible disorder configurations. By symmetry, the configurations with flipped disorder or with left-right mirrored disorder yield equivalent results. This reduces the number of inequivalent configurations to in the case of open boundary conditions (OBC). For periodic boundary conditions (PBC), the shift symmetry leads to a further reduction to . For instance for , the dimension of the Hilbert space is , and there are (OBC) and (PBC) unique disorder configurations, resp. For and (PBC) we typically perform complete disorder averages; for (OBC) we sample inequivalent random configurations. This is in contrast with the LCRG algorithm which works in the much larger Hilbert space of spins and ancillas and produces the complete disorder average in a single run AndraschkoEnssSirker ().

#### ED Level statistics

For each disorder configuration we define the level spacing between adjacent energy eigenvalues . In order to normalize the energy scale we consider the ratios which lie between and . The level distribution is then averaged over all binary disorder configurations. In the presence of an extensive set of local conserved charges the level spacing is Poisson distributed with and average value . In the ergodic phase, instead, a Wigner-Dyson distribution (GOE) of favors larger ratios with .

Level spectra for open boundary conditions show less degeneracies as compared to those for periodic boundary conditions and are therefore better suited to determine the phase boundary. In Fig. 9 we show the values for different anisotropies as a function of disorder at fixed . The points where these curves cross the intermediate determine the phase boundary shown in Fig. 1(b) in the main text.

#### ED magnetization

**Time evolution.** The time evolution of the staggered magnetization from an initial
Néel state is computed as the disorder
average of the quantum evolution

In the eigenbasis for each disorder configuration, one can write

Fig. 6 shows the time evolution of the staggered magnetization for a Heisenberg chain with small disorder . The LCRG results are exact for an infinite system and extend to finite times . They provide strong evidence that and that the system is therefore in the MBL phase in accordance with the phase diagram Fig. 1(a) in the main text. The ED time evolution for can be computed for arbitrarily long times but deviates from the LCRG result already for short times due to finite-size effects. The localization length just beyond the MBL transition is much larger than any system size accessible by ED so that ED can only capture the short-time dynamics correctly, making an extrapolation of time averaged data to lengths impossible.

For larger disorder shown in Fig. 10, the ED results for periodic boundary conditions are much closer to the LCRG result and differ visibly only for ( PBC). This is likely due to the proliferation of small localized clusters which are well captured by ED and which dominate the dynamics well inside the MBL phase. In contrast, the ED time evolution for open boundary conditions (upper set of curves) is far from the result even for (OBC) and converges only slowly with increasing .

**Time averaged magnetization.** At long times, the staggered magnetization oscillates around the
average magnetization

(5) |

The contributions with unequal energy dephase and do not contribute to the time average, such that only the energy diagonal terms remain. Note that matrix elements of between different degenerate states vanish, such that the sum in (5) reduces to a single sum with .

The dependence of the average magnetization on disorder and system size is shown in Fig. 11. A nonzero value is obtained for any disordered system with . Note that decreases with increasing for small , while it increases for larger . While a crossing point exists, it does not agree with the phase transition point found by a finite-size scaling analysis in Fig. 5(b) nor with the phase boundary obtained in LCRG for the infinite chain.

A finite-size scaling analysis for the average magnetization from ED with periodic boundary conditions (PBC) is shown in Fig. 12 for the isotropic point. In the ergodic phase we expect the magnetization to decay for an infinite system as where is the critical exponent. Since the total magnetization is conserved, , we expect that spin transport occurs as a random walk similar to energy transport leading to a scaling .VoskHuse (); PotterVasseur () This scaling argument would suggest that . For we find that the scaling of the magnetization in Fig. 12 appears to follow a power law with exponent , or . This seems to further support our findings from the analysis of the entanglement entropy for infinite chains presented in the main text that the dynamics at intermediate times (intermediate lengths) in the MBL phase close to the transition is diffusive. For larger disorder the average magnetization saturates to a finite value. The apparent position of the MBL phase transition with PBC is consistent with, but slightly larger than, the phase boundary obtained for open boundary conditions as shown in Fig. 1(b) in the main text.

**Magnetization spectra.** Using ED we can calculate a time averaged magnetization,
Eq. (5), for each disorder configuration. In
Fig. 13 the values as a function of
the disorder configuration are exemplarily shown for and
.
The two magnetization spectra are qualitatively very different. While
the spectrum shown in Fig.13(b) for large disorder
(deep inside the MBL phase) shows a gap, there is no gap for
(near the phase transition) visible, see
Fig.13(a). For fixed we find that the gap for
is about a factor larger than the gap for . This
provides an estimate for the phase transition which is consistent with
the estimate based on the level spectra for the same system size.

### References

- J. Z. Imbrie, Phys. Rev. Lett. 117, 027201 (2016).
- R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
- E. Altman and R. Vosk, Annu. Rev. Condens. Matter Phys. 6, 383 (2015).
- A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
- D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
- D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 93, 060201 (2016).
- Y. Bar Lev, G. Cohen, and D. R. Reichman, Phys. Rev. Lett. 114, 100601 (2015).
- K. Agarwal, S. Gopalakrishnan, M. Knap, M. Müller, and E. Demler, Phys. Rev. Lett. 114, 160401 (2015).
- J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
- M. Žnidarič, T. Prosen, and P. Prelovšek, Phys. Rev. B 77, 064426 (2008).
- R. Nandkishore and A. C. Potter, Phys. Rev. B 90, 195115 (2014).
- A. Chandran, C. R. Laumann, and V. Oganesyan (2015), arXiv:1509.04285.
- R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X 5, 031032 (2015).
- A. C. Potter, R. Vasseur, and S. A. Parameswaran, Phys. Rev. X 5, 031033 (2015).
- F. Andraschko, T. Enss, and J. Sirker, Phys. Rev. Lett. 113, 217201 (2014).
- B. Tang, D. Iyer, and M. Rigol, Phys. Rev. B 91, 161109 (2015).
- Y. Zhao, F. Andraschko, and J. Sirker, Phys. Rev. B 93, 205146 (2016).
- B. Paredes, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 95, 140501 (2005).
- M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).
- J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe (2016), Nat. Phys. 12, 907 (2016).
- P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider, Phys. Rev. Lett. 116, 140401 (2016).
- T. Enss and J. Sirker, New J. Phys. 14, 023008 (2012).
- P. Barmettler, M. Punk, V. Gritsev, E. Demler, and E. Altman, Phys. Rev. Lett. 102, 130603 (2009).
- P. Barmettler, M. Punk, V. Gritsev, E. Demler, and E. Altman, New J. Phys. 12, 055017 (2010).
- G. Roósz, U. Divakaran, H. Rieger, and F. Iglói, Phys. Rev. B 90, 184202 (2014).
- R. Vosk and E. Altman, Phys. Rev. Lett. 110, 067204 (2013).
- M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 110, 260601 (2013).
- M. Serbyn, Z. Papić, and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013).
- A. Nanduri, H. Kim, and D. A. Huse, Phys. Rev. B 90, 064201 (2014).
- D. A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
- K. Fabricius and B. M. McCoy, Phys. Rev. B 57, 8340 (1998).
- J. Sirker, Phys. Rev. B 73, 224424 (2006).
- C. Karrasch, R. G. Pereira, and J. Sirker, New J. Phys. 17, 103003 (2015).
- V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).