Are many-body localized systems stable in the presence of a small bath?

TL;DR

This paper investigates the stability of many-body localization in large disordered quantum systems, finding that localization may be more robust than previously thought despite the presence of ergodic regions.

Contribution

The study introduces a numerical approach to assess many-body localization stability in the thermodynamic limit considering low-disorder regions.

Findings

Many-body localization appears more stable than recent theories suggest.

Low-disorder regions do not necessarily lead to delocalization in large systems.

Numerical results support the persistence of localization despite ergodic grains.

Abstract

When pushed out of equilibrium, generic interacting quantum systems equilibrate locally and are expected to evolve towards a locally thermal description despite their unitary time evolution. Systems in which disorder competes with interactions and transport can violate this expectation by exhibiting many-body localization. The strength of the disorder with respect to the other two parameters drives a transition from a thermalizing system towards a non-thermalizing one. The existence of this transition is well established both in experimental and numerical studies for finite systems. However, the stability of many-body localization in the thermodynamic limit is largely unclear. With increasing system size, a generic disordered system will contain with high probability areas of low disorder variation. If large and frequent enough, those areas constitute ergodic grains which can hybridize and thus compete with localization. While the details of this process are not yet settled, it is conceivable that if such regions appear sufficiently often, they might be powerful enough to restore thermalization. We set out to shed light on this problem by constructing potential landscapes with low disorder regions and numerically investigating their localization behavior in the Heisenberg model. Our findings suggest that many-body localization may be much more stable than anticipated in other recent theoretical works.

Equations

$$H$$

$$+\sum _{i=1}^sϵh_i{\sigma }_i^z+\sum _{i=s+1}^L\Delta h_i{\sigma }_i^z,$$

$$U_D{\mathcal{Z}}_i{U}_D^{†}={\alpha }_i{\sigma }_i^z+\sum _{\mu ,\nu \in C}{\beta }_i(\mu ,\nu ){\Sigma }^z(\mu ){\Sigma }^x(\nu ),$$

Full text

I Introduction

One of the long-standing puzzles of physics is how the postulates of quantum statistical mechanics and thermodynamics can be made compatible with the unitary time evolution of quantum systems. It is increasingly becoming clear that generic interacting quantum systems – once pushed out of equilibrium – are expected to dynamically evolve towards a locally thermal description again Rigol et al. (2008); Polkovnikov et al. (2011); Eisert et al. (2015); Gogolin and Eisert (2016); Langen et al. (2015). This constitutes an interesting state of affairs, since it reconciles the apparent contradiction between the description of time-evolving states and of equilibrium ensembles. Such an interpretive scheme also provides a picture in which interacting systems can be described by only a small number of parameters for almost all time intervals, thus avoiding the curse of dimensionality. A few exceptions are known to exist but these are fine-tuned integrable systems featuring local constants of motion that prohibit a general description in terms of thermal ensembles.

At the turn of the millennium, a new class of quasi-integrable systems emerged which are fail to thermalize over a wide range of parameters. As this effect is caused by the subtle interplay of transport, interactions and disorder, it has been dubbed many-body localization (MBL) Gornyi et al. (2005); Basko et al. (2006); Oganesyan and Huse (2007). A many-body localized system does not exhibit transport of particle-like quantities and therefore can be effectively described by an extensive set of quasi-locally conserved constants of motion (qLCOMs) Serbyn et al. (2013); Huse et al. (2014). This leads to local memory of particle configurations Schreiber et al. (2015). Note that unlike the non-interacting Anderson insulator, systems exhibiting MBL may well transport information-like quantities such as correlations between particles. This is reflected e.g. in the logarithmic growth of the entanglement entropy in time Znidaric et al. (2008); Bardarson et al. (2012).

While many of the properties ascribed to MBL have been observed either experimentally or numerically in finite systems, the question of the stability of MBL in the thermodynamic limit is as of yet unresolved. The expected leading sources of instability are rare regions of low disorder that might conceivably thermalize the rest of the system Ponte et al. (2017); De Roeck and Huveneers (2017); Luitz et al. (2017); Thiery et al. (2018). Questions regarding the effect of local regions that are localized in the otherwise ergodic phase or ergodic regions in the otherwise localized phase have been investigated in the field of Griffiths effects Nandkishore et al. (2014); Gopalakrishnan et al. (2015); Agarwal et al. (2015); Luitz and Lev (2017); Nandkishore and Gopalakrishnan (2017); Gopalakrishnan and Huse (2019) (for a review, see Ref. Agarwal et al. (2017)). The question of whether a closed MBL system is stable to uncharacteristic disorder potentials is hence part of this sub-field.

In this work, we deliberately construct potential landscapes with regions of improbably low disorder and study their influence on a localized chain surrounding them. A previous study has reported that a constant size thermal region is able to thermalize a localized chain coupled to it independent of its size if the constants of motion do not decay sufficiently strongly Luitz et al. (2017). The study has been conducted within an effective description of MBL Serbyn et al. (2013); Huse et al. (2014), one in which the localized part has been modelled by constants of motion and the bath by a suitable random matrix. The postulated coupling between bath and localized part is characterized by a decay length which is related to the localization length of the constants of motion.

Here, instead we use its real space equivalent; the disordered Heisenberg chain. Employing exact diagonalization, we analyze the statistics of local expectation values of the eigenstates which allow us to draw conclusions about the locality of the constants of motion of the system. Our results are not entirely compatible with those of Ref. Luitz et al. (2017), as we find MBL is not compromised by the presence of low disorder regions of constant size. When the size of the region is allowed to scale with system size however, we find that localization vanishes when extrapolating our results to large systems.

II Setting

II.1 Hamiltonian model

We consider the “drosophila” of MBL, the disordered spin-$1/2$-Heisenberg chain. In order to study the effect of small subregions of low disorder we consider systems where the potential on the sites $1,\dots,s$ is close to zero. For a system of $L$ sites our Hamiltonian then reads

[TABLE]

where the $h_{i}$ are drawn uniformly and independently from the interval $[-1,1]$, $\epsilon$ and $\Delta$ denote the disorder strength in the low and high disorder region and $\sigma^{a}_{i}$ is the Pauli-$a$ operator acting on the $i$-th site. The fully disordered model (without ergodic subregion) is expected to undergo a localization transition at a critical disorder strength $\Delta_{c}\approx 7$ (note that since we use Pauli instead of spin operators, the transition is shifted with respect to other literature). Moreover, we use periodic boundary conditions and work in the zero magnetization sector. We set $\epsilon=10^{-6}$ independent of the disorder strength $\Delta$. This creates the situation that the local flatness of the potential may thermalize the subsystem, which is expected to compete with localization effects in the remaining system. In the following, we will hence refer to the first $s$ sites as the thermal sites and refer to the remaining $L-s$ sites as the disordered sites or part of the system. We consider two different scenarios: In the first, the number of thermal sites is independent of the system size. In the second setting, we take a fixed fraction of the full size to be thermal.

II.2 Distribution of disorder

Within this model, we investigate below the effect of a local cluster with uncharacteristically small disorder on its environment. The formation of such a cluster should be understood as a rare instance of a generic MBL model and is specifically relevant when considering the stability of an interacting localized phase in the thermodynamic limit, where the effect of these rare regions is not fully understood yet. For large systems, the probability of having small thermal regions increases and poses the question if these are able to compromise localization of the full chain as well. As a prerequisite, these must first thermalize their neighbours which can also be checked in small systems and which is precisely the context of this work. Note that if one tries to draw conclusions for larger systems from our results, these can only hold if the reductions of the eigenstates resemble the eigenstates of the small system when probed locally. Global quantities such as energy levels commonly used to detect MBL Oganesyan and Huse (2007); Luitz et al. (2015) will likely not be faithfully recovered in our small system and are hence not considered here.

III Measure of localization

The proposed effective integrability of MBL comes about due to the presence of extensively many quasi-local conserved operators $\mathcal{Z}_{i}$ Serbyn et al. (2013); Huse et al. (2014). This is seen most readily in the infinite disorder limit ($\Delta\rightarrow\infty$) where these are given by the local fields $\mathcal{Z}_{i}=\sigma^{z}_{i}$ and hence act on a single site only. Moreover, they allow to label the eigenstates by the occupation of the $\mathcal{Z}_{i}$. When one moves away from the infinite disorder case, energy space and real space no longer coincide and the above identity becomes more intricate. The Hamiltonian formulated in real space is then diagonalized by a unitary $U_{D}$ which relates energy and real space. The real space representation of the $\mathcal{Z}_{i}$ is given by a decomposition in which all possible Pauli operator combinations can appear

[TABLE]

where $\Sigma^{z},\Sigma^{x}$ are Pauli words consisting of local Pauli-z and Pauli-x operators respectively and $\mu,\nu\in\{0,1\}^{L}$ indicate the position of the Pauli operators in the chain, e.g. $\Sigma^{a}(\mu)=\otimes_{i}(\sigma^{a})^{\mu_{i}}$. We deliberately singled out the weight on $\sigma^{z}_{i}$ corresponding to the infinite disorder conserved operator. All remaining weights are subsumed by the index set $C$ which contains all possible $\mu,\nu$ configurations except the one which would yield $\Sigma^{z}(\mu)\Sigma^{x}(\nu)=\sigma^{z}_{i}$. In practice, $\alpha_{i},\beta_{i}(\mu,\nu)$ can be calculated using the Hilbert-Schmidt scalar product $\beta_{i}(\mu,\nu)=\text{Tr}(U_{D}\mathcal{Z}_{i}U_{D}^{\dagger}\sigma^{z}(\mu)\sigma^{x}(\nu))/2^{L}$. This representation might seem cumbersome at first glance, but also allows one to understand how a transition between strictly local constants of motion in the infinite disorder case to non-local constants of motion in the ergodic regime can be captured formally in terms of the structure of the constants of motion. In fact this is nothing but a formalization of the ”dressing” process Serbyn et al. (2013); Huse et al. (2014). The above decomposition is completely general in the sense that the constants of motion of any system of qubits can be expanded as in Eq. (2) and can hence be applied to the localized as well as to ergodic phase of an MBL system.

The behavior of the conserved operators can be tracked by analyzing the statistics of $\alpha_{i}$. In the localized regime the $\mathcal{Z}_{i}$ are expected to be quasi-local such that the corresponding weights need to decay strongly for operators that have support outside of the localization length, an intuition that has been confirmed numerically Chandran et al. (2015); He and Lu (2018); Rademaker and Ortuño (2016); Mierzejewski et al. (2018); Thomson and Schiró (2018); Kulshreshtha et al. (2018); Goihl et al. (2018, 2019). The largest weight will hence still be given by $\alpha_{i}\sim 1$ and all other weights should be significantly smaller. In the ergodic regime, however, $\mathcal{Z}_{i}$ are not local at all causing the weights to smear out over many different operators and therefore $\alpha_{i}$ will be essentially random and small.

To access the statistics of $\alpha_{i}$, our main numerical tool is constituted by local magnetizations of eigenstates, i.e.$\left\langle E\right|\sigma^{z}_{i}\left|E\right\rangle$ Pal and Huse (2010); Luitz et al. (2017). As the eigenstates can be expressed in terms of the projectors $(\mathds{1}\pm\mathcal{Z}_{i})/2$ onto the (un)occupied sectors of the qLCOM, calculating the expectation value with $\sigma^{z}_{i}$ yields $\pm\alpha_{i}$ depending on the occupation. Therefore, we analyze the histograms over all eigenstates, expecting two very distinct regimes. A localized phase is expected to exhibit a bimodal distribution peaked at $\pm 1$. This holds whenever the constants of motion are quasi-local and hence most of the operator weight is still on the onsite magnetization implying $\alpha_{i}\sim 1$.

For ergodic systems however, the distribution of the $\alpha_{i}$ values should feature a close-to-Gaussian shape with zero mean. Here, the constant of motion will be spread out over many different operators and therefore the weight on $\sigma^{z}_{i}$ is essentially random. Note that in principle, one could also check the other weights $\beta_{i}(\mu,\nu)$ which certainly yield further insights into the decomposition of the qLCOMs. There are however $4^{L}-1$ many of these and without prior knowledge about which one should be sampled, this task is computationally infeasible.

IV Results

In this section, we show and discuss the obtained results. We worked with system sizes $L\in\{8,10,12,14,16\}$ and various disorder strengths. Each point is an average over at least 2000 realizations. Errors have been calculated either using the standard deviation or bootstrapping which amounts to resampling from the obtained data to obtain a distribution of the quantity of interest. In all plots, the error bars are smaller than the symbols used. We either set three or $L/2$ sites to be thermal to cover the system size independent and dependent case. We first discuss the case of the constant size thermal region.

IV.1 Constant size thermal region

In the model considered in this section, we set $s=3$, i.e., three sites are thermal independent of the system size. The two plots on the left in Fig. 1 show the histogram of the local expectation values of magnetization operators $\left\langle E\right|\sigma^{z}_{i}\left|E\right\rangle$ for different sites $i$ encoded in color. We use a system of size $L=16$. As pointed out above, this amounts to sampling the weight of the constants of motion on the $\sigma^{z}_{i}$ operator and hence is an indirect measure for their locality. For low disorder $\Delta=3$, the distributions show a distinct peak at zero and feature strongly decaying tails consistent with the predictions for the thermal region. When comparing the thermal sites to the ones that experience the full disorder strength, we find that the decay for sites inside the thermal region is stronger – as one would intuitively expect. This implies that the constants of motion are not very dominantly supported on the respective $\sigma^{z}_{i}$ operator but presumably rather spread over all remaining operators and hence non-local.

The picture changes drastically with increasing disorder. For $\Delta=9$, the distributions originating from sites in the disordered part of the chain show a distinct bimodality. This shows that the $\sigma^{z}_{i}$ are very close the exact constants of motion as they inherit their occupation statistics which is a sign of quasi-locality. For sites in the thermal region, we find that their distributions feature a peak around zero and small peaks in their heavy tails as well. Since the disorder that these sites experience directly is of order $\epsilon=10^{-6}$, the peaks and the heavy tails should be ascribed to the proximity to the localized chain.

We will now turn to a system size scaling of the observed locality behavior of the constants of motion. Since the standard deviation measures the width of a distribution, it is a genuine measure to compare the distributions for different system sizes Luitz et al. (2017). This is shown in the two plots on the right in Fig. 1, where we plot the standard deviation for all available system sizes encoded in color as a function of disorder strength in the Heisenberg model with three thermal sites. As a guide to the eye, the ergodic region is hatched in red. For both, weak and strong disorder, we find in Fig. 1 that the distributions of local expectation values of energy eigenstates are narrower for sites in the thermal region than in the disordered one. In the ergodic regime for $\Delta=3$, we see that in both regions the standard deviation decreases on all sites with increasing system size. This is an indicator that the system tends towards thermalization as even in the disordered part of the system the distributions become more narrow upon increasing the system size, which is to be expected for $\Delta<\Delta_{c}$ However, for $\Delta=9$, we find that the standard deviation shows no systematic shift with increasing system size, instead a saturation seems to be the most accurate description. This suggests that localization is not compromised by a constant size thermal region.

To detail this observation, we show the standard deviation of a site in the middle of the disordered part of the system at $i=L-\lfloor(L-3)/2\rfloor$ as a function of disorder strength in Fig. 2. In the ergodic regime $\Delta<\Delta_{c}$, the standard deviation of the middle site decays with increasing system size. This implies that the constants of motion become less local when increasing the system size. Upon increasing the disorder strength, we observe a crossing behavior of the standard deviation for disorder values higher than the critical disorder $\Delta>\Delta_{c}$. This shows that in the localized phase there is a(n admittedly weak) tendency towards higher standard deviation for increasing system sizes indicating that the model is driven towards localization upon increasing system size in the sense that the constants motion become more and more localized. These results indicate, that the presence of a fixed size region of low disorder acting as ergodic grain does not alter the qualitative behavior of the full system — depending on the disorder present we find it to be in the ergodic or localized phase without a shift of the transition point compared to the standard disordered Heisenberg chain.

In the next section, we will increase the size of the thermal region with the system size and carry out the same analysis.

IV.2 Constant fraction thermal region

Let us now consider the case $s=L/2$ such that the thermal region covers half of the system and is thus extensive. The two plots on the left in Fig. 3 show the histogram of the local expectation values of magnetization operators $\left\langle E\right|\sigma^{z}_{i}\left|E\right\rangle$ for different sites $i$ in a system of size $L=16$. For low disorder $\Delta=3$ the distribution essentially has the same shape as for the three site thermal region with the only exception being that now the eight thermal sites show the strong decay and the eight disordered sites decay less strongly. Again, this implies that the constants of motion are not very dominantly supported on the respective $\sigma^{z}_{i}$ operator. For $\Delta=9$ the distributions originating from sites in the disordered part of the chain show the expected bimodality but furthermore also feature a peak centered around zero indicating the proximity of the thermal region. For the distribution of the thermal sites, we also notice the difference to the fix size setting discussed above that they do not show a bimodality anymore but only heavier tails than in the case of low disorder.

To investigate these effects in a system size scaling, we again analyze the standard deviation of the distributions in the two plots on the right of Fig. 3. Here, we find that the two regimes show the same qualitative behavior namely a decaying standard deviation with increasing system size. As this is also the case for both parts of the chain, it strongly suggests that the system tends towards narrower distributions and hence delocalization upon increasing the system size.

Let us emphasize this last observation by showing the standard deviation of a site in the middle of the disordered part of the chain at $i=L-\lfloor L/4\rfloor$ as a function of disorder strength in Fig. 4. We find that the gap between standard deviations for different system sizes is diminished at higher disorder but there is no crossover as in the case of a constant size thermal region. For finite systems this implies that significantly stronger disorder scaling with the system size is needed in order to localize even the disordered part of the system. Furthermore, the extrapolation from the available system sizes leads to the conclusion that localization vanishes in the large system limit for this model as constants of motion become more non-local even in the regime of high disorder.

Let us summarize the results. We find that constant thermal regions which are independent of the system size seemingly can not hinder localization and their influence is suppressed by increasing the system size. An ergodic region that scales with the system size, however, changes the physics of the model and apparently delocalizes the system. It is important to note that an ergodic bubble of the order of the system size is exponentially unlikely to appear. It is, however, unclear, if already a weaker scaling of the size of the thermal region with the system size would be sufficient to delocalize the system which is not reasonable to test with the system sizes available to exact diagonalization based schemes.

Let us relate our findings to results obtained in Ref. Luitz et al. (2017) where the authors predict and numerically show a parameter region in which a localized chain can be thermalized by a constant thermal bath. There the authors work in the effective description of MBL and present theoretical arguments for a possible mechanism of instability of MBL phases. Note however, that the disorder strength in their model is only implicitly defined. In Ref. Luitz et al. (2017) the authors use an effective model of the disordered part of the system in terms of quasi-local constants of motions. The disorder strength as used in our work here controls the locality of these constants of motion and by this changes the decay length of the coupling of the bath to different qLCOMs. Upon changing this decay length the authors of Ref. Luitz et al. (2017) identify a critical value whereupon a small thermal grain thermalizes the full system. Based on the fact that we only find such an instability for thermal regions which scale with the system size, we present two explanations for this apparent contradiction.

The first one being that the two Hamiltonians do not exhibit similar physical signatures. Both display features of MBL, but the structure of the bath and coupling to the disordered part of the system might be incompatible. On the level of effective models and not the Hamiltonians as such, it may be that the two effective models in terms of qLCOMs may be incompatible. In Ref. Luitz et al. (2017) the authors employ off-diagonal terms which couple bath and MBL chain; we instead suggest that the important ingredient governing the localization effects is solely the support of the qLCOMs and hence the diagonalizing unitary $U_{D}$. In the setting, in which the bath is constant, the unitary $U_{D}$ should still qualify to be quasi-local with an increased localization length inside the bath, whereas for the extensive bath the unitary will not be local at all in sufficiently large systems.

The second explanation is that the critical coupling decay length at the transition found in Ref. Luitz et al. (2017) might actually correspond to the critical disorder strength separating the ergodic and MBL phase, giving rise to the detected instability. Since no comprehensive theory for this transition exists as well, it might be fruitful to combine our findings with the arguments of Ref.  Luitz et al. (2017) in order to possibly establish a better understanding for the ergodic to MBL transition.

V Conclusions

In this work, we have investigated the fate of the qLCOMs in the disordered Heisenberg model in the presence of a small thermal region modelled by improbably low disorder. We have examined the influence of these regions on the localization behavior. As such regions occur with high probability in large systems this allows us insights into the stability of the MBL phase in the large system limit. As a measure of the locality of the qLCOMs we have employed the expectation value of local magnetization operators obtained from exact diagonalization. In the system sizes accessible to us, we find that the qLCOMs and thereby the MBL phase is stable when coupled to a finite bath. If the thermal region is allowed to scale with the system size, however, our findings suggest that MBL will vanish when approaching the thermodynamic limit.

These results suggest that an isolated MBL system is stable upon increasing the system size. If, in contrast, coupled to an external bath, the question of the stability of MBL may depend on more subtle details as the precise size and shape of the bath may have a strong influence. An interesting further research direction would be to carry out a similar analysis for two dimensional systems and devise local probes for delocalization which could then be used in state of the art experimental realizations of MBL.

VI Acknowledgements

MG is grateful for the feedback on the manuscript by Nicolas Tarantino as well as numerous discussion with the participants of the Anderson Localization and Interactions workshop at the MPIPKS – specifically with Wojciech De Roeck, Alan Morningstar, Anna Goremykina and Maksym Serbyn as well as Antonio Rubio Abadal and Jun Rui at MPQ. Moreover, we would like to thank Henrik Wilming for many discussions at earlier stages of this work. This work has been supported by the ERC (TAQ), the DFG (FOR 2724, CRC 183, EI 519/14-1), and the Templeton Foundation. This work has also received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 817482 (PASQuanS).