Slow dynamics and strong finite-size effects in many-body localization with random and quasi-periodic potential

TL;DR

This study explores slow dynamics and finite-size effects in many-body localization within disordered and quasi-periodic quantum wires, proposing a two-phase scenario with distinct subphases and challenging existing Griffiths effect explanations.

Contribution

It introduces a novel two-subphase framework for MBL, based on the divergence behavior of the density correlation width, and discusses finite-size effects and transition characteristics.

Findings

Finite-size effects limit the growth of the correlation width.

Evidence suggests a two-subphase structure within the MBL regime.

Transition length scale diverges with an essential singularity.

Abstract

We investigate charge relaxation in disordered and quasi-periodic quantum-wires of spin-less fermions ($t{-}V$-model) at different inhomogeneity strength $W$ in the localized and nearly-localized regime. Our observable is the time-dependent density correlation function, $\Phi(x,t)$, at infinite temperature. We find that disordered and quasi-periodic models behave qualitatively similar: Although even at longest observation times the width $\Delta x(t)$ of $\Phi(x,t)$ does not exceed significantly the non-interacting localization length, $\xi_0$, strong finite-size effects are encountered. Our findings appear difficult to reconcile with the rare-region physics (Griffiths effects) that often is invoked as an explanation for the slow dynamics observed by us and earlier computational studies. As a relatively reliable indicator for the boundary towards the many-body localized (MBL) regime even under these conditions, we consider the exponent function $\beta(t) {=} d\ln \Delta x(t) / d\ln t$. Motivated by our numerical data for $\beta$, we discuss a scenario in which the MBL-phase splits into two subphases: in MBL$_\text{A}$ $\Delta x(t)$ diverges slower than any power, while it converges towards a finite value in MBL$_\text{B}$. Within the scenario the transition between MBL$_\text{A}$ and the ergodic phase is characterized by a length scale, $\xi$, that exhibits an essential singularity $\ln \xi \sim 1/|W-W_\text{c}|$. Relations to earlier numerics and proposals of two-phase scenarios will be discussed.

Figures

b of this Ref. Khemani *et al.*, [2017b](#bib.bib52). A pair of traces of $r_{n}(W,L)$ for $L{=}12,14$ intersects near $W{\approx}5.5$ while the pair $L{=}16,18$ intersects near $W{\approx}7$. The true scaling limit would announce itself by all pairs intersecting at the same point $W{=}W_{\text{c}}$.

Equations

$$\beta (t)=\frac{d\ln \Delta x(t)}{d\ln t}.$$

$$\hat{\mathcal{H}}=$$

$$+V\sum _{x=1}^{L-1}\left({\hat{n}}_x-\frac{1}{2}\right)\left({\hat{n}}_{x+1}-\frac{1}{2}\right),$$

$$\Phi (x,t)$$

$$\beta [\Delta x]=\frac{\mathrm{b}}{1+{\mathrm{b}}_1\mathcal{K}({\xi }_0/\Delta x)}.$$

$$\mathcal{K}(y)\approx \tilde{\mathcal{K}}(y{)}^{W-W_{{\text{c}}_1}},\tilde{\mathcal{K}}(y):=e^{{\frac{d\ln \mathcal{K}(y)}{dW}|}_{W_{{\text{c}}_1}}}.$$

$$\xi \sim {\xi }_0{\mathrm{b}}_1^{\frac{1}{\kappa (W_{{\text{c}}_1}-W)}}.$$

Full text

Many body localization

I Introduction

The interplay of quantum interference and interactions exhibits some of the most fascinating phenomena to be encountered in many-body systems. A typical representative of concern to us here is many-body localization Basko et al. (2006); Gornyi et al. (2005) (MBL), which has attracted considerable attention over the last decade Žnidarič et al. (2008); Pal and Huse (2010); Bardarson et al. (2012); Nandkishore and Huse (2015); Altman and Vosk (2015); Luitz et al. (2015); Bera et al. (2015); Žnidarič et al. (2016); Ann. Phys. (2017) (Berlin); Alet and Laflorencie (2018); Abanin et al. (2018). At its heart is the basic observation that Anderson localization prevalent in low-dimensional, disordered systems of non-interacting fermions can be robust against interaction effects at finite temperature despite of dephasing as is illustrated in Fig. 1.

Starting from the MBL phase, with decreasing disorder a quantum-critical point is reached separating the MBL region from an ergodic one. Much progress has been made in recent years with respect to experimental investigations of the corresponding localization-delocalization transition in diverse setups including cold atoms Schreiber et al. (2015); Choi et al. (2016); Lüschen et al. (2017); Bordia et al. (2017); Kohlert et al. (2018); Rispoli et al. (2018), nuclear spin chains Wei et al. (2018), Yb ion system Smith et al. (2016), dipolar spins Kucsko et al. (2018), or in InO thin films Ovadia et al. (2015). Importantly, in several experiments a very pronounced slow-down of the relaxation dynamics has indeed been reported, when the (quasi-)disorder strength exceeds a threshold value Lüschen et al. (2017); Bordia et al. (2017); Rispoli et al. (2018).

The theoretical investigation of MBL and the related effects poses great challenges. On the one hand side the phenomenon escapes controlled analytical investigations. On the other hand, the computationally available window of observation times, $t_{\text{obs}}$, and system sizes often is too small in order to predict with confidence the thermodynamic, long time asymptotics. In fact as shown recently, due to dynamical slowing down, even in computational studies of the ergodic regime the width of the diffusion propagator, ${\Delta x}(t)$, is hardly ever seen to exceed the non-interacting localization length, $\xi_{0}$, significantly Bera et al. (2017). The resulting difficulties to safely extrapolate plague virtually all computational studies of dynamical behavior near and in the MBL phase. In fact, even present days experimental studies do not (yet) allow to go very deeply into the asymptotic regime.

It is noteworthy that in computations in addition to achieving sufficiently large system sizes and long observation times other difficulties arise; we list two prominent examples: (i) Computer simulations carry nonphysical numerical parameters that need to be converged, such as the bond dimension, time increments and number of Chebyshev polynomials. Convergence can be difficult to achieve and will, in general, be sensitive to the calculated observable and the physical parameter regime. (ii) Quenches may produce transients or even asymptotic time traces that reflect properties of the initializing wavefunction more than the thermodynamic equilibrium. In this case the results of studies using different initialization, e.g. random versus staggered states, may be difficult to compare.111For example, Luitz et al. (2016) report values for the imbalance exponent $\zeta$ that are smaller than the values detected by Doggen et al. (2018) by a factor of two or more, e.g., near $W{=}2$. Both works have tested convergence with respect to numerical parameters, while they differ in the choice of the initial state: random Luitz et al. (2016) versus Néel Doggen et al. (2018).

Summarizing, a possibility to get stuck with transient or untypical behavior has to be accounted for. Nonwithstanding these difficulties, from combining computational studies employing a variety of numerical techniques and analytical arguments based, e.g., on phenomenological renormalization-group treatments, a certain picture has emerged. It has been surveyed in several recent reviews Abanin et al. (2018); Vosk et al. (2015); Potter et al. (2015); Agarwal et al. (2017); Luitz and Lev (2017); Prelovšek et al. (2017); Alet and Laflorencie (2018).

In this work we present a computational study of the MBL region employing the $t{-}V$-model with fully random and quasi-periodic (Andre-Aubry-type) potentials Iyer et al. (2013); Lee et al. (2017); Khemani et al. (2017a); Setiawan et al. (2017); Doggen and Mirlin (2019). We invoke as observable the density-density correlation function, $\Phi(x,t)$. As a first characterization its variance, $\Delta x(t)$, and the associated exponent function

[TABLE]

is being analyzed. Even though the diffusion propagator $\Phi(x,t)$ has been studied before, Bar Lev et al. (2015); Lev and Reichman (2016); Luitz and Lev (2017); Prelovšek et al. (2017) it is not the usual object of investigation. When the exponent $\beta$ has been considered, it has mostly been derived from closely related observables, such as the frequency dependent conductivity at zero wavenumber calculated in finite size systems Agarwal et al. (2015); Gopalakrishnan et al. (2015); Khait et al. (2016); Steinigeweg et al. (2016); Barišić et al. (2016). The density response at wavenumber $q{=}\pi$ (imbalance) and the corresponding dynamical exponent $\tilde{\beta}$ has been observed more frequently Luitz and Lev (2017); Doggen et al. (2018), though its relation to $\beta$ is not well understood, in general. Also the exponent of the return probability, $\alpha$, has been employed, although the relation to $\beta$ involves an extra exponent, which is not known with good accuracy Bera et al. (2017).

The interest of previous authors has been in the long time limit, $\beta_{\infty}(W){\equiv}\beta(t{\to}\infty)$, or closely related exponents, where $\beta$ can serve as an effective “order parameter”; in the ergodic region prevailing at disorder $W{<}W_{\text{c}}$ one has $\beta_{\infty}{>}0$, while the MBL-phase, $W{>}W_{\text{c}}$, has been defined by $\beta_{\infty}{=}0$. However, the data we here report leaves open a possibility for an intermediate window, $W_{\text{c}_{1}}{<}W{<}W_{\text{c}_{2}}$, which constitutes within the MBL-phase a subphase MBL${}_{\text{A}}$ ; It exhibits a width $\Delta x(t)$ diverging in time but with a growth weaker than any power, see Fig. 2 for a sketch and Fig. 3 for illustrating data. At very large disorder, $W{>}W_{\text{c}_{2}}$, one enters the subphase, MBL${}_{\text{B}}$ , in which $\beta(t){\to}0$ upon ${\Delta x}(t)$ approaching the localization length $\xi(W)$.

In view of the difficulty accessing the asymptotic regime of $\beta(t)$, we consider as a pragmatic indicator in addition to $\beta$ also its slope, $\dot{\beta}$, taken at the computationally still accessible (i.e. relatively short) times $t_{\xi}$, where $\Delta x(t)$ has passed the largest relevant microscopic length. 222We include $\dot{\beta}(t)$ as an indicator with the idea that $t_{\xi}$ is the only relevant time scale. Correspondingly, at $t\gtrsim t_{\xi}$ the sign of $\beta$ is converged, while it numerical value may not be, yet. Following this idea, entering the MBL${}_{\text{A}}$ -phase from the ergodic side at $W{=}W_{\text{c}_{1}}$ is signalized by a change of sign in $\dot{\beta}(t)$ at $t\gtrsim t_{\xi}$. While the ergodic phase is characterized by $\beta(t){>}0,\dot{\beta}(t){>}0$, we have $\beta(t){>}0,\dot{\beta}(t){<}0$ in MBL${}_{\text{A}}$ . A natural candidate for the microscopic length appearing in $t_{\xi}$ is the localization length, $\xi_{0}(W)$, of the non-interacting reference system, possibly having undergone an interaction-mediated (static) renormalization. Hence, we define $t_{\xi}$ implicitly via $\Delta x(t_{\xi})\approx\xi_{0}.$

We list the salient aspects of our numerics and the scenario displayed schematically in Fig. 2: (i) In the ergodic phase, $\Delta x(t)\sim t^{\beta}$, where $\beta(t)$ converges to a constant with $\dot{\beta}(t)$ approaching zero from above; correspondingly, the slope $\dot{\beta}(t_{\xi})$ seen in the numerical data should be positive when converged with respect to computational parameters. In MBL${}_{\text{A}}$ the asymptotic growth of $\Delta x(t)$ is also unbound, but $\beta(t){\to}0$ with $\dot{\beta}(t)$ approaching zero from below; correspondingly, $\dot{\beta}(t_{\xi})$ seen in the numerical data is negative when converged. Finally, in MBL${}_{\text{B}}$ $\Delta x(t)$ is bounded from above and $\beta$ extrapolates to zero at a finite ${\Delta x}$. We interpret the numerical data shown in Fig. 4 as reflecting the evolution with increasing disorder strength $W$ described here and in Fig. 2. The existence of an intermediate phase has been proposed beforeDumitrescu et al. (2019); the relation to our scenario is discussed below.

(ii) We have applied our analysis method to the $t{-}V$ model with fully disordered and quasi-periodic potential. Importantly, comparing both cases we cannot detect an essential difference in the evolution of ${\Delta x}(t)$ with the inhomogeneity strength. Therefore, our results seem hard to reconcile with Griffiths-type scenarios that rely on rare-regions to explain the origin of slow dynamics Vosk et al. (2015); Agarwal et al. (2017); Luitz and Lev (2017); Lev et al. (2017).

(iii) To formulate our scenario we adopt the spirit of mode-coupling theory extrapolating into the asymptotic regime. It suggests that the correlation length $\xi$, when entering MBL${}_{\text{A}}$ from the ergodic phase, exhibits an essential singularity rather than a power-law, which is consistent with the Harris-Chayes criterion Chayes et al. (1986).

II Model and Computational Method

We consider the $t{-}V$-model

[TABLE]

where hopping is over $L$ sites with amplitude $t_{\text{h}}$. We work at a half filling and with hard-wall boundary conditions. When employing disordered potentials, the onsite-energies, $\epsilon_{x}$, are drawn from a box distribution $[-W,W]$. 333 When working with quasi-periodic potential (Aubry-André model) we choose, $\epsilon_{x}{=}W\cos(2\pi\alpha x+\phi)$ with $\alpha=2/(\sqrt{5}+1)$ and $\phi$ is the random phase distributed uniformly between $[0,2\pi]$ In subsequent calculations the interaction is taken twice larger, $V{=}2t_{\text{h}}$, than in previous studies Bera et al. (2017) to emphasize interaction effects. At the disorder values shown in Fig. 3, the spectral statistics is close to Poissonian and in this sense the system is many-body localized, see Sec. IV.3.

Our observable is the density-density correlator,

[TABLE]

with primary focus on the second moment, $\Delta x(t)=\left[\sum_{x=-L/2}^{L/2}x^{2}\Phi(x,t)-(\sum_{x=-L/2}^{L/2}x\Phi(x,t))^{2}\right]^{1/2}.$ The disorder average is indicated by an overline; the thermal average is performed at infinite temperature and denoted by angular brackets, $\langle\ldots\rangle$; corresponding traces are performed stochastically. 444Here, we have used that in the infinite temperature limit the density is homogeneous and equal to $1/2$. For the time evolution, Eq. (3), we employ a standard Chebyshev-polynomial propagation Weiße et al. (2006); further details are given in Ref. Bera et al., 2017.

III Results.

III.1 Simulations

We systematically analyze traces of the kind shown in Fig. 3 in terms of the logarithmic derivative (1); the results are shown in Fig. 4. The top-row displays the evolution of $\Phi(x,t)$ with system size and disorder at $t{\mathrel{\mathop{:}}=}t_{\text{obs}}{\sim}400t^{-1}_{\text{h}}$. Fully consistent with our earlier findings, $\Phi(x,t)$ exhibits a pronounced non-Gaussian, nearly exponential shape – here seen in a wide range of $W$-values. Not surprisingly, the interaction effects on $\Phi(x,t)$ are strongest in the tails; they emerge gradually, with ever increasing system size $L$ and observation time. Note that at larger disorder the interacting and non-interacting traces are seen to almost coincide for smaller distances at the time observed, see, e.g., Fig. 4(d), $W{>}5$. At these (possibly transient) times one has ${\Delta x}(t)\approx{\Delta x}^{(0)}$. Therefore, interaction-related static renormalization effects of the disorder-strength are small at $W\gtrsim V$.

The lower row of Fig. 4(e-h) is the basis of the schematic plot Fig. 2. It shows the exponent function $\beta(t)$ that highlights the important effects related to finite $L$. Consistent with our previous findings Bera et al. (2017), at weak disorder $W{\lesssim}3$ and ${\Delta x}(t){\sim}{\Delta x}_{0}$ there is a system size convergence with respect to the sign of the slope of $\beta(t)$; there is no notion of convergence at ${\Delta x}(t){\gtrsim}{\Delta x}_{0}$ beyond the statement that $\beta({\Delta x})$ is curved to the right. At larger disorder, $W{\gtrsim}5$, a qualitative change is observed: the curvature seen in the inset of Fig. 4(g) is about to turn left at ${\Delta x}(t){>}{\Delta x}_{0}$. The MBL${}_{\text{A}}$ scenario advocated in Fig. 2 is based on the observation that the curvature of $\beta({\Delta x})$ changes sign with growing $W$, so what is a right-curvature in Fig. 4 turns into a left-curvature at later times; also we assume that there is no turn-around in $\beta({\Delta x})$ at even longer times, but rather a smooth decay towards the horizontal axis at ${\Delta x}{\to}\infty$. (A discussion of the pronounced short-time feature of the data Fig. 4(e-h) lower row, namely the spike growing in $\beta({\Delta x})$ with increasing $W$ near ${\Delta x}_{0}$, is given in the appendix.)

III.2 Phenomenological modeling.

The data for the $\beta$-function displayed in Fig. 4 condenses into the phenomenological ansatz

[TABLE]

Similar to our previous work Bera et al. (2017), we here consider $\beta$ as a functional of the slowest dynamical mode, which for simplicity we parametrize with ${\Delta x}(t)$; the parameters $\mathrm{b},\mathrm{b}_{1}$ and the kernel $\mathscr{K}$ are evolving with disorder $W$. The ansatz helps formulating a scenario for the qualitative properties of the long-time limit, which are just barely seen to emerge in Fig. 4.

Embarking on (4) we describe a passage from the ergodic regime through MBL${}_{\text{A}}$ into MBL${}_{\text{B}}$ with increasing $W$. The kernel $\mathscr{K}(y)$ equals unity at $y{=}1$ and otherwise is assumed to be a monotonous function in the interval $1{\gtrsim}y{\gtrsim}0$. Its precise form is evolving with the disorder strength, $W$: With decreasing $y$, the kernel decreases in the ergodic phase, increases in MBL${}_{\text{A}}$ and is seen to diverge at a finite $y$ upon entering MBL${}_{\text{B}}$ .

Ergodic regime $W{<}W_{\text{c}_{1}}$. The power-law growth of $\Delta x(t)$ characteristic for the ergodic regime is captured by (4) with $\mathrm{b}{>}0,\mathrm{b}_{1}{>}0$ and taking $\mathscr{K}(y)$ as decreasing to zero in the limit $y\to 0$. Namely, with the variance $\Delta x(t)$ growing in time, the kernel $\mathscr{K}(\xi_{0}/\Delta x(t))$ vanishes and so the denominator converges to unity. At this point, $\mathrm{b}$ should be identified with the diffusion exponent: $\mathrm{b}=\beta_{\infty}$

Transition into MBL${}_{\text{A}}$ . The transition is associated with a vanishing slope, $\dot{\beta}{=}0$. Within in the setting (4) this requires the existence of a particular value of the disorder $W_{\text{c}_{1}}$, such that the contribution of the term containing $\mathrm{b}_{1}\mathscr{K}(\lambda/\Delta x)$ is independent of $\Delta x$. This situation could be realized in several ways; the most plausible scenario is $\left.\mathscr{K}(y)\right|_{W_{\text{c}_{1}}}=1$. 555 Alternative possibilities could be (i) $\mathrm{b}_{1}(W_{\text{c}_{1}}){=}0$ or (ii) $\mathrm{b}_{1}(W_{\text{c}_{1}}){\to}\infty$. However, in our scenario we do not follow these directions.

We thus parameterize $\mathscr{K}(y)$ in the vicinity of $W\approx W_{\text{c}_{1}}$ as follows:

[TABLE]

This shape has consequences for the correlation lengths, $1{\approx}\mathrm{b}_{1}\ \mathscr{K}(\xi_{0}/\xi).$ which we here define implicitly; we have $\tilde{\mathscr{K}}(\xi_{0}/\xi)\approx\mathrm{b}_{1}^{\frac{1}{(W_{\text{c}_{1}}-W)}}.$ To illustrate the impact of this statement, we adopt a power-law shape $\tilde{\mathscr{K}}(\xi_{0}/\xi)\sim(\xi/\xi_{0})^{\kappa}$, $\kappa{>}0$, so that

[TABLE]

It is thus seen that the scenario formulated in (4) tends to generate length scales that are diverging at the transition in an exponential rather than in a power-law fashion.

MBL${}_{\text{A}}$ -phase. According to (5), in MBL${}_{\text{A}}$ the kernel $\mathscr{K}(y)$ diverges with $y{\to}0$ keeping $\mathrm{b}{>}0$ and $\mathrm{b}_{1}{>}0$. With growing time, the second term in the denominator grows without bound and eventually $\beta(\Delta x)\approx\frac{\mathrm{b}}{\mathrm{b}_{1}}\ \mathscr{K}\left(\xi_{0}/\Delta x\right)^{-1}.$ The behavior is consistent with the data Fig. 4 as seen, e.g., at $W{=}5$. It predicts for the $\beta$-function a decay approaching zero upon ${\Delta x}(t){\to}\infty$.

Transition into the MBL${}_{\text{B}}$ -phase. Upon approaching $W_{\text{c}_{2}}$ and MBL${}_{\text{B}}$ the kernel $\mathscr{K}(y)$ diverges at a nonvanishing value $y_{c}$. For illustration we adopt the generic form $\mathscr{K}(y)=\frac{y^{\alpha}}{1+\mathrm{b}_{2}y^{1/\nu}}$ with $\alpha,\nu>0$. In MBL${}_{\text{A}}$ , we have $\mathrm{b}_{2}{>}0$ and so the kernel diverges in a power-law fashion, $\mathscr{K}(y)\sim y^{\alpha-1/\nu}$, as required. When $W$ crosses $W_{\text{c}_{2}}$ from below, the coefficient $\mathrm{b}_{2}(W)$ changes sign from positive to negative, so in MBL${}_{\text{B}}$ the kernel diverges at $y^{1/\nu}=|\mathrm{b}_{2}|^{-1}.$ we conclude that the localization length diverges $\xi\approx\lambda|W_{\text{c}_{2}}-W|^{-\nu}$.

III.3 Discussion

The model (4) has implications, which can be tested further. We mention the following two: (i) When entering MBL${}_{\text{A}}$ from the ergodic side, $\beta_{\infty}(W_{\text{c}_{1}})$ is taking a non-vanishing, critical value. In contrast, at the critical point separating MBL${}_{\text{B}}$ and MBL${}_{\text{A}}$ $\beta_{\infty}(W_{\text{c}_{2}}){=}0$ is expected. (ii) Previous computational studies have attempted a scaling analysis near $W{\approx}W_{\text{c}_{1}}$ in order to extract a localization length exponent based on a power-law dependency $\xi{\sim}|W-W_{\text{c}_{1}}|^{-\nu}$. The value of $\nu$ has not consolidated, yet: Luitz et al. (2015) obtained $\nu{\approx}0.7$ from the system size scaling of the magnetization density, violating the Harris-Chayes bound $\nu{>}1$ in one dimension Chayes et al. (1986). This led to a subsequent discussion Khemani et al. (2017b) attempting at rationalizing this result. On the other hand, recent exact diagonalization studies of the scaling of the Schmidt gap Gray et al. (2018) or local temperature fluctuations Lenarčič et al. (2018) found exponents exceeding $2$, which is consistent with the Harris bound. We point out that within the scenario here lined out the transition that has been investigated in these studies would corresponds to the entry into MBL${}_{\text{A}}$ , so the correlation length diverges stronger than any power. Hence, the Harris-Chayes bound is satisfied automatically.

Saying this we emphasize once again the exploratory character of the scenario developed with Eq. (4). Due to creep, we do not see at present a reliable way of telling, whether the system actually localizes or not. We stress that creep is a phenomenon that gradually emerges at largest system sizes, longest time scales and only in the tails of the distribution (4). For observing or ruling out creep with frequently employed observables, such as the level statistics, the return probability or the sublattice imbalance, deviations have to be identified from the non-interacting background-behavior; such deviations manifest as small signals towards delocalization that slowly grow with increasing time and system size and thus are challenging to consolidate.

IV Relation to earlier studies

IV.1 RG-studies.

We compare our scenario with recent renormalization group (RG) calculations Morningstar and Huse (2019); Dumitrescu et al. (2019); Goremykina et al. (2019); Zhang et al. (2016); Thiery et al. (2018) and a numerical study Herviou et al. (2018). DumitrescuDumitrescu et al. (2019) et al. foresee a Kosterlitz-Thouless-type scenario that also exhibits an intermediate MBL phase; also in this case an exponentially diverging length scale was predicted at the transition into the ergodic phase  (implying $\nu{=}\infty$), however with an exponent $|W-W_{\text{c}_{1}}|^{-1/2}$ (rather than $|W-W_{\text{c}_{1}}|^{-1}$). This intermediate phase is characterized by power-law distributed rare thermal blocks with an exponent dependent on disorder strength, taking a universal value only at the critical point Goremykina et al. (2019); Dumitrescu et al. (2019); Herviou et al. (2018); at strong disorder the distribution becomes exponential.

In most recent workMorningstar and Huse (2019) a modification of the RG scheme underlying Ref. Goremykina et al., 2019 was proposed, introducing a second length scale. The intermediate phase of the original model appears to be unstable under this modification and, in fact, a power-law distribution occurs only at the critical point. Therefore, from the perspective of phenomenological RGs, the existence of an intermediate phase is still under debate. In order to clarify if creep relates to rare disorder fluctuations, we have been comparing dynamics in random and quasi-periodic potentials in Fig. 2, see also the appendix. Remarkably, we observe a qualitatively similar slow dynamics in both cases. This points towards a possibility that rare events are not the only viable explanation for creep as observed, e.g., in Fig. 3.

Our assertion of slow dynamics receives (tentative) support from most recent experiments in optical lattices Rispoli et al. (2018). The Bose-Hubbard model was simulated with quasi-periodic disorder potential; a slow particle transport has been seen at finite system sizes ($L=8,12$) above the critical disorder strength. We would like to interpret this observation as a possible experimental manifestation of a creep-type dynamics.

IV.2 Numerical studies.

The exponent function $\beta(t)$, Eq. (1), was considered before by Bar Lev, Cohen and Reichman.Bar Lev et al. (2015) The fitted exponents together with selected raw data are shown in Fig. 3 of their work. Even though system sizes and observation times are considerably smaller than the ones used in our paper, the onset of creep at large disorder is clearly visible in their raw data, e.g. at $W{\approx}7$.

One would expect that manifestations of the “creep phenomenon” exist in the time-dependency of generic correlation functions and we believe that this indeed is the case. For instance, Doggen et al. Doggen et al. (2018) have simulated the quench of an initial state with maximum sublattice imbalance. In Fig. 2 of their work they observe a power-law decay of the imbalance with an exponent $\tilde{\beta}(W)$. Specifically, they report an exponential decrease that roughly parametrizes as $\tilde{\beta}(W)\approx\tilde{\beta}(2)\exp(-(W-2))$ in $2<W<5$; at larger $W$ the exponent has become so small that it can no longer be resolved within the error bars. We interpret this result as broadly consistent with our observations; in particular, it implies that $W_{\text{c}}>5$ for the model considered and the localized phase proper, strictly speaking, is outside of the computational observation window. 666With respect to $W_{\text{c}}$ a similar conclusion can also be drawn from the scaling of the Schmidt-gap displayed in Fig. 5 of the same work. Doggen et al. (2018)

Gopalakrishnan et. al. Gopalakrishnan et al. (2015) investigated the low-frequency asymptotics of the conductivity. They reported $\sigma(\omega)\sim|\omega|^{\alpha}$ with an exponent $\alpha$ varying continuously within a region in the MBL-phase. At stronger disorder the agreement reflects the fact that interaction effects are almost invisible in the numerical data, so the non-interacting result $\alpha{\sim}2$ is reached trivially at large $W$ within the available range of computational parameters. At weaker disorder near the border to the thermal phase, Gopalakrishnan et al. (2015) give a smaller value $\alpha{\sim}1$ corresponding to a logarithmic growth of $\Delta x(t)$. This is consistent with the behavior we observe in MBL${}_{\text{A}}$ (see Fig. 2). In fact, a more recent study Steinigeweg et al. (2016) on significantly larger system sizes reports $\alpha\sim 1$ in a broad window of disorder strengths near the transition. We mention that an exponent varying continuously across the transition, $\alpha(W)$, has also been found in a self-consistent perturbation theory in disorder Prelovšek and Herbrych (2017).

Similarly, also Barišić et. al. Barišić et al. (2016) report values $\alpha{\sim}1$ in the MBL-phase, employing system sizes up to $L{=}28$. In their fitting, the authors are allowing for a non-vanishing residual dc-conductivity, however, and indicate that finite size corrections are not fully under control. So, we understand that the numerical value $\alpha{=}1$ is subject to a significant uncertainty. The same group later extended their studies to related correlation functions Mierzejewski et al. (2016); Prelovšek et al. (2017). Reminiscent of creep, also in these correlators a very slow dynamics with significant finite-size effects can be identified, e.g., in Fig. 2 (c,d) of Mierzejewski et al. (2016).

We mention that the sparsity of many-body states in Fock-space can give rise to a kind of multifractality, which exhibits a non-universal and, in general, basis dependent spectrum Atas and Bogomolny (2014); Gornyi et al. (2017). Most recent numerical work has confirmed this earlier result at least within a certain window of system sizes computationally available Macé et al. (2018). If and how multifractality of the many-body wavefunction manifests in the evolution of $\alpha(W)$ remains to be clarified.

Finally, also the study by Serbyn, Papic and Abanin exhibits discernible finite size effects, which may be consistent with creep Serbyn et al. (2015). Specifically, these authors study expectation values of energy-normalized overlap matrix elements of local operators. Fig. 3 of this work shows that with increasing system size, $L$, the flow of the average turns around from insultating and redirects towards thermalizing behavior. The system size of turn-around, $L^{*}(W)$, is seen to grow with disorder. Since in this work only system sizes up to $L{=}16$ have been studied, it remains open if a finite critical value of disorder, $W_{\text{c}}$, exists at which $L^{*}(W)$ diverges, and what numerical value it takes. 777We express our gratitude to Ehud Altman for drawing our attention to this work.

IV.3 Further discussion: Spectral statistics and creep

Level spacing distribution $P(s)$.

We comment on the spectral statistics of the random potential and the AAH model. Figure 5 shows the evolution of the level spacing distribution $P(s)$ from weak (Wigner-Dyson statistics) to strong (Poisson distribution) disorder as it is observed in a relatively small system of size $L{=}16$. The data is qualitatively similar to the one presented, e.g., in Ref. Bertrand and García-García, 2016. It suggests that at disorder values exceeding $W\approx 3$ the $P(s)$ distribution is nearly Poisson and therefore samples are in a certain spectral sense close to the many-body localized phase. As is also seen from the Fig. 5 that statement holds irrespective of whether the potential is fully random or quasi-periodic.

Importantly, the distribution $P(s)$ undergoes an evolution with increasing system size $L$, see e.g. Fig. 7 in Ref. Bertrand and García-García, 2016 and Fig. 9 in Ref. Sierant and Zakrzewski, 2019. At stronger disorder values and with $L$ increasing from small values, e.g. $L{\sim}12$, the initial flow is towards a Poisson shape. This observation is fully consistent with - and in fact an indispensable requirement for - the applicability of the concept of local integrals of motion that underlies the current understanding of many-body localization.

Notice, however, that at present there is no proof that indeed a critical value of the disorder, $W_{\text{c}}$, exists, such that the asymptotic flow of $P(s)$ for the $t{-}V$-model or Heisenberg model is indeed towards Poisson at disorder values exceeding $W_{\text{c}}$. The observation of “creep” indicates that the critical disorder strength, $W_{\text{c}}$, may be very large; if it is finite at all, then it must correspond to a non-interacting localisation length that is considerably smaller than the lattice spacing.

A conclusive numerical study of $P(s)$ to this end is lacking as of today. However, inspecting the data for $P(s)$ shown in Fig. 9 of Ref. Sierant and Zakrzewski, 2019 at very small $s$ and disorder, $W{=}4$, the flow towards Poisson indeed appears to stop already near $L\approx 18$. Conceivably, at larger $W$ a similar termination can be observed, if only at larger system sizes $L(W)$. Therefore, the present evidence concerning $P(s)$ being close to Poisson at smaller system sizes is not inconsistent with the creep-scenario.

Level statistics ratio $r_{n}$.

A simplified diagnostics of the level statistics is provided by the ratio $r_{n}(W,L){\coloneqq}\text{min}(\Delta_{n},\Delta_{n+1})/\text{max}(\Delta_{n},\Delta_{n+1})$, where $\Delta_{n}{\coloneqq}E_{n}-E_{n+1}$ is the spacing between eigenvalues. For the Gaussian-orthogonal ensemble (nearly chaotic situation, $W{\ll}W_{\text{c}}$) one has $r_{n}{\simeq}0.53$ while for the Poisson case (localized situation, $W{\ll}W_{\text{c}}$) $r_{n}{\simeq}0.39$. In analogy to conventional second-order phase transitions, it is tempting to test for single-parameter scaling near the critical point, i.e. $r_{n}(\xi/L)$, defining the localization length, $\xi\propto|W-W_{\text{c}}|^{-\nu}$, with the corresponding exponent $\nu$. Indeed, the earlier attempts to determine $W_{\text{c}}$ and $\nu$ did employ such an analysis Luitz et al. (2015).

It has been pointed out by Khemani et al. Khemani et al. (2017b) already that due to the small system sizes available, – typically not exceeding $L{=}20$, – one is stuck in a pre-asymptotic regime. A clear illustration indicating how large the actual distance to criticality really is can be drawn, e.g., from Fig. 2b of this Ref. Khemani et al., 2017b. A pair of traces of $r_{n}(W,L)$ for $L{=}12,14$ intersects near $W{\approx}5.5$ while the pair $L{=}16,18$ intersects near $W{\approx}7$. The true scaling limit would announce itself by all pairs intersecting at the same point $W{=}W_{\text{c}}$.

Attempting in this situation a scaling analysis and fully ignoring finite size corrections yields $\nu{\approx}1$. The result does not violate the Harris criterion, because the validity of the latter is restricted to the critical point, while the scaling analysis has been done in a preasymptotic regime. Following the same type of analysis, one expects with increasing system sizes a drift of $\nu$ and $W_{c}$ to ever larger values. Again, the existing data seems consistent with the creep scenario.

We mention that a drift in pair-wise intersection points is very common in computational studies of Anderson transitions. For a recent example how to include finite-size corrections in a case relevant to MBL, i.e. random-regular graphs, we refer to a recent study by Ref. Tikhonov et al., 2016.

V Conclusion

In this work we have studied the density autocorrelation function of spinless fermions in the $t-V$-model at infinite temperature. We detect strong finite-size effects and a very slow dynamics - “creep” - with a tendency towards thermalization even at strong disorder, where the level statistics is close to Poissonian. Our observations are very similar for truly random and quasi-periodic potentials, so that rare-region effects seem discouraged as a likely origin of subdiffusion and creep. We believe that our results imply one out of two possibilities:

(i) “Creep” is a transient phenomenon, perhaps indicative of the critical fan Khemani et al. (2017b) (similarly for quasi-periodic potential Setiawan et al. (2017)). In this case, the width of a wavepacket, $\Delta x(t)$, either grows in a power law manner, $\Delta x(t){\sim}t^{\beta}$, in a thermalizing phase, or exhibits a logarithmic (or weaker) growth window in time after which $\Delta x(t)$ eventually saturates at a finite value, $\xi(W)$, in the MBL-phase.

If this is the case, an important conclusion of our study is that due to creep the asymptotic regime, $\Delta x(t)\approx\xi(W)$, in simulations is very difficult to reach in a controlled manner. We find that even at disorder values such that the non-interacting localization length is considerably smaller than the lattice spacing, $\Delta x(t)$ does not saturate within the simulation window of $1000t_{\text{h}}^{-1}$. Hence, a numerical study of the true asymptotics including $\xi(W)$

• especially with respect to criticality - may have been out of reach of simulations of physical observables thus far, and, presumably, even in experiments Abanin et al. (2018).

(ii) An alternative possibility is that “creep” may indicate the existence of marginally thermalizing phases, such as the subphase MBL${}_{\text{A}}$ discussed in this work. The characterizing feature is that $\Delta x(t)$ keeps growing in time weaker than any power, eventually crossing any finite upper bound $\xi$. In this case, the numerical observation of creep-type phenomena would be closer to the asymptotic regime, but the present understanding of the MBL-phase would be incomplete. Conceivably, such phases can be described by complementing the concept of local integrals of motionSerbyn et al. (2013); Huse et al. (2014), which has been developed for the case $W{>}W_{\text{c}}$, with a weak integrability breaking mechanism.

Acknowledgements.

Acknowledgments. We thank E. Altman, I. Gornyi, V. Khemani, Z. Lenarcic, J. Moore, A. Mirlin and D. Polyakov for inspiring discussion. FE expresses his sincere gratitude to Ehud Altman and the team at UC Berkeley for their generous hospitality . SB would like to thank G. De Tomasi for earlier collaboration on a similar topic and for many insightful discussions. SB acknowledges support from DST, India, through Ramanujan Fellowship Grant No. SB/S2/RJN-128/2016 and ECR/2018/000876; FE was supported from the DFG under Grants No. EV30/11-1 and EV30/12-1. We acknowledge computing time on the supercomputer ForHLR funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Education and Research.

Appendix A Time dependence of $\Phi(x,t)$

In the localized phase ${\Delta x}(t)$ increases slowly (as shown in Fig. 2) with time above the non-interacting saturation within our observation time ($\sim 10^{3}$ in unit of inverse $t_{\text{h}}$). In the distribution function $\Phi(x,t)$ this propagation is reflected via the lifting of the tails with time as seen in Fig. 6. The tail-lifting eventually produces a stretched-exponential shape of $\Phi(x,t)$.

Appendix B Oscillations in $\beta$

We analyze the spike observed in $\beta({\Delta x})$ (see Fig. 3), which becomes more pronounced at increasing disorder strenth $W$. It occurs at $\Delta x\lesssim 1$ (cf. Fig. 7), and therefore we interpret it as a signature of the lattice spacing, which is is taken as the unit of length. In essence, the spike is a manifestation of strong oscillations that induce a non-mononous growth of $\Delta x(t)$ in time. Such oscillatory behaviour arises from neighboring pairs of sites that have a common local potential, which deviates strongly from the environment and thus is forming a ”trap”.

The trap constitutes a two-site system with an associate energy doublet that is split by the hopping term. The dynamics of the associated charge density exhibits a characteristic frequency $\omega\approx t_{\text{h}}=1$. The corresponding oscillations in the correlator $\Phi(x,t)$ are illustrated in Fig. 7 for a given disorder configuration with trap near $x{=}0$. Strong oscillations are also seen in Fig. 1 differing from Fig. 8 by rescaling the on-site potential to smaller values with factor $3/5$ (blue) and $1.5/5$ (red).

We discuss the effect of interactions on time-dependent fluctuations in $\Delta x(t)$: In the non-interacting limit, $V{=}0$, pronounced oscillations are seen large disorder (corresponding to strength $W{=}2$) in Fig. 8; they persist up to very long times. With decreasing disorder, Fig. 1, lower frequencies are typically mixing in because resonances between more distant sites become accessible, converting oscillations to mesoscopic fluctuations. This type of oscillation survives the ensemble average, because return times tend to be integer multiples of $2\pi/t_{\text{h}}$.

With finite interaction the random potential will, in general, undergo a static renormalization; it shifts fluctuation frequencies without affecting the qualitative physics. New effects can be expected from dynamical renormalizations, i.e., from dephasing; it is a characteristic feature of many-body systems expressing the possibility of particles to scatter off each other and thus exchanging momentum and energy. In the correlator $\Phi(x,t)$ dephasing manifests itself in a damping of the mesoscopic fluctuations, see Fig. 8, bottom and 1. The defining characteristics of the mbl-phenomenon is that despite mesoscopic fluctuations washing out very quickly with increasing $V$, the average return probability is seen to keep values of order unity that hardly decay at all, see Fig. 8 and 1.

Appendix C Evolution of $\Delta x(t)$ in a quasi-periodic (AAH) potential

In the main text we have analysed the evolution of the correlator $\Phi(x,t)$ for the case of fully disordered potentials, Fig. 3. We here display the analogous data obtained for the quasi-periodic case.

Thermalizing regime.

Before we discuss the main data set, Fig. 10, we report an effect that we observe for quasi-periodic potentials, which is not seen in the fully random case. In Fig. 9 we show the correlation function $\Phi(x,t)$ at very large times. Despite the fact that the dynamics has explored the full system size, oscillations are seen that do not show a tendency to equilibrate, see also Fig. 10 for other disorder values. While a detailed study of this effect is beyond the scope of this work, a few comments are in place.

A preliminary analysis suggests that the amplitude of the oscillations decays with increasing system size in a power law-fashion; therefore, we tentatively consider them as a finite size effect. Presumably, the oscillations are related to the correlations in the underlying quasi-periodic potential; hence, one would expect that their shape is sensitive to the choice of the parameter $\alpha$ in the definition of $\epsilon_{x}$ (see 888When working with quasi-periodic potential (Aubry-André model) we choose, $\epsilon_{x}{=}W\cos(2\pi\alpha x+\phi)$ with $\alpha=2/(\sqrt{5}+1)$ and $\phi$ is the random phase distributed uniformly between $[0,2\pi]$ for further details). Furthermore our data Fig. 10 indicates that due to these oscillations at moderate disorder $W{=}1.5-3.0$, i.e. in a regime that is believed to be thermal, finite size effects can be strong. We believe that they should be observable, in principle, also in typical cold-atom experiments Lüschen et al. (2017). Finally, in recent computational work it was found that the critical disorder strength could also be dependent on $\alpha$: $W_{\text{c}}(\alpha)$ Doggen and Mirlin (2019). Conceivably, this finding is related to effects such as displayed in Fig. 9.

The analysis of $\Phi(x,t)$ for the quasi-periodic case in terms of the effective exponent $\beta(t)$ is done in analogy to the fully random case. The middle panel of Fig. 10 displays the evolution of ${\Delta x}(t)$ at system sizes $L{=}12,16,20,24$ for interaction $V/t_{\text{h}}=2.0$. At small disorder, ${\Delta x}(t)$ is seen to reach a plateau value $\approx L/3.5$ indicating saturation due to the finite system size. While the saturation width is only slowly decreasing with disorder, $W$, the saturation time is clearly seen to be a strongly growing function of $W$ exceeding the computational time window already at moderate disorder $W{=}2.5$.

The lower panel shows the corresponding evolution of the logarithmic derivative of ${\Delta x}(t)$, i.e. the exponent function $\beta(t)$. As announced before, the evolution for the quasi-periodic case is analogous to the random disorder case (see Fig. 3 lower panel). In particular, extremely strong finite size effects hamper a reliable extraction of the asymptotic exponent. The overall trend is that the short-time exponent at weak disorder takes larger values which are not far from the diffusive regime, $\beta{=}1/2$.

Strong disorder. At stronger disorder the level statistics is seen to approach the Poissionian limit in (Fig. 8). We thus associate with disorder values $W\gtrsim 5$ a regime close or in the many-body localized phase. Correspondingly, with increasing disorder the dynamics as displayed by ${\Delta x}(t)$ slows down considerably, see Fig. 11; in the quasi-periodic and the fully random situations it shows the same qualitative behavior. Since effects of rare regions are not expected in the quasi-periodic situation, our observations suggest that the typical strong disorder effects may not be the origin of slow dynamics and creep.