Slow Dynamics in Translation-Invariant Quantum Lattice Models

TL;DR

This paper investigates the potential for slow, localized-like dynamics in translation-invariant quantum lattice models, revealing regimes where systems exhibit slow delocalization despite lacking disorder, challenging traditional views on many-body localization.

Contribution

The authors develop a degenerate perturbation theory approach to identify and estimate slow dynamical regimes and delocalization times in translation-invariant quantum systems, highlighting robustness of integrability.

Findings

Identification of slow dynamical regimes in translation-invariant models

Construction of models with transient localization-like behavior

Estimation of delocalization time scales in these models

Abstract

Many-body quantum systems typically display fast dynamics and ballistic spreading of information. Here we address the open problem of how slow the dynamics can be after a generic breaking of integrability by local interactions. We develop a method based on degenerate perturbation theory that reveals slow dynamical regimes and delocalization processes in general translation invariant models, along with accurate estimates of their delocalization time scales. Our results shed light on the fundamental questions of robustness of quantum integrable systems and the possibility of many-body localization without disorder. As an example, we construct a large class of one-dimensional lattice models where, despite the absence of asymptotic localization, the transient dynamics is exceptionally slow, i.e., the dynamics is indistinguishable from that of many-body localized systems for the system sizes and time scales accessible in experiment and numerical simulations.

Tables

Table 1. Table 1: Behaviour of the saturation plateau of 1st and 2nd order of DPT as L → ∞ → 𝐿 L\to\infty for various models. ↗ ↗ \nearrow and ↘ ↘ \searrow illustrate the fact that D ​ ( t → ∞ ) 𝐷 → 𝑡 D(t\to\infty) increases or decreases as a function of the system size.

Model	1st order	2nd order
XXZ	$\searrow$	$\searrow$
Hopping+NNN	$\searrow$	$\searrow$
Hopping+3-body	$\searrow$	$\searrow$
Hopping+NNN+3-body	$\searrow$	$\searrow$
XXZ+NNN	$\nearrow$	$\searrow$
XXZ+3-body	$\nearrow$	$\searrow$
XXZ+NNN+3-body	$\nearrow$	$\searrow$
XXZ+NNN+3-body+4-body	$\nearrow$	$\searrow$
XXZ+ all up to range-4	$\nearrow$	$\nearrow$

Equations

$$H_0$$

$$V$$

$$S=-{\mathrm{tr}}_A\left({\rho }_A\log {\rho }_A\right),$$

$$|\psi ⟩=\underset{j}{⨂}\left(\cos \frac{{\theta }_j}{2}|{\downarrow }_j⟩+e^{i{\varphi }_j}\sin \frac{{\theta }_j}{2}|{\uparrow }_j⟩\right).$$

$${\cos }^r{\theta }_j={\xi }_j.$$

$$H=H_0+V,$$

$$[H_0,T_{\mathbf{n}}]=J_{\mathbf{n}}{T}_{\mathbf{n}},$$

$$T_{\mathbf{n}}=\sum _{b=1}^M{\lambda }_b\sum _{i=1}^L{F}_{\mathbf{n}}^{ib}.$$

$$H^{[1]}=e^{S_1}H{e}^{-S_1}=H+[S_1,H]+\dots$$

$$=H_0+T_0+O({\lambda }^2/J).$$

$$S_1=\sum _{\mathbf{n}\ne 0}\frac{T_{\mathbf{n}}}{J_{\mathbf{n}}}.$$

$$T_0^{\text{3-body}}$$

$$T_0^{\text{XXZ}}$$

$$H^{[2]}=e^{S_2}{e}^{S_1}H{e}^{-S_1}{e}^{-S_2}=H+[S_1+S_2,H]+\dots$$

$$=H_0+T_0+\sum _{\mathbf{n}\ne 0}\frac{1}{2J_{\mathbf{n}}}[T_{-\mathbf{n}},T_{\mathbf{n}}]+O({\lambda }^3/J^2),$$

$$S_2=\sum _{\mathbf{n}\ne 0}\frac{1}{J_{\mathbf{n}}^2}[T_{\mathbf{n}},T_0]-\sum _{\begin{array}{c}\mathbf{n}\ne 0\\ {\mathbf{n}}^{\prime }\ne \{-\mathbf{n},0\}\end{array}}\frac{1}{2J_{\mathbf{n}}{J}_{\mathbf{n}+{\mathbf{n}}^{\prime }}}[T_{\mathbf{n}},T_{{\mathbf{n}}^{\prime }}].$$

$$H^{[3]}=H^{[2]}+\sum _{\begin{array}{c}\{\mathbf{n},{\mathbf{n}}^{\prime \prime }\}\ne 0\\ {\mathbf{n}}^{\prime }\end{array}}\frac{1}{2J_{\mathbf{n}}{J}_{{\mathbf{n}}^{\prime \prime }}}[T_{\mathbf{n}},[T_{{\mathbf{n}}^{\prime \prime }},T_{{\mathbf{n}}^{\prime }}]]+$$

$$-\sum _{\{\mathbf{n},{\mathbf{n}}^{\prime },{\mathbf{n}}^{\prime \prime }\}\ne 0}\frac{1}{6J_{\mathbf{n}}{J}_{{\mathbf{n}}^{\prime }}}[T_{\mathbf{n}},[T_{{\mathbf{n}}^{\prime }},T_{{\mathbf{n}}^{\prime \prime }}]]+O({\lambda }^4/J^3).$$

$$D(t,k)=\frac{1}{Z}\text{tr}\Big{(}e^{iHt}{\tilde{\sigma}}^{(3)}_{k}e^{-iHt}{\tilde{\sigma}}^{(3)}_{k}\Big{)},$$

$${\tilde{\sigma }}_k^{(3)}=\frac{1}{\sqrt{L}}\sum _j{\sigma }_j^{(3)}\exp (i2\pi jk/L),$$

$$Z=\text{tr}({\tilde{\sigma }}_k^{(3)†}{\tilde{\sigma }}_k^{(3)}),$$

$$e^{-iHt}={\mathcal{U}}^{[m]†}{e}^{-i{H}^{[m]}t}{\mathcal{U}}^{[m]},$$

$$\text{tr}\left({\tau }_k^{(3)\ast }(t){\tau }_k^{(3)}\right),$$

$$\displaystyle\text{tr}\Bigl{\{}e^{iH^{[m]}t}{\tilde{\sigma}}^{(3)\dagger}_{k}e^{-iH^{[m]}t}{\tilde{\sigma}}^{(3)}_{k}\Bigr{\}}+O(\lambda^{2}/J^{2}).$$

$$\forall |{\psi }_1⟩,|{\psi }_2⟩\in \mathcal{S}{H}_0|{\psi }_1⟩-H_0|{\psi }_2⟩=0.$$

$$D(t,k)\sim \exp (-A{t}^b),$$

$$D(t,k)\sim \exp (-A\ln (t)),$$

$$F_{-400}^i$$

$$F_{-44-4}^i$$

$$F_{-444}^i$$

$$F_{-480}^i$$

$$F_{0-4-4}^i$$

$$F_{0-44}^i$$

$$F_{000}^i$$

$$F_{04-4}^i$$

$$F_{044}^i$$

$$F_{4-80}^i$$

Full text

I Introduction

One of the central questions of quantum statistical physics is how isolated many-body systems reach thermal equilibrium. The process of thermalization results from the spreading of quantum information into non-local degrees of freedom during the system’s unitary evolution. In ergodic systems, this spreading is fast (ballistic), since the individual eigenstates of the system are highly-entangled thermal states Deutsch (1991); Srednicki (1994); Rigol et al. (2008). On the other hand, there is growing interest in non-ergodic systems, which include integrable models Baxter (1989) and many-body localized (MBL) systems Anderson (1958); Basko et al. (2006); Gornyi et al. (2005). In the latter, strong disorder significantly constrains quantum dynamics due to the emergence of an extensive number of Local Integrals of Motion (LIOMs) Serbyn et al. (2013); Huse et al. (2014), which cause the information to spread very slowly, i.e., logarithmically in time Znidaric et al. (2008); Bardarson et al. (2012).

In this paper we investigate the possibility of slow dynamics in quantum lattice models with local interactions and in the absence of disorder. This is motivated by the open question of integrability breaking in quantum systems: what are the constraints on quantum dynamics following a weak but generic breaking of integrability? Does the integrability breakdown in large, translation-invariant lattice systems proceed in a smooth, classical-like KAM style Kolmogorov (1954), where non-ergodic regions remain for finite integrability breaking, or is ergodicity immediately restored at asymptotic time scales Prosen (1998, 1999)? We answer this question by providing a theoretical formalism that explains the appearance of slow ergodic dynamics on very long time scales.

For simplicity, we limit ourselves to the case of spin systems described by the Hamiltonian $H=H_{0}+\lambda V$, where $H_{0}$ is a classical potential energy and $V$ is a quantum hopping (tunneling) term. The eigenstates of $H_{0}$ are classical product states of spins, while $V$ introduces quantum dynamics. Such models have recently been investigated as potential analogs of MBL in translation-invariant systems when $\lambda$ is small Carleo et al. (2012); De Roeck and Huveneers (2014a); Schiulaz and Müller (2014); Grover and Fisher (2014); Schiulaz et al. (2015); Yao et al. (2016); De Roeck and Huveneers (2014b); van Horssen et al. (2015); Hickey et al. (2016); Kim and Haah (2016); He and Weng (2016); Bols and De Roeck (2016); Mondaini and Cai (2017); Yarloo et al. (2018); Schecter et al. (2017); Lan et al. (2017); Schiulaz et al. (2017). Despite much effort, the understanding of models for small but finite $\lambda$ remains less complete than that of strongly disordered systems, partly due to more pronounced finite size effects Papić et al. (2015). For example, while some signatures of MBL-like dynamics have been observed in such models Carleo et al. (2012); Schiulaz and Müller (2014); Schiulaz et al. (2015); van Horssen et al. (2015), the relevance of non-perturbative effects for delocalization has also been pointed out De Roeck and Huveneers (2014b). Furthermore, Ref. Yao et al., 2016 argued that ergodicity is restored in the thermodynamic limit, resulting in a “quasi-MBL” phase.

In this paper we develop a general formalism based on degenerate perturbation theory (DPT) that accurately describes the long-time dynamics of systems where $H_{0}$ consists of $k\geq 2$-body interactions between particles and $V$ is a single-particle hopping term. ($H_{0}$ is taken to be diagonal in the computational basis and can be viewed as a model of an integrable system.) In contrast to the related perturbative arguments De Roeck and Huveneers (2014b); Schiulaz et al. (2015), previously used for qualitative analysis of particular models, the DPT below is shown to yield a quantitatively accurate description of the time evolution of the initial inhomogeneity, thus serving as a general diagnostic of the possible delocalization processes.

The physical picture resulting from our study is that of a slow dynamical regime at intermediate time, which is exponentially long in the range of interaction terms in $H_{0}$. Thus, even though the studied models are asymptotically not localized, the time scales that reveal delocalization can be very large. To illustrate this slow transient dynamics, we introduce an example of a clean 1D lattice model with $k=3$-body interactions whose dynamics is shown to be indistinguishable from MBL systems at the experimentally accessible time scales. In contrast to previous work Schiulaz and Müller (2014); Yao et al. (2016); Grover and Fisher (2014); Schiulaz et al. (2015), the model studied here does not require very small hopping energy scales to exhibit MBL-like features, and moreover it includes only one particle species, which makes it more amenable to numerical simulations. This is demonstrated by matrix product state simulations on large systems ($L\lesssim 200$ particles), which provide strong numerical evidence that the 3-body model displays a clear log-like growth of entanglement entropy, in stark contrast with the integrable XXZ model. However, using DPT, we also access much longer time scales and show that the 3-body model does not display true MBL as it delocalizes in the 2nd order in DPT.

The remainder of this paper is organized as follows. In Sec. II we introduce the model and numerically demonstrate that it features slow dynamics. In Sec. III we present the general formalism of DPT. In Sec. IV we apply DPT to study the relaxation of the 3-body model, while Sec. V contains a generalization of our results to other types of local models. Our conclusions are presented in Sec. VI. Appendices contain further discussion on the finite-size effects of numerical simulations and details of DPT for the 3-body model.

II A model with slow dynamics

Consider the following 1D open chain of spins $\frac{1}{2}$ with length $L$:

[TABLE]

where $\{\sigma_{j}^{(\alpha)}\}$ is the Pauli basis on site $j$ and $\sigma_{j}^{\pm}=(\sigma_{j}^{(1)}\pm i\sigma_{j}^{(2)})/2$. Interaction amplitudes $J_{k}$ are taken to be of the same order, $J$, and irrational to avoid commensurate terms in the DPT below. For numerical demonstrations we choose $J_{1}=\sqrt{2}/4$, $J_{2}=\sqrt{3}/4$, $J_{3}=\sqrt{5}/8$, but our results are not sensitive to these precise values. Moreover, as we show below, our results are insensitive to the precise value of $\lambda$, as long as $\lambda\ll J_{i}$. Models like Eq. (II) physically arise in the large-$U$ limit of the Bose-Hubbard model Pachos and Plenio (2004) and in polar molecules Buchler et al. (2007).

We are interested in the dynamics for weak breaking of integrability, $\lambda\ll J$. We characterize dynamics by the entanglement entropy,

[TABLE]

where the reduced density matrix $\rho_{A}={\rm tr}_{B}\left|\psi\right\rangle\left\langle\psi\right|$ is obtained by tracing out the degrees of freedom of the subsystem $B$ for a bipartition $A\cup B$ of the entire system. We perform a global quench from the product state

[TABLE]

Here $\varphi_{j}$ is a uniform random phase, while $\theta_{j}$ is obtained from a random uniform variable $\xi_{j}\in[-1,1]$ via the transform

[TABLE]

Parameter $r$ biases the orientation of each spin on the Bloch sphere. For $r=1$ one has spin-$1/2$ states that are random uniform on the Bloch sphere. For large $r$, on the other hand, the distribution is biased towards the poles of the Bloch sphere, with the width scaling as $\approx 1/\sqrt{r}$. In the limit $r\to\infty$ one recovers random computational states, i.e., states where each spin is either $\left|\uparrow\right\rangle$ or $\left|\downarrow\right\rangle$. Our results are independent of the choice of $r$ (see Appendix A) and we fix $r=11$ for optimal balance between state-to-state fluctuations and the magnitude of $S$, allowing for longest simulation times.

In Fig. 1 we show the representative dynamics of $S$ in the model (II) starting from a single product state. Time evolution is carried out using time-evolving block-decimation (TEBD) algorithm (Vidal, 2003). We consider a large chain of $L=64$ sites and evolve the system using bond dimensions up to $\chi=350$. We see that for all values of $\lambda$, except $\lambda=1$, there is a clear difference between the XXZ model ($J_{2}=J_{3}=0$) and the 3-body model in Eq. (II). In particular, even for $\lambda$ as large as 0.4, we find a growth of entropy in the 3-body case which is much slower than linear.

In Fig. 2(a) we show the entropy growth for $\lambda=0.2$ for the 3-body and XXZ model, in both cases starting from the same initial state in Eq. (3). On the time scales $t\lesssim 200$, we observe a clear difference between the 3-body model in Eq. (II) (red) and the XXZ model (blue). In the XXZ case, $S\propto t$, while in the 3-body case the data is consistent with $S(t)\propto\log t$. Phenomenologically, the linear growth of entropy in the XXZ model is due to the propagation of coherent quasiparticles with a velocity $\sim\lambda$ Calabrese and Cardy (2016, 2005). Reducing the hopping only affects the slope $S(t)\propto c(\lambda)t$, where $c(\lambda)\propto\lambda$, and the XXZ chain remains delocalized for arbitrarily small $\lambda$. Tiny deviation from linear growth in the XXZ model is likely due to finite $L=64$ (see Appendix A). In contrast, the spreading of entropy in the 3-body model for small $\lambda$ appears logarithmic, even in a very large system, which is reminiscent of MBL physics Znidaric et al. (2008); Bardarson et al. (2012). We emphasize, however, that our numerical result in Fig. 2 does not rule out the possibility of a power growth with a small exponent, even though the logarithmic dependence appears to give a better fit, as discussed in Appendix A.

To understand the mechanism of the slow dynamics, in Figs. 2(b),(c) we examine the snapshots of entropy, evaluated at all bonds $j,j+1$, and the local magnetization $\langle\sigma_{j}^{(3)}\rangle$ at different times. For the given initial configuration, we observe a “blocking region” (shaded), which does not decay on the accessible time scales in the 3-body model, but decays by the time $t\sim 50$ in the XXZ model.

In order to understand the long-time behaviour of the model, one must go beyond the TEBD simulations (which are limited to short times) and exact diagonalization (which is susceptible to finite size effects). We therefore introduce an analytic method based on perturbation theory for the dynamics. This method will show that the model in Eq. (II) delocalizes at long (but finite) times corresponding to the 2nd order in $\lambda$. More generally, this method will allow us to understand the nature of the delocalizing processes order by order, and to quantify the role of finite-size effects.

III Degenerate perturbation theory for the dynamics

When $\lambda\ll J$, the Hamiltonian of Eq. (II) separates into the unperturbed part $H_{0}(\{\sigma^{(3)}\})$ and the perturbation $V$. Without $V$, the system is integrable – each computational state is an eigenstate. In a semi-classical picture of small $\lambda$, the interactions still tend to localize domain walls, resulting in slow dynamics. We want to know how slow this dynamics is, and, specifically, to treat also finite values of $\lambda$. We use the Schrieffer-Wolff transformation Schrieffer and Wolff (1966) (well-known in the context of ground state physics) as a framework for degenerate perturbation theory (DPT) to systematically keep track of corrections in orders of $\lambda/J$. Note that the first few orders of DPT may not capture the eigenstates at arbitrary energy. However, we show that it is possible to access dynamics in the time scales in orders of $1/\lambda$ at all energies, thereby accurately revealing the breakdown of integrability.

In this Section, we outline the DPT formalism for a general Hamiltonian of the form

[TABLE]

where we assume that $H_{0}$ is a $K$-local operator and $V$ is $R$-local, with $K>R$. We aim to find all different energy variations on the unperturbed eigenstates after applying the perturbation. $H_{0}$ may contain $N$-many terms with different amplitudes $\{J_{i}\}$, while $V$ contains $M$ hopping terms with amplitudes $\{\lambda_{i}\}$, and we assume $J_{1},\ldots,J_{N}\gg\lambda_{1},\ldots,\lambda_{M}$. In addition, $M$ must be small enough for the perturbation to remain a correction to the unperturbed Hamiltonian.

We work in the computational $\{\sigma^{(3)}\}$ eigenbasis in which $H_{0}$ is diagonal. The idea of DPT is to find such a unitary transformation that will generate a block-diagonal form while eliminating higher orders in $\lambda$. To get an intuitive feeling, let us look at the effect of $V$ on an unperturbed eigenstate $\left|\psi\right\rangle$, $H_{0}\left|\psi\right\rangle=E_{0}\left|\psi\right\rangle$. In the ordinary 1st order DPT, one would diagonalize $V$ in each of the subspaces ${\cal S}_{0}$ corresponding to a given $E_{0}$. As it turns out, $V$ can still have a block-diagonal structure on ${\cal S}_{0}$, i.e., ${\cal S}_{0}=\oplus_{k}{\cal S}_{0}^{(k)}$, where subspaces ${\cal S}_{0}^{(k)}$ consist of degenerate states $\left|\psi\right\rangle$ connected by a single application of $V$. We call them a path-connected degenerate subspace (PCDS), see Fig 3. In 1st order their energy $E_{0}$ will be split by $O(\lambda)$, but most importantly, in a given $S_{0}^{(k)}$ some spins have a fixed orientation for all states. This means that the dynamics happening in 1st order (on a timescale $\sim 1/\lambda$) will affect only certain spins, while others remain frozen. Those spins form “blocking” regions, like in Fig. 2, and are responsible for the slow dynamics. We now formalize this reasoning and systematically extend it to higher orders.

We start by defining operators $T$ such that $T_{{\bf n}}\subseteq V$, $\sum_{{\bf n}}T_{{\bf n}}=V$ and

[TABLE]

where $J_{{\bf n}}=\left(\sum_{a=1}^{N}n_{a}J_{a}\right)$. To ensure the unitarity of the transformations at any order, $T$’s have to satisfy $T^{\dagger}_{{\bf n}}=T_{-{\bf n}}$. The operator $T_{{\bf 0}}$ (if it exists) commutes with $H_{0}$ and thus spans a degenerate subspace at 1st order. Simply put, $T_{{\bf 0}}$ is a projection of $V$ to the degenerate subspaces of $H_{0}$ (the block-diagonal part of $V$), while $T_{{\bf n}\neq{\bf 0}}$ denote the corresponding off-diagonal blocks of $V$. The number of different operators $T_{{\bf n}}$ is system dependent. Each operator $T$ is translation invariant and as such it can be decomposed as

[TABLE]

Operators $F$ are at most ($2K+R-2$)-local and are the starting points of the expansion that we describe below. Each $F^{ib}_{{\bf n}}$ transforms a basis state into another basis state, i.e., it just flips certain spins, and the vector index ${\bf n}\equiv\{n_{a}\}$ labels the difference in the unperturbed energy of a basis state $\left|\psi\right\rangle$ and the flipped one, $F^{ib}_{{\bf n}}\left|\psi\right\rangle$. Index $b$ labels different perturbations. In this work, $b$ can be omitted since we have a single perturbation and $\lambda_{1}\equiv\lambda/2$. Once $\{F\}$ are known, any order in perturbation theory can in principle be computed. The main issue is that the support of the operators increases with respect to the order of the expansion since nested commutators are involved, thus the calculation beyond the first few orders is limited by computational resources. In Appendix B we explicitly evaluate $F^{i}_{{\bf n}}$ for the 3-body model in Eq. (II), and show they are 6-local operators.

To find the 1st order expansion of the Hamiltonian, a unitary transformation $\mathcal{U}^{[1]}=e^{S_{1}}$ rotates the system to the subspace where the perturbing terms that change the unperturbed energy (${\bf n}\neq 0$) are removed, i.e., every process is resonant:

[TABLE]

For the last expression we used Eq. (6) to pick the correct transformation

[TABLE]

For example, as shown in Appendix B [see Eq. B], for the 3-body and XXZ models:

[TABLE]

In 2nd order DPT, the expansion is calculated in the same spirit. The rotation removes all perturbative terms of order $O(\lambda^{2}/J)$. This can be calculated iteratively from the 1st order,

[TABLE]

where

[TABLE]

The generator of the unitary transformation of the $i$th order expansion using the iterative method is $S_{i}\sim O(\lambda^{i}/J^{i})$. For example, the 3rd order hamiltonian is

[TABLE]

Note that the subscripts in Eq. (14) also obey ${\bf n}+{\bf n}^{\prime}+{\bf n}^{\prime\prime}={\bf 0}$ in order to keep the unperturbed energy constant.

We see that the 1st order allows dynamics only within each PCDS [Fig. 3(a)], while the 2nd order connects different PCDS through a single virtual hop [Fig. 3(b)]. The $m$th order allows connections through $m-1$ virtual hops. The rotation of the basis consists of operators that jump between different energies and in most cases generate dephasing without transport. We note that $[H^{[m]},H_{0}]=0$ for any order $m$ automatically follows from Eq. (6). When we apply the DPT below, we numerically diagonalize the Hamiltonian at each order, which is a simple way to account for the splitting of the degenerate levels.

We note that certain models, e.g., the one introduced in Ref. van Horssen et al., 2015, do not have a degenerate subspace in 1st order, i.e., $T_{{\bf 0}}=\{\}$. This means that the first non-trivial order in such models is the 2nd order. In the case of Ref. van Horssen et al., 2015 we only have $T_{{\bf n}_{1}},T_{-{\bf n}_{1}}$, ${\bf n}_{1}=1$. Consequently, not only the 1st, but also all odd orders do not generate any new terms because odd orders include nested commutators of an odd number of $T$’s, which require that a sum of an odd number of ${\bf n}$’s must equal zero. Since the two choices are $({\bf n}_{1},-{\bf n}_{1})$ this is impossible. Such models, which usually result from imposing classical kinetic constraints on the Hamiltonian, are expected to show the absence of relaxation for longer times due to the vanishing of odd orders of perturbation theory.

IV Polarization decay

We now focus on a general dynamical probe of relaxation Pal and Huse (2010); Yao et al. (2016): we prepare an initial inhomogeneity in the spin magnetization and monitor its decay as a function of time,

[TABLE]

where we have introduced the Fourier transform of the Pauli operator

[TABLE]

assuming periodic boundary conditions. The normalization of $D(t,k)$ is

[TABLE]

where the trace is taken over zero magnetization sector.

The interpretation of $D(t,k)$ is: throw a particle of momentum $k$ into the system; after some time remove the particle and measure the state overlap with the initial state. If the particle scatters, the memory of the initial state gets lost and $D(t,k)\ll 1$; if the particle does not scatter, by removing it one returns to the original state and $D(t,k)=1$. For momentum $k\approx 1$, scattering will only take place if eigenstates are extensive, thus we interpret $D(t,1)$ as a probe of delocalization of the system. Due to translation invariance, in a finite system the polarization always vanishes at $t\to\infty$. In the thermodynamic limit, if the system is in a quasi-MBL phase, one expects a time scale for the decay of $D(t,1)$ that exponentially diverges with the system size Yao et al. (2016).

We now apply DPT to Eq. (15). The denominator of Eq. (15) is invariant under unitary basis rotations. The time evolution operator in the numerator is transformed to the $m$th order as

[TABLE]

where $\mathcal{U}^{[m]}=\prod_{i=1}^{m}e^{S_{m}}$. Using the cyclic property of the trace, the numerator of Eq. (15) is written as

[TABLE]

where $\tau^{(3)}_{k}=\mathcal{U}^{[m]}{\tilde{\sigma}}^{(3)}_{k}\mathcal{U}^{[m]\dagger}$, $\tau^{(3)*}_{k}=\mathcal{U}^{[m]}{\tilde{\sigma}}^{(3)\dagger}_{k}\mathcal{U}^{[m]\dagger}$ are the effective quasiparticles of the rotated picture. The corrections to ${\tilde{\sigma}}^{(3)}$ due to the rotation are given in orders of $O(\lambda/J)$. The expansion of the unitary in Eq. (19) results in

[TABLE]

The terms of order $O(\lambda/J)$ vanish since the operators inside the trace are purely off-diagonal. To understand that, assume that the basis used to evaluate the trace is the unperturbed eigenbasis. The Hamiltonians in DPT are block-diagonal in the computational basis at any order. Every block ${\cal S}$ has some basis $\{\left|\psi\right\rangle\}$ spanned by vectors of equal unperturbed energy,

[TABLE]

By default $e^{-iH^{[m]}t}$, $e^{+iH^{[m]}t}$ have the same block-diagonal structure and thus map states from ${\cal S}\rightarrow{\cal S}$. Operators ${\tilde{\sigma}}^{(3)\dagger}_{k}$, ${\tilde{\sigma}}^{(3)}_{k}$ have trivial action as they are diagonal in the unperturbed eigenbasis, so they trivially map states from ${\cal S}\rightarrow{\cal S}$. On the other hand, $[S_{1},{\tilde{\sigma}}^{(3)}_{k}]$, $[S_{1},{\tilde{\sigma}}^{(3)\dagger}_{k}]$ always map to states outside the block, which follows trivially from the definition of $S_{1}$ in Eq. (9). Thus an operator product which contains block-conserving and only one of $[S_{1},{\tilde{\sigma}}^{(3)}_{k}]$, $[S_{1},{\tilde{\sigma}}^{(3)\dagger}_{k}]$ can only have vanishing diagonal elements. This means that magnetization decay does not feature first order basis corrections.

IV.1 Plateaus in polarization decay

Using Eqs. (8)-(12), we numerically compute the magnetization $D(t,k)$ in DPT and contrast it against exact time evolution in Fig. 4. As explained above, in this particular calculation we can ignore basis rotations up to the time-scale to which the given order is accurate and thus model the dynamics according to Eq. (19). (By contrast, the calculation of entanglement entropy would be sensitive to the basis rotations.)

The comparison between the first three orders of DPT and exact evolution is shown in Fig.4. Evidently, DPT is practically exact up to the relevant breakdown time scale $t\sim J^{m}/\lambda^{m+1}$ for each order $m$. We note that a small value of $\lambda$ is chosen to resolve the plateaus in $D$, which correspond to different orders in the DPT. However, the values of the 1st and 2nd order plateau are independent of $\lambda$ and depend only on the size and number of disconnected subspaces. Moreover, even at larger $\lambda$, when the plateaus are no longer separated, we find excellent agreement between DPT and exact evolution.

IV.2 Comparison of 3-body with XXZ model

After successfully benchmarking DPT against exact time evolution, we now use DPT to compare the relaxation in the 3-body model against the well-studied example of the XXZ chain which shows fast relaxation. Fig.5 shows $D(t,1)$ plotted for different sizes $L$ and orders of DPT. Notice that we do not terminate the evolution at the time scale that would be relevant to each order, but we allow the system to evolve until it reaches a saturation plateau. This method allows us to measure how much the system relaxes in each order of DPT. It is obvious that including additional orders only lowers the values of saturation plateaus as each order contains the previous order plus some extra terms whose value is independent of the strength of the perturbation.

Interestingly, the 1st order already reveals a clear difference between the 3-body and XXZ model, Fig.5(a). In the latter case, $D$ quickly decays to a small value ($\lesssim 0.1$), which further decreases with $L$ (inset). In the 3-body case, the plateau is close to 1 and grows with the system size $L$. This means the system does not relax on the time scales where only the 1st order is relevant, and is a direct consequence of the PCDSs in Fig.3. As the system becomes larger, the more extended modes do not scatter, indicating that the fraction of extended non-local 1st order eigenstates vanishes.

In the second order, Fig.5(b), we observe a completely different behaviour of the 3-body plateau which now decreases with $L$. Finite size scaling (inset) suggests that as one increases $L$, the system becomes progressively more delocalized in the 2nd order of DPT. We note, however, that delocalization in our DPT does not automatically imply ergodicity, since any order in DPT will have disconnected subspaces whose support is a vanishing fraction of the total Hilbert space. Once the model is in a delocalized regime, the entire Hilbert space can become connected by the action of unitary rotations, which generate subleading corrections to $D$.

IV.3 Relaxation time

Using DPT we can furthermore scrutinize the possible similarity of our model in Eq. (II) with models showing “quasi-MBL” behaviour Yao et al. (2016). This can be done by investigating the finite-size scaling of the saturation time in DPT.

If the system delocalizes at a finite order in DPT, it is natural to expect that $D(t,k)$ obeys the following asymptotic behavior in time:

[TABLE]

where $A$ depends on momentum (and, therefore, system size $L$). On the other hand, if the system is “quasi-MBL” Yao et al. (2016), $D(t,k)$ should have the asymptotic form

[TABLE]

which would yield a time scale for the relaxation of the smallest Fourier mode ($k=1$) that diverges exponentially with the size of the system.

In Fig. 6 we show the delocalization time (defined as the time it takes for magnetization to drop to 5% of its initial value) as a function of system size $L$. This plot was obtained by exact diagonalization of the 3-body model for the fixed hopping $\lambda=0.2$. Despite small system sizes, Fig. 6 suggests that Eq. (21) is a much better fit to the data, suggesting that our 3-body model is in a different class from “quasi-MBL” models.

As an alternative way to probe the difference between our model in Eq. (II) and quasi-MBL models, we have also considered the standard quantifiers of integrability breaking used in random matrix theory, such as the statistics of energy level spacing. We found that the 3-body model is described by the Wigner-Dyson statistics already for $\lambda$ as small as 0.1. For this value of $\lambda$, we find the average value of the level statistics “$r$” parameter Pal and Huse (2010) to be $\langle r\rangle\approx 0.53$ (this value is obtained for $L=20$ spins, after resolving translation and discrete symmetries of the model, and it was averaged over all eigenstates, which corresponds to the infinite temperature). The obtained value is very close to the Wigner-Dyson value $\langle r\rangle\approx 0.53$, and clearly inconsistent with the value expected for Poisson statistics, $\langle r\rangle\approx 0.39$. Since the level statistics probes the smallest energy scale in the system or, equivalently, the asymptotically long time scales, this confirms our claim in Sec. IV.2 that the 3-body model delocalizes at asymptotically long times, even though it displays slow dynamics over surprisingly long intermediate timescales.

V Generalization to other models

Finally, we discuss some generalizations of the model in Eq. (II) in order to establish a more general understanding of the possible types of slow dynamics due to interaction constraints. We compare the polarization decay between different models from the point of view of 1st and 2nd order of DPT. The perturbation $V$ is always assumed to be the nearest neighbor (NN) hopping. The unperturbed Hamiltonian is chosen to have a combination of different terms, denoted by the following abbreviations: NNN stands for $\sum_{i}\sigma^{(3)}_{i}\sigma^{(3)}_{i+2}$, 3-body for $\sum_{i}\sigma^{(3)}_{i}\sigma^{(3)}_{i+1}\sigma^{(3)}_{i+2}$, 4-body for $\sum_{i}\sigma^{(3)}_{i}\ldots\sigma^{(3)}_{i+3}$. “Range-4” is used to denote all possible range-4 interactions.

By combining these interaction terms (with irrational coefficients, as mentioned in Sec. II), various models can be constructed, and their dynamical behaviour (according to the behaviour of the plateau in 1st and 2nd order DPT) is summarized in Table I. For example, we can observe that relaxation of the system is suppressed up to order $m$ if $K-R\geq m$, where $K$,$R$ are the ranges of $H_{0}$ and $V$, respectively. However, this is a necessary condition but it is not sufficient. For example, by just adding 4-body interactions to the 3-body Hamiltonian, the system still delocalizes at 2nd order DPT. We believe that this condition becomes sufficient only if $H_{0}$ contains all possible interactions up to that range.

In order to corroborate the previous statement, we show that it is possible to prevent relaxation in the second order of DPT, corresponding to the 2nd order plateau increasing with system size. For this, we require a Hamiltonian with 4-body interaction terms. Fig. 7 illustrates the saturation plateau of the 2nd order Hamiltonian for two different models where $K=4$, $R=2$. In the first case, a 4-body term ($\sum_{i}\sigma^{(3)}_{i}\ldots\sigma^{(3)}_{i+3}$) is added to the 3-body Hamiltonian. In the second case, the most generic range-4 unperturbed Hamiltonian is chosen by adding all possible range-4 combinations of 2,3,4-body interactions. We observe that 4-body interactions in themselves are not enough to prevent relaxation of the system in the 2nd order as the saturation plateau decreases. On the other hand, the most generic range-4 interaction, obtained by taking the 3-body model and adding to it terms such as $\sum_{i}\sigma^{(3)}_{i}\ldots\sigma^{(3)}_{i+3}$, $\sum_{i}\sigma^{(3)}_{i}\sigma^{(3)}_{i+1}\sigma^{(3)}_{i+3}$, $\sum_{i}\sigma^{(3)}_{i}\sigma^{(3)}_{i+2}\sigma^{(3)}_{i+3}$, $\sum_{i}\sigma^{(3)}_{i}\sigma^{(3)}_{i+3}$, does indeed prevent relaxation in the 2nd order.

A general picture which emerges is that higher range terms in the diagonal part $H_{0}$ inhibit transport. More precisely, the previous results support our conjecture that a generic, range-$K$ translation-invariant interaction leads to the absence of relaxation up to order $m=K-R$, where $R$ is the range of the hopping term $V$. This is in line with the situation in MBL systems, where $H_{0}$ is expressed in terms of Local Integrals of Motion (LIOMs) Serbyn et al. (2013); Huse et al. (2014) (LIOMs) and contains terms of arbitrary range (with decaying strengths). The LIOMs are expected to be robust to adding a small $V$, thus the system should stay localized in all orders of DPT. The DPT picture therefore presents a general framework which allows to understand truly localized systems, like disordered MBL models, as well as (local) translation-invariant systems that may display localization-like features only up to large but finite times.

VI Conclusion

We have introduced a general formalism to characterize slow dynamics in a broad class of systems with finite-range interactions and bounded local Hilbert space in any dimension. We illustrated the formalism and slow dynamics in a particular 1D model by demonstrating the plateaus in the decay of spin polarization and the log-like spreading of entanglement entropy. These results are insensitive to the choice of the parameters of the Hamiltonian as long as $\lambda\ll J$, i.e., they depend solely on the structure of the DPT expansion.

We showed that the dynamics can be significantly inhibited by changing the range of the diagonal term $H_{0}$. More precisely, the order $m$ plateau in the DPT approaches $1$ as $L\to\infty$ if all interaction terms with range $\leq m+2$ and with incommensurate amplitudes are included in $H_{0}$. Nevertheless, the system delocalizes by the order $m+1$ of DPT. Our 3-body model is an explicit example of this: it has a robust $m=1$ plateau, but delocalizes in $m=2$ order of DPT. Higher order plateaus, e.g., $m=2$, can be stabilized at the expense of including all interaction terms of range $\leq 4$.

The general scenario of the absence of relaxation up to a finite time in translation-invariant systems should be contrasted with disordered MBL systems. In the latter case, $H_{0}$ is given in terms of LIOMs and contains interactions of arbitrary range with a decaying strength Serbyn et al. (2013); Huse et al. (2014). For small nonzero $\lambda$, the LIOMs are redefined and relaxation would be absent in all orders in DPT. We also note that local models without a degenerate subspace exist (van Horssen et al., 2015), where odd orders in DPT do not contribute and thus delay the onset of delocalization. Finally, it would be of interest to extend the DPT method to two component models Schiulaz and Müller (2014); Yao et al. (2016); Grover and Fisher (2014); Garrison et al. (2017); Veness et al. (2017); Smith et al. (2017), which generally become non-local when one particle species is integrated out.

VII Acknowledgements

We thank François Huveneers for useful discussions. Z.P. and A.M. acknowledge support by EPSRC grant EP/P009409/1 and and the Royal Society Research Grant RG160635. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data. D.A. acknowledges support by the Swiss National Science Foundation. M.Z., M.M. and T.P. acknowledge Grants J1-7279 (M.Z.) and N1-0025 (M.M. and T.P.) of Slovenian Research Agency, and ERC grant No. 694544 - OMNES (T.P.).

Appendix A Convergence with system size and the choice of the initial states

Our results in the main text hold for generic initial product states. Two important special cases of such states are random computational states, i.e., states for which each spin is pointing either up or down (with equal probability), and random states in the sense of the Haar measure, i.e., states drawn uniformly on the Bloch sphere. In order to be able to smoothly interpolate between these two cases, we introduced a biased Bloch ensemble in Eq. (3) parametrized by $r$ in Eq. (4).

State-to-state fluctuations in the entanglement entropy $S(t)$ are largest for $r=\infty$ and smallest for $r=1$ (uniform Bloch), while the average value of $S(t)$ is the largest for $r=1$ and smallest for $r=\infty$. The maximal time that can be simulated by the TEBD algorithm depends on $S(t)$, and therefore large $r$ would be preferred. However, large $r$ would also necessitate a large ensemble size, in order to suppress state-to-state fluctuations. Therefore, some intermediate choice of $r$ would be optimal in practice. In the main text we have used $r=11$. We emphasize, however, that this choice is just for numerical convenience, and qualitatively similar results are obtained for other choices which we now demonstrate.

In Fig. 8 we show data similar to Fig. 2 in the main text, but here the initial state is unform on the Bloch sphere ($r=1$). One can see that even though the initial spins are not fully polarized, similar slow dynamics emerges as for the states with $r=11$. The difference is only in the prefactor of the log-like entropy growth, which is larger here, and thus only shorter times can be reached in the simulation. Spatial profile of the entanglement entropy $S(j)$ for all cuts as well as magnetization profiles $\langle\psi(t)|\sigma_{j}^{(3)}|\psi(t)\rangle$ are again qualitatively different for the 3-body model compared to the XXZ chain.

The initial state used in the main text is an instance of a state with $r=11$. This is the reason why initially spins are almost fully polarized in the $\pm z$-direction. In Fig. 9 we show the average entanglement entropy, where the averaging is done over the ensemble with $r=11$. The results shown in Fig. 9 demonstrate that for the sizes used, the behavior becomes system-size independent.

Finally, we also illustrate the difficulty in numerically distinguishing the logarithmic dependence from a power law with a small power. In Fig. 10 we demonstrate that, while we can not exclude the possibility of a power-law growth of entanglement entropy in the 3-body model, logarithmic dependence appears to give better fit to the data.

Appendix B Application of degenerate perturbation theory to the 3-body model

In this section we apply the general formalism of DPT in Sec. III to the 3-body model in Eq. (II). This Hamiltonian is 3-local, leading to 6-local $F^{i}_{n_{1}n_{2}n_{3}}$’s:

[TABLE]

where $n_{1}$, $n_{2}$, $n_{3}$ are associated with the nearest-neighbour (NN), next-nearest-neighbour (NNN), and 3-body interactions, respectively, while $p_{i}$, $q_{i}$ are the projectors to the $\uparrow$, $\downarrow$ spin on the $i$th site, i.e., $p_{i}=\text{diag}(1,0)$, $q_{i}=\text{diag}(0,1)$. The operator $\Pi_{i}$ performs a reflection of an operator around the bond $i$, i.e., $\Pi_{i}(O_{i}O_{i+2})=O_{i+1}O_{i-1}$. The reflection symmetry of $F$’s is a consequence of the reflection symmetry of the full Hamiltonian.