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

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.


I. INTRODUCTION
Technological advances in synthetic quantum systems [1,2] have started an era where nonequilibrium 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) [3] and Bethe-ansatz integrable systems [4], may host nontrivial 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 toward extending quantum coherence in many-body systems.
Progress toward extending quantum coherence is tied with developing a more complete understanding of quantum many-body chaos and thermalization [5]. In this direction, recent studies [6][7][8] consider a model of "typical" quantum dynamics generated by applying random unitary operators in a quantum circuit, which allow one to obtain results for the dynamics of entanglement and other physical observables. In such models, the thermalization rate should not depend strongly on the initial state. Recent experiments, however, reveal that the dynamics in physical many-body systems may be significantly more complex [9,10]. In particular, for certain initial states, the chain of interacting Rydberg atoms [10] is found to exhibit much slower thermalization or even its absence on the experimentally accessible timescales. The 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 are explained by "quantum many-body scars" [11][12][13]-a subset of anomalous eigenstates with strongly nonthermal properties [11]. 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 [14]. Atypical eigenstates had previously been constructed analytically in the Afflect-Kennedy-Lieb-Tasaki spin chain [15,16], and it has been suggested that some of these nonthermalizing states may be closely related to the ones in the Rydberg atom chain [17,18]. Additionally, various types of nonthermalizing behaviors have been reported in a number of other models [17,[19][20][21][22][23][24][25][26].
In the context of few-body quantum chaos [27], the strong dependence of the 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 [28,29]. According to the Kolmogorov-Arnold-Moser (KAM) theorem [30], 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 [31]. The semiclassical limit allows one to introduce a notion of regular and chaotic eigenstates that are dominated by the corresponding regions of phase space [28,29]. 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 [27].
In contrast to few-body quantum systems, in the manybody case one expects that conditions for the 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) [32], 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 (MPSs) [33] with a fixed number of degrees of freedom (d.o.f.) [10,13], as well as higher-dimensional systems where we use the infinite tensor tree states (TTSs) [34]. 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 three-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 infinite-time evolving block decimation (iTEBD) [35].
Similar to the known examples of few-body systems [29], we expect that regular regions of mixed phase space give rise to nonthermal eigenstates in the many-body quantum system. Thus, mixed phase space provides a more general mechanism for slow thermalization compared to quantum scars. We establish TDVP as a practical method for finding models with nonergodic quantum dynamics in arbitrary spatial dimensions, which complements its recent applications to quantum thermalization [36,37]. Our approach provides a potential direction to generalize 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 approach is distinct from other approaches that rely on broken Bethe-ansatz integrability [38] or the absence of quantum resonances in the MBL phase [39][40][41].
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 [10,42,43]. We demonstrate that TDVP allows one to identify a stable trajectory that gives rise to new quantum revivals beyond those that were probed in recent experiments [10]. 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.

A. A brief overview of TDVP
TDVP equations of motion are obtained by extremizing the action [44], Such action yields a set of nonlinear classical equations of motion [44], where the variational parameters are real valued, x i ∈ 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 the Appendix A for a detailed derivation. TDVP has been successfully applied to the cases where the variational manifold is spanned by MPS [33,45] as well as finite TTSs [46]. Below, we introduce the basics of MPS parametrization, which is 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 is used in Sec. III to capture the dynamics of the Rydberg model generalized to higher dimensions.

B. 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)]: where A σ ðxÞ is χ × χ matrix (χ is the bond dimension) that depends on N real on-site variational parameters, x ¼ ðx 1 ; …; x N Þ. The physical indices σ label the basis of on-site Hilbert space, which we take to be a spin-1=2 d.o.f. with σ ¼ ↑; ↓. The index m ¼ 1; …; L=K labels the unit cells, and in what follows we take the thermodynamic limit L → ∞. 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; …; NK, to fully fix the state. Below, we 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 a larger number of parameters, a similar analysis could be performed numerically.
As another extension, a close relative of the MPS Ansatz are the TTSs with higher connectivity-see Fig. 1(b). On the one hand, TTSs allow one 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 TTSs which approximate bipartite lattices with a natural choice of a two-site unit cell. The manifold of TTSs with two inequivalent sites, for a system of size L and connectivity C, is defined by generalizing Eq. (3) to the case where A σ ðxÞ is a rank-C tensor, with bond dimension χ. The product of matrices in Eq. (3) is extended to contraction of indices on all bonds connecting tensors according to Fig. 1(b). While TTSs with a fixed bond dimension cannot capture the loop effects of higherdimensional lattices, they can still describe the local physics of weakly entangled states. In Sec. III, we find the TTS Ansatz to be very accurate in capturing the  1. (a) The state of a quantum system that is short-range entangled and translationally invariant with a unit cell of size K can be described by an MPS with the same size of unit cell. (b) The wave function of the square lattice is approximated by a TTS in the form of a Cayley tree with connectivity C ¼ 4. The state with K ¼ 2 unit cell is described by a tree where two different tensors alternate.
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.

A. PXP model
In this section, we focus on the PXP model [10], describing a 1D chain of Rydberg atoms defined by the Hamiltonian where σ x is the standard Pauli matrix and operator P j ¼ ð1 − σ z j Þ=2 projects to the ↓ 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 → ∞, and restrict the Hilbert space to spin configurations without two adjacent up spins, j…↑↑…i. This space is the largest connected component of the full Hilbert space for the Hamiltonian in Eq. (4). The PXP model is interacting and nonintegrable [11], yet its relaxation strongly depends on the initial conditions [10]. For instance, there is fast thermalization when the system is prepared in a down-polarized state, j↓↓…i, while other initial states, such as the jZ 2 i ¼ j↑↓↑↓…i state, exhibit revivals of local observables [10], entanglement entropy [11], and even the many-body wave function [12].
Below, we apply the Ansatz in Eq. (3) with K ¼ 3 to the PXP model. We parametrize the matrices A σ by two angles, The matrix A ↑ ∝ σ þ satisfies the condition A ↑ ðθ; ϕÞA ↑ ðθ 0 ; ϕ 0 Þ ¼ 0, thus effectively imposing the constraint that no two adjacent spins are in j↑i state.
The Hamiltonian in Eq. (4) has particle-hole symmetry in the spectrum: E → −E, since it anticommutes with C ¼ Q i σ z i . In this case, the system's equations of motion have a class of solutions with phase angles being stationary, ϕ 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 particlehole symmetry and time-reversal invariance of the PXP Hamiltonian in Eq. (4); see Appendix A. Note that hHi ¼ 0 vanishes identically when ϕ i ¼ 0, thus imposing no additional constraint. Restricting to ϕ i ¼ 0 subspace for a K ¼ 2 unit cell, Ref. [13] finds 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 jZ 2 i ¼ j↓↑↓↑…i.
However, a two-dimensional phase space is nongeneric 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 three or more parameters, which have regular and chaotic regions coexisting in phase space. We show below that such behavior, known as "mixed phase space," is robust and persists for various deformations of the PXP model. We then investigate the implications for quantum dynamics, finding new initial states that exhibit revivals and slow thermalization.

B. Revivals in a three-site unit cell
Higher-dimensional variational phase space naturally arises when the unit-cell size is increased, while particlehole symmetry is kept intact. This space 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 jZ 3 i ¼ j↑↓↓…i, which can be experimentally prepared in a Rydberg system with a larger blockade radius [10]. The resulting equations of motion for θ i , i ¼ 1, 2, 3, obtained from the TDVP projection of the PXP model can be found in Appendix A 2 c.
In order to analyze the dynamical flow, we consider its Poincaré section. The flow generates a discrete mapping known as the Poincaré map [47] that maps a given point ðθ 1 ; θ Ã 2 ; θ 3 Þ on the chosen hyperplane θ 2 ¼ θ Ã 2 to a position ðθ 0 1 ; θ Ã 2 ; θ 0 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). Figure 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 shortestperiod orbit denoted by a star in Fig. 2(a) and located at ðθ 1 ; θ 2 ; θ 3 Þ ¼ ð0.8356π; 0; 0.1923πÞ. 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). Figure 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 result, note that, when the entanglement for one cut, say, S 12 , vanishes, the entanglement for a different cut is nonzero, 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 jZ 3 i ¼ j↓↑↓…i. 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, we show the negative logarithm of fidelity 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 a higher return amplitude) when the system is initialized on the TDVP trajectory compared to the jZ 3 i product state. The bottom in 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.

C. Effect of regular regions of phase space
Above, we consider 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 jψðcÞi, where the parameter c ∈ ½0; 1 interpolates along the dashed line in Fig. 2(a) between the periodic orbit and the point θ 1 ¼ 0; θ 3 ¼ π.
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 ≈ 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 nonergodic behavior, 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 determination smoothens out any abrupt transition between the nonergodic 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 (IPRs) of local density matrices, measured at the time when the dynamics . 2. (a) Poincaré sections at θ 2 ¼ 0 for a three-site unit cell reveal mixed phase space with large regular regions. A stable periodic orbit that gives rise to robust revivals is denoted by star, while diamond denotes a longer orbit that does not give rise to revivals in exact dynamics. The dashed line represents a particular cut through the regular region, discussed in Sec. II C. (b) Variational TDVP prediction for the local density of the Rydberg excitation (top) and the entanglement entropy (bottom) for the shortest-period trajectory labeled by star on the Poincaré section in (a). The inset shows the unit-cell choice and entanglement cuts. (c) Quantum dynamics for jZ 3 i initial state and two entangled initial states marked in (a) obtained from iTEBD with χ ¼ 900. Top: The decay of fidelity per spin has regular minima that correspond to revivals of fidelity in a finite-size system. Bottom: Entanglement entropy. The stable "star" trajectory features a better fidelity revival, as well as strongly suppressed entanglement growth. The growth of entanglement limits the iTEBD simulation to shorter times for jZ 3 i, diamond states.
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 in 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 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 focus on a 2D square lattice and compare the predictions of the TTS dynamics against exact diagonalization. We also discuss generalizations of these results to any bipartite lattice in an arbitrary number of spatial dimensions.

A. PXP model in higher dimensions
We start by introducing a constrained model on the square lattice, inspired by the 1D PXP model: where the pair of indices i, j is used to label lattice sites. In addition to the PXP term, we also include a possible perturbation in terms of chemical potential μ z , which couples to the Rydberg excitation density operator n ¼ ð1 þ σ 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 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 σ i 1 ;…;i 4 , which forbid excitations on adjacent sites. Relegating the detailed discussion to Appendix B, here we list only the nonzero elements of the tensor A σ that is parametrized by two angles (θ; ϕ): where i 1;…;4 are virtual indices and I ¼ P 3 j¼1 i j counts the number of "excitations" on the three legs of the tensor. The fourth leg of the tensor produces an "output 1-bit" on the last virtual index i 4 if the corresponding site is excited and 0 if it is in the ground state. In 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 two-site unit cell results in the dynamical system of four variables that come in two canonically conjugate pairs, ðθ 1 ; ϕ 1 Þ and ðθ 2 ; ϕ 2 Þ. The general form of equations of motion reads where the square lattice case (ignoring loops) corresponds to C ¼ 4. The equations for _ θ 2 ; _ ϕ 2 can be obtained by substitution 1 ↔ 2. In addition, the energy density The entropy of a region of size L at fixed time t 0 ¼ 15 remains low when the quantum dynamics is launched near the center of the regular region (we use the iTEBD algorithm to simulate quantum dynamics). When the border of the regular region is approached, the value of entanglement rapidly increases with changing the initial state, defined by the value of c. When c approaches one, the dynamics becomes insensitive to the changes in the initial state due to rapid thermalization. (c) The derivative of entropy along the trajectory is most sensitive around the boundary of the regular region.
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.

B. Generalizing scars to arbitrary bipartite lattices
When μ z ¼ 0 and one restricts to a flow-invariant subspace ϕ 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. [13] 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 ↑ 1 and all atoms on the second sublattice are in the ground state ↓ 2 . In Fig. 4(a), we show that quantum dynamics of a local observable hσ z i on the 4 × 4 square lattice agrees well with TDVP dynamics based on C ¼ 4 TTS. Moreover, we compare the TDVP dynamics of a wave function satisfying the kinetic constraints and having the same weights as the tree parametrization [see Eq. (B4) in Appendix 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 the 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 → ∞. In this limit, the sites that are in their ground state cannot flip unless all their neighbors are also in their ground state. The periodic trajectory can be understood by studying the C → ∞ limit of Eq. (9). At the first quarter of the period, the excited sites on the first sublattice relax ↑ 1 → ↓ 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 ↓ 2 → ↑ 2 . The half period over which two sublattices exchange their state is given by π in our notations, so the total period of the orbit is T C→∞ ¼ 2π. Figure 4(b) shows that period of the orbit approaches its C → ∞ limit approximately as

C. Phase space and revivals in higher dimensions
When detuning μ z ≠ 0, the dynamics generated by Eq. (9) takes place on constant energy surfaces in fourdimensional space. In order to visualize such dynamics, we fix θ 1 variable and show the resulting Poincaré sections in Fig. 5(a) [the precise value of θ 1 is not important, and we choose ðθ Ã 1 ¼ 1.25π; _ θ 1 < 0Þ]. Energy conservation results in complicated surfaces in the space ðϕ 1 ; θ 2 ; ϕ 2 Þ that cannot be globally projected. Therefore, we show a small part of the Poincaré section that can be projected onto the ðϕ 1 ; θ 2 Þ plane. We find that the point ðϕ 1 ; θ 2 Þ ≈ ð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. Figure 5(b) shows this comparison for local observables hσ x;z i i, for the system initialized on the periodic trajectory (indicated by ⋆ 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 ≲ 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 d.o.f.) and the true many-body wave function on a two-dimensional lattice (see Appendix B 2). In order to visualize the variational (classical) dynamics, Fig. 5(c) shows the evolution of the trajectory of an individual spin ½hσ x 1 ðtÞi; hσ y 1 ðtÞi; hσ z 1 ðtÞi on the Bloch sphere as perturbation μ z ≠ 0 is turned on. The second spin in the unit cell performs similar oscillations shifted in phase by π.
For μ z ¼ 0, the trajectory gives oscillations of the spin between ↑ and ↓ states, and the expectation value hσ x ðtÞi is zero at all times. In agreement with physical intuition, the nonzero 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 μ ¼ 0.
In the entire range of μ z ∈ ð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 μ z suggests that the most stable revivals occur for initial states that are different from the Neél product state Z 2 by having some amount of entanglement.

IV. CHARACTERIZING TRAJECTORIES BY LEAKAGE
Variational dynamics of the PXP model and its deformations, studied above, reveal mixed phase space with multiple stable orbits. The Lyapunov exponent, which is often used to quantify chaos in TDVP [37] and semiclassical [48] 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.

A. Rate of leaving the variational manifold
The TDVP approach accurately approximates short-time quantum evolution of initial states in the chosen MPS manifold, but at long times the errors 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 a low bond dimension can represent only 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 [13,33]: 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 γ assumes a finite value in the thermodynamic limit.
In the previous section, we demonstrate 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, jΨi ¼ jψ½x ð0Þ a ðtÞi. The disagreement between quantum dynamics and TDVP generally precludes such perfect revivals (except for models with perfect scars [49]). 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 → ∞ limit, and using the large system size limit of Eq. (11), we posit the following upper bound on f T : 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 nontrivial, 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.

B. 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: 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 PXP þ H μ 3 allows us to use the same Ansatz as in Sec. II B with three variational parameters θ 1;2;3 . The resulting TDVP equations of motion are again cumbersome and can be found in Appendix A 2 c. Analysis of these Poincaré sections reveals a pronounced asymmetry between the effect of this deformation for μ 3 ¼ AE0.3. The deformation with μ 3 < 0 increases the size of regular regions in the phase space; see Fig. 10. In contrast, positive μ 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 μ 3 : For μ 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 γ½fx a g to calculate the integrated leakage Γ T , defined in Eqs. (11) and (12), along the trajectory. Figure 6 shows that the leakage attains a minimum at a small negative value μ TDVP 3 ≈ −0.11. For larger negative values of μ 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 Γ T . The f T obtained from exact quantum dynamics achieves a minimum at μ ED 3 ≈ −0.06, which corresponds to the best fidelity revivals. The optimal value of the deformation, μ ED 3 , is approximately twice smaller compared to the TDVP prediction μ TDVP

3
. Nevertheless, TDVP may be used to find the best operators suitable for stabilizing quantum many-body revivals [49]. Moreover, in contrast to the revivals from the jZ 2 i state [13], our results show that generally one must account for the effect of the deformation on the initial state.
Eventually, for large absolute values of μ 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 result to fast thermalization at such deformation parameters. At these values, quantum dynamics no longer exhibits many-body revivals, and simultaneously Γ T is on the order of unity. Thus, we propose as a tentative criterion for when a periodic trajectory ceases to result in revivals. This criterion was previously heuristically conjectured in Ref. [13] in the context of dynamics of the 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 μ 3 ¼ 0.3, which shows that the majority of trajectories have leakage that is larger than one, thus not giving rise to quantum revivals.

C. Floquet exponents and quantum leakage
We demonstrate 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 d.o.f.-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 reparametrizationinvariant measure [50]; see Appendix 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: Figure 6 (inset) shows how a nonzero Floquet exponent λ emerges for μ 3 ≥ 0.255. However, the nonzero value of λ does not have any influence on the leakage, and The data illustrate that f T is system size independent and therefore is well defined for μ 3 ∈ ½−0.6; 0.4. The inset shows that the Floquet exponent λ becomes nonzero only for μ 3 > 0.255, suggesting that there is no direct relation between λ and leakage. Finally, for μ 3 ≥ 0.45, the periodic trajectory disappears. this unstable trajectory still gives rise to fidelity revivals for μ 3 ∈ ð0.255; 0.4. This example also suggests that mixed phase space provides a more generic mechanism for weak ergodicity breaking. The stable trajectory in the center of a regular region of mixed phase space that gives rise to "regular eigenstates" [28,29] may become unstable upon a deformation of the Hamiltonian, leading to the appearance of quantum scars [11,13], which are quantum eigenstates affected by the unstable trajectories.

V. NONUNIVERSAL THERMALIZATION
Above, we focus 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 [49]). 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 the order of Oð1Þ. Nevertheless, we show that mixed phase space still leaves an imprint on quantum dynamics, leading to nonuniversal thermalization.

A. 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 In order to capture the dynamics, we use the following MPS Ansatz with bond dimension χ ¼ 2: that depends on four parameters, ðθ; ϕ; ξ; χÞ, per spin. Such a choice of the Ansatz is inspired by "trotterization" of the evolution operator e −iH TFIM δt ≈ e −iH 0 δt e −iH 1 δt , where H 0 and H 1 are the transverse-field term and remaining terms in H TFIM , respectively. We restrict ourselves to translationinvariant states, setting the unit-cell size K ¼ 1. We note that 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.

B. 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 hHi=L ¼ 0.19, which corresponds to initial states at a 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: Γ 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 the order of t ∼ 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 ent ∈ ½ð2=3ÞS 0 ; ð4=3ÞS 0 , where S 0 ∼ 0.08 is the minimal value of entanglement entropy of the given periodic orbit. Since all these states have the same energy density and nearly the same initial value of the entanglement and are translationally invariant, we expect that, after a short time sufficient to equilibrate local observables, the entanglement dynamics should be independent of the initial state. Figure 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 eigenstate thermalization hypothesis (ETH); hence, such dynamics cannot be explained by the existence of anomalous eigenstates. We note that a similar phenomenon is observed in a Floquet version of TFIM with a somewhat weaker integrability-breaking field [52].
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 a 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 are well known in the context of integrable models [53][54][55]. 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 nonintegrable 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 [56] remains an intriguing open question.

VI. DISCUSSION
In summary, we examine the relation between exact quantum dynamics and an effective classical nonlinear 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 χ ¼ 1 MPS) and semiclassical treatments used in few-body quantum chaos in that it incorporates short-range entanglement. We demonstrate the relevance of mixed phase space, identified in the TDVP dynamics, for the exact dynamics of quantum manybody systems. We use quantum leakage to distinguish two qualitatively different situations.
(i) In the small-leakage regime, stable periodic trajectories lead to long-time revivals of the manybody fidelity and oscillations of local observables. This result provides a more general mechanism of weak ergodicity breaking compared to quantum many-body scars in an arbitrary number of spatial dimensions. (ii) In the strong-leakage regime, exact quantum 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 tool 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 space is shown to have an 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 demonstrate that the revivals in the PXP model are much more robust to deformations of the Hamiltonian than was previously reported in the literature [11,12]. 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 manybody system onto a classical dynamical system; see Fig. 1. While some TDVP Ansätze may find unstable trajectories, thus motivating the term quantum scars [13], more generic embeddings may result in mixed phase space. Hence, it remains an open question which of the atypical eigenstates observed in different models [11][12][13]15,16,[22][23][24][25][26][57][58][59] 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 mention in the introduction, mixed phase space is known to influence the semiclassical limit of few-body quantum systems [29]. In this case, an intermediate spectral statistics is known to arise [60,61], and the system has an increased stability with respect to noise when it is prepared in the regular region of phase space [27]. 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 nontrivial, 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 manybody 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 ≪ 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 [62] and MPS-based dynamical mean-field theory [63].
In a complementary direction, TDVP was recently proposed as a method to capture thermalization. In particular, Ref. [36] demonstrates that physical properties saturate in the case of TFIM when bond dimension χ > 2. In addition, Ref. [37] studies 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 nonchaotic behavior in the TDVP dynamics, limiting the utility of the Lyapunov exponent. On the one hand, one may expect that mixed phase space does not persist for large bond dimensions, as increasing χ 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 [64] 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 the small-χ 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 establish numerous examples of anomalous thermalization. In particular, these include slow thermalization emerging from confinement [19,20], the existence of "slow operators" in thermalizing models [56], dynamical phase transitions [65], dependence of dynamics on the initial state in Floquet systems [52], mixed phase space in quantum maps [66], 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 analog of the KAM theorem. The robustness of quasilocal integrals of motion in MBL phases [39][40][41] provides an established particular case of the quantum KAM [3]. In a different setting, Ref. [38] demonstrates the existence of quasiconserved 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 toward its quantum counterpart. Indeed, recent work [49] conjectures 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. 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.

General framework a. 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 A. A. MICHAILIDIS et al.
PHYS. REV. X 10, 011055 (2020) 011055-12 TDVP in MPS manifolds is obtained in Ref. [33], here we aim to provide a simpler recipe which can be used to obtain the analytical form of equations of motion for states with a small number of variational parameters. We consider translationally invariant 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 σ αβ ðxÞ has one physical index, σ ¼ ↑; ↓ for spin-1=2 d.o.f., and two virtual indices α; β ¼ 1; …; χ, where χ is the (fixed) bond dimension. This tensor depends on a set of N real variational parameters, x ∈ 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 fx i g with i ¼ 1; …; K, where i labels sites within the unit cell. In order to simplify the notations, in the remainder of this Appendix we suppress the dependence of A on x i , A i ≡ Aðx i Þ.
In order for the state jψi [Eq. (3)] to have a system-sizeindependent norm, its unit-cell transfer matrix must be nondegenerate with the largest eigenvalue λ max ¼ 1. The unit-cell transfer matrix is obtained from a product of individual one-site transfer matrices: ðA1Þ over the entire unit cell: In the above notations, all transfer matrices operate in the

b. Action and gauge choice
The TDVP equations are obtained by extremizing the action (1) [33]. This symmetric choice of the Lagrangian is invariant under time-dependent unitary transformation jψi ¼ Ujψi, 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 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. (C2)-from both the equations of motion and the error [Eqs. (2) and (C1)].
In the case when the MPS is sparsely parametrized and satisfies the normalization condition, it is often possible to calculate analytically hψj _ ψi and hψjHjψi. In such cases, the calculation of equations of motion is greatly simplified compared to Ref. [33]. Below, we concentrate on such a case and derive the general form of equations of motion.

c. Hamilton's equations of motion
In the case when the MPS manifold is parametrized by complex parameters, the complex conjugate pairs z i andz i can be interpreted as momenta and coordinates. Since we use an explicitly real parametrization, we also separate the NK parameters fx i g into NK=2 coordinate variables fθ i g and NK=2 phase variables fϕ i g, where N is the number of real on-site parameters that are assumed to be even. The amplitude variables satisfy the condition of vanishing overlap between the wave function and tangential vector in any of θ direction: In contrast, the phase variables have a nonvanishing overlap characterized by the set of real functions f κ ðfθ i gÞ that depend solely on amplitudes. It is convenient to define the matrix which can be understood as the imaginary part of the manifold metric g ij ¼ h∂ x i ψj∂ x j ψi.
Extremizing the action [Eq. (A3)] over phase variables fϕ i ðtÞg gives equations of motion for fθg: while extremizing the action over amplitude variables yields equations of motion for fϕg: thus making apparent their symplectic structure. The equations of motion (EOMs) obtained above are determined by the expectation value of the Hamiltonian and the set of functions f κ ðfθ i gÞ. Below, we discuss the calculation of these ingredients in the specific case of PXP and TFIM.

TDVP for the PXP model
The MPS Ansatz used for the PXP model [Eq. (5)] has a left-canonical form, since the matrix A satisfies the condition In such a case, all single-site transfer matrices T i have largest eigenvalue λ max ¼ 1, with dominant left eigenvector ðL i j ¼ ð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. EOMs for PXP with K = 2 Ansatz
For the two-site MPS Ansatz, we get the following transfer matrix: cos θ 1 cos 2 θ 2 sin θ 1 0 0 cos θ 1 sin θ 1 cos θ 1 cos 2 θ 2 sin θ 1 0 0 cos θ 1 sin θ 1 cos 2 θ 2 sin 2 θ 1 0 0 and the corresponding right dominant eigenvector: jRÞ ¼ cos 2 θ 2 tan θ 1 ðcos −2 θ 2 tan −1 θ 1 ; 1; 1; tan θ 1 Þ: ðA11Þ Using the expression for the left dominant eigenvector, we calculate their overlap: that enters as a normalization factor in all expressions. The functions f 1;2 ðfθgÞ are calculated as follows. We explicitly use translational invariance of the MPS state to obtain ðA13Þ where the tensor ∂ ϕ 1 A 1 ≡ ∂ ϕ 1 A σ ðθ 1 ; ϕ 1 Þ can be obtained from Eq. (5). Replacing the environments by the dominant left and right vectors, we get the following expression: where we define the matrix T ∂ ϕ 1 as Using the explicit form of the matrix T ∂ ϕ 1 , we obtain the following expression for f 1 : Repeating the calculation for f 2 gives us 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 calculation results in matrices operating in the double virtual space: The matrix H 2;1;2 PXP where the operator σ x acts on the A 2 site can be obtained from Eq. (A18) by replacing θ 1 ↔ θ 2 and ϕ 1 ↔ ϕ 2 . The case of the chemical potential can be treated in a similar way, but the resulting matrices H 1;2 μ z are independent of variational parameters: 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: Summing these expressions, we obtain Eq. (10) in the main text. For C ¼ 2, the resulting equations of motion are obtained by substituting the specific form of the functions f 1;2 and the Hamiltonian density into Eqs. (A7) and (A8). 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.
The emergence of the flow-invariant subspace ϕ 1;2 ¼ 0 is not restricted to the 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 timereversal 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 σ y matrices in the Hamiltonian). Indeed, particle-hole symmetry requires that the operator C ¼ Q σ z i anticommutes with H, or in other words CHC ¼ −H. In addition, if we use the fact that the Hamiltonian does not change under complex conjugation, we may write hψjHjψi ¼ −hψjCHCjψi Ã ¼ −ðhψjCÞ Ã HðCjψiÞ Ã . Now we notice that ðCjψiÞ Ã is equivalent to the MPS state jψi with ϕ i → −ϕ i . Hence, when ϕ i ¼ 0 we obtain hψjHjψi ¼ −hψjHjψi ¼ 0 for any values of θ i . This result is sufficient to give a flow-invariant subspace _ ϕ i ¼ 0 when ϕ i ¼ 0 according to Eq. (A8). 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 σ x matrices. This result motivates the choice of the deformation H μ 3 considered below and also in the main text for the K ¼ 3 unit cell. In passing, we note that all terms in the quasilocal deformation conjectured by Ref. [49] to give perfect quantum scars satisfy both of these symmetries, thus leaving the flowinvariant subspace intact.
Finally, we note that, despite the angles ϕ 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 (afterward, they can be set to zero). Indeed, if one sets the phases ϕ i to zero before taking the derivatives in Eq. (A7), 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. (C1) below] with respect to θ 1;2 .

c. 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 ðθ 1 ; θ 2 ; θ 3 Þ. As we discuss above, such a derivation is more conveniently done with the help of conjugate variables ðϕ 1 ; ϕ 2 ; ϕ 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 the 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 jRÞ ¼ rð1=r; − cos θ 1 sin θ 1 ; − cos θ 1 sin θ 1 ; −sin 2 θ 1 Þ; where r ¼ cos 2 θ 2 cos 2 θ 3 þ sin 2 θ 3 cos 2 θ 2 sin 2 θ 1 − 1 : Using the form of the right dominant vector, we calculate the normalization factor The functions f i read SLOW QUANTUM THERMALIZATION AND MANY-BODY … PHYS. REV. X 10, 011055 (2020) The local terms from the PXP Hamiltonian and deformation H μ 3 [Eq. (13)], contracted with the environment, result in the following matrices: where the indices i are understood to be mod 3 and we use the following shorthand notations: Substituting these matrices in an expression analogous to Eq. (A20), 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 ϕ i ¼ 0. Plugging hψjH PXP þ H μ 3 jψi and expressions for f i into Eq. (A7), we obtain the equations of motion for variables θ i . Since we are interested in the flow-invariant subspace, we set ϕ i ¼ 0, resulting in the following equations of motion: where the factors that multiply the derivatives read

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 where the approximation works for small values of δt. The operator on the right-hand side can be written as a translationally invariant matrix product operator (MPO) [67] with bond dimension χ ¼ 2. This MPO has two physical indices, σ 1 , σ 2 ¼ ↑; ↓, and can be written as Values of the angles are found to be ðθ;φ;χÞ ¼ ðh x δt; h z δt; 2Jδ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 the j↑↑…i state. This choice results in the MPS with χ ¼ 2 that is specified by matrices: Now we allow the angles to be independent parameters, obtaining an MPS Ansatz with three variational parameters. Two of these parameters are phase variables ðφ;χÞ, and one is an amplitude variableθ. Hence, it is natural to extend this Ansatz by adding yet another amplitude variable ξ that is conjugate toχ. This extension can be done by replacing e −iχ=2 → cosξe −iχ=2 and e iχ=2 → sinξe iχ=2 . After such an extension, supplemented by a multiplication with the phase e iφ and a shift in parameters Eq. (A32) turns into the MPS Ansatz (16) used in the main text. The MPS Ansatz (16) is obtained as the twist of a product state by a χ ¼ 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 TFIM . Indeed, one can explicitly check that by restricting the angles to values ðχ; ξÞ ¼ ð0; π=4Þ we obtain a generic translationally invariant product state parametrized by (θ; ϕ).
After justifying the form of the Ansatz, we discuss its equations of motion. The unit cell consists of the K ¼ 1 site; thus, the transfer matrix can be obtained straightforwardly. Its largest eigenvalue reads 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 nondegenerate. Crucially, the largest eigenvalue is generally different from one, λ max ≠ 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 nonnormalized states [44], we prefer to normalize each tensor in the numerical implementationsÃ → A= ffiffiffiffiffiffiffiffi ffi λ max p , so that λ 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 evaluation 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. [68].

APPENDIX B: TDVP FOR TREE TENSOR STATES
In this section, we introduce the TTS Ansatz for trees of arbitrary connectivity. Afterward, 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 reduction allows one to use the method described in Appendix 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 [69].

Tensor tree state Ansatz
The basic building block of the TTS Ansatz is the tensor A σ z i 1;…;i C with C virtual indices i 1;…;C and one physical index σ z ; see Fig. 8(a). Out of C ¼ 3 virtual indices shown in Fig. 8(a), one is designated as an "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: 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, σ z ¼ ↓. In turn, the quantity I ¼ P C−1 j¼1 i j in Eq. (B1) 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 ¼ Á Á Á ¼ i C−1 ¼ 0, can the given site be in the excited state, thus implementing the constraint of no adjacent ↑ states.
We select 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. (B1) is invariant under permutations of all remaining input legs: where P ∈ 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. (B1) respects the Rydberg blockade constraints. This constraint imposes restrictions on the nonzero elements of the tensor A. The particular form of these elements in Eq. (B1) is dictated by the normalization condition. For the bipartite tree, one can check that this Ansatz generates a wave function 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. (B9). Numbers n 1;2 count the number of ↑ states in the corresponding sublattice, which are assigned the following weights: 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 ↓ i ¼ cos θ i and c ↑ i ¼ −ie −iϕ i sin θ i . However, in the presence of a constraint, these mean-field weights lead to a wave function with the norm vanishing as kjψik ∝ 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 for arbitrary values of θ; ϕ. This condition is enough to allow for normalization in the thermodynamic limit. Moreover, it greatly simplifies the calculation of expectation values, as we discuss below.

Mapping the tensor tree to a one-dimensional lattice
As explained in Appendix A 1 c, the main ingredient of TDVP are one-point correlation functions of the form hψj∂ x k jψi (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
where the indices inside parentheses are vectorized. In the first step, we show a part of the tree which includes the site k, where the local operator acts. The gray colored vertices imply contracted tensors in those sites. In the second step, the gauge symmetry (b) is used, and the whole tree is reduced to a semi-infinite one-dimensional tensor network (the green line is a guide to the eye). In the final step, the infinite product of local transfer matrices is replaced by the right dominant eigenvector of the matrix; see Eq. (B8). identity ðLj a;b ¼ δ a;b . We also define the generalized transfer matrix that includes the local operator [see Fig. 8(d)]: 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 where ðL l L m j ≡ ðL l j ⊗ ðL m j. The subscript indicates which input bonds are traced by the corresponding vector, while the transfer matrix T C is applied to the output bonds (i C , i 0 C ); see Fig. 8(d). Equation (B8) 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. (A13)] 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. (B8).
The normalization factor N of Eq. (B3) is derived by calculating hψjψi, which is a trivial modification of Eq. (B8) and corresponds to

Comparing variational wave function on a lattice and TTS
In the previous section, we explain 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 wave function for the square lattice is also the form of Eq. (B3) but with a different normalization factor. Such a wave function 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 [69].
The calculation of one-point correlation functions hψjO k jψi in PEPS is known to be a sharp P-hard problem for generic PEPS states [70] 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 dynamics on the true variational wave function for the lattice (PEPS) and the TTS. We perform the comparison by parametrizing the PEPS Ansatz as Eq. (B3). We use a finite system with periodic boundaries (4 × 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 a detailed justification of the fidelity bound. Afterward, we apply quantum leakage to characterize the different periodic trajectories observed in the deformed PXP model.

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 disagreement 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: 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 the explicit calculation of two-body correlators, it is not possible in the FIG. 9. We derive the general bound on the norm of the vector jφðtÞi that represents the mismatch between the exact wave function and its TDVP "image" after the TDVP evolution returns to its initial state. We note that the projection of the exact evolution onto the TDVP manifold does not coincide with the TDVP trajectory at finite times. present case. Nevertheless, the use of the gauge of Eq. (A3) allows one to cancel all disconnected correlators. This cancellation can be accounted by replacing hÃi → hÃi c in Eq. (C1), where the connected correlators are defined as Upon substitution of the connected correlators, the error scales as Λ 2 ∝ 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 γ 2 ¼ Λ 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 Ref. [33] 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 is performed in Ref. [13] for the flow-invariant subspace (ϕ i ¼ 0, hHi ¼ 0) of the PXP model.
We are interested in calculating the correlators defined in Eq. (C2). We show how to perform the resummation of two-point correlators of a uniform MPS, initially displayed in Ref. [33]. Extensions to larger unit cells follow straightforwardly by using the transfer matrix of the corresponding unit cell. For two parameters a and b, we get ðC5Þ where n and m are lattice site labels and in the second line we replace the summation over n and m by a sum over n and q ¼ jn − m − 1j. From such expressions, it is evident that one has to resum geometric series of the form P q ½T u:c: q , which is possible if the operator has spectral radius ρðT 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: where we define the projector onto the dominant subspace P ¼ jRÞðLj=ðLjRÞ, its complement Q ¼ 1 − P, and introduce the matrix T as Note that the same resummation formula holds for any size of the unit cell. Substituting Eq. (C6) into Eq. (C5), we get The first line contains the resummation of disconnected correlator, and these terms are canceled 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 hψjH 2 jψi and h∂ b ψjHjψi are also calculated using Eq. (C6). For the Hamiltonian which can be written as a sum of local operators H ¼ P i h i , with h i having support on an finite number of sites, chosen to be two for simplicity, we find These definitions are valid under the assumption that the Hamiltonian has only two-site terms, but generalization to longer-range terms in the Hamiltonian is straightforward. The disconnected correlators in Eqs. (C2)-(C4) are calculated as Collecting together these expressions with Eqs. (C8)-(C10), we observe that terms in the connected correlators that scale as L 2 vanish. Thus, the instantaneous error scales linearly with the system size, Λ 2 ∝ L.

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 jψ½fx a ðtÞgi and exact unitary evolution e −iHt jψ½fx a ð0Þgi 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. Because of nonzero leakage, the quantum and TDVP evolution would disagree after a single period.
In what follows, we obtain a lower bound on the manybody fidelity [Eq. (6)] after the time T, i.e., the period of the trajectory.
To obtain the bound, we represent the exact wave function jΨi that follows the exact evolution according to the Schrödinger equation i∂ t jΨðtÞi ¼ HjΨðtÞi as a sum of two terms: where jψ½fx a ðtÞgi belongs to the variational manifold and is obtained from the TDVP equations of motion, while jφðtÞi is the error vector defined as the difference between these two wave functions; see Fig. 9.
Expressing the error vector via jΨðtÞi and the TDVP wave function, we obtain the following equation of motion: where we explicitly use the Schrödinger equation for jΨðtÞi. The first term in the last line transports the error forward in time and does not change the norm of jφðtÞi.
The second term describes the change of the error due to mismatch between exact evolution jΨðtÞi and its TDVP projection jψ½fx a ðtÞgi, and its norm coincides with the quantum leakage defined in Eq. (C1). Calculating the change in the norm of vector jφi from Eq. (C17), we obtain where we omit time dependence of φ to simplify notations.
Using the triangle and Cauchy-Schwartz inequalities, we bound the growth of the norm of the error vector as where kφk ¼ ffiffiffiffiffiffiffiffiffiffiffi hφjφi p is the norm of the error vector. Using this notation, we obtain an upper bound on the rate of increase of the norm of jφi: where the quantity ΛðtÞ is the instantaneous geometric error that was already defined in Eq. (C1). 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: by the integrated geometric error I t .
The Fubini-Study metric or quantum angle can be defined as SLOW QUANTUM THERMALIZATION AND MANY-BODY … PHYS. REV. X 10, 011055 (2020) 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 γ can be visualized as an angle between the two lines. The three vectors jΨðtÞi, jψ½fx a ðtÞgi, and jφðtÞi 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, Assuming that we initialize 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: where F Ψ ¼ jhψ½fx a ðTÞgjΨðtÞij 2 ¼ jhΨje −iHT jΨij 2 and we use the fact that the quantum system is initialized on the TDVP periodic trajectory, jΨð0Þi ¼ jψ½fx a ð0Þgi ¼ jψ½fx a ðTÞgi. 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. (C24), this bound would appear useless, since I 2 t ∝ 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 ≪ 1 and the bound is still effective. In such a case, assuming the scaling form of the fidelity F Ψ ¼ e −f T L , with f T ≪ 1, we can expand ffiffiffiffiffiffi 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 jψ½fx a ðTÞgi, we can limit the system size L needed to be effectively in the thermodynamic limit, thus justifying the expansion used in deriving Eq. (C25). 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.

a. 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. (C8)-(C10). While in the case of the nondeformed 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 [68].

b. Calculation of Floquet exponent
In order to characterize the periodic orbits observed in the TDVP dynamics, we use the integrated leakage that is 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 [50]. The latter quantity is used to characterize chaotic dynamical flow in the case of large bond dimension TDVP [37]. In contrast, here we are dealing with periodic orbits. In this case, the Floquet exponent allows one 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 characteristic of the periodic orbit.
Since we do not use the Lyapunov exponent in this work, we keep the notation λ for the Floquet exponent. It is defined as where Λ max is the largest-norm Floquet multiplier of the orbit in the space of TDVP parameters, x ð0Þ a ðtÞ, and t ∈ ½0; T, where T is the orbit period. Intuitively, Λ max characterizes how the unit volume, chosen at some initial point on the periodic orbit, transforms after one traversal of the orbit. More specifically, Λ max is the largest-norm eigenvalue of the orbit Floquet matrix, ½J T ab ¼ ∂x a ðTÞ= ∂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 [50]:  C27) is calculated by solving a system of nonautonomous first-order differential equations. Since we have three variational parameters θ 1;2;3 , J T is 3 × 3 matrix. As expected, we find that the spectrum of J T always has an eigenvalue Λ ¼ 1 [50], and the remaining two eigenvalues satisfy jΛ max Λ min j ¼ 1, reflecting the symplectic (volumepreserving) nature of the flow.
If an orbit is stable and surrounded by a KAM torus, we find that all eigenvalues of the Floquet matrix have jΛj ¼ 1. Thus, the Floquet exponent vanishes, λ ¼ 0, signaling the stability of the orbit to small deformations. In contrast, for unstable periodic orbits, we have jΛ max j > 1 and jΛ min j < 1. This result implies that a unit volume after one period is elongated in one direction and compresses in a different direction.

c. 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 μ 3 ¼ −0.3. 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 μ 3 ¼ 0. Nevertheless, despite an increase in the size of the regular regions, the fidelity revivals are degraded by the deformation with large negative values of μ 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 λ ¼ 0 throughout the entire range of negative values of μ 3 .
Next, we turn to the case of μ 3 ¼ 0.3. The Poincaré section in Fig. 11(a) shows that the deformation with μ 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 are able to find using Newton's search method [47,50] (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 Γ P and the Floquet exponent. In particular, star 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 a finite value of λ. However, there are still many stable orbits present in the Poincaré section (e.g., the ones denoted by circle and square). These orbits have much stronger leakage, despite the vanishing Floquet exponent λ. This result confirms our conclusion that the Floquet  _ θ 2 < 0Þ reveals a smaller size of regular islands in phase space. Black symbols star, diamond, up-pointing triangle, circle, square, and down-pointing triangle show the locations of some periodic orbits. For each of the orbits, in (b) we show the value of the Floquet exponent, quantum leakage Γ P , and the dynamics of local observables n 1;2;3 . The physical period P is defined as the time after which the local observables return to their initial values and can be smaller than the period of the orbit for the MPS parameters. exponent is not obviously associated with quantum leakage.
On the other hand, the diamond 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 diamond 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.

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 ðξ ¼ 0.9; _ ξ < 0Þ. The periodic orbit with period T ¼ 2.09 crosses this plane at the point ðχ; ξ; ϕ; θÞ ¼ ð0.2607; 0.9; 4.888; 0.4308Þ, shown by "star" in Fig. 12. This orbit is stable surrounded by a small KAM torus, but nevertheless it does not give rise to many-body fidelity revivals. This result may be attributed to large leakage, Γ T ¼ 0.57.
In the main text, we illustrate 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),   13. Dynamics of entanglement of subsystems of size l at one boundary, S l ¼ SðlÞ=2, in TFIM for "fast" (solid lines) and "slow" (dashed lines) initial states defined in the text. Different colors correspond to different subsystem sizes that range from l ¼ 1 (magenta) to l ¼ 9 (red). Black lines correspond to the entanglement entropy of a half-infinite subsystem. The data are obtained with iTEBD with truncation error ϵ < 10 −5 . The inset shows the same data for l ¼ 5; …; 9 with rescaled axes, highlighting that the difference in entanglement spreading cannot be explained by the constant time shift and persists for largest subsystems. with parameters ðχ; ξ; ϕ; θÞ slow ¼ ð0.1; 0.9; 1.17; 3.86Þ and ðχ; ξ; ϕ; θÞ fast ¼ ð−0.32; 0.9; 1.6; 4.28Þ. These states can also be identified by local expectation values ðσ x ; σ y ; σ z Þ slow ¼ ð0.38; 0.85; 0.072Þ and ðσ x ; σ y ; σ z Þ 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 result 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 system 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.

Entanglement dynamics and ETH indicators in the deformed PXP model
We consider the following deformation of the PXP model, and values of the couplings are fixed to μ 2 ¼ 0.24 and μ z ¼ 0.4. This deformation is chosen in such a way that both terms in δ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 þ δ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ðδ i ; δ iþ1 Þ=maxðδ i ; δ iþ1 Þ yields hri ¼ 0.528, a value that is close to the Gaussian orthogonal ensemble prediction 0.5307 [71] (averaging is done over the middle 1=3 of the full many-body spectrum).
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 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 μ z ¼ 0.4 and slowly ramp up μ 2 at fixed energy density hHi=L ¼ 0. In practice, we can choose the initial state to be any phase space point of the periodic trajectory. We observe 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 ðθ 0 1 ; ϕ 0 1 ; θ 0 2 ; ϕ 0 2 Þ ¼ ð3.926; 0.289; 2.987; 0.121Þ. This state has bipartite entanglement entropy S ≈ 0.036. In the chosen MPS Ansatz, the bipartite entanglement entropy Sðθ 1 ; θ 2 Þ is independent of the phases. Thus, we can easily pick different initial states with exactly the same entropy by varying ϕ 2 while ϕ 1 is fixed by the energy density constraint. The initial states in Fig. 16 are generated by splitting the domain ½ϕ 0 2 − π=4; ϕ 0 2 þ π=4 to 200 points and then discarding the points for which the energy density constraint cannot be satisfied for any value of ϕ 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 the 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. ] strongly depends on initial conditions. All initial conditions are chosen to have the same energy density hHi=L¼0 and the same entanglement entropy. The light blue line corresponds to the periodic trajectory. The inset shows the correlation between the value of entanglement at late time and quantum leakage Γ t f =t 2 f , t f ¼ 10.