Slow quantum thermalization and many-body revivals from mixed phase space

TL;DR

This paper explores how mixed phase space in semiclassical dynamics influences quantum thermalization and revivals, revealing new quantum scars and slow thermalization in many-body systems through variational classical approximations.

Contribution

It demonstrates that mixed phase space structures in classical variational dynamics predict anomalous quantum behaviors like revivals and slow thermalization in many-body systems.

Findings

Persistent many-body quantum revivals linked to regular phase space regions.

Identification of new quantum many-body scars causing long-time oscillations.

Slowdown of thermalization in systems initialized in regular phase space regions.

Abstract

Relaxation of few-body quantum systems can strongly depend on the initial state when the system's semiclassical phase space is mixed, i.e., regions of chaotic motion coexist with regular islands. In recent years, there has been much effort to understand the process of thermalization in strongly interacting quantum systems that often lack an obvious semiclassical limit. Time-dependent variational principle (TDVP) allows to systematically derive an effective classical (nonlinear) dynamical system by projecting unitary many-body dynamics onto a manifold of weakly-entangled variational states. We demonstrate that such dynamical systems generally possess mixed phase space. When TDVP errors are small, the mixed phase space leaves a footprint on the exact dynamics of the quantum model. For example, when the system is initialized in a state belonging to a stable periodic orbit or the surrounding regular region, it exhibits persistent many-body quantum revivals. As a proof of principle, we identify new types of "quantum many-body scars", i.e., initial states that lead to long-time oscillations in a model of interacting Rydberg atoms in one and two dimensions. Intriguingly, the initial states that give rise to most robust revivals are typically entangled states. On the other hand, even when TDVP errors are large, as in the thermalizing tilted-field Ising model, initializing the system in a regular region of phase space leads to slowdown of thermalization. Our work establishes TDVP as a method for identifying interacting quantum systems with anomalous dynamics in arbitrary dimensions. Moreover, the mixed-phase space classical variational equations allow to find slowly-thermalizing initial conditions in interacting models. Our results shed light on a link between classical and quantum chaos, pointing towards possible extensions of classical Kolmogorov-Arnold-Moser theorem to quantum systems.

Figures

Equations

$$S=\int dt\mathcal{L},\mathcal{L}=\frac{i}{2}(⟨\psi ⟩\dot{\psi }-⟨\dot{\psi }⟩\psi )-⟨\psi ⟩H\psi .$$

$$\sum _a{\dot{x}}_a\text{Im}⟨{\partial }_{x_b}\psi ⟩{\partial }_{x_a}\psi =-\frac{1}{2}{\partial }_{x_b}⟨\psi ⟩H\psi ,$$

$$\ket{\psi(\{{\bf x}_{i}\})}=\sum_{\{\sigma\}}\Big{(}V^{L}\prod_{m=0}^{L/K-1}\left[A^{\sigma_{1+mK}}({\bf x}_{{1}})A^{\sigma_{2+mK}}({\bf x}_{2})\ldots\right.\\ \left.A^{\sigma_{K+mK}}({\bf x}_{K})\right]V^{R}\Big{)}\ket{\sigma_{1}}\ket{\sigma_{2}}\ldots\ket{\sigma_{L}},$$

$$H_{\text{PXP}}=\sum _{i=1}^L{P}_{i-1}{\sigma }_i^x{P}_{i+1},$$

$$A^{\uparrow }({\theta }_i,{\varphi }_i)=\left(\begin{array}{cc}0& i{e}^{-i{\varphi }_i}\\ 0& 0\end{array}\right),A^{\downarrow }({\theta }_i,{\varphi }_i)=\left(\begin{array}{cc}\cos {\theta }_i& 0\\ \sin {\theta }_i& 0\end{array}\right).$$

$$F(t)=|⟨{\psi }_0|\exp (-iHt)|{\psi }_0⟩{|}^2,$$

$$H_{\text{PXP}}=\sum _{i,j}{P}_{i-1,j}{P}_{i,j-1}{\sigma }_{i,j}^x{P}_{i+1,j}{P}_{i,j+1}+{\mu }_z{n}_{i,j},$$

$$\begin{array}{rl}{A}_{i_1,i_2,i_3,0}^{\downarrow }& ={\cos }^I\theta {\sin }^{3-I}\theta ,\\ A_{0,\dots 0,1}^{\uparrow }& =i{e}^{-i\varphi },\end{array}$$

$${\dot{\theta }}_1$$

$${\dot{\varphi }}_1$$

$$\frac{⟨H⟩}{L}=\frac{{\mu }_z{\sin }^2{\theta }_2+{\cos }^{C-1}{\theta }_1\sin 2{\theta }_2\sin {\varphi }_1}{2+2{\cos }^2{\theta }_2{\tan }^2{\theta }_1}+1\leftrightarrow 2$$

$${\gamma }^2(\{x_a\})=\underset{L\to \infty }{\lim }\frac{1}{L}||(iH+{\dot{x}}_b{\partial }_{x_b})|\psi (\{x_a\})⟩|{|}^2.$$

$$f_T\le {\Gamma }_T,\text{with}{\Gamma }_T={\left[{\int }_0^{T}dt\gamma (\{x_a^{(0)}(t)\})\right]}^2,$$

$$H_{{\mu }_3}={\mu }_3\sum _i{P}_{i-2}\left({\sigma }_{i-1}^{+}{\sigma }_i^{-}{\sigma }_{i+1}^{+}+{\sigma }_{i-1}^{-}{\sigma }_i^{+}{\sigma }_{i+1}^{-}\right)P_{i+2},$$

$${\Gamma }_T\sim 1$$

$$H_{\text{TFIM}}=\sum _i{J}_z{\sigma }_i^z{\sigma }_{i+1}^z+h_z{\sigma }_i^z+h_x{\sigma }_i^x,$$

$$A^{\uparrow }$$

$$A^{\downarrow }$$

$$(T_{i})_{\alpha\gamma,\beta\delta}=\sum_{\sigma}\ \raisebox{-20.95781pt}{\includegraphics[width=32.52017pt]{mps_1.pdf}}\ =\sum_{\sigma}(A^{\dagger}_{i})^{\sigma}_{\alpha\beta}(A_{i}^{\sigma})_{\gamma\delta},$$

$$T_{\text{u.c.}}=T_1\cdot T_2\cdot \dots \cdot T_K.$$

$$\begin{array}{rl}⟨\tilde{\psi }⟩\dot{\tilde{\psi }}& =-i⟨\tilde{\psi }⟩H\tilde{\psi }\\ \Rightarrow \frac{da}{dt}& =-⟨\psi ⟩H\psi +i⟨\psi ⟩\dot{\psi },\end{array}$$

$$⟨\psi ⟩{\partial }_{{\theta }_{\kappa }}\psi =0,\forall \kappa \in 1,\dots ,\frac{NK}{2}.$$

$$-i⟨\psi ⟩{\partial }_{{\varphi }_{\kappa }}\psi =f_{\kappa }(\{{\mathbf{\theta }}_i\}),$$

$${\eta }_{\kappa {\kappa }^{\prime }}={\partial }_{{\theta }_{{\kappa }^{\prime }}}{f}_{\kappa }=2\text{Im}⟨{\partial }_{{\theta }_{{\kappa }^{\prime }}}\psi ⟩{\partial }_{{\varphi }_{\kappa }}\psi ,$$

$$\sum _{{\kappa }^{\prime }}{\eta }_{\kappa {\kappa }^{\prime }}{\dot{\theta }}_{{\kappa }^{\prime }}={\partial }_{{\varphi }_{\kappa }}⟨\psi ⟩H\psi ,$$

$$\sum _{{\kappa }^{\prime }}{\eta }_{\kappa {\kappa }^{\prime }}^T{\dot{\varphi }}_{{\kappa }^{\prime }}=-{\partial }_{{\theta }_{\kappa }}⟨\psi ⟩H\psi .$$

$$\{g,\rho \}={\partial }_{\mathbf{\theta }}g{\hat{\eta }}^{-1}{\partial }_{\mathbf{\varphi }}\rho -{\partial }_{\mathbf{\varphi }}\rho {\hat{\eta }}^{-1}{\partial }_{\mathbf{\theta }}g,$$

$$T=\left(\begin{array}{cccc}{\cos }^2{\theta }_1{\cos }^2{\theta }_2+{\sin }^2{\theta }_2& 0& 0& {\cos }^2{\theta }_1\\ \cos {\theta }_1{\cos }^2{\theta }_2\sin {\theta }_1& 0& 0& \cos {\theta }_1\sin {\theta }_1\\ \cos {\theta }_1{\cos }^2{\theta }_2\sin {\theta }_1& 0& 0& \cos {\theta }_1\sin {\theta }_1\\ {\cos }^2{\theta }_2{\sin }^2{\theta }_1& 0& 0& {\sin }^2{\theta }_1\end{array}\right),$$

$$|R)={\cos }^2{\theta }_2\tan {\theta }_1({\cos }^{-2}{\theta }_2{\tan }^{-1}{\theta }_1,1,1,\tan {\theta }_1).$$

$$(L|R)=1+(\cos {\theta }_2\tan {\theta }_1{)}^2$$

$$f_{1}=-i\langle{\psi|\partial_{\phi_{1}}\psi}\rangle\\ =-i\frac{L}{2}\ \raisebox{-14.45377pt}{\includegraphics[width=151.76633pt]{mps_4.pdf}}\ \ ,$$

Full text

I Introduction

Technological advances in synthetic quantum systems Polkovnikov et al. (2011); Bloch et al. (2008) have started an era where non-equilibrium dynamics of isolated quantum matter can be experimentally probed. The process of intrinsic thermalization that may occur in isolated many-body systems results in featureless thermal states and scrambling of quantum information. Systems which avoid thermalization due to an extensive number of integrals of motion, such as many-body localized (MBL) Abanin et al. (2018) and Bethe-ansatz integrable systems Sutherland (2004), may host non-trivial quantum-coherent states, which have been fruitfully studied in recent years. However, these systems require either the presence of quenched disorder (MBL) or fine tuning (Bethe-ansatz integrability). Therefore, it is desirable to find other routes towards extending quantum coherence in many-body systems.

Progress towards extending quantum coherence is tied with developing a more complete understanding of quantum many-body chaos and thermalization D’Alessio et al. (2016). In this direction, recent studies Nahum et al. (2018); Chan et al. (2018); Kos et al. (2018) considered a model of “typical” quantum dynamics generated by applying random unitary operators in a quantum circuit, which allowed to obtain results for the dynamics of entanglement and other physical observables. In such models, thermalization rate should not depend strongly on the initial state. Recent experiments, however, revealed that the dynamics in physical many-body systems may be significantly more complex Greiner et al. (2002); Bernien et al. (2017). In particular, for certain initial states, the chain of interacting Rydberg atoms Bernien et al. (2017) was found to exhibit much slower thermalization or even its absence on the experimentally accessible time scales. Dynamics following from such initial states features persistent periodic revivals of local observables, such as the density of domain walls in the initial state as a function of time.

The observed periodic revivals were explained by “quantum many-body scars” Turner et al. (2018a, b); Ho et al. (2019) – a subset of anomalous eigenstates with strongly non-thermal properties Turner et al. (2018a). Many-body scarring is a generalization of the well-known phenomenon of quantum scars in stadium billiards, where anomalous eigenstates exhibit an increased concentration of probability density around unstable classical periodic orbits Heller (1984). Atypical eigenstates had previously been constructed analytically in the AKLT spin chain Moudgalya et al. (2017, 2018), and it has been suggested that some of these non-thermalizing states may be closely related to the ones in the Rydberg atom chain Lin and Motrunich (2017); Shiraishi (2019). Additionally, various types of non-thermalizing behaviors have since been reported in a number of other models Kormos et al. (2016); Lin and Motrunich (2017); James et al. (2018); Shiraishi and Mori (2017); *Rigol18; Ok et al. (2019); Pai and Pretko (2019); Khemani and Nandkishore (2019); Sala et al. (2019); Bull et al. (2019).

In the context of few-body quantum chaos Haake (2013), the strong dependence of quantum relaxation rate on the initial state more often emerges not from unstable periodic orbits (as in quantum scars), but from the phenomenon of mixed phase space Percival (1973); Bohigas et al. (1993). According to Kolmogorov-Arnold-Moser (KAM) theorem Kolmogorov (1954); *Arnold, weak deformations of an integrable classical system destroy the regular structure of its phase space only in some regions, leading to a coexistence of regular islands with regions of chaotic motion Zaslavsky (2007). The semiclassical limit allows to introduce a notion of regular and chaotic eigenstates that are dominated by the corresponding regions of phase space Percival (1973); Bohigas et al. (1993). The regular eigenstates strongly affect the dynamics. When a quantum system is initialized in a wave packet residing predominantly in the mixed region of the phase space, it exhibits much slower relaxation compared to the case when the quantum evolution begins in the chaotic region of the phase space Haake (2013).

In contrast to few-body quantum systems, in the many-body case one expects that conditions for KAM theorem become very stringent, potentially leading to a quick disappearance of regular regions of phase space already at very small integrability-breaking perturbations. Thus, one may naively expect that quantum many-body systems should display “typical” relaxation, irrespective of the initial conditions. In this work we demonstrate that this intuition is incomplete. A mixed classical phase space leaves an imprint on quantum dynamics in interacting many-body systems, giving rise to slow, atypical thermalization for certain initial conditions. Our starting point is the time-dependent variational principle (TDVP) Dirac (1930), which we use to project quantum dynamics onto a classical nonlinear dynamical system. In what follows, we consider both 1D systems, where we choose the variational manifold to consist of translation-invariant matrix-product states (MPS) Haegeman et al. (2011) with a fixed number of degrees of freedom Bernien et al. (2017); Ho et al. (2019), as well as higher dimensional systems where we use the infinite tensor tree states (TTS) Nagaj et al. (2008). We demonstrate that the resulting dynamical systems generally have mixed phase space.

Regular regions of phase space are governed by stable periodic trajectories, such as the one in Fig. 2 below. These periodic orbits have a vanishing Lyapunov exponent, but may be classified according to their “quantum leakage”, i.e., a measure of discrepancy between TDVP and exact quantum dynamics, which intuitively corresponds to the irreversible entanglement growth. On the one hand, short periodic trajectories, with low quantum leakage, generate many-body revivals in quantum dynamics. As a proof of principle, we identify a trajectory with a 3-site unit cell that gives rise to robust revivals, and also generalize the notion of scars to a 2D square lattice of Rydberg atoms. On the other hand, even the trajectories with high leakage out of the variational manifold leave a measurable signature on the quantum system: initializing the quantum system in the vicinity of these trajectories significantly slows down thermalization compared to generic initial conditions. We independently confirm these findings by numerical simulations based on exact diagonalization and iTEBD Vidal (2007).

Similar to the known examples of few-body systems Bohigas et al. (1993), we expect that regular regions of mixed phase space give rise to non-thermal eigenstates in the many-body quantum system. Thus, mixed phase space would provide a more general mechanism for slow thermalization compared to quantum scars. We establish TDVP as a practical method for finding models with non-ergodic quantum dynamics in arbitrary spatial dimensions, which complements its recent applications to quantum thermalization Leviatan et al. (2017); Hallam et al. (2018). Our results also open the door to generalizing the results on few-body chaos to many-body systems, and for approaching the KAM theorem in quantum systems by utilizing the classical KAM for the TDVP equations of motion. This is distinct from other approaches that rely on broken Bethe-ansatz integrability Brandino et al. (2015) or the absence of quantum resonances in the MBL phase Serbyn et al. (2013); Huse et al. (2014); Imbrie (2016).

The rest of the paper is organized as follows. In the remainder of this section, we briefly introduce the TDVP and the tensor network states of interest. In Sec. II, we show how mixed phase space emerges in the TDVP equations that describe the dynamics of Rydberg atom chains, recently realized experimentally Schauß et al. (2012); Labuhn et al. (2016); Bernien et al. (2017). We demonstrate that TDVP allows to identify a stable trajectory that gives rise to new quantum revivals beyond those that were probed in recent experiments Bernien et al. (2017). In Sec. III we show that higher dimensional bipartite lattices of Rydberg atoms also display mixed phase space and the associated quantum revivals. In Sec. IV we show that quantum leakage can be used to distinguish trajectories that lead to quantum revivals, and to find deformations of the model that improve the revivals originating from a given periodic orbit. In Sec. V, we study the transverse-field Ising model (TFIM) in a longitudinal field, which is a typical example of a thermalizing system. In this case we show that mixed phase space does not give rise to many-body revivals, but leads to a state-dependent thermalization rate. Finally, we conclude with the discussion of open directions and outlook in Sec. VI.

I.1 A brief overview of TDVP

TDVP equations of motion are obtained by extremizing the action Kramer and Saraceno (1981),

[TABLE]

Such action yields a set of nonlinear classical equations of motion Kramer and Saraceno (1981),

[TABLE]

where the variational parameters are real valued, $x_{i}\in\mathbb{R}$. The equations of motion can also be obtained by minimizing the discrepancy between exact quantum dynamics and its projection onto the variational manifold, see Appendix A for a detailed derivation.

TDVP has been succesfully applied to the cases where the variational manifold is spanned by MPS Ueda et al. (2006); Haegeman et al. (2011) as well as finite tensor tree states (TTS) Bauernfeind and Aichhorn (2019). Below we introduce the basics of MPS parametrization, which will be used to describe the dynamics of the Rydberg atom model in Sec. II and the Ising model in Sec.V. In addition, we discuss TTS parametrization that will be used in Sec. III to capture the dynamics of the Rydberg model generalized to higher dimensions.

I.2 Matrix product states and tensor tree states

The simplest tensor network state is the matrix product state. The manifold of translationally-invariant MPS for a system of size $L$ with a unit cell of fixed size $K$ is defined as follows [see Fig. 1(a)]:

[TABLE]

where $A^{\sigma}(\bf x)$ is $\chi\times\chi$ matrix ($\chi$ is the bond dimension) that depends on $N$ real on-site variational parameters, ${\bf x}=(x_{1},\ldots,x_{N})$. The physical indices $\sigma$ label the basis of on-site Hilbert space, which we take to be a spin-1/2 degree of freedom with $\sigma=\uparrow,\downarrow$. The index $m=1,\ldots,L/K$ labels the unit cells, and in what follows we take the thermodynamic limit $L\to\infty$. The state in Eq. (3) is assumed to be normalized, and $V^{L,R}$ denote some boundary vectors the choice of which is not important for our purposes.

The MPS state (3) has $N$ parameters for each of the $K$ sites within its unit cell. Hence, one needs to specify overall $NK$ real parameters, $x_{a}$ with $a=1,\ldots,NK$, to fully fix the state. Below we will use the MPS ansatz restricted to a unit cell of small size and with a small number of variational parameters. These restrictions allow us to obtain analytical equations of motion, making the analysis more transparent. For larger number of parameters, similar analysis could be performed numerically.

As another extension, a close relative of the MPS ansatz are the TTS with higher connectivity – see Fig. 1(b). On the one hand, TTS allow to mimic the topology of two-dimensional and higher-dimensional lattices. On the other hand, the absence of loops still allows for an analytical derivation of equations of motion. In this work we restrict ourselves to TTS which approximate bipartite lattices with a natural choice of two-site unit cell. The manifold of TTS with two inequivalent sites, for a system of size $L$ and connectivity $C$, is defined by generalizing Eq. (3) to the case where $A^{\sigma}(\bf x)$ is a rank-$C$ tensor, with bond dimension $\chi$. The product of matrices in Eq. (3) is extended to contraction of indices on all bonds connecting tensors according to Fig. 1(b). While TTS with fixed bond dimension cannot capture the loop effects of higher dimensional lattices, they can still describe the local physics of weakly entangled states. In Sec. III we find TTS ansatz to be very accurate in capturing the quantum dynamics of certain initial states in dimensions $d>1$.

II Mixed phase space in TDVP dynamics in one dimension

In this Section, we investigate the phase portrait of the classical nonlinear TDVP equations for an MPS manifold of low bond dimension. Below we introduce the model and variational state, and then demonstrate the mixed nature of phase space and its effects on quantum dynamics.

II.1 PXP model

In this section we focus on the PXP model Bernien et al. (2017), describing a 1D chain of Rydberg atoms defined by the Hamiltonian:

[TABLE]

where $\sigma_{x}$ is the standard Pauli matrix, and operator $P_{j}=(1-\sigma_{j}^{z})/2$ projects to $\downarrow$ state on site $j$. The projectors encode the kinetic constraint due to the Rydberg blockade: two neighboring atoms are not allowed to be simultaneously excited. Unless explicitly stated otherwise, we work in the thermodynamic limit, $L\to\infty$, and restrict the Hilbert space to spin configurations without two adjacent up spins, $|\ldots{\uparrow}{\uparrow}\ldots\rangle$. This is the largest connected component of the full Hilbert space for the Hamiltonian in Eq. (4). The PXP model is interacting and non-integrable Turner et al. (2018a), yet its relaxation strongly depends on the initial conditions Bernien et al. (2017). For instance, there is fast thermalization when the system is prepared in a down-polarized state, $\ket{\downarrow\downarrow\ldots}$, while other initial states, such as $|\mathbb{Z}_{2}\rangle=\ket{\uparrow\downarrow\uparrow\downarrow\ldots}$ state, exhibit revivals of local observables Bernien et al. (2017), entanglement entropy Turner et al. (2018a), and even the many-body wave function Turner et al. (2018b).

Below, we apply the ansatz in Eq. (3) with $K=3$ to PXP model. We parametrize the matrices $A^{\sigma}$ by two angles, ${\bf x}_{i}=(\theta_{i},\phi_{i})$:

[TABLE]

The matrix $A^{\uparrow}\propto\sigma^{+}$ satisfies the condition $A^{\uparrow}(\theta,\phi)A^{\uparrow}(\theta^{\prime},\phi^{\prime})=0$, thus effectively imposing the constraint that no two adjacent spins are in $\ket{\uparrow}$ state.

The Hamiltonian in Eq. (7) has particle-hole symmetry in the spectrum: $E\to-E$, since it anti-commutes with ${\cal C}=\prod_{i}\sigma^{z}_{i}$. In this case the system’s equations of motion have a class of solutions with phase angles being stationary, $\phi_{i}=0$, which corresponds to dynamics restricted to a flow-invariant subspace in the language of dynamical systems. This subspace arises from the presence of particle-hole symmetry and time-reversal invariance of the PXP Hamiltonian in Eq. (4), see Appendix A. Note that $\langle{H}\rangle=0$ vanishes identically when $\phi_{i}=0$, thus imposing no additional constraint. Restricting to $\phi_{i}=0$ subspace for $K=2$ unit cell, Ref. Ho et al. (2019) found an unstable periodic trajectory in a two-dimensional flow invariant subspace. This trajectory is intimately connected with the revivals from the initial Neél state $\ket{\mathbb{Z}_{2}}=\ket{\downarrow\uparrow\downarrow\uparrow\ldots}$.

However, a two-dimensional phase space is non-generic and does not allow for chaos: in two dimensions any periodic trajectory fragments the phase space into dynamically isolated regions. Thus, we focus on ansätze with 3 or more parameters, which have regular and chaotic regions coexisting in phase space. We will show below that such behavior, known as “mixed phase space”, is robust and persists for various deformations of the PXP model. We will then investigate the implications for quantum dynamics, finding new initial states that exhibit revivals and slow thermalization.

II.2 Revivals in a three-site unit cell

Higher-dimensional variational phase space naturally arises when unit cell size is increased, while particle-hole symmetry is kept intact. This describes the dynamics of initial states that are periodic in space, with period $K=3$ (or more) lattice sites. An experimentally relevant example of such initial states is $|\mathbb{Z}_{3}\rangle=|{\uparrow\downarrow\downarrow}\ldots\rangle$, which can be experimentally prepared in a Rydberg system with larger blockade radius Bernien et al. (2017). The resulting equations of motion for $\theta_{i}$, $i=1,2,3$, obtained from the TDVP projection of the PXP model can be found in Appendix A.2.3.

In order to analyze the dynamical flow, we consider its Poincaré section. The flow generates a discrete mapping known as the Poincaré map Strogatz (2014) that maps a given point $(\theta_{1},\theta^{*}_{2},\theta_{3})$ on the chosen hyperplane $\theta_{2}=\theta_{2}^{*}$ to a position $(\theta^{\prime}_{1},\theta^{*}_{2},\theta^{\prime}_{3})$ where the trajectory intersects this plane again. Periodic trajectories correspond to stationary points of the Poincaré map, while in the case of chaotic behavior the system returns to the same plane at a location that is generally far away from the previous encounter. Successively iterating the Poincaré map for many different initial conditions yields the Poincaré section, shown in Fig. 2(a).

Fig. 2(a) reveals the phase portrait characteristic of dynamical systems with mixed phase space. Cross-sections of at least four large stable KAM tori are evident, surrounded by a number of smaller tori and chaotic regions. The analysis of periodic orbits reveals a stable shortest-period orbit denoted by $\star$ in Fig. 2(a) and located at $(\theta_{1},\theta_{2},\theta_{3})=(0.8356\pi,0,0.1923\pi)$. Other large tori surround the orbits related to this particular orbit via symmetry transformations.

Next, we analyze the local observables and entanglement entropy of the $\star$ periodic orbit from Fig. 2(a). Fig. 2(b) shows the time evolution of the entanglement entropy for the variational solution, which indicates that the dynamics on this orbit never passes through a product state. To see this, note that when the entanglement for one cut, say $S_{12}$, vanishes, the entanglement for a different cut is non-zero, e.g., $S_{23}>0$. Interestingly, this short-range entanglement improves the fidelity of revivals in the corresponding quantum system. Indeed, the closest product state to this trajectory is $|\mathbb{Z}_{3}\rangle=|{\downarrow\uparrow\downarrow}\ldots\rangle$. However, preparing the quantum system in a slightly entangled state positioned on the periodic TDVP orbit gives rise to more robust revivals. The improvement of revivals is illustrated in Fig. 2(c). In the top panel, we show the negative logarithm of fidelity,

[TABLE]

normalized by the system size, $-(\ln F)/L$, that quantifies the fidelity decay per spin. The minima in this quantity correspond to the best fidelity revivals of all finite-size subsystems. These minima display a regular pattern, with smaller values (corresponding to higher return amplitude) when the system is initialized on the TDVP trajectory compared to $|\mathbb{Z}_{3}\rangle$ product state. The bottom panel of Fig. 2(c) shows that the growth of entanglement is notably slower for the initial state positioned on the TDVP orbit, consistent with better fidelity revivals.

II.3 Effect of regular regions of phase space

Above we considered the dynamics of a quantum system initialized on the periodic orbit. However, any such classically stable orbit is surrounded by a regular region of phase space where TDVP dynamics is constrained to a torus. In order to investigate the effect of this region on quantum dynamics, we consider a family of initial states $\ket{\psi(c)}$, where the parameter $c\in[0,1]$ interpolates along the dashed line in Fig. 2(a) between the periodic orbit and the point $\theta_{1}=0,\theta_{3}=\pi$.

Launching classical dynamics for different initial conditions, we clearly see the termination of the regular region in the Fourier spectrum of local observables, Fig. 3(a). Inside the regular region, the Fourier spectrum is strongly peaked. In contrast, upon entering the chaotic sea when $c\approx 0.42$, one observes a rapid drop in the amplitude of the maximal Fourier harmonic component that signals a transition to chaotic dynamics with the continuous Fourier spectrum. For the same set of initial conditions, quantum dynamics shows a crossover between non-ergodic behaviour, with slow entanglement growth, and thermalizing dynamics, as witnessed by the behavior of entanglement entropy, see Fig. 3(b). The changes in the behavior of the quantum system are much more smooth: intuitively, one may think of a quantum system initialized in a finite region of the phase space, with the size of this region being determined by quantum fluctuations. This smoothens out any abrupt transition between the non-ergodic and thermalizing behaviors. Nevertheless, Fig. 3(c) shows that the maximal sensitivity of quantum dynamics to a change in initial conditions is achieved near the border of the regular region.

In addition to entanglement entropy, we observed similar behavior of inverse participation ratios (IPR) of local density matrices, measured at the time when the dynamics has saturated (not shown). The IPRs are found to be low in the regular region, consistent with the expectation that the system does not explore the full configuration space when launched there. In contrast, the participation ratios rapidly increase near the border of the regular region and attain thermal values outside of it.

The classical dynamical system also has many longer orbits surrounded by “thinner” tori – for example, see the orbit marked by the blue diamond in Fig. 2(a). We observe that these orbits do not give rise to long-lived oscillations in exact dynamics of the PXP model, and they are characterized by rapid entanglement growth – see the dashed line in the bottom panel of Fig. 2(c). Since numerous stable orbits appear in the Poincaré sections, we need to understand which of these orbits give rise to strong quantum revivals in exact quantum dynamics. We will return to the question of how to distinguish orbits in Sec. IV, after first considering the dynamics and mixed phase space on higher-dimensional lattices in the next section.

III Revivals and mixed phase space in higher dimensions

In this section we introduce the PXP model in higher dimensions and apply the TTS ansatz to capture its dynamics. We will focus on a 2D square lattice, and compare the predictions of the TTS dynamics against exact diagonalization. We will also discuss generalizations of these results to any bipartite lattice in an arbitrary number of spatial dimensions.

III.1 PXP model in higher dimensions

We start by introducing a constrained model on the square lattice, inspired by the 1D PXP model:

[TABLE]

where the pair of indices $i,j$ is used to label lattice sites. In addition to the PXP term, we have also include a possible perturbation in terms of chemical potential $\mu_{z}$, which couples to the Rydberg excitation density operator $n=(1+\sigma^{z})/2$. The presence of four projectors allows a given atom to get excited only if all neighboring atoms are in the ground state. Similarly to the PXP model in 1D, we will restrict our attention to the largest connected subspace of the Hilbert space, where no adjacent atoms are excited.

After generalizing the quantum model to the 2D case, we turn to the variational state that can capture the quantum dynamics. Inspired by the ansätze considered above for the 1D lattice, we introduce tree tensors $A^{\sigma}_{i_{1},\ldots,i_{4}}$, which forbid excitations on adjacent sites. Relegating the detailed discussion to Appendix B, here we list only the non-zero elements of the tensor $A^{\sigma}$ that is parametrized by two angles, $(\theta,\phi)$:

[TABLE]

where $i_{1,\ldots 4}$ are virtual indices, and $I=\sum^{3}_{j=1}i_{j}$ counts the number of “excitations” on the three legs of the tensor. The fourth leg of the tensor tensor produces an “output $1$-bit” on the last virtual index, $i_{4}$, if the corresponding site is excited and [math] if it’s in the ground state. In the Appendix B we give the form of the ansatz for general connectivity $C$ of the tree. When $C=2$, which corresponds to the 1D limit, this ansatz reproduces Eq. (5).

Projection of quantum dynamics generated by Eq. (7) onto the TTS manifold with a 2-site unit cell results in the dynamical system of four variables that come in two canonically conjugate pairs, $(\theta_{1},\phi_{1})$ and $(\theta_{2},\phi_{2})$. The general form of equations of motion reads:

[TABLE]

where the square lattice case (ignoring loops) corresponds to $C=4$. The equations for $\dot{\theta}_{2},\dot{\phi}_{2}$ can be obtained by substitution $1\leftrightarrow 2$. In addition, the energy density,

[TABLE]

is a conserved quantity, provided the system satisfies the equations of motion, Eqs. (9). The existence of a conserved quantity effectively reduces the dimensionality of the phase space to three dimensions, as dynamics is restricted to constant energy surfaces.

III.2 Generalizing scars to arbitrary bipartite lattices

When $\mu_{z}=0$ and one restricts to a flow-invariant subspace $\phi_{1,2}=0$, Eqs. (9) reduce to a two-dimensional flow. This flow has an isolated periodic trajectory that generalizes the trajectory found in Ref. Ho et al. (2019) to an arbitrary bipartite lattice. This trajectory is responsible for the revivals of a state in which all atoms on the first sublattice are in the Rydberg state, $\uparrow_{1}$ and all atoms on the second sublattice are in the ground state, $\downarrow_{2}$. In Fig. 4(a) we show that quantum dynamics of a local observable $\braket{\sigma^{z}}$ on the $4\times 4$ square lattice agrees well with TDVP dynamics based on $C=4$ TTS. Moreover, we compare the TDVP dynamics of a wavefunction satisfying the kinetic constraints and having the same weights as the tree parametrization, see Eq. (54) in App. B.1, to the tree dynamics and find perfect agreement (not shown). The fact that TDVP dynamics that takes into account loops agrees with TTS TDVP demonstrates that presence of loops is not important for local observables.

As connectivity increases, the isolated periodic orbit of the TTS manifold is expected to become more accurate, and effectively exact at $C\rightarrow\infty$. In this limit, the sites that are in their ground state cannot flip unless all their neighbours are also in their ground state. The periodic trajectory can be understood by studying the $C\rightarrow\infty$ limit of Eq. (9). At the first quarter of the period, the excited sites on the first sublattice relax $\uparrow_{1}\to\downarrow_{1}$ while their neighbors remain in the ground state due to infinite connectivity. At the second quarter of the period, the second sublattice undergoes the transition $\downarrow_{2}\to\uparrow_{2}$. The half period over which two sublattices exchange their state is given by $\pi$ in our notations, so the total period of the orbit is $T_{C\to\infty}=2\pi$. Figure 4(b) shows that period of the orbit approaches its $C\to\infty$ limit approximately as $T_{C}=2\pi-2.4/\sqrt{C}$.

III.3 Phase space and revivals in higher dimensions

When detuning $\mu_{z}\neq 0$, the dynamics generated by Eq. (9) takes place on constant energy surfaces in four-dimensional space. In order to visualize such dynamics, we fix $\theta_{1}$ variable and show the resulting Poincaré sections in Fig. 5(a) [the precise value of $\theta_{1}$ is not important, and we have chosen $(\theta^{*}_{1}=1.25\pi,\dot{\theta}_{1}<0)$]. Energy conservation results in complicated surfaces in the space $(\phi_{1},\theta_{2},\phi_{2})$ that cannot be globally projected. Therefore, we show a small part of the Poincaré section that can be projected onto $(\phi_{1},\theta_{2})$ plane. We find that the point $(\phi_{1},\theta_{2})\approx(2.921,2.981)$ is a stationary point of the Poincaré map that corresponds to a stable periodic trajectory. We observe that this periodic trajectory is surrounded by circle-shaped contours which are the intersections of KAM tori by our section plane.

A key question concerns the correspondence between regular variational dynamics and exact quantum dynamics. Fig. 5(b) shows this comparison for local observables $\langle\sigma_{i}^{x,z}\rangle$, for the system initialized on the periodic trajectory (indicated by $\star$ in Fig. 5). It demonstrates that the classical periodic trajectory gives rise to pronounced quantum revivals; moreover, the exact quantum dynamics agrees with the TDVP predictions up to $t\lesssim 10$. We note that the slight disagreement at $t=0$ is caused by the difference between the structure of a TTS (that cuts loops and introduces redundant degrees of freedom) and true many-body wave function on a two-dimensional lattice (see App. B.2). In order to visualize the variational (classical) dynamics, Fig. 5(c) shows the evolution of the trajectory of an individual spin $(\langle{\sigma^{x}_{1}(t)}\rangle,\langle{\sigma^{y}_{1}(t)}\rangle,\langle{\sigma^{z}_{1}(t)}\rangle)$ on the Bloch sphere as perturbation $\mu_{z}\neq 0$ is turned on. The second spin in the unit cell performs similar oscillations shifted in phase by $\pi$.

For $\mu_{z}=0$ the trajectory gives oscillations of the spin between $\uparrow$ and $\downarrow$ states, and the expectation value $\langle{\sigma^{x}(t)}\rangle$ is zero at all times. In agreement with physical intuition, non-zero chemical potential acts as a “magnetic field” in the $z$-direction, tilting the plane in which the spin rotates. Surprisingly, in addition to the tilt, the trajectory slightly winds around the south pole of the Bloch sphere. This winding reflects the fact that we are dealing with an interacting system rather than with a precession of independent spins, which shows that this orbit cannot be obtained as a smooth deformation of the unstable periodic orbit found in the case of $\mu=0$.

In the entire range of $\mu_{z}\in(0,0.65]$ the regular regions remain pronounced in the phase space of the system in Eqs. (9), thus demonstrating the persistence of mixed phase space. The deformation of the trajectory away from the poles upon increasing $\mu_{z}$ suggests that the most stable revivals occur for initial states that are different from the Neél product state, $\mathbb{Z}_{2}$, by having some amount of entanglement.

IV Characterizing trajectories by leakage

Variational dynamics of the PXP model and its deformations, studied above, revealed mixed phase space with multiple stable orbits. The Lyapunov exponent, which is often used to quantify chaos in TDVP Hallam et al. (2018) and semiclassical Scaffidi and Altman (2017) dynamics, vanishes for all these orbits and thus cannot be used to distinguish them. In this section we show that quantum leakage can be used instead as a heuristic for distinguishing trajectories that give rise to revivals in exact quantum dynamics.

IV.1 Rate of leaving the variational manifold

TDVP approach accurately approximates short-time quantum evolution of initial states in the chosen MPS manifold, but at long times the errors will grow because exact quantum dynamics generally brings the system out of the variational manifold. Intuitively, the TDVP error for a given initial state is linked to the rate of entanglement growth for that state, since MPS with low bond dimension can only represent weakly entangled states. When the entanglement grows significantly, the state requires an MPS description with a larger bond dimension that lies outside the variational manifold. Making this connection mathematically precise is challenging and is beyond the scope of this work. Instead, here we characterize the error by quantum leakage defined as the instantaneous rate at which the exact quantum wave function leaves the variational manifold Haegeman et al. (2011); Ho et al. (2019):

[TABLE]

In Appendix C we discuss the computation of this quantity in the TDVP framework. This calculation involves the square of the Hamiltonian operator, thus effectively it contains information that goes beyond the TDVP equations of motion. The normalization factor is chosen such that $\gamma$ assumes a finite value in the thermodynamic limit.

In the previous section we demonstrated that TDVP equations of motion may have periodic trajectories, ${x}_{a}(t)={x}^{(0)}_{a}(t)$, with a period $T$. If the quantum system were following the TDVP dynamics exactly, that would imply persistent oscillations in local observables and perfect revivals of the many-body quantum fidelity, i.e., $F(t=nT)=1$ in Eq. (6), where the initial state can be any MPS wave function that belongs to the periodic trajectory, $|\Psi\rangle=|\psi({x}^{(0)}_{a}(t))\rangle$. The disagreement between quantum dynamics and TDVP will generally preclude such perfect revivals (except for models with perfect scars Choi et al. (2018)). However, we conjecture that one can obtain a lower bound on the many-body fidelity revival using quantum leakage.

Assuming an extensive scaling of the fidelity, $F(T)=e^{-f_{T}L}$ at some fixed short time and $L\to\infty$ limit, and using the large system size limit of Eq. (11), we posit the following upper bound on $f_{T}$,

[TABLE]

that translates to a lower bound on the fidelity. This bound can be justified in the limit of small leakage for a finite size system, see Appendix C. However, its extension to the thermodynamic limit is non-trivial and at present it remains a conjecture. Below we present a specific example which demonstrates the usefulness of quantum leakage and provides a test of the bound.

IV.2 Leakage and revivals criterion

We study the evolution of the Poincaré sections and corresponding trajectory under the deformation of the PXP model that makes it more thermalizing. Specifically, we consider the following deformation,

[TABLE]

which is the lowest order term that induces spin flips (expected to facilitate thermalization) while maintaining particle-hole symmetry.

The presence of particle-hole symmetry in $H=H_{\text{PXP}}+H_{\mu_{3}}$ allows us to use the same ansatz as in Sec. II.2 with three variational parameters $\theta_{1,2,3}$. The resulting TDVP equations of motion are again cumbersome and can be found in Appendix A.2.3. Analysis of these Poincaré sections reveals a pronounced asymmetry between the effect of this deformation for $\mu_{3}=\pm 0.25$. The deformation with $\mu_{3}<0$ increases the size of regular regions in the phase space, see Fig. 10. In contrast, positive $\mu_{3}$ causes the regular regions of the phase space to shrink in size and the phase portrait looks progressively more chaotic (Fig. 11). Exact diagonalization shows a similar asymmetry between positive and negative $\mu_{3}$: for $\mu_{3}>0$ thermalization, diagnosed via level repulsion and average eigenstate entanglement entropy, becomes more pronounced.

In order to quantify the effect of the deformation, we use the expression from Appendix C for $\gamma(\{x_{a}\})$ to calculate the integrated leakage, $\Gamma_{T}$, defined in Eqs. (11)-(12), along the trajectory. Fig. 6 shows that the leakage attains a minimum at a small negative value $\mu^{\text{TDVP}}_{3}\approx-0.11$. For larger negative values of $\mu_{3}$ the leakage starts to increase again. Although the leakage overestimates the suppression of the fidelity revivals, the exact diagonalization results for $f_{T}=-\ln F(T)/L$ follow qualitatively the same trend as $\Gamma_{T}$. The $f_{T}$ obtained from exact quantum dynamics achieves a minimum at $\mu_{3}^{\text{ED}}\approx-0.06$, which corresponds to the best fidelity revivals. The optimal value of the deformation, $\mu_{3}^{\text{ED}}$, is approximately twice smaller compared to the TDVP prediction, $\mu^{\text{TDVP}}_{3}$. Nevertheless, TDVP may be used to find the best operators suitable for stabilizing quantum many-body revivals Choi et al. (2018). Moreover, in contrast to the revivals from $|\mathbb{Z}_{2}\rangle$ state Ho et al. (2019), our results show that generally one must account for the effect of the deformation on the initial state.

Eventually, for large absolute values of $\mu_{3}$, we observe that the quantity $f_{T}$ is no longer well defined, since the data for different system sizes do not collapse onto each other in Fig. 6. We attribute this to fast thermalization at such deformation parameters. At these values, quantum dynamics no longer exhibits many-body revivals, and simultaneously $\Gamma_{T}$ is on the order of unity. Thus, we propose

[TABLE]

as a tentative criterion for when a periodic trajectory ceases to result in revivals. This criterion was previously heuristically conjectured in Ref. Ho et al. (2019) in the context of dynamics of $Z_{2}$ state. With this criterion, it is natural to test other periodic trajectories visible in the Poincaré sections of the TDVP dynamics. In Appendix C.3 we present such an analysis for $\mu_{3}=0.25$, which shows that the majority of trajectories have leakage that is larger than one, thus not giving rise to quantum revivals.

IV.3 Floquet exponents and quantum leakage

We demonstrated the utility of quantum leakage in classifying trajectories and predicting which trajectories give rise to fidelity revivals. The leakage gives a particular measure of instability of the dynamics with respect to external degrees of freedom – it quantifies how quickly quantum evolution leaves the variational manifold. At the same time, the TDVP equations of motion can be classified by an internal measure of instability of the TDVP dynamics. To that end, in the case of chaotic motion, one uses the Lyapunov exponent, whereas in the case of periodic orbits the Floquet exponents provide a conceptually similar but reparametrization-invariant measure Cvitanović et al. (2016), see App. C.3 for details.

It is natural to ask if the Floquet exponent and quantum leakage are directly related to each other. The answer is negative: Fig. 6(inset) shows how a non-zero Floquet exponent $\lambda$ emerges for $\mu_{3}\geq 0.255$. However, the non-zero value of $\lambda$ does not have any influence on the leakage and this unstable trajectory still gives rise to fidelity revivals for $\mu_{3}\in(0.255,0.4]$. This example also suggests that mixed phase space provides a more generic mechanism for weak ergodicity breaking. Stable trajectory in the center of a regular region of mixed phase space that gives rise to “regular eigenstates” Percival (1973); Bohigas et al. (1993) may become unstable upon a deformation of the Hamiltonian, leading to the appearance of quantum scars Turner et al. (2018a); Ho et al. (2019) which are quantum eigenstates affected by the unstable trajectories.

V Non-universal thermalization

Above we focused on periodic trajectories that arise in classical TDVP equations, with small quantum leakage. These trajectories give rise to spectacular revivals of the fidelity and slow thermalization (or even its absence, as is the case for the deformed PXP model Choi et al. (2018)). In this Section we consider the opposite situation where quantum leakage is strong. In this regime, TDVP fails to capture the dynamics of local observables beyond times of order $O(1)$. Nevertheless, we show that mixed phase space still leaves an imprint on quantum dynamics, leading to non-universal thermalization.

V.1 Transverse-field Ising model

We study the quantum Ising model with transverse and longitudinal fields (TFIM), a popular model for investigating eigenstate thermalization, defined by a Hamiltonian,

[TABLE]

with fixed values of parameters $(J_{z},h_{z},h_{x})=(1,0.4,1)$. Thermalization in this model was established for somewhat different parameters $(J_{z},h_{z},h_{x})=(1,0.8090,0.9045)$ in Ref. Kim et al. (2014); we choose a smaller value of $h_{z}$ to facilitate the TDVP analysis. Very similar behavior is also observed for other values of couplings, $h_{x}=0.9$ and $1.1$.

In order to capture the dynamics, we use the following MPS ansatz with bond dimension $\chi=2$,

[TABLE]

that depends on four parameters, $(\theta,\phi,\xi,\chi)$, per spin. Such a choice of the ansatz is inspired by “trotterization” of the evolution operator $e^{-iH_{\text{TFIM}}\delta t}\approx e^{-iH_{0}\delta t}e^{-iH_{1}\delta t}$, where $H_{0}$ and $H_{1}$ are the transverse-field term and remaining terms in $H_{\text{TFIM}}$, respectively. We restrict ourselves to translation-invariant states, setting the unit cell size $K=1$. We note the such a choice of tensors does not result in a normalized wave function. Hence, the derivation of the equations of motion is more complicated and we follow the general procedure outlined in Appendix A.3. We perform the numerical integration of the equations of motion and obtain the phase portraits of the TDVP dynamics.

V.2 Dependence of dynamics on the initial state

We consider the TDVP and exact quantum dynamics for initial states specified by the ansatz in Eq. (16), with fixed energy density $\langle{H}\rangle/L=0.19$. This corresponds to initial states at high but finite temperature in the thermodynamic limit. In this case, the TDVP dynamics reveals the presence of mixed phase space (see Fig. 12 in Appendix D.1). However, the integrated leakage for the shortest stable trajectory with the period $T=2.097$ is close to one, $\Gamma_{T}=0.57$. Fully consistent with such values, the trajectory does not result in oscillations of the fidelity and all local observables relax by the time of order $t\sim 2$.

In order to study the effect of the periodic trajectory, we consider an ensemble of initial states specified by the ansatz in Eq. (16), with the same energy density and a small initial value of the entanglement entropy $S_{\text{ent}}\in[(2/3)S_{0},(4/3)S_{0}]$, where $S_{0}\sim 0.08$ is the minimal value of entanglement entropy of the given periodic orbit. Since all these states have the same energy density, nearly the same initial value of the entanglement, and are translationally invariant, we expect that after short time sufficient to equilibrate local observables, the entanglement dynamics should be independent of the initial state. Fig. 7 shows that these expectations are not correct, and the entanglement at late times is still strongly influenced by the vicinity to the periodic orbit. At the same time, all these initial conditions result in the same saturation value of entanglement for finite subsystems, see Fig. 13.

Thus we conclude that these initial conditions not only have different short-time behavior of the entanglement, but are also characterized by a different velocity of the entanglement spreading. In Appendix D.1 we demonstrate that, for the present choice of parameters, the eigenstates of this model have properties fully consistent with ETH, hence such dynamics cannot be explained by the existence of anomalous eigenstates. We note that a similar phenomenon was observed in a Floquet version of TFIM with somewhat weaker integrability-breaking field Bertini et al. (2018).

The correlation between TDVP leakage integrated over a short time and the entanglement spreading rate – see Fig. 7(inset) – provides further support for the relevance of TDVP dynamics for entanglement spreading even at late times. The statistical analysis gives a correlation coefficient of 0.53 with $p$-value of 0.0007, suggesting that correlations are statistically significant. Appendix D.2 reveals similar phenomenology in the case of a strongly deformed PXP model, where the different velocity of entanglement spreading and correlation between leakage and entanglement is even more apparent. The existence of such a correlation suggests that it may be possible to use the leakage as an upper bound on the value of entanglement at long times.

Finally, we note that different rates of entanglement spreading observed in a thermalizing system is well-known in the context of integrable models Calabrese and Cardy (2005); Alba and Calabrese (2017, 2018). In the latter case the dependence of entanglement dynamics on the initial state stems from the presence of multiple conserved quantities. In the present case we are dealing with a non-integrable TFIM; however, it is tempting to conjecture the existence of one (or a few) approximately conserved quantities related to regular regions of the phase space. Also, a possible relation between such entanglement spreading and the existence of slowly growing operators found in the TFIM Kim et al. (2015) remains an intriguing open question.

VI Discussion

To summarize, we examined the relation between exact quantum dynamics and an effective classical non-linear system, obtained by projecting quantum dynamics via TDVP onto the restricted MPS manifold. This approach is different from the time-dependent mean field (which can be viewed as a particular case of $\chi=1$ MPS) and semiclassical treatments used in few-body quantum chaos in that it incorporates short-range entanglement. We demonstrated the relevance of mixed phase space, identified in the TDVP dynamics, for the exact dynamics of quantum many-body systems. We used quantum leakage to distinguish two qualitatively different situations:

1. (i)

In the small-leakage regime, stable periodic trajectories lead to long-time revivals of the many-body fidelity and oscillations of local observables. This provides a more general mechanism of weak ergodicity breaking compared to quantum many-body scars in an arbitrary number of spatial dimensions. 2. (ii)

In the strong-leakage regime, exactquantum dynamics does not follow the TDVP predictions, but the mixed nature of phase space influences the rate of entanglement growth.

In the weak leakage case, (i), we establish the TDVP method as an indispensable tool in searching for new stable periodic trajectories leading to fidelity revivals. This is demonstrated in Secs. II and III by finding revivals in one-dimensional PXP model and its higher-dimensional generalizations. In addition to revivals, the entire regular region of the phase phase space was shown to have influence on quantum dynamics. Moreover, this regular region is expected to give rise to special regular eigenstates in the many-body spectrum (we use the term “regular eigenstates” to distinguish the atypical eigenstates that originate from the regular regions in the phase space from “scars” that come from unstable periodic trajectories). The detailed investigation of these states is left to future work. Next, in Sec. IV we demonstrated that the revivals in the PXP model are much more robust to deformations of the Hamiltonian than was previously reported in the literature Turner et al. (2018a, b). However, in the course of the deformation one has to follow the periodic trajectory that is influenced by the deformation, which can be inferred from TDVP.

In the case of strong leakage, (ii), considered in Sec. V, the thermalizing TFIM is also shown to have mixed phase space. Initializing the system in the regular region of the phase space near the stable periodic trajectory yields a slower rate of entanglement growth compared to other states with the same energy density. These results show that the low bond dimension TDVP, despite its failure to capture rapid entanglement growth, can be useful in searching for slowly thermalizing initial conditions in interacting quantum systems.

It is important to keep in mind that TDVP, in general, gives a multitude of possible ways to map a quantum many-body system onto a classical dynamical system, see Fig. 1. While some TDVP ansätze may find unstable trajectories, thus motivating the term “quantum scars” Ho et al. (2019), more generic embeddings may result in mixed phase space. Hence it remains an open question which of the atypical eigenstates observed in different models Turner et al. (2018a); Ho et al. (2019); Turner et al. (2018b); Khemani et al. (2019); Iadecola et al. (2019); Lin and Motrunich (2018); Moudgalya et al. (2017, 2018); Ok et al. (2019); Pai and Pretko (2019); Khemani and Nandkishore (2019); Sala et al. (2019); Bull et al. (2019) correspond to many-body scars, and which to regular eigenstates.

In addition to establishing the methodological utility of low bond dimension TDVP as a diagnostic of anomalous thermalization, our results pose questions related to the consequences of mixed phase space for the spectrum, level statistics and other indicators of thermalization in quantum many-body systems. As we mentioned in the introduction, mixed phase space is known to influence the semiclassical limit of few-body quantum systems Bohigas et al. (1993). In this case an intermediate spectral statistics is known to arise Berry and Robnik (1984); Prosen and Robnik (1993), and the system has an increased stability with respect to noise when it is prepared in the regular region of phase space Haake (2013). Our work suggests that TDVP can be used to generalize these results to the case of quantum many-body chaos. However, we note that this generalization may be highly non-trivial as quantum many-body eigenstates in the zero momentum sector may see imprints from mixed phase space in TDVP ansätze with different unit cell choices, as in Fig. 1. These studies are of practical importance, as they can result in ways of delaying thermalization and increasing the stability of a quantum many-body system to noise. Thus it is interesting to explore different mechanisms that can lead to low-leakage trajectories. Two prospective directions include weakly perturbed mean-field models (e.g. TFIM with small $J\ll h_{x},h_{z}$) and models with constrained Hilbert spaces. In addition, it is interesting to check for the effect of mixed phase space in other methods used to treat interacting systems, such as Gutzwiller projected dynamics Schiró and Fabrizio (2010) and MPS-based DMFT Wolf et al. (2014).

In a complementary direction, TDVP was recently proposed as a method to capture thermalization. In particular, Ref. Leviatan et al. (2017) demonstrated that physical properties saturate in the case of TFIM when bond dimension $\chi>2$. In addition, Ref. Hallam et al. (2018) studied the spectrum of Lyapunov exponents in the case of high bond dimension TDVP. Our work shows that in the case of low bond dimension one may encounter non-chaotic behavior in the TDVP dynamics, limiting the utility of the Lyapunov exponent in such cases. On the one hand, one may expect that mixed phase space does not persist for large bond dimensions, as increasing $\chi$ increases the dimension of the phase space where the TDVP dynamics occurs, making it more susceptible to chaos. On the other hand, for local Hamiltonians the Lieb-Robinson bound Lieb and Robinson (1972) suggests that a small bond dimension is sufficient to capture quantum dynamics at early times. Thus, it is important to understand how the mixed phase space evolves upon including additional MPS parameters that describe longer-range entanglement. Specifically, one can embed one of small-$\chi$ MPS ansätze considered here into an MPS with a larger bond dimension and track the behavior of the variational parameters that correspond to longer-range entanglement.

At the more phenomenological level, recent works established numerous examples of anomalous thermalization. In particular, these include slow thermalization emerging from confinement Kormos et al. (2016); James et al. (2018), the existence of “slow operators” in thermalizing models Kim et al. (2015), dynamical phase transitions Heyl et al. (2013), dependence of dynamics on the initial state in Floquet systems Bertini et al. (2018), mixed phase space in quantum maps Moudgalya et al. (2019), etc. It would be interesting to understand if all or some of these phenomena can be interpreted using the framework developed here and its straightforward extension to Floquet systems.

Finally, TDVP formulated in the MPS basis could provide a novel pathway to a quantum-mechanical analogue of the KAM theorem. The robustness of quasi-local integrals of motion in MBL phases Serbyn et al. (2013); Huse et al. (2014); Imbrie (2016) provides an established particular case of the quantum KAM Abanin et al. (2018). In a different setting, Ref. Brandino et al. (2015) demonstrated the existence of quasi-conserved quantities in Bethe-ansatz integrable systems deformed by an integrability-breaking perturbation. These works argue for a quantum KAM without making reference to its classical counterpart. By contrast, our results suggest that one may use classical KAM as a route towards its quantum counterpart. Indeed, recent work Choi et al. (2018) conjectured that a weak deformation of the PXP model may lead to perfect many-body revivals in the thermodynamic limit. It remains an open question if such a deformation of the PXP Hamiltonian can allow quantum dynamics to exactly follow the TDVP trajectory. More broadly, searching for deformations of quantum models that reduce quantum leakage in the TDVP dynamics or yield integrable TDVP dynamics may provide a specific route to quantum KAM that is expressed in terms of the KAM for the corresponding classical system.

Acknowledgements.

We acknowledge useful discussions with E. Altman, B. Budanur, P. Calabrese, A. Green, F. Pollmann, and T. Prosen. The work of D.A. is supported by the Swiss National Science Foundation. C.J.T. and Z.P. acknowledge support by EPSRC grants EP/P009409/1, EP/R020612/1 and EP/M50807X/1. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data.

Appendix A Variational ansätze, equations of motion, and leakage

In this Appendix we provide details on the derivation of TDVP equations of motion used in the main text. We begin by setting up the general framework that gives an efficient way of obtaining the equations of motion. Next, we apply this framework to derive equations of motion for the (deformed) PXP model and TFIM which are analyzed in the main text.

A.1 General framework

A.1.1 Notations

In what follows we use TDVP formulated for a wave function that is represented in MPS form. We note that this formulation also includes dynamical mean-field theory as a particular case if one restricts to the product-state form for the variational wave function. While the general theory of TDVP in MPS manifolds was obtained in Ref. Haegeman et al. (2011), here we aim to provide a simpler recipe which can be used to obtain analytical form of equations of motion for states with a small number of variational parameters.

We consider translationally-invariant matrix-product state (MPS) ansatz for the wave function, with the size of the unit cell being fixed at $k$, see Eq. (3) and Fig. 1. Tensor $A^{\sigma}_{\alpha\beta}({\bf x})$ has one physical index, $\sigma=\uparrow,\downarrow$ for spin-1/2 degrees of freedom, and two virtual indices $\alpha,\beta=1,\ldots\chi$, where $\chi$ is the (fixed) bond dimension. This tensor depends on a set of $N$ real variational parameters, ${\bf x}\in\mathbb{R}^{N}$. While our Hamiltonian is translationally invariant, we consider initial states that break translational symmetry. More specifically, we allow for initial states with a unit cell of size $K$. Thus, the complete set of variational parameters consists of a set $\{{\bf x}_{i}\}$ with $i=1,\ldots K$, where $i$ labels sites within the unit cell. In order to simplify the notations, in the remainder of Appendix we suppress the dependence of $A$ on ${\bf x}_{i}$, $A_{i}\equiv A({\bf x}_{i})$.

In order for the state $\ket{\psi}$, Eq. (3), to have system-size independent norm, its unit cell transfer matrix must be non-degenerate with the largest eigenvalue $\lambda_{\text{max}}=1$. The unit cell transfer matrix is obtained from a product of individual one-site transfer matrices,

[TABLE]

over the entire unit cell:

[TABLE]

In the above notations all transfer matrices operate in the $\chi^{2}$-dimensional “double” virtual space labeled by a pair of indices [e.g., $\alpha\gamma$ in Eq. (17)]. In what follows, vectors in the physical Hilbert space are denoted as $\ket{*}$ and vectors of the double virtual space by $|*)$. The dominant left/right eigenvectors of the unit-cell transfer matrix $T_{\text{u.c.}}$ that correspond to an eigenvalue with maximal norm, $\lambda_{\text{max}}=1$ are labelled as $|L)$ and $|R)$. The boundary tensors $v,u$ in Eq. (3) do not affect the calculations in the thermodynamic limit as long as $(v^{{\dagger}}v|R)\neq 0$, $(L|uu^{{\dagger}})\neq 0$.

A.1.2 Action and gauge choice

The TDVP equations are obtained by extremizing the action (1) Haegeman et al. (2011). The symmetric choice of the Lagrangian is invariant under time-dependent unitary transformation $\ket{\tilde{\psi}}=U\ket{\psi}$, where $U=e^{ia(t)}$ is a pure phase. Such a phase is fixed by requiring the system to satisfy the “mean” Schrödinger equation,

[TABLE]

where the phase $a(t)$ is real since the norm of the wave function is time-independent. In what follows we show that this choice of global phase effectively removes all disconnected correlations — see Eq. (61) — from both the equations of motion and the error, Eqs. (2) and (60).

In the case when the MPS is sparsely parametrized and satisfies the normalization condition, it is often possible to calculate analytically $\braket{\psi}{\dot{\psi}}$ and $\braket{\psi}{H}{\psi}$. In such cases the calculation of equations of motions is greatly simplified compared to Ref. Haegeman et al. (2011). Below we concentrate on such a case, and derive the general form of equations of motion.

A.1.3 Hamilton’s Equations of Motion

In the case when the MPS manifold is parametrized by complex parameters, the complex conjugate pairs, $z_{i}$ and $\bar{z}_{i}$ can be interpreted as momenta and coordinates. Since we use an explicitly real parametrization, we also separate the $NK$ parameters $\{{\bf x}_{i}\}$ into $NK/2$ coordinate variables $\{\bm{\theta}_{i}\}$ and $NK/2$ phase variables $\{\bm{\phi}_{i}\}$, where $N$ is the number of real on-site parameters that is assumed to be even. The amplitude variables satisfy the condition of vanishing overlap between the wave function and tangential vector in any of $\bm{\theta}$ direction:

[TABLE]

In contrast, the phase variables have a non-vanishing overlap,

[TABLE]

characterized by the set of real functions $f_{\kappa}(\{\bm{\theta}_{i}\})\in\mathbb{R}$ that depend solely on amplitudes. It is convenient to define the matrix,

[TABLE]

which can be understood as the imaginary part of the manifold metric, $g_{ij}=\braket{\partial_{x_{i}}\psi}{\partial_{x_{j}}\psi}$.

Extremizing the action, Eq. (1), over phase variables $\{\bm{\phi}_{i}(t)\}$ gives equations of motion for $\{\bm{\theta}\}$,

[TABLE]

while extremizing the action over amplitude variables yields equations of motion for $\{\bm{\phi}\}$,

[TABLE]

The equations (23)-(24) are Hamilton’s equations in a manifold with the Poisson bracket,

[TABLE]

thus making apparent their symplectic structure.

The equations of motion obtained above are determined by the expectation value of the Hamiltonian and the set of functions $f_{\kappa}(\{\bm{\theta}_{i}\})$. Below we discuss the calculation of these ingredients in the specific case of PXP and TFIM.

A.2 TDVP for PXP model

The MPS ansatz used for PXP model, Eq. (5), has a left-canonical form since the matrix $A$ satisfies the condition $\sum_{\sigma}A_{i}^{\sigma{\dagger}}A_{i}^{\sigma}=\mathds{1}$ for any values of $(\theta_{i},\phi_{i})$. In such a case, all single-site transfer matrices $T_{i}$ have largest eigenvalue $\lambda_{\text{max}}=1$, with dominant left eigenvector $(L^{i}|=(1,0,0,1)$. On the other hand, the right dominant eigenvector depends on the size of the unit cell. It should be noted that the left-canonical form of the MPS ansatz is not crucial for the calculation, but it greatly simplifies the analytical expressions.

A.2.1 EOMS for PXP with $K=2$ ansatz

For the two-site MPS ansatz we get the following transfer matrix,

[TABLE]

and the corresponding right dominant eigenvector,

[TABLE]

Using the expression for the left dominant eigenvector, we calculate their overlap,

[TABLE]

that enters as a normalization factor in all expressions.

The functions $f_{1,2}(\{\bm{\theta}\})$ are calculated as follows. We explicitly use translational invariance of the MPS state to obtain

[TABLE]

where the tensor $\partial_{\phi_{1}}A_{1}\equiv\partial_{\phi_{1}}A^{\sigma}(\theta_{1},\phi_{1})$ can be obtained from Eq. (5). Replacing the environments by the dominant left and right vectors, we get the following expression

[TABLE]

where we defined the matrix $T^{\partial_{\phi_{1}}}$ as

[TABLE]

Using the explicit form of the matrix $T^{\partial_{\phi_{1}}}$ we obtain the following expression for $f_{1}$:

[TABLE]

Repeating the calculation for $f_{2}$ gives us:

[TABLE]

The expectation value of the PXP Hamiltonian has two different contributions, depending on which sites the Hamiltonian acts. First we calculate the contraction of the local Hamiltonian term with the immediate environment. This results in matrices operating in the double virtual space,

[TABLE]

The matrix $\mathcal{H}^{2,1,2}_{PXP}$ where the operator $\sigma^{x}$ operates on the $A_{2}$ site can be obtained from Eq. (34) by replacing $\theta_{1}\leftrightarrow\theta_{2}$ and $\phi_{1}\leftrightarrow\phi_{2}$. The case of the chemical potential can be treated in a similar way, but the resulting matrices $\mathcal{H}^{1,2}_{\mu_{z}}$ are independent of variational parameters,

[TABLE]

Now by contracting these expressions with the left and right dominant vectors and taking into account the normalization, we obtain the expectation value of the respective terms in the Hamiltonian,

[TABLE]

[TABLE]

Summing these expressions we obtain Eq. (10) in the main text. The resulting equations of motion are obtained by substituting the specific form of the functions $f_{1,2}$ and the Hamiltonian density into Eqs. (23)-(24). We note that the system size $L$ enters both the functions $f_{i}$ and the expectation value of the Hamiltonian. Hence, it gets canceled in the equations of motion. The equations are the same as Eq. (9) for $C=2$.

A.2.2 Flow invariant subspace

As we noted in the main text, Eqs. (9) possess an additional symmetry when $\mu_{z}=0$. Namely, $\langle{\psi|H_{PXP}|\psi}\rangle$ in Eq. (36) vanishes for any values of $\theta_{1,2}$ when $\phi_{1,2}=0$. This leads to the emergence of a flow-invariant subspace: if we start the dynamics from any point with $\phi_{1,2}=0$, the phase variables will remain zero throughout the flow. In other words, the hyperplane $\phi_{1,2}=0$ is preserved under the flow generated by equations of motion (9) with $\mu_{z}=0$ since the time derivative $\dot{\phi}_{1,2}=0$ when $\phi_{1,2}=0$, for any values of $\theta_{1,2}$.

The emergence of the flow invariant subspace $\phi_{1,2}=0$ is not restricted to $K=2$ unit cell. In facts, it is a generic feature of the MPS ansatz (5) with any unit cell, provided the Hamiltonian is invariant under particle-hole and time-reversal symmetries (by time-reversal symmetry we mean the invariance of $H$ under complex conjugation, or equivalently the absence of any terms with an odd number of $\sigma^{y}$ matrices in the Hamiltonian). Indeed, particle-hole symmetry requires that the operator ${\cal C}=\prod\sigma^{z}_{i}$ anticommutes with $H$, or in other words ${\cal C}H{\cal C}=-H$. In addition, if we use the fact that the Hamiltonian does not change under complex conjugation, we may write $\langle{\psi|H|\psi}\rangle=-\langle{\psi|{\cal C}H{\cal C}|\psi}\rangle^{*}=-(\langle{\psi|{\cal C})^{*}H({\cal C}|\psi}\rangle)^{*}$. Now we notice that $({\cal C}|\psi\rangle)^{*}$ is equivalent to the MPS state $|\psi\rangle$ with $\phi_{i}\to-\phi_{i}$. Hence, when $\phi_{i}=0$ we obtain $\langle{\psi|H|\psi}\rangle=-\langle{\psi|H|\psi}\rangle=0$ for any values of $\theta_{i}$. This is sufficient to give a flow invariant subspace, $\dot{\phi}_{i}=0$ when $\phi_{i}=0$ according to Eq. (24).

The above conditions of particle-hole symmetry and time reversal are met not only by the pure Hamiltonian of the PXP model, $H_{PXP}$, but also by any Hamiltonian which contains an odd number of $\sigma^{x}$ matrices. This motivates the choice of the deformation $H_{\mu_{3}}$ considered below and also in the main text for $K=3$ unit cell. In passing, we also note that all terms in the quasi-local deformation conjectured by Ref. Choi et al. (2018) to give perfect quantum scars satisfy both of these symmetries, thus leaving the flow-invariant subspace intact.

Finally, we note that despite the angles $\phi_{i}$ being redundant when one restricts to the flow-invariant subspace, as we do in the next section, their presence is crucial for deriving the equations of motion (afterwards they can be set to zero). Indeed, if one sets the phases $\phi_{i}$ to zero before taking the derivatives in Eq. (23), one cannot obtain the equations of motion for amplitude variables. One possible route to the equations of motion, if one wants to set the phase variables to zero from the very start, would be to directly minimize the quantum leakage (see Eq. (60) below) with respect to $\theta_{1,2}$.

A.2.3 EOMS for PXP with $K=3$ ansatz

In this section we derive the equations of motion for the three-site MPS ansatz with three variational parameters $(\theta_{1},\theta_{2},\theta_{3})$. However, as we discussed above, such a derivation is more conveniently done with the help of conjugate variables, $(\phi_{1},\phi_{2},\phi_{3})$. Hence, in what follows we keep the phase variables and set them to zero at the end of the calculation.

All calculations are done analogously to the previous case of $K=2$ site unit cell. Hence, we omit the details of the calculations, listing the resulting expressions. The right dominant vector, is calculated to be

[TABLE]

Using the form of the right dominant vector, we calculate the normalization factor,

[TABLE]

The functions $f_{i}$ read:

[TABLE]

[TABLE]

[TABLE]

The local terms from the PXP Hamiltonian and deformation $H_{\mu_{3}}$, Eq. (13), contracted with the environment, result in the following matrices:

[TABLE]

where the indices $i$ are understood to be $\mathop{\rm mod}3$ and we use the following short-hand notations:

[TABLE]

Substituting these matrices in an expression analogous to Eq. (36), one can obtain the expectation value of the Hamiltonian. It has a cumbersome form and we do not need it here, since it vanishes when $\phi_{i}=0$. Plugging $\langle{\psi|H_{PXP}+H_{\mu_{3}}|\psi}\rangle$ and expressions for $f_{i}$ into Eq. (23), we obtain the equations of motion for variables $\theta_{i}$. Since we are interested in the flow-invariant subspace, we set $\phi_{i}=0$ resulting in the following equations of motion:

[TABLE]

where the factors that multiply the derivatives read:

[TABLE]

A.3 TDVP for TFIM

We begin by motivating the form of the TDVP ansatz used in the Ising model, Eq. (16). Its form can be obtained by considering the approximate “trotterized” unitary evolution,

[TABLE]

where the approximation works for small values of $\delta t$. The operator on the right-hand side can be written as a translationally-invariant matrix product operator (MPO) Pirvu et al. (2010) with bond-dimension $\chi=2$. This MPO has two physical indices, $\sigma_{1},\sigma_{2}=\uparrow,\downarrow$ and can be written as:

[TABLE]

Values of the angles are found to be $(\tilde{\theta},\tilde{\phi},\tilde{\chi})=(h_{x}\delta t,h_{z}\delta t,J\delta t)$.

We can obtain an MPS from this MPO by applying it to a reference state. For convenience, we choose the reference state to be $\ket{\uparrow\uparrow\ldots}$ state. This results in the MPS with $\chi=2$ that is specified by matrices:

[TABLE]

Now we allow the angles to be independent parameters, obtaining an MPS ansatz with three variational parameters. Two of these parameters are phase variables $(\tilde{\phi},\tilde{\chi})$ and one is an amplitude variable $\tilde{\theta}$. Hence, it is natural to extend this ansatz by adding yet another amplitude variable $\xi$ that is conjugate to $\tilde{\chi}$. This can be done by replacing $e^{-i\tilde{\chi}}\rightarrow\cos{\tilde{\xi}}e^{-i\tilde{\chi}}$ and $e^{i\tilde{\chi}}\rightarrow\sin{\tilde{\xi}}e^{i\tilde{\chi}}$. After such an extension, supplemented by a multiplication with the phase $e^{i\tilde{\phi}}$ and a shift in parameters

[TABLE]

Eq. (48) turns into the MPS ansatz (16) used in the main text.

The MPS ansatz (16) was obtained as the twist of a product state by a $\chi=2$ MPO operator that approximates the unitary evolution with the Hamiltonian of transverse-field Ising model. Thus, this state can be understood as the generalization of the product-state mean-field ansatz that incorporates short-range entanglement generated by unitary evolution with $H_{\text{TFIM}}$. Indeed, one can explicitly check that by restricting the angles to values $(\chi,\xi)=(0,\pi/4)$ we obtain a generic translationally invariant product state parametrized by $(\theta,\phi)$.

After justifying the form of the ansatz, we discuss its equations of motion. The unit cell consists of $K=1$ site, thus the transfer matrix can be obtained straightforwardly. Its largest eigenvalue reads

[TABLE]

This eigenvalue is doubly degenerate when the expression under the square root vanishes. In what follows we ignore the degenerate surface and assume that the density matrix is non-degenerate. Crucially, the largest eigenvalue is generally different from one, $\lambda_{\text{max}}\neq 1$. Hence, this ansatz does not result in the normalized wave function and we have to modify the derivation of our equations of motion.

Even though TDVP can be reformulated to apply to non-normalized states Kramer and Saraceno (1981), we prefer to normalize each tensor in the numerical implementations $\tilde{A}\rightarrow A/\sqrt{\lambda_{\text{max}}}$, so that $\tilde{\lambda}_{\text{max}}=1$. The rest of the calculation is similar to the previous section. The main difference is that due to the increased complexity of the tensors and the presence of normalization in tensors, analytical calculations are more complicated even if the method is exactly the same as in the PXP model. The simplest way to avoid the complexity is to evaluate all expressions including the derivatives numerically, however this may result in instabilities in the integration of the EOMS. In our implementation we evaluate all derivatives analytically but substitute numeric values for the variables for which no more derivatives are to be taken. We find that using such approach is computationally efficient, while providing numerically stable time integration. Resulting dynamical system can be found in Ref. SOM .

Appendix B TDVP for tree tensor states

In this section, we introduce the TTS ansatz for trees of arbitrary connectivity. Afterwards we discuss how to reduce the calculation of local correlation functions of a TTS to a calculation for an effective one dimensional tensor network. This allows to use the method described in Section (A.2) for calculating equations of motion. In Fig. (8) we give a graphical illustration of the algebraic calculations presented in this section. Finally, we discuss the relation between the TTS ansatz and the projected entangled pair states (PEPS) ansatz Verstraete and Cirac (2004).

B.1 Tensor tree state ansatz

The basic building block of the TTS ansatz is the tensor $A^{\sigma^{z}}_{i_{1,\ldots,i_{C}}}$ with $C$ virtual indices $i_{1,\ldots C}$ and one physical index, $\sigma^{z}$ see Fig. 8(a). Out of $C=3$ virtual indices shown in Fig. 8(a) one is designated as “output leg” with an outgoing arrow; the other indices are input legs. In order to obtain the normalized wave function that obeys the Rydberg blockade constraint, one may choose the following form of tensors:

[TABLE]

with all other elements of $A$ being zero. These expressions generalize Eq. (8) for arbitrary connectivity, $C$. Such parametrization produces an “output” 1-bit, $i_{C}=1$ when the corresponding site is excited, and 0-bit, $i_{C}=0$ when the given site is in the ground state, $\sigma^{z}=\downarrow$. In turn, the quantity $I=\sum^{C-1}_{j=1}i_{j}$ in Eq. (51) counts the number of “excitations” on $C-1$ adjacent sites. Only in the case when all of these $C-1$ adjacent sites do not have an excitation, $i_{1}=\ldots=i_{C-1}=0$, can the given site be in the excited state, thus implementing the constraint of no adjacent $\uparrow$ states.

We selected out one of the virtual indices by designating it as “output”. However, the particular choice of this leg is not important, and the TTS ansatz of Eq. (51) is invariant under permutations of all remaining input legs,

[TABLE]

where $\mathcal{P}\in S_{C-1}$ and $S_{C-1}$ denotes the symmetric group of $C-1$ symbols. This property allows us to swap input legs at will, and it also holds for TTS with two alternating tensors that we use to imitate states with $K=2$ sites unit cell on a lattice.

As we discussed, the TTS ansatz with tensors in Eq. (51) respects the Rydberg blockade constraints. This constraint imposes restrictions on the non-zero elements of the tensor $A$. The particular form of these elements in Eq. (51) is dictated by the normalization condition. For the bipartite tree, one can check that this ansatz generates a wave function

[TABLE]

where the sum runs over all product states that obey the constraint, and $N$ ensures the normalization in the thermodynamic limit and is calculated in Eq. (59). Numbers $n_{1,2}$ count the number of $\uparrow$ states in the corresponding sublattice, which are assigned the following weights:

[TABLE]

Note that the weights are the same for all connectivities. This property is reminiscent of the mean field product state ansatz where the weights would be $c_{\downarrow_{i}}=\cos{\theta_{i}}$, $c_{\uparrow_{i}}=-ie^{-i\phi_{i}}\sin{\theta_{i}}$. However, in the presence of a constraint these mean-field weights would lead to a wave function with the norm vanishing as $||\ket{\psi}||\propto e^{-aL}$ in the thermodynamic limit.

In addition to being normalized, one can check that the tensor $A$ obeys a generalized canonical gauge, see Fig. 8(b),

[TABLE]

for arbitrary values of $\theta,\phi$. This condition is enough to allow for normalization in thermodynamic limit. Moreover, it greatly simplifies the calculation of expectation values as we discuss below.

B.2 Mapping the tensor tree to a one-dimensional lattice

As explained in Section A.1.3, the main ingredient of TDVP are one-point correlation functions of the form $\braket{\psi}{\partial_{x_{k}}}{\psi}$ (or similar expressions in which the local operator is the Hamiltonian density), where $k$ denotes the site of the lattice. For clarity we work with a one site unit cell $K=1$, however the results also hold for $K=2$. We define the generalized transfer matrix, see Fig. 8(c), as:

[TABLE]

where the indices inside parenthesis are vectorized. For $C=2$, $\mathbb{T}$ is the same as the single site transfer matrix defined in Eq. (18). We note that tracing any set of $C-2$ input indices is equivalent due to the permutation symmetry of the input legs, Eq. (52). Due to the canonical gauge of Eq. (55), the left eigenvector of the transfer matrix is just the identity, ${(L|_{a,b}=\delta_{a,b}}$. We also define the generalized transfer matrix that includes local operator [see Fig. 8(d)],

[TABLE]

Using the canonical gauge [Fig. 8(b)] and property of the right eigenvector [Fig. 8(e)], we write the one-point correlation functions in Fig. 8(f) as:

[TABLE]

where $(L_{l}L_{m}|\equiv(L_{l}|\otimes(L_{m}|$. The subscript indicates which input bonds are traced by the corresponding vector, while the transfer matrix $\mathbb{T}_{C}$ is applied to the output bonds ($i_{C},i_{C}^{\prime}$), see Fig. 8(d). Equation (B.2) has the remarkable feature of being able to trace the whole (infinite) lattice by the repeated application of a transfer matrix of system size independent dimensions. Such a feature is reminiscent of the one-point function calculation using MPS, see Eq. (29), and reflects the loop-free structure of TTS. The calculation of one-point functions of local operators which are supported in more than a single site is a straightforward generalization of Eq. (B.2).

The normalization factor $N$, of Eq. (53) is derived by calculating $\braket{\psi}{\psi}$, which is a trivial modification of Eq. (B.2) and corresponds to

[TABLE]

B.3 Comparing variational wave function on a lattice and TTS

In the previous section we explained the main steps in the analytical calculation of the TDVP equations of motion on a TTS. However, our goal is to approximate the dynamics on higher dimensional lattices of the same connectivity but different topology, e.g. square lattice. The true variational wavefunction for the square lattice is also the form of Eq. (53), but with a different normalization factor. Such a wavefunction is a tensor network which has the same geometry as the underlying lattice, so for the lattices of interest it belongs to the PEPS family Verstraete and Cirac (2004).

The calculation of one point correlation functions, $\braket{\psi}{O_{k}}{\psi}$ in PEPS is known to be a sharp P-Hard problem for generic PEPS states Schuch et al. (2007) and is associated with the presence of loops in the lattice. To circumvent this problem, numerical algorithms typically employ some form of truncation in the calculation of correlation functions. However, the analytical calculation of the equations of motion is practically impossible for an infinite PEPS. In the main text we mention the comparison of the the dynamics on the true variational wavefunction for the lattice (PEPS) and the TTS. We performed the comparison by parametrizing the PEPS ansatz as Eq. (53). We use a finite system with periodic boundaries (4 by 4 lattice), in order to be able to calculate the correlation functions analytically. The TDVP equations of motion are calculated for the unnormalized state (not including the factor $N$), to avoid dealing with the complicated normalization factor.

Appendix C Classifying trajectories by quantum leakage

In this Appendix we provide additional details on how to calculate the rate at which the quantum system is leaving the MPS variational manifold – which we call quantum leakage. We start with the general framework and provide detailed justification of the fidelity bound. Afterwards, we apply quantum leakage to characterize the different periodic trajectories observed in the deformed PXP model.

C.1 General framework: quantum leakage

If one initializes the quantum state in the form of Eq. (3), at $t=0$ the values of all local observables and their time derivative are captured by TDVP equations of motion. At later times the TDVP evolution begins to disagree with the exact unitary dynamics, see Fig. 9. Intuitively, this is visualized by “leakage” of the exact quantum wave function from the variational manifold. The instantaneous leakage rate is given by the disagreement between quantum evolution and TDVP dynamics,

[TABLE]

We note that this quantity goes beyond TDVP equations of motion, as it contains the square of the Hamiltonian operator, $H^{2}$, which depends on quantum commutation rules and operator algebra. While we were able to obtain the equations of motion, Eq. (2), without explicit calculation of two-body correlators, this is not possible in the present case. Nevertheless, the use of the gauge of Eq. (19) allows to cancel all disconnected correlators. This can be accounted by replacing $\langle{*}\rangle\to\langle{*}\rangle_{c}$ in Eq. (60), where the connected correlators are defined as

[TABLE]

Upon substitution of the connected correlators the error scales as $\Lambda^{2}\propto L$. As we show in the next section, the disconnected component that is proportional to $L^{2}$ vanishes. In the thermodynamic limit the quantity of interest is the error density, $\gamma^{2}=\Lambda^{2}/L$, as shown below.

In contrast to EOMS, the error calculation is more complicated because it involves two-point correlators. We briefly describe the method presented in Haegeman et al. (2011) on how to resum two-point correlators in the thermodynamic limit and present an analytical recipe used to calculate the leakage in this work. We note that a similar calculation was performed in Ref. Ho et al. (2019) for the flow invariant subspace (${\phi}_{i}=0$, $\braket{H}=0$) of the 2-site PXP model.

We are interested in calculating the correlators defined in Eq. (61). We show how to perform the resummation of two-point correlators of a uniform MPS, initially displayed in Ref. Haegeman et al. (2011). Extensions to larger unit cells follow straightforwardly by using the transfer matrix of the corresponding unit cell. For two parameters $a,b$ we get,

[TABLE]

where $n,m$ are lattice site labels and in the second line we replaced the summation over $n,m$ by a sum over $n$ and $q=|n-m-1|$. From such expressions it is evident that one has to resum geometric series of the form $\sum_{q}[T_{u.c.}]^{q}$, which is possible if the operator has spectral radius $\rho(T_{u.c.})<1$. Since the transfer matrix has a single largest unit eigenvalue, the dominant subspace has to be projected out and resummed separately,

[TABLE]

where we have defined the projector onto the dominant subspace, $\mathcal{P}=|R)(L|/(L|R)$, its complement, $\mathcal{Q}=\mathds{1}-\mathcal{P}$, and introduced the matrix $\cal T$ as

[TABLE]

Note that the same resummation formula holds for any size of the unit cell. Substituting Eq. (65) into Eq. (64) we get,

[TABLE]

The first line contains the resummation of disconnected correlator, and these terms are cancelled by the proper gauge choice. The terms in the second line correspond to cases when both derivatives are taken in the same unit cell and long range resummations of correlations for cases when $n>m$ and $m<n$.

The expressions for $\braket{\psi}{H^{2}}{\psi}$ and $\braket{\partial_{b}\psi}{H}{\psi}$ are also calculated using Eq. (65). For the Hamiltonian which can be written as a sum of local operators $H=\sum_{i}h_{i}$, with $h_{i}$ having support on an finite number of sites, chosen to be two for simplicity, we find:

[TABLE]

[TABLE]

where $\mathcal{H}$ is the contraction of the local Hamiltonian density with environment, and the matrices $\mathcal{H}_{\partial_{b}A}$, $\mathcal{H}^{(2)}$ are defined as

[TABLE]

These definitions are valid under the assumption that the Hamiltonian has only two-site terms, but generalization to longer-range terms in the Hamiltonian are straightforward.

The disconnected correlators in Eqs. (61)-(63) are calculated as:

[TABLE]

Collecting together these expressions with Eqs. (67)-(69) we observe that terms in the connected correlators that scale as $L^{2}$ vanish. Thus, the instantaneous error scales linearly with the system size, $\Lambda^{2}\propto L$.

C.2 Fidelity bound

At long times, a generic quantum many-body system is eventually expected to develop entanglement that cannot be captured by the MPS ansatz in Eq. (3). At this point, the TDVP dynamics $\ket{\psi(\{x_{a}(t)\})}$ and exact unitary evolution, $e^{-iHt}\ket{\psi(\{x_{a}(0)\})}$, are expected to strongly disagree. However, as we demonstrate in the main text, the TDVP equations of motion are generically expected to have short periodic trajectories, see the blue line in Fig. 9. If the quantum system were following the TDVP dynamics exactly, that would imply persistent oscillations in local observables and revivals of the many-body quantum fidelity to unity. Due to non-zero leakage, the quantum and TDVP evolution would disagree after a single period. In what follows we obtain a lower bound on the many-body fidelity, Eq. (6), after the time $T$, i.e., the period of the trajectory.

To obtain the bound we represent the exact wave function, $\ket{\Psi}$, that follows the exact evolution according to the Schrödinger equation, $i\partial_{t}\ket{\Psi(t)}=H\ket{\Psi(t)}$, as a sum of two terms,

[TABLE]

where $\ket{\psi(\{x_{a}(t)\})}$ belongs to the variational manifold and is obtained from the TDVP equations of motion, while $\ket{\varphi(t)}$ is the error vector defined as the difference between these two wave functions, see Fig. 9.

Expressing the error vector via $\ket{\Psi(t)}$ and the TDVP wave function, we obtain the following equation of motion:

[TABLE]

where we explicitly used the Schrödinger equation for $\ket{\Psi(t)}$. The first term in the last line transports the error forwards in time and does not change the norm of $\ket{\varphi(t)}$. The second term describes the change of the error due to mismatch between exact evolution $\ket{\Psi(t)}$ and its TDVP projection, $\ket{\psi(\{x_{a}(t)\})}$, and its norm coincides with the quantum leakage defined in Eq. (60). Calculating the change in the norm of vector $\ket{\varphi}$ from Eq. (76) we obtain

[TABLE]

where we omit time dependence of $\varphi$ to simplify notations. Using the triangle and Cauchy-Schwartz inequalities we bound the growth of the norm of the error vector as

[TABLE]

where $\|\varphi\|=\sqrt{\langle\varphi|\varphi\rangle}$ is the norm of the error vector. Using this notation, we obtain an upper bound on the rate of increase of the norm of $\ket{\varphi}$,

[TABLE]

where the quantity $\Lambda(t)$ is the instantaneous geometric error that was already defined in Eq. (60). Now, the triangle inequality can be used to bound the norm of the error vector accumulated during evolution over the finite period of time $t$,

[TABLE]

by the integrated geometric error $I_{t}$.

The Fubini-Study metric or quantum angle can be defined as

[TABLE]

It is a natural metric on projective space, where pure states are represented by lines through the origin of Hilbert space, reflecting the fact that global phases are not physically meaningful. It is a quotient of the standard Euclidean metric restricted to the unit sphere. The distance $\gamma$ can be visualized as an angle between the two lines.

The three vectors $\ket{\Psi(t)}$, $\ket{\psi(\{x_{a}(t)\})}$ and $\ket{\varphi(t)}$ form sides of a triangle with the angle between the exact quantum wave function and its TDVP approximation equal to the Fubini-Study distance between those two states. Using the law of cosines,

[TABLE]

Assuming that we initialized the system on the periodic TDVP trajectory, and taking the time $t=T$ to be the period of this trajectory, we arrive at the following bound on fidelity,

[TABLE]

where $F_{\Psi}=|\langle\psi(\{x_{a}(T)\})|\Psi(t)\rangle|^{2}=|\langle\Psi|e^{-iHT}|\Psi\rangle|^{2}$ and we used the fact that the quantum system is initialized on the TDVP periodic trajectory, $\ket{\Psi(0)}=\ket{\psi(\{x_{a}(0)\})}=\ket{\psi(\{x_{a}(T)\})}$.

Until now, we did not discuss the scaling of fidelity and error with the system size. Naïvely if we take the thermodynamic limit in Eq. (83), this bound would appear useless since $I^{2}_{t}\propto L$ increases linearly with the system size and eventually becomes larger than one. However, assuming that the entanglement growth is not too fast (the weak leakage case), one may find an intermediate system size such that $I^{2}_{t}\ll 1$ and the bound is still effective. In such a case, assuming the scaling form of the fidelity $F_{\Psi}=e^{-f_{T}L}$, with $f_{T}\ll 1$, we can expand $\sqrt{F_{\Psi}}=1-f_{T}L/2$. Plugging this expansion in Eq. (83) along with the definition of $\Gamma_{T}=[\int_{0}^{T}dt\gamma(t)]^{2}=I^{2}_{t}/L$, we obtain,

[TABLE]

recovering the upper bound on $f_{T}$ that governs the fidelity decay, see Eq. (12) in the main text.

We note that the above argument cannot be regarded as a rigorous proof of the fidelity bound in the many-body case. However, using the limited velocity of entanglement growth from the initial weakly-entangled state $\ket{\psi(\{x_{a}(T)\})}$, we can limit the system size $L$ needed to be effectively in the thermodynamic limit, thus justifying the expansion used in deriving Eq. (84). It is an interesting open question if a similar bound can be rigorously proven under some assumptions limiting the entanglement growth, and if a tighter bound could be obtained.

C.3 Trajectories in the deformed PXP model with $K=3$ ansatz

C.3.1 Calculation of leakage

In what follows we use the leakage function calculated for the case of the Hamiltonian given by Eqs. (4) and (13). To calculate the leakage we use Mathematica to perform the required contractions and substitute them into Eqs. (67)-(69). While in the case of the non-deformed PXP model one can obtain simple analytical expressions for the leakage, for the deformed PXP model these expressions become too lengthy to be presented here SOM .

C.3.2 Calculation of Floquet exponent

In order to characterize the periodic orbits observed in the TDVP dynamics we use the integrated leakage that was discussed above. In addition, we study the stability of periodic orbits by calculating the Floquet exponent. The Floquet exponent characterizes the stability of a given periodic orbit, in a manner similar to the Lyapunov exponent. However, despite qualitative similarities, the Floquet exponent quantitatively differs from the Lyapunov exponent Cvitanović et al. (2016). The latter quantity was used to characterize chaotic dynamical flow in the case of large bond dimension TDVP Hallam et al. (2018). In contrast, here we are dealing with periodic orbits. In this case, the Floquet exponent allows to quantify the orbit stability in a way that is invariant under all local smooth nonlinear coordinate transformations. Thus, the Floquet exponent is an intrinsic characteristics of the periodic orbit.

Since we do not use the Lyapunov exponent in this work, we keep the notation $\lambda$ for the Floquet exponent. It is defined as

[TABLE]

where $\Lambda_{\text{max}}$ is the largest-norm Floquet multiplier of the orbit in the space of TDVP parameters, $x_{a}^{(0)}(t)$, and $t\in[0,T]$ where $T$ is the orbit period. Intuitively, $\Lambda_{\text{max}}$ characterizes how the unit volume, chosen at some initial point on the periodic orbit, transforms after one traversal of the orbit. More specifically, $\Lambda_{\text{max}}$ is the largest-norm eigenvalue of the orbit Floquet matrix, $[J_{T}]_{ab}=\partial x_{a}(T)/\partial x_{b}(0)$, that quantifies the effect of small perturbation at time $t=0$ after one period of the orbit, $T$. The orbit Floquet matrix is calculated as the time-ordered integral of the instantaneous Jacobian matrix along the periodic orbit Cvitanović et al. (2016),

[TABLE]

where we assumed that $v_{a}(x)$ is the vector of velocities that generates the flow [in other words, the equations of motion are written as $\dot{x}_{a}=v_{a}(x)$].

Numerically, the time-ordered integral in Eq. (86) is calculated by solving a system of non-autonomous first-order differential equations. Since we have three variational parameters $\theta_{1,2,3}$, $J_{T}$ is $3\times 3$ matrix. As expected, we find that the spectrum of $J_{T}$ always has an eigenvalue $\Lambda=1$ Cvitanović et al. (2016), and the remaining two eigenvalues satisfy $|\Lambda_{\text{max}}\Lambda_{\text{min}}|=1$, reflecting the symplectic (volume-preserving) nature of the flow.

If an orbit is stable and surrounded by KAM torus, we find that all eigenvalues of the Floquet matrix have $|\Lambda|=1$. Thus, the Floquet exponent vanishes, $\lambda=0$, signaling the stability of the orbit to small deformations. In contrast, for unstable periodic orbits we have $|\Lambda_{\text{max}}|>1$ and $|\Lambda_{\text{min}}|<1$. This implies that a unit volume after one period is elongated in one direction and compresses in a different direction.

C.3.3 Trajectory analysis

With the expressions for quantum leakage and Floquet exponent at hand, we turn to the analysis of different trajectories in the deformed PXP model. First, we consider the deformation with a large negative value of $\mu_{3}=-0.25$. Figure 10 shows that this deformation leads to a drastic increase of the size of regular islands compared to Fig. 2(a), which shows the Poincaré section for $\mu_{3}=0$. Nevertheless, despite an increase in size of the regular regions, the fidelity revivals are degraded by the deformation with large negative values of $\mu_{3}$, see Fig. 6. The leakage follows the same trend, as is apparent from the same figure. However, the orbit remains stable, hence its Floquet exponent $\lambda=0$ throughout the entire range of negative values of $\mu_{3}$.

Next, we turn to the case of $\mu_{3}=0.25$. The Poincaré section in Fig. 11(a) shows that the deformation with $\mu_{3}>0$ decreases the size of regular regions and makes the TDVP dynamics more chaotic. Different symbols in Fig. 11(a) show the location of several periodic orbits with short periods that we were able to find using Newton’s search method Strogatz (2014); Cvitanović et al. (2016) (we show only a subset of short orbits and omit those which are symmetry-related).

Figure 11(b) visualizes the TDVP dynamics of local observables in these different orbits and also shows the integrated leakage, $\Gamma_{P}$, and the Floquet exponent. In particular, $\star$ symbol shows the location of an orbit related by symmetry to the one studied in the main text. We observe that this orbit is not surrounded by a torus any more, consistent with finite value of $\lambda$. However, there are still many stable orbits present in the Poincaré section (e.g., the ones denoted by $\bullet$ and ${\small\blacksquare}$ symbols). These orbits have much stronger leakage, despite the vanishing Floquet exponent $\lambda$. This confirms our conclusion that the Floquet exponent is not obviously associated with quantum leakage.

On the other hand, the $\blacklozenge$ orbit in Fig. 11(b) has a slightly smaller leakage compared to the $\star$ orbit. By comparing quantum dynamics we indeed find that these two orbits have similar behavior of the fidelity revivals. However, due to relatively large leakage, the fidelity revivals and oscillations of local observables are damped. Hence, while the orbit $\blacklozenge$ in principle represents another type of oscillations, it remains an open question if the quality of these oscillations can be improved, e.g. if they can be stabilized by a deformation of the Hamiltonian.

Appendix D Thermalization in the strong leakage regime

In this Appendix we present additional results for the cases when the system has no low-leakage trajectories. We start with additional data for the TFIM. In addition, we consider a strongly deformed PXP model with broken particle-hole symmetry and demonstrate that it still has entanglement dynamics that depends on the initial conditions.

D.1 Entanglement dynamics and ETH indicators in TFIM

We begin with illustrating the periodic orbit found in the TDVP dynamics of a TFIM using Newton’s search algorithm. Figure 12 illustrates the Poincaré section defined by $(\xi=0.9,\dot{\xi}<0)$. The periodic orbit with period $T=2.09$ crosses this plane at the point $(\chi,\xi,\phi,\theta)=(0.2607,0.9,4.888,0.4308)$, shown by “$\star$” symbol in Fig. 12. This is a stable orbit surrounded by a small KAM torus, but nevertheless it does not give rise to many-body fidelity revivals. This may be attributed to large leakage, $\Gamma_{T}=0.57$.

In the main text we illustrated the strong influence of initial conditions on entanglement growth. Here, we compare the initial states with slowest and fastest entanglement dynamics. These states are defined by the MPS in Eq. (16), with parameters $(\chi,\xi,\phi,\theta)_{\text{slow}}=(0.1,0.9,1.17,3.86)$ and $(\chi,\xi,\phi,\theta)_{\text{fast}}=(-0.32,0.9,1.6,4.28)$. These states can also be identified by local expectation values, $(\sigma^{x},\sigma^{y},\sigma^{z})_{\text{slow}}=(0.38,0.85,0.072)$ and $(\sigma^{x},\sigma^{y},\sigma^{z})_{\text{fast}}=(0.27,0.74,-0.49)$. Figure 13 shows the dynamics of entanglement of small subsystems for slow and fast initial states. Both states reach identical values of the saturated entanglement entropy for all subsystems. This is expected, since both initial states have the same energy density per site. As energy is the only conserved quantity, these states are expected to have the same entanglement saturation value. However, the time when the entanglement saturates is different, proving that these two initial states have different entanglement velocity.

Finally, Fig. 14 shows a test of ETH. For the TFIM for the considered values of parameters, i.e., $(J_{z},h_{z},h_{x})=(1,0.4,1)$. Fig. 14(a) illustrates the expectation values of local observables and confirms that average differences of local magnetization between adjacent eigenstates at energy density $E/L=0.19$ decay with systems size according to the ETH prediction. In Fig. 14(b) we observe that the entanglement of eigenstates essentially has no outliers near the middle of the spectrum.

D.2 Entanglement dynamics and ETH indicators in the deformed PXP model

We consider the following deformation of the PXP model, $\delta H=H_{\mu_{z}}+H_{\mu_{2}}$ where

[TABLE]

and values of the couplings are fixed to $\mu_{2}=0.24$ and $\mu_{z}=0.4$. This deformation is chosen in a such way that both terms in $\delta H$ break the particle-hole symmetry of the model. In addition, the correlated spin flip term is expected to enhance thermalization.

Indeed, the entanglement entropy of many-body eigenstates of the Hamiltonian $H_{PXP}+\delta H$ shown in Fig. 15 shows that there are no “outliers”, i.e., states with anomalously low entanglement. Likewise, the level statistics indicator $r={\min(\delta_{i},\delta_{i+1})}/{\max(\delta_{i},\delta_{i+1})}$ yields $\langle{r}\rangle=0.528$, a value that is close to the GOE prediction $0.5307$ Atas et al. (2013) (averaging is done over the middle 1/3 of the full many-body spectrum).

The TDVP equations of motion for $K=2$ MPS ansatz in presence of this deformation can be obtained using the framework laid out in Appendix A. Using the expectation value of $H_{\mu_{2}}$ that is given by

[TABLE]

we obtain additional terms that have to be added to the equations of motion (9)

[TABLE]

The EOMS for $\delta\dot{\theta}_{2}$, $\delta\dot{\phi}_{2}$ are obtained by swapping the subscripts $1\leftrightarrow 2$.

Similarly to the case of TFIM discussed in the main text, we study the dependence of entanglement dynamics on the initial condition within the MPS manifold. Initially, we find a periodic trajectory at zero energy density. The trajectory is found by starting from the $\mathbb{Z}_{2}$ trajectory of the PXP model in a similar fashion to the trajectory with the chemical potential discussed in the main text. In this case, we start from the periodic trajectory at $\mu_{z}=0.4$, and slowly ramp up $\mu_{2}$ at fixed energy density, $\braket{H}/L=0$. In practice, we can choose the initial state to be any phase space point of the periodic trajectory. We have observed that even though the entropy of states in this periodic trajectory fluctuates, it does not considerably affect the long time entropy growth in the exact quench dynamics. Therefore, we fix the initial point to be $(\theta^{\prime}_{1},\phi^{\prime}_{1},\theta^{\prime}_{2},\phi^{\prime}_{2})=(3.926,0.289,2.987,0.121)$. This state has bipartite entanglement entropy $S\approx 0.036$. In the chosen MPS ansatz the bipartite entanglement entropy $S(\theta_{1},\theta_{2})$ is independent of the phases. Thus, we can easily pick different initial states with exactly the same entropy by varying $\phi_{2}$ while $\phi_{1}$ is fixed by the energy density constraint. The initial states of Fig. 16 are generated by splitting the domain $[\phi^{\prime}_{2}-\pi/4,\phi^{\prime}_{2}+\pi/4]$ to $200$ points, and then discarding the points for which the energy density constraint cannot be satisfied for any value of $\phi_{1}$. The number of remaining states are $15$, including the periodic trajectory.

In Fig. 16 we observe that the entropy growth of a state is closely related to the leakage of the state out of the MPS manifold. Compared to the case of TFIM studied in the main text, the PXP model features stronger correlations between leakage and entanglement with correlation coefficient being 0.95 (compared to 0.53 for TFIM). In this case there are two sets of states which have different velocities of entropy growth, and all the states in the “slow” set have considerably smaller leakage than any of the states in the “fast” set. The periodic trajectory belongs to the set of slow states but it is not the slowest one.