{'fraction_non_alphanumeric': 0.038416806660636155, 'fraction_numerical': 0.012091786815337185, 'mean_word_length': 4.57902054688707, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weak ergodicity breaking through the lens of quantum entanglement

Recent studies of interacting systems of quantum spins, ultracold atoms and correlated fermions have shed a new light on how isolated many-body systems can avoid rapid equilibration to their thermal state. It has been shown that many such systems can"weakly"break ergodicity: they possess a small number of non-thermalising eigenstates and/or display slow relaxation from certain initial conditions, while majority of other initial states equilibrate fast, like in conventional thermalising systems. In this chapter, we provide a pedagogical introduction to weak ergodicity breaking phenomena, including Hilbert space fragmentation and quantum many-body scars. Central to these developments have been the tools based on quantum entanglement, in particular matrix product states and tangent space techniques, which have allowed to analytically construct non-thermal eigenstates in various non-integrable quantum models, and to explore semiclassical quantisation of such systems in the absence of a large-N or mean field limit. We also discuss recent experimental realisations of weak ergodicity breaking phenomena in systems of Rydberg atoms and tilted optical lattices.


I. INTRODUCTION
Experimental progress in cold atoms, trapped ions and superconducting circuits 1,2 has generated a flurry of interest into foundational questions of many-body quantum mechanics, such as how a statistical-mechanics description emerges in isolated quantum systems comprising many degrees of freedom. In experiments with quantum simulators, the process of thermalisation can be conveniently probed in real time by quenching the system: one prepares a non-equilibrium initial state |ψ(0) , typically a product state of atoms, and observes its fate after time t, see Fig. 1(a).
Well-isolated systems can be assumed to evolve according to the Schrödinger equation for the system's Hamiltonian H. Thermalisation can be characterised by the time evolution of local observables, O(t) , or entanglement entropy, S(t)=−trρ A (t) ln ρ A (t). The latter is defined as the von Neumann entropy for the reduced density matrix, ρ A , of the subsystem A. Here we assume that the entire system is bipartitioned into A and B subsystems, and ρ A (t)=tr B |ψ(t) ψ(t)| is obtained by tracing out the degrees of freedom belonging to B. Entanglement entropy quantifies the spreading of quantum correlations between spatial regions as the entire system remains in a pure state.
The results of measurements described above are schematically illustrated in Fig. 1(b), which contrasts the behaviour of two large classes of physical systems: quantum-ergodic systems 3,4 and many-body localised (MBL) systems [5][6][7] . In the first case, parts of the system act as heat reservoirs for other parts, and an initial non-equilibrium state relaxes to a thermal equilibrium, with some well-defined effective temperature. This behaviour is reminiscent of classical chaotic systems, which effectively FIG. 1. Strong vs. weak breakdown of thermalisation. (a) Experiments probe thermalisation of isolated many-body systems using a quantum quench: the system is prepared in a simple initial state |ψ(0) and the dynamics of some local observable O and entanglement entropy S are measured during unitary evolution. (b) Strong breakdown of ergodicity caused by (many-body) localisation leads to qualitatively different dynamical behaviour compared to an ergodic system. This difference persists for broad classes of initial states. (c) In contrast to (b), weak ergodicity breaking results in strikingly diffferent dynamical behaviour when different initial states evolve under the same thermalising Hamiltonian. Quantum scarred systems, discussed in Sec. VI, are an important class of systems where such behaviour has been experimentally observed. (d): Manifestations of weak ergodicity breaking in entanglement and expectation values of local observables in the system's eigenstates. Red dots correspond to eigenstates of a weakly non-ergodic system, where enhanced fluctuations and outlier states (red diamonds) are visible, compared to typical eigenstates of conventional thermalising systems (black dots) whose properties vary smoothly with energy E. This chapter focuses on the mechanisms and physical realisations of the phenomenology summarised in (c)-(d).
arXiv:2108.03460v1 [cond-mat.quant-gas] 7 Aug 2021 "forget" their initial condition in the course of time evolution. In contrast, in MBL systems local observables reach some stationary value which is non-thermal, retaining the memory of the initial state. As a result, systems with strong ergodicity breaking can sustain much longer quantum coherence. For example, the information stored in certain local observables would not decay in MBL systems, while it would rapidly decohere in an ergodic system [ Fig. 1(b), bottom panel]. The behaviour of local observables is mirrored by that of quantum entanglement: in thermalising systems, entanglement typically grows linearly with time, S(t)∼t, corresponding to quasiparticles moving at a finite speed 8 . In MBL systems, quasiparticles are localised, with residual, exponentially-decaying interactions that lead to slow, logarithmic increase of entanglement, S(t)∝ log t 9-11 , see the top panel of Fig. 1(b).
The generic behaviour of quantum-chaotic systems sketched in Fig. 1(b) is expected to hold irrespective of the chosen initial state. This is a consequence of the Eigenstate Thermalisation Hypothesis (ETH) 3,4 , a powerful conjecture governing the behaviour of quantum-ergodic systems -see Box 1. In contrast, the present chapter focuses on recent theoretical works and experiments which pointed to the existence of a new class of behaviour, which can be loosely thought of as "intermediate" between thermalisation in fully chaotic systems and strong ETH breakdown like in MBL systems. In such weakly non-ergodic systems, certain initial states show relaxation to thermal ensembles, yet other states exhibit non-stationary dynamics including persistent oscillations, illustrated in Fig. 1(c). Moreover, such systems also weakly violate the ETH in their eigenstate properties: e.g., they exhibit a few "outlier" states with anomalously low entropy and local-observable matrix elements which do not vary smoothly with energy, see Fig. 1(d).
This new regime of ergodicity breaking attracted broader attention after the experimental observation of dynamical revivals in large-scale Rydberg atom quantum simulators 12 . Weak ergodicity breaking observed in such systems was subsequently named quantum many-body scarring 13,14 , highlighting its analogy with chaotic stadium billiards, which had been known to host nonergodic eigenfunctions bearing "scars" of classical periodic orbits 15 . The discovery of many-body quantum scars fuelled a broader quest to understand physical consequences of weak ergodicity breaking. Importantly, during this quest, it came to light that some examples of weak ergodicity breaking had previously been known. For example, non-integrable models such as the Affleck-Kennedy-Lieb-Tasaki (AKLT) spin chain, had been rigorously proven to possess non-thermal eigenstates 16,17 , which are now understood to share a similar algebraic structure with the non-thermal eigenstates in Rydberg atom systems. Other, more general, mechanisms of weak ergodicity breaking, such as Hilbert space fragmentation 18,19 and "embedding" constructions 20 have also come into light and are being experimentally probed 21 .
This chapter provides a pedagogical introduction to weak ergodicity breaking phenomena, expanding upon a recent short overview in Ref. 22. As will become clear from the many physical examples presented below, the phenomenology of weak ergodicity breaking is quite rich and still rapidly evolving as this chapter is being written. Instead of giving an exhaustive review of all these developments, the aim is to highlight common themes between currently known examples of non-thermal dynamics, in particular focusing on the insights obtained by studying quantum entanglement in such systems. In Sec. II we start by introducing the relevant methodology based on matrix product states, that has successfully been used in recent works to analytically construct non-thermal eigenstates and to define the semiclassical limit of non-integrable many-body systems. In Sec. III we discuss in detail three main mechanisms of weak ergodicity breaking and highlight their realisations in physical systems. Sec. IV introduces the so-called PXP model, which has played one of the key roles in understanding weak ergodicity breaking in Rydberg-atom experiments. In Sec. V we relate the semiclassical quantisation of many-body systems with the timedependent variational principle applied to manifolds spanned by tensor network states. In Sec. VI we explain how these ideas allow to explore parallels between weak ergodicity breaking in many-body systems and quantum scars in few-body systems. Finally, in Sec. VII we discuss two experimental platforms -Rydberg atoms and cold atoms in a tilted optical lattice -that have recently observed signatures of weak ergodicity breaking. Conclusions and some open questions are presented in Sec. VIII. 3. For a large finite subsystem A of an infinite ETH system, the reduced density matrix ρ A is equal to the thermodynamic density matrix at the effective temperature set by the energy of the corresponding eigenstate.
The equality between thermal and reduced density matrices (3) implies that the entanglement entropy scales as the volume of region A. For example, in one spatial dimension, volume law implies S∝L A . This reflects the fact that ergodic eigenstates are highly entangled, and agrees with the intuition that the ETH eigenstates are similar to random vectors. In some physical systems considered below, we will encounter global kinetic constraints, similar to those occuring in classical glasses 26 . In such systems, the Hilbert space is globally constrained, and a "random" vector is understood to be one compatible with the constraints. Consequently, the entropy of such constrained random vectors can differ 27 from the so-called Page value attained by random qubit states 28 .
In systems which obey the so-called "strong" ETH 29 , the above properties are expected to hold for all states 30 . On the other end, maximal violation of ETH is known to occur in integrable systems 31 and MBL systems 32 , both of which have macroscopic numbers of conservation laws. Between these extremes lie the weak ergodicity breaking phenomena discussed in this chapter, where the majority of eigenstates follow the predictions of the ETH while a smaller number of eigenstates (e.g., polynomially many in system size) violate these properties. Similarly, out-of-equilibrium dynamics from most initial conditions results in fast relaxation, consistent with ETH expectations, while special initial conditions can lead to non-stationary dynamics at relatively long times.

II. MATRIX PRODUCT STATE METHODS
The investigation of weak ergodicity breaking in non-integrable quantum systems has greatly benefited from matrix product states (MPS) 33 , a formalism designed to compactly represent and perform algebraic manipulations on a class of quantum states that are weakly entangled. These methods are naturally suited to capture aspects of weak ergodicity breaking, as the main signatures of the latter are the suppression of entanglement or its growth rate compared to conventional thermalising systems. Recent works have demonstrated the utility of MPS methods in two new settings: (i) the MPS have allowed to exactly construct highly-excited eigenstates (i.e., eigenstates at finite energy density above the system's ground state) of a wide class of physical systems; (ii) the MPS have been used to effectively define the system's semiclassical dynamics by projecting the Schroödinger time evolution into the manifold of MPS states whilst conserving the total energy of the system. In this section, we provide a brief overview of this methodology which yielded much of the physical insights into weak ergodicity breaking phenomena presented in subsequent sections.

A. Towers of quasiparticles
Consider a one-dimensional quantum chain with a d-dimensional Hilbert space on each of the L sites. The many-body basis of the system is formed by tensor products of single-site Hilbert spaces, |σ 1 , σ 2 , · · · , σ L . Any state can be decomposed in this basis as |ψ = σ1,σ2,...,σ L c σ1σ2...σ L |σ 1 , σ 2 , . . . , σ L (1) by specifying its d L coefficients c σ1σ2...σ L . This illustrates the "exponential barrier" to studying many-body systems: one must specify ∼ exp(L) coefficients to fully describe a generic state of a system of size L. The idea of MPS is that each of the coefficients in Eq. (1) can be viewed as resulting from a product of matrices A The matrices A are defined over the so-called auxiliary space with dimension χ, but they also depend on the physical degrees of freedom σ i . Moreover, for open boundary conditions, the matrices also vary from site to site in general, hence they carry a label i. The boundary χ-dimensional vectors, b l and b r , determine the boundary conditions for the wavefunction. Thus, we have reexpressed a d L -dimensional tensor c σ1σ2...σ L in terms of d × χ × χ tensors A i . If the system is translation-invariant, the matrices are the same on every site (or between different unit cells, more generally) and the boundary vectors are replaced with a trace over the auxiliary space, Eqs.
(2)-(3) furnish an MPS representation for the wave function in Eq. (1). Intuitively, χ controls the degree of entanglement in the wave function. When χ=1, we have a mean-field description where degrees of freedom on different sites are independent of each other; increasing χ builds in quantum correlations between different sites. For a general state, rewriting the coefficients c σ1σ2...σ L in MPS form would require χ to be exponentially large in L. However, for states in one spatial dimension that obey the "area law" 34 , i.e., whose entanglement entropy obeys S ≤ const, χ is also bounded by a constant. This results in a major simplification when describing such states using the MPS. Moreover, the same language also offers a complete algebraic framework to manipulate the MPS, for example we can efficiently add two MPS states or calculate expectation values of local observables sandwiched between MPS states. The latter is achieved by expressing local operators in a similar representation known as "matrix product operator" (MPO) 35 .
Ground states of many important condensed-matter systems are known to possess elegant MPS representations. A notable example is the Affleck-Kennedy-Lieb-Tasaki (AKLT) 36 quantum spin-1 model: where S are the spin-1 operators and we assume periodic boundary condition, i.e., L + 1 ≡ 1. The ground state of this model has an exact MPS representation, Eq. (3), with where the labels +, 0, and − denote the S z = + 1, 0, and −1 spin-1 basis states, respectively, and σ ± , σ z are the spin-1/2 Pauli operators. Thus, for the AKLT state, the physical dimension d=3 and the bond dimension is only χ=2. This state and the model in Eq. (4) played an important role in establishing the existence of the Haldane gap in integer-spin chains 37 .
Beyond examples such as AKLT, where the MPS exactly describe certain states, it has been established more generally that ground-state wave functions of gapped local Hamiltonians can be approximated by an MPS with a small bond dimension χ 38 . This result can be intuitively understood from the fact that the ground states of gapped local Hamiltonians are necessarily weakly entangled 39 and hence well-represented by the MPS. However, this intuition does not immediately extend to highly excited states at a finite energy density above the ground state. As previously discussed in Sec. I, in a generic quantum many-body system, finite-energy-density states are expected to be governed by the ETH, hence their entanglement is expected to be high (i.e., scaling with the volume of the subsystem), so their MPS representation would not be efficient (as mentioned above, it would require χ ∝ exp(L)).
Recent work by Moudgalya et al. 17 on the AKLT model has shown that the MPS can nevertheless capture some ETH-violating eigenstates of a non-integrable model, regardless of how high in the energy spectrum such states occur. The technique used in these works in the construction of quasiparticle excitations above the ground state, originally introduced in works on tangent space methods [40][41][42] . Following Ref. 17, a single-site quasiparticle excitation with momentum k on top of a general MPS state, |ψ A , is given by where B [σj ] is a χ × χ matrix with physical dimension d and k denotes the momentum, see Fig. 2(a). In the framework of the Single-Mode Approximation 43 , the quasiparticles are usually described in terms of a single-site "quasiparticle creation operator" , which we denote in shorthand as |B = O |A . Note that this is a special case, as for example O could act on several neighbouring sites -such states have been considered in Ref. 44. To illustrate Eqs. (6), some examples of low-lying excited eigenstates of the AKLT include the Arovas "A" and "B" states 45 : where |ψ 0 is the AKLT ground state, {...} denotes the anticommutator, and we have omitted the normalisation factors. Another exact excited state, as shown by Moudgalya et al. 17 , is the spin-2 magnon state with momentum π,  (6) and (10). (b) Time-dependent variational principle captures the optimal projection of quantum dynamics onto a given variational manifold M, here parametrised by MPS. Limiting the bond dimension of the MPS to be low, this approach defines a semiclassical limit for quantum dynamics. The corresponding classical system is formed by equations of motions for the coefficients xa, which parametrise the projection of the wave function into the MPS manifold. The equations of motions follow from minimising the leakage outside the manifold, i.e., the deviation of −iH |ψ(t) fromẋa∂x a |ψ(t) , see Eq. (14).
All of these states can be written in the form of Eq. (6). For example, the state in Eq. (9) can be expressed in such a form using the ground state matrices A [σ] in Eq. (5) and the B [σ] matrices given by is the only non-trivial matrix, a direct consequence of the (S + ) 2 operator acting on spin-1.
In addition to single quasiparticles, multiple quasiparticle states can be described in the MPS formalism using multiple tensors. For example, a state with two quasiparticles described by tensor B with momentum k is given by 17 Expanding the square, we see that we can have two cases: when Bs act on the same site or different sites. Since in the AKLT chain we have O = (S + ) 2 , we will haveÔ 2 |ψ 0 = 0 and the expression simplifies Extending this formally to more quasiparticles appears straightforward. However, the difficulty that arises is that the successive applications of the quasiparticle operator lead to an increase in the bond dimension of the resulting state. Thus, one needs to finds a way to compress this state and compute its properties. These advances have been made in Ref. 17, allowing to evaluate the entanglement properties of highly excited states in the AKLT model, obtained by acting withÔ on the ground state an extensive number of times n∼L/2. The physical consequences of this for weak ergodicity breaking will be discussed further in Sec. III A below. Finally, while the emphasis in this section has been on analytical constructions of exact eigenstates containing finite density of quasiparticle excitations, one could envision employing similar methods to variationally construct quasiparticle towers in other models which are less tractable than the AKLT. Indeed, such studies have been used to numerically characterise dynamical properties of quantum spin systems such as the XXZ model, including quasiparticle dispersion relations and dynamical structure factor, albeit at much lower energies than the phenomena discussed here 46 .

B. Time-dependent variational principle
In the previous subsection, we mentioned that some non-integrable quantum spin chains have highly-excited energy eigenstates with a particularly simple structure, encoded in a local MPO acting on the ground state written as MPS. "Simple" structure means that such eigenstates have sub-volume law entanglement, i.e., they are much less entangled than random vectors in the Hilbert space (cf. Box 1). Now we ask a complementary question: if the system undergoes anomalously slow unitary dynamics, e.g., the spreading of correlations is inhibited compared to a thermalising system, could we compactly describe such dynamics using the MPS?
The natural framework to address the previous question is the time-dependent variational principle (TDVP), originally formulated by Dirac 47 , which yields the optimal projection of quantum dynamics onto a given variational manifold, M = span{|ψ(x) }, parametrized by some parameters x ≡ {x a } ∈ R N , Fig. 2(b). TDVP equations of motion are obtained by extremizing the action dt L, with an effective Lagrangian 48 In the present case, we assume that M is spanned by the MPS of some given bond dimension, i.e., by matrices A σ (x) like in Eq. (3) with bond dimension χ and the matrices depend on variational parameters, x∈R N . These parameters can be given an interpretation in terms of pairs of coordinate and momenta (we assume N is even). In Box 2, we motivate this approach using a very simple example of a single spin. Extremising the action in Eq. (12) results in Euler-Lagrange equations of motion for the x variables 48,49 : This set of equations constitutes an effective mapping of quantum dynamics onto a classical non-linear dynamical system. While TDVP has been successfully applied to manifolds spanned by the MPS 49 and finite tensor tree states 50 , in cases when it is possible to analytically calculate ψ|ψ and ψ|H|ψ , further insights can be obtained from studying the non-linear system using the tools of classical dynamical systems 51 and few-body chaos 52 . In practice, such studies are naturally limited to sufficiently simple manifolds which have low bond dimensions. Nevertheless, this approach proved extremely valuable in establishing an analogy between non-thermal eigenstates in Rydberg atom chains and their counterparts in few-body systems, justifying the name "quantum many-body scars" -these developments will be the subject of Secs. V and VI.
We note that there are some caveats to mapping quantum dynamics onto low-dimensional manifolds M. At late times, the error of the TDVP approximation necessarily grows because the exact quantum dynamics generally brings the system out of the variational manifold -see Fig. 2(b). Intuitively, the TDVP error for a given initial state is linked to the rate of entanglement growth for that state, since the MPS with low bond dimension can only successfully capture weakly entangled states. When the entanglement grows significantly, the state starts to require an MPS description with a bond dimension that extends beyond the variational manifold. This discrepancy is the error of the TDVP approximation, also known as quantum leakage 14,49 : This is the instantaneous rate at which the exact quantum wave function leaves the variational manifold. Equations of motion can be obtained by minimising the discrepancy between exact quantum dynamics and its projection onto M. Explicit computation of γ 2 is possible within the TDVP framework, but since it involves the square of the Hamiltonian operator, it contains information that goes beyond the TDVP equations of motion. Understanding γ 2 is key to understanding the relation between quantum dynamics and its classical counterpart. We postpone a more extensive discussion of the subtleties of the TDVP approach to many-body systems to Sec. V.

Box 2 | TDVP for a single spin
To illustrate the TDVP way of thinking, consider the simplest example of a quantum spin prepared in a product state |(ψ(θ, φ)) = cos θ |↑ + sin θe iφ |↓ , parametrised by angles θ and φ on the Bloch sphere. For concreteness, we can think of the spin Hamiltonian as H= σ x , corresponding to the magnetic field along the x-direction. Our aim is to find the equations of motion for the angles, θ(t) and φ(t), as the spin evolves under the influence of the magnetic field. The expectation value of the Hamiltonian in the given state is H = sin(2θ) cos φ. Moreover, we have ψ|ψ = −iφ sin 2 θ. Plugging this into Eq. (12), our Lagrangian takes the form L =φ sin 2 θ − H . We see that we can identify the canonical "position" and "momentum" variables as φ → q and sin 2 θ → p. Thus, we have mapped the quantum spin dynamics onto a classical dynamical system described by the equations of motionφ = 2 cot(2θ) cos φ andθ = − sin φ. Notice that the line φ=π/2 is stationary -this corresponds to the simple precession dynamics,θ = − , which in this case is easy to derive directly by solving the Schrödinger equation.
This example of a single spin can also be viewed as a particularly simple mean-field picture for the dynamics. This approach is generally insufficient to describe dynamics in more complex many-body systems and, as discussed in Sec. V, recent works have been exploring generalisations to variational manifolds beyond the mean-field using MPS states.

III. MECHANISMS OF WEAK ERGODICITY BREAKING
ETH-violating eigenstates have recently been theoretically identified in a variety of non-integrable quantum models, revealing a complex landscape of weak ergodicity breaking phenomena. The common feature of these models is the emergence of a decoupled subspace within the many-body Hilbert space, in general without any underlying symmetry. In this section, we focus on idealised models where such a subspace is perfectly decoupled from the rest of the spectrum, resulting in a decomposition of the Hamiltonian where H non−ETH is the non-thermalising subspace, exactly decoupled from the thermalising bulk of the spectrum, H thermal . The eigenstates that inhabit the subspace H non−ETH violate the ETH and have different properties compared to the majority of thermal eigenstates residing in H thermal . Below we elucidate three commonly encountered mechanisms that produce such decoupled subspaces. We would like to point out that in much of the recent literature, the non-thermalising eigenstates that span H non−ETH are commonly referred to as "quantum many-body scars". We will discuss the precise meaning of quantum many-body scarring in Sec. VI.

A. Spectrum generating algebra
Using the MPS techniques introduced in Sec. II, Moudgalya at al. 17 , building on the early work by Arovas 45 , analytically constructed a tower of exact eigenstates in the AKLT model in Eq. (4), which were shown to have sub-volume entanglement entropy, thus providing the first rigorous demonstration of the strong-ETH violation. Later on, Schecter and Iadecola 53 introduced a tower of exact eigenstates in a family of non-integrable spin-1 XY models as well in a a particular non-integrable spin-1/2 model that conserves the number of domain walls 54 . In these and many other similar examples that followed, the non-thermal eigenstates form a so-called "spectrum generating algebra" (SGA), first introduced in the context of high energy physics 55 and subsequently applied to the Hubbard model 56,57 .
To define the SGA, we follow the presentation given by Mark, Lin and Motrunich 58 . Suppose we have a Hamiltonian H; a linear subspace W ; a state |Ψ 0 ∈ W , which is an eigenstate of H with energy E 0 ; and an operator Q † such that Q † W ⊂ W and Then it can be easily proven that the family of states as long as they are non-vanishing vectors, are eigenstates of H with eigenvalues E 0 +nω. The non-thermalising nature of these states becomes non-trivial when W is not the entire Hilbert space and Q † is not associated with a symmetry of the Hamiltonian. The latter crucially distinguishes these non-thermal states from somewhat similar "eta-pairing" states in the Hubbard model 56,57 , a point discussed further at the end of this subsection.
As an example, in the AKLT model one of the possible SGA operators is a spin-2 magnon excitation that we have previously seen in Eq. (9), with This generates the tower of states in Eq. (17) for n = 0, ..., L/2, which extend from the ground state (n=0) up to the ferromagnetic state |1, 1..., 1 (n=L/2), and the state in Eq. (9) is a member of this family. It is also worth noting that these states are not unique, as they have total spin s = 2n, hence we can obtain equivalent spin-rotated versions of these states by applying the SU(2) spin-lowering operator. Generally, it is clear that states of the form in Eq. (17) can be non-thermal, provided |Ψ 0 is sufficiently "simple" (e.g., the ground state of H if the latter is gapped) andQ † is a local operator. The AKLT tower of states in Eqs. (17)-(18) is the first example of a rigorous construction of a non-thermal family of states in a non-integrable model that are not protected by a global symmetry. Any Hamiltonian with the SGA property explicitly decomposes as in Eq. (15), where H non−ETH contains the tower of states obtained by the action ofQ † . Note that this algebraic structure is not powerful enough to fully diagonalise H; indeed, H could well be a non-integrable Hamiltonian with an energy spectrum obeying the Wigner-Dyson statistics.
Physically,Q † creates a wave packet corresponding to a "quasiparticle" excitation (e.g., a magnon), and repeated applications ofQ † create a condensate of such quasiparticles. In a class of frustration-free models which include the AKLT 44 , the quasiparticle condensates are non-thermal, e.g., their entanglement entropy scales logarithmically with the subsystem size: Another important clue when looking for non-thermalising SGA states is that they appear at energies which are integer multiples of ω. Thus, SGA eigenstates can be detected as regularly spaced entropy-outliers, with entanglement much lower than the ETH volume law, as depicted in Fig. 1(d). Perhaps more practically, the SGA states also have very low entanglement rank, i.e., many eigenvalues of ρ A strictly vanish, which is how these states were originally identified numerically in Ref. 17. We mentioned that the SGA traces back to the so-called "eta-pairing states" in the Hubbard model 56,57 . For eta-pairing states, the subspace W is the entire Hilbert space and Q † essentially corresponds to a symmetry of the Hamiltonian 59 . However, it is possible to perturb the Hubbard model 60,61 in a way which destroys the aforementioned symmetry and makes these states similar to the SGA states in the AKLT model discussed above. Note that the high degree of symmetry in Hubbard-type models can also lead to disconnected subspaces which host free-particle eigenstates with ballistic dynamics 62 . This "fragmentation" of the Hilbert space is a much more general phenomenon beyond just the Hubbard model, as we discuss in the following section. Finally, we note that generalisations of the SGA construction can be found in a number of models in various dimensions 54,58,63-68 (see Box 3). It has also been pointed out that the SGA can arise in open quantum systems in the presence of dissipation or driving 69,70 .

Box 3 | Constructions of non-thermal towers of eigenstates
Towers of non-thermal eigenstates can be systematically constructed using the "tunnels-to-towers" approach in Ref. 68 based on generators of the Lie algebra of a symmetry group G. For simplicity, consider the case when G is SU(2) and we have a model defined by an SU(2)-symmetric Hamiltonian, H sym . The generators of the symmetry {Q + , Q − , Q z } are associated with the corresponding su(2) algebra. The spectrum of H sym is organised into 'tunnels' of degenerate eigenstates, with the same eigenvalue for the Casimir Q 2 but different eigenvalues for Q z . One can 'move' between states in a tunnel using Q ± .
Now imagine perturbing the model by adding H SGA ∝ Q z . This perturbation preserves the eigenstates, but breaks the degeneracy of the tunnels. Instead, states in each tunnel get promoted to 'towers' and acquire an evenly spaced harmonic spectrum because of the SGA, [Q z , Q ± ] = ±Q ± . Finally, we can further add a symmetry-breaking term H SB , such that it annihilates a specific tower of states but generically breaks all symmetries and mixes between the other states so as to make the rest of the spectrum thermal. In the full model, H = H sym + H SGA + H SB , our chosen tower of states are a collection of non-thermalising eigenstates, evenly distributed throughout the spectrum but not protected by a global symmetry, hence we arrive at a similar phenomenology to the AKLT model discussed in the main text. Various extension of this construction to non-Abelian groups and their q-deformed versions are possible 68 . Figure B3: Tunnels-to-towers scheme for constructing non-thermal eigenstates from Ref. 68. To illustrate the tunnels-to-towers construction, consider the spin-1 XY-like model introduced in Ref. 53: A tower of non-thermalising states is built by the action of the same Q † as in the AKLT model 58 : where |Ω is the fully polarised down state |Ω = | − − − · · · − . Note that Q + (and the corresponding Q − ), together with form an su(2) algebra but they are distinct from the spin-su(2) operators. The first term ∝ J breaks Q-SU(2) symmetry and annihilates the tower, the term hS z acts as H SGA and gives energy to the states in the tower, while the term ∝ D commutes with Q z and Q + . The third neighbour term is added to break a non-local SU(2) symmetry for which the states belonging to the tower are the only states in their symmetry sector.

B. Hilbert space fragmentation
Spectrum generating algebra relies on the tower operator Q † to construct the non-thermalising subspace. We now describe a related mechanism which produces exactly embedded subspaces, but does not require a priori knowledge of Q † . For a Hamiltonian H and some arbitrary vector in the Hilbert space, |ψ 0 , the Krylov subspace, S, is defined as the set of all vectors obtained by repeated action of H on |ψ 0 , Readers familiar with numerical linear algebra will recall the same subspace S is used in iterative methods for finding extremal eigenvalues of large matrices, such as the Arnoldi and Lanczos algorithms. By definition, S is closed under the action of H. While |ψ 0 in Eq. (22) can in principle be an arbitrary state, in physics one is primarily interested in initial product states (also called the "root" states), which are more easily preparable in experiment. For a generic non-integrable Hamiltonian H without any symmetries, one expects that S (H, |ψ 0 ) for any initial product state |ψ 0 is the full Hilbert space of the system. If H has some symmetry (and assuming |ψ 0 an eigenstate of the symmetry), S (H, |ψ 0 ) is expected to span all states with the same symmetry quantum number as |ψ 0 . Surprisingly, it has been shown that in many cases the system can exhibit fragmentation i.e., even after resolving the symmetries, S (H, |ψ 0 ) does not span all states with the same symmetry quantum numbers as |ψ 0 18,19,71,72 . Following the notations in Ref. 72, we can formally state this as where s labels the distinct symmetry quantum numbers, # (s) denotes the number of disjoint Krylov subspaces generated from product states with the same symmetry quantum numbers, and |ψ (s) i are the root states generating the Krylov subspaces. Note that the root states in Eq. (23) are chosen such that they generate distinct disconnected Krylov subspaces (the same subspace can be generated by different root states). If one (or more) Krylov sectors are non-thermalising, we recognise Eq. (23) is of the same form as our previous Eq. (15).
Fragmentation as in Eq. (23), where the total number of Krylov subspaces is exponentially large in the system size, was recently shown to always exist in Hamiltonians and random-circuit-models with conservation of dipole moment (see Box 4). Before illustrating this for a particular model, we note that, generally, one can distinguish between "strong" and "weak" fragmentation, depending respectively on whether or not the ratio of the largest Krylov subspace to the Hilbert space within a given global symmetry sector vanishes in the thermodynamic limit. Strong (resp. weak) fragmentation is associated with the violation of weak (resp. strong) ETH with respect to the full Hilbert space. Moreover, different fragments can exhibit vastly different dynamical properties. For example, some fragments (even though exponentially large in system size) may be integrable, while others may be non-integrable. Among the non-integrable ones, more subtle ETH-breaking phenomena are also possible 73 .
To illustrate fragmentation, consider the following model of fermions hopping on an open 1D chain while preserving their centre of mass 72 : where c † j , c j are the standard fermion creation/annihilation operators on site j. This model has a variety of physical realisations, which we discuss at the end of this section. We now demonstrate the dynamical fragmentation in this model, closely following the presentation in Ref. 72.
At filling factor ν = 1/2, i.e., with half as many fermions as sites in the chain L, it is convenient to split the chain into 2-site cells assuming L is even. For each cell, define the new degrees of freedom: It is fruitful to name these composite degrees of freedom: |+ , |− are called "fractons", |+− , |−+ are "dipoles", and |↑ , |↓ are "spins". By inspecting the possible action of the Hamiltonian in Eq. (24) on these states, we find several types of allowed processes. For example, the following processes can be interpreted as free propagation of dipoles when separated by spins: Similarly, the following processes imply that a fracton can only move through the emission or absorption of a dipole: Here, Eqs. (26)-(29) resemble the rules restricting the mobility of fracton phases of matter 74 . However, in contrast to usual fracton phenomenology, here the movement of fractons is also sensitive to the background spin configuration. For example, the fracton in the configuration |· · · ↓ + ↑ · · · can move by emitting a dipole (see Eq. (28)) while that in the configuration |· · · ↑ + ↓ · · · cannot. The model in Eq. (24) has several symmetries, most notably the total charge (number of fermions) and the total dipole moment are conserved, with Q j ≡n 2j−1 +n 2j − 1, where j is the unit cell index and 2j − 1, 2j are the site indices of the original configuration. Exponentially many of Krylov subspaces in the model in Eq. (24) are one-dimensional frozen configurations. For instance, the Hamiltonian vanishes on any product state that does not contain the patterns " · · · 0110 · · · " or " · · · 1001 · · · ", since those are the only configurations on which terms of H act non-trivially. The state |1111000011110000 . . . . . . 1111000011110000 is one example of a static configuration that is an eigenstate. In terms of the composite degrees of freedom we can equivalently consider configurations with only +, −, and no spins, such as with a pattern that alternates between + and − with 'domain walls' that are at least 2 sites apart. Once again, it is easy to see that all terms of the Hamiltonian vanish on these configurations: since there are exponentially many such patterns, there are equally many one-dimensional Krylov subspaces.
In Ref. 72 it was found that large fragments (with dimension exponential in system size L) can be integrable (mappable to spin-1/2 XX spin model) as well as non-integrable. See Fig. 3 for an illustration of some of the disconnected sectors of the Hilbert space. The dynamics initialised in any of the states, e.g., in Fig. 3(a) cannot reach any of the states in plots (b) and (c).
We previously announced that the pair-hopping model in Eq. (24) arises in a number of different physical contexts: (i) it arises as the dominant hopping process for electrons in the regime of the fractional quantum Hall effect in a quasi-1D limit 73 ; (ii) it arises in the Wannier-Stark problem, i.e., spinless fermions hopping on a one-dimensional lattice, subject to a large electric field 75 ; (iii) it can be mapped to the following spin-1 fractonic model studied in Ref. 19 : (iv) at filling factor ν=1/3, one of the Krylov sectors can be mapped to the so-called PXP model, introduced in Eq. (35) below, which describes a chain of strongly-interacting Rydberg atoms. Note that in realisations (i), (ii) there are usually additional diagonal terms in the Hamiltonian that are of the same order as the hopping. While this does not affect the fragmentation, it may significantly impact other dynamical properties within the fragments. The realisation (ii) has recently been investigated as a platform for many-body localisation without disorder 76,77 .
Signatures of fragmentation have been observed in other models including, e.g., the Fermi-Hubbard model and its cousins 59,62 , various constrained models [78][79][80] , and bosons in optical lattices 81,82 . In the latter case, the Krylov subspaces are only approximately exact in the sense that there exist non-zero matrix elements which connect different Krylov subspaces but their magnitude is much smaller than the matrix elements within a given Krylov subspace.

Box 4 | Fragmentation in quantum circuits
Beyond Hamiltonian systems, fragmentation generally arises in models of 1D random unitary circuit dynamics, constrained to conserve both U(1) charge and its dipole moment 18,71 . Consider a model from Ref. 71 with a chain of S=1 quantum spins, with the local z-basis |+ , |− , |0 , and unitary gates which locally conserve charge Q = j S z j and dipole moment D = j jS z j . These intertwined conservation laws greatly restrict the allowed movement of charges, e.g., a single + or − charge on site x has dipole moment D= ± x. Such a charge cannot hop to the left or right, because this would change the net dipole moment by one unit. On the other hand, bound states of charges or "dipoles" of the form (−+) have net charge zero and net dipole moment D= ± 1 independent of position, and these can move freely. Additionally, dipoles can enable the movement of charges, because a charge can move if it simultaneously emits a dipole to keep D unchanged: The above phenomenology is very similar to the pair hopping model discussed in the main text, but following Ref. 18 , we realise it using Floquet circuits composed from -site unitary gates, which locally conserve Q and D. These -site gates are therefore block-diagonal, and they are applied periodically in time. Assuming for simplicity =3, we tile the chain by 3-site unitaries, staggered across three layers, i.e., B4 below (and similarly for U 2 , U 3 in red and green colour, respectively). The complete circuit consists of repeatedly applying U F some number of times. Note that U 1 , U 2 and U 3 can be chosen at random for a given realisation, but then remain fixed throughout the circuit.  Ref. 18 showed that local fractonic circuits of this type must have exponentially many Krylov sectors. To see this, note that any pattern that alternates between locally 'all plus' and locally 'all minus,' with domain walls in between at least sites apart, must be inert. These are states of the form | + + + + − − − − − + + + + · · · . The argument is similar to the one given in the main text for the Hamiltonian system. One can then lower bound the size of the disconnected subspace by dividing the system up into blocks of length , and allowing each block to be either 'all plus' or 'all minus.' This yields an inert subspace of dimension at least 2 L/ = c L , where c = 2 1/ . This bound is not tight 18 , but it proves the existence of an exponentially large, localised subspace for any finite gate size. We emphasise that the exponentially large number of sectors goes beyond simple symmetry considerations, which would only predict ∝ L 3 sectors, since the allowed values of quantum numbers are −L ≤ Q ≤ L and −L(L − 1)/2 ≤ D ≤ L(L − 1)/2, where L is the length of the chain.

C. Projector embedding
At the highest level of abstraction, we can ask if one could embed an arbitrary subspace into the spectrum of a thermalising system? For concreteness, assume we are given an arbitrary set of states |ψ i that span our target non-thermalising subspace H non−ETH = span{|ψ i }, which we wish to embed into a thermalising Hamiltonian H as in Eq. (15). This can be achieved via the "projector embedding" construction first introduced by Shiraishi and Mori 20 .
Our target states |ψ i are assumed to be non-thermal, hence we can furthermore assume there exists a set of local projectors, P i , which annihilate these states, P i |ψ j = 0, for any i ranging over lattice sites 1, 2, . . . , L. Next, consider a lattice Hamiltonian of the form where h i are arbitrary operators that have support on a finite number of sites around i, and [H , P i ] = 0 for all i. It follows thus [H, i P i ] = 0. Therefore, H takes the desired block diagonal form in Eq. (15). The target eigenstates can in principle be embedded at arbitrarily high energies for a suitably chosen H , which also ensures that the model is overall non-integrable 20 . However, there is no guarantee that the embedded states must be equidistant in energy and they may even be degenerate, such that this scheme could clearly result in a model without an SGA property discussed in Sec. III A. Physical applications of projector embedding include various topologically ordered systems 83,84 and models with lattice supersymmetry 85 , which are themselves defined in terms of local projectors (see Box 5). While the projector construction allows to embed general classes of states into a thermalising spectrum, in practice one is often more interested in the reverse question: given a Hamiltonian H belonging to some known physical system, can one identify non-thermal embedded states |ψ i ? Although this task is obviously much harder, for some physical models such as the PXP model, introduced in Sec. IV below, this question has been answered affirmatively for a single target state at zero energy in the middle of the spectrum, which has been identified with the AKLT ground state in Eq. (5) 86 . However, for the AKLT non-thermal eigenstates constructed via the SGA in Sec. III A, in particular for the tower of states in Eqs. (17)-(18), at present it is not known whether it is possible to find a set of local projectors and cast the entire tower of AKLT states in the Shiraishi-Mori form. This suggests there may be broader classes of models with non-thermal eigenstates that extend beyond the form in Eq. (31). Indeed, in a very recent work 87 a scheme based on a double copy of the system has been proposed for embedding eigenstates with high entanglement.

Box 5 | Physical realisations of projector embedding
To illustrate constructions of non-thermal eigenstates via projector embedding, we consider a class of topological lattice models defined in terms of local projectors, following Ref. 83. A typical Hamiltonian of such models has the form where p labels, e.g., elementary plaquettes of a lattice, such as in toric code 88 . The operators Q p (β) are Hermitian, positivesemidefinite, and local. They contain only sums of products of operators defined within the bounded region labeled by p, e.g., at the Rokhsar-Kivelson point of the quantum dimer model on the square lattice 89 , Q p (β) are projectors that encode both the potential and kinetic (plaquette flip) terms. The dimensionless parameter β is used to deform solvable models and break integrability.
The operators Q p (β) are built so as to share a common null state |Ψ(β) , i.e., for every p we have Q p (β) |Ψ(β) = 0. If all the couplings α p are positive, the state |Ψ(β) is the ground state of H(β). If, instead, α p take both positive and negative values, then it is not guaranteed that |Ψ(β) is a ground state. Nevertheless, |Ψ(β) is still an eigenstate with energy E=0. When this state is a high-energy eigenstate of H(β), it is an atypical state in that it displays area law entanglement entropy, since it is also a ground state of a different local Hamiltonian: Hence |Ψ(β) is a non-thermal state, if H(β) is nonintegrable. By deforming exactly solvable models -the toric code, for instance -one can break integrability while retaining the E=0 state 83 . The construction outlined here can be shown to be equivalent to the embedding construction in Eq. (31) by diagonalising the operators Q p (β) and separating the terms belonging to zero eigenvalues (which must exist by construction) from other eigenvalues.

IV. PXP MODEL
In the previous section, we introduced several mechanisms of weak ergodicity breaking, realised by diverse physical systems. While these mechanisms are mutually independent and different systems may only exhibit some of them, we next show how these mechanisms "co-exist" within a paradigmatic model of weak ergodicity breaking known as the PXP model. This model arises as an effective description of Rydberg atoms in the regime of a strong Rydberg blockade, and its experimental realisation will be the topic of Sec. VII. As we explain below, the PXP model realises much of the weak ergodicity breaking phenomenology discussed previously, in particular its Hilbert space is fragmented and it contains non-thermal eigenstates which form an approximate SGA.

A. The model
The PXP model describes a chain of coupled two-level systems, where each system can be in one of the two possible states, |• and |• . In the Rydberg atom realisation, these states correspond to an atom being in its ground state or in the excited Rydberg state, respectively. However, for present purposes, we can view these two states as ↓, ↑ projections of a spin-1/2 degree on the given site. We consider a 1D chain of such two-level systems coupled according to the "PXP" Hamiltonian 90 , where is the projector on the ground state at site i. Equivalently, the projectors can be defined as Without projectors, the Hamiltonian in Eq. (35) is that of a free paramagnet: each spin would independently precess with the Rabi frequency set to 1. P introduces a kinetic constraint: it allows a spin to flip only if both of its nearest neighbours are in • state. For example, the process · · · • • • · · · ↔ · · · • • • · · · is allowed, while · · · • • • · · · ↔ · · · • • • · · · is forbidden. This makes the model intrinsically interacting, as it is no longer possible to describe the state of each spin independently of other spins. Numerical simulations based on exact diagonalisation of the PXP model with up to L=32 spins have demonstrated that the statistics of its energy-level spacings approached the prediction of random matrix theory 13 as the system size was increased. Hence, despite its very simple form, we expect the PXP model cannot be fully "solved" using the known integrability techniques 31 .
Another consequence of projectors is the fragmentation of the PXP Hamiltonian: H PXP splits into sectors corresponding to different numbers of adjacent spin excitations. The largest connected component of the Hilbert space is one that excludes any configurations with adjacent excitations, . . . •• . . .. The number of classical configurations that satisfy such a constraint is still exponentially large -more precisely, it scales asymptotically where ϕ is the golden ratio. This unusual scaling is a manifestation of the global constraint. Similar type of constraints can arise due to emergent gauge fields and have been used to model interactions between anyon excitations in topological phases of matter 91 . Apart from this largest sector, there are further sectors which contain some number of nearest neighbour spin flips; however, in the remainder of this section we will focus on the largest connected component of the Hilbert space, which already displays non-trivial weak ergodicity breaking phenomenology.

B. Ergodicity breaking in the PXP model
While the PXP model is quantum-chaotic, numerical simulations of its quench dynamics 13,92-94 have revealed surprising nonergodic behaviour. For example, the return probability in global quenches with the PXP Hamiltonian is shown in Fig. 4(a) 13,94 . 'Global quench' means that the system is prepared in a highly non-equilibrium initial state, |ψ 0 , at time t=0, and the system is subsequently evolved with the many-body Hamiltonian, H PXP in Eq. (35). Since the PXP Hamiltonian is purely off-diagonal in the standard z-basis, any classical product state has an average energy equal to zero and extensive energy variance, thus effectively playing the role of an "infinite temperature" ensemble. In Fig. 4(a), three choices of density-wave states were considered for |ψ 0 : |Z 2 ≡ |· · · •••• · · · , |Z 3 ≡ |· · · •••••• · · · and |Z 4 ≡ |· · · •••••••• · · · . We note that these states can be prepared in Rydberg atom experiments by modulating the so-called detuning term introduced below in Eq. (45).
As the system evolves following the quench, one can characterise its behaviour by the return probability, also known as many-body fidelity, which quantifies the probability of observing the initial state after unitary dynamics, Intuitively, thermalising dynamics leads to a quick spreading of the many-body wave function over the full Hilbert space. According to Fig. 1(b), fidelity is expected to rapidly decrease to an exponentially small value, F ∝ 1/D L , and remain near that value at late times. The numerical result in Fig. 4(a) defies this expectation for certain initial states. In particular, the Z 2 and Z 3 density wave states show pronounced revivals at certain times, F (nT ) ∼ O(1). In contrast, quench dynamics from Z 4 (and many other initial product states not shown) displays fast relaxation without revivals. We note that the dynamics of local  observables, such as the density of domain walls nucleated in Z 2 , i.e., the number of · · · •• · · · patterns, is very similar to the fidelity dynamics. In particular, the frequency of revivals in local observables is the same as that of fidelity for the initial state |Z 2 . Below we focus on understanding the dynamics for this initial state.
In order to understand the origin of the atypical dynamical behaviour in the PXP model, we write the return amplitude as where E n , |E n denote the eigenenergies and eigenvectors of H PXP , respectively. Thus, quantum evolution is fully determined by the eigenenergies, {E n }, and overlaps of energy eigenstates with the initial state, {| E n |Z 2 | 2 }. Any special features in the dynamics translate into atypical properties of these overlaps. Indeed, Fig. 4(b) reveals a set of eigenstates that have strongly enhanced overlaps with Z 2 state, | E|Z 2 | 2 . These eigenstates violate the ETH 3,4,30 , since we expect that individual eigenstates of thermalising systems should behave like thermal ensembles. By contrast, Fig. 4(b) demonstrates that special eigenstates have strongly enhanced overlaps with a particular product state, thus they are anomalously concentrated in the Hilbert space and do not resemble random vectors. Note that violations of the ETH are traditionally probed by studying the matrix elements of local observables 95 . The overlaps with |Z 2 state indeed probe the matrix elements of the alternating magnetic field operator, j (−1) j σ z j , for which |Z 2 is the highest weight state. The higher the overlap with |Z 2 , the stronger the violation of the ETH -thus, the special eigenstates are the "most" ETH violating states in the spectrum of the PXP model. Intriguingly, the highest overlaps with |Z 2 are achieved in the middle of the spectrum, i.e., at the highest effective temperature in this system.
The total number of special states in Fig. 4(b) scales with the number of spins as L+1 94 . Moreover, the special states are equidistant in energy, especially near the middle of the spectrum. Along with their large overlap with |Z 2 state, this explains the existence of fidelity revivals. However, as seen from the coloured density of data points in Fig. 4(b), while the majority of other states in the spectrum have vanishingly small overlap with |Z 2 , there is still a considerable number of states forming "towers" which cluster around the energy of the special eigenstates. Towers of scarred eigenstates are also found in the PXP model in the presence of periodically driven detuning 96 , and in a two-dimensional PXP model 97 , where the overlap is computed with a charge-density wave state having a checkerboard pattern 98,99 .
Furthermore, the same L+1 special eigenstates show anomalously low entanglement entropy as seen in Fig. 4(c). For a thermalising eigenstate, the bipartite entanglement entropy of a subsystem of length L A is expected to scale extensively with the volume of the subsystem, S ∝ L A . In contrast, the entanglement entropy of special eigenstates is found to scale approximately with the logarithm of the subsystem size, as in Eq. (19) 94 . This is another evidence of ETH violation as it shows that special eigenstates are not uniformly spread over the Hilbert space.
Recall that logarithmic scaling of entanglement entropy is commonly found in other models with SGA, Eq. (19). Indeed, approximations to special eigenstates can be constructed starting from an approximate ground state of the PXP model and creating spin wave excitations on top of it 100 . Although the accuracy of this scheme deteriorates for special eigenstates near the middle of the spectrum, the scheme provides an intuitive picture of special eigenstates as condensates of weakly-interacting magnons, similar to models with an SGA. An alternative approach in Ref. 101 focused on the middle of the energy spectrum of the PXP model, where a few exact eigenstates can be constructed, thus rigorously demonstrating the ETH violation (see Box 6). While exact MPS constructions similar to Ref. 101 can be generalised to other models, e.g., a transverse Ising ladder 102 and the Floquet versions of the PXP model 103,104 , they remain limited to a small number of states in the middle of the spectrum, and in particular they do not capture the top states in the towers, which are not simple MPS due to having larger-than-area-law entanglement.
Finally, we note that the model in Eq. (35) with the addition of a few diagonal terms can be made integrable 91,105,106 and it has a very rich phase diagram in its the ground state. However, these terms that render the PXP model integrable have magnitude O(1), i.e., they are strong perturbations and are not expected to be relevant for explaining the physics presented in this section. Furthermore, several studies explored the effect of (generally, weaker) perturbations on the PXP model 94,107,108 . While nonthermal eigenstates are generally not robust to perturbations, some aspects of atypical dynamical behaviour are found to persist, including in the presence of disorder 109 . Finally, we mention that similar phenomenology of weak ergodicity breaking has also been found in the two-dimensional PXP model 98,99 , and in related models of transverse Ising ladders 102 and the periodically driven PXP model 96,103,104 .

Box 6 | Ergodicity breaking in null spaces
Among the special properties of the PXP model is the existence of an exponentially large number of states which are exactly annihilated by the PXP Hamiltonian in Eq. (35). This null space contains states with energy E=0 that reside in the middle of the spectrum 13,94,110 . This exponential zero-energy degeneracy is a consequence of the combined action of particle-hole transformation C = i σ z i , which anticommutes with the PXP Hamiltonian, and spatial inversion which maps site j to L−j+1 111,112 . Among the exponentially many E=0 eigenstates in the PXP model, there are also a few weakly entangled states which have been analytically constructed by Lin and Motrunich 101 . Such states are compactly written as MPS in Eq. (3) with χ=2: The local Hilbert space on which these matrices are defined is a 2-site block of the chain; due to the constraint, this block can be in three states, ( In fact, from the state defined in Eq. (39) we can construct its partner translated by one site, also at E= 0. These two states thus display translational symmetry breaking with a period-2 bond-centered pattern, despite being in 1D and at an infinite temperature. They also manifestly violate the ETH via their entanglement entropy, which must obey the area law due to their explicit MPS form. We note that Ref. 113 has recently shown that MPS states such as the one in Eq. (39) can arise more generally in random quantum networks due to connectivity bottlenecks.
In Ref. 114, a more general exploration of null spaces was performed by applying a "subspace disentangling" algorithm in order to construct the least entangled state belonging to E=0 subspace. It was numerically demonstrated on several examples that this least entangled null state obeys area-law entanglement scaling, and as such, violates strong ETH. Thus, although protected by symmetry, null spaces provide another non-trivial mechanism for realising non-thermal states.

C. The origin of non-thermal eigenstates and quantum revivals
In quantum physics, the simplest system that exhibits non-trivial dynamics and revivals is an elementary spin: a spin pointing along the z-direction precesses when the magnetic field is turned on in the x-direction, causing the spin to periodically return to its initial orientation as time passes. The revivals in the PXP model can also be understood as precession of a spin with magnitude s=L/2. The latter spin, however, is a collective degree of freedom representing the system of L atoms 115 . More precisely, the L+1 scarred eigenstates identified in Fig. 4 form an approximate representation of an su(2) algebra for spin-L/2. This perspective brings the PXP model in line with other models discussed in Sec. III A, with the main difference being that the SGA in the PXP model is only approximate.
Specifically, the spin picture follows from the mapping of the PXP model to a tight-binding chain -see Box 7. For the initial |Z 2 state, the "big spin" raising operator excites an atom anywhere on the even sublattice and deexcites an atom on the odd sublattice, whereσ ± i = P i−1 σ ± i P i+1 are the on-site raising and lowering operators that respect the constraint 13 . Similarly, the spin lowering operator H − performs the same process with the sublattices exchanged.
The reason for this choice of H ± is that their commutator defines the z-projection of spin, H z ≡ 1 2 [H + , H − ], for which |Z 2 plays the role of the extremal weight state. Now the analogy with spin precession is almost complete because the PXP Hamiltonian in Eq. (35) is given by the sum H PXP = H + + H − , i.e., it plays the role of an x-component of spin. Thus, preparing the atoms in |Z 2 state is equivalent to initialising the spin along the z-axis, and the state revives because the PXP Hamiltonian acts like a transverse magnetic field.
If the above spin picture were exact, the revivals in the PXP model would be perfect, with fidelity in Eq. (37) reaching 1 at certain late times. This is not seen, either in experiments discussed below or in the numerical simulations of the PXP modelrecall Fig. 4(a). The reason is that the mentioned su(2) spin algebra is only approximate, [H z , H ± ] ≈ ±H ±115,116 , where " ≈ " means there are additional operators on the right-hand side with smaller numerical prefactors from the leading H ± .
It has been realised that the structure of the su(2) algebra and the robustness of the revivals can be significantly improved by small deformations of the PXP model 115,117 . The inclusion of deformations completely arrests the entanglement growth, resulting in the band of scarred eigenstates nearly fully separated from the rest of the spectrum. We note that a similar procedure can be used to enhance revivals from other initial states, such as |Z 3 and |Z 4 , for suitably redefined H ± operators 116 . Thus, the PXP model could be deformed to stabilise different embedded su(2) algebras. The deformations of the PXP model that enhance the su(2) algebra and improve the |Z 2 revival fidelity suggest the existence of an idealised parent model that hosts "perfect" many-body scars. Unfortunately, the parent model is defined by a long-ranged Hamiltonian with exponentially-decaying tails 115 , and it remains unknown whether it can be expressed in a more compact (short-range) form.

Box 7 | Forward scattering approximation
A free paramagnet is an example of a perfectly reviving system: any product state of spins pointing along the z-axis exhibits periodic dynamics under the magnetic field along the x-direction. Alternatively, the free paramagnet can be viewed as a hypercube graph of dimension L, where each of the 2 L spin product states represents a vertex, while the edges connect vertices that can be reached by flipping one spin. This picture helps to understand the PXP model, which is a partial cube, i.e., a hypercube where all the vertices violating the constraint have been removed -see Fig. B7. The resulting PXP graph contains two smaller hypercubes of dimension L/2, with |Z 2 and |Z 2 states as extremal vertices, and the polarised state |•• . . . at the intersection point of the two hypercubes. Additionally, there are also "bridges" that connect the two hypercubes, e.g., in Fig  Quantum walk on the PXP graph starting in the |Z 2 vertex can be accurately modelled by assuming the wave function spreads only in the forward direction, a scheme known as the Forward Scattering Approximation (FSA) 13 . The FSA compresses the partial cube down to a one-dimensional chain, where each site |n is a superposition of states with the same number of excitations relative to |Z 2 vertex. Thus, for L atoms, the chain contains L+1 sites. The FSA interprets each site |n as an eigenstate of a big spin with magnitude s=L/2 and pointing at an angle πn/L with respect to z-axis. As mentioned in the text, the spin-raising and lowering operators directly follow from a decomposition of the PXP Hamiltonian 13 . The FSA not only allows to construct remarkably accurate approximations to tthe exact eigenstates of the PXP and several other models 13,73,94,118 , but it also provides the foundation for the elegant interpretation of the revivals as precession of an emergent spin.

V. SEMICLASSICAL DYNAMICS
Quantum-classical correspondence has been at the heart of quantum mechanics since the early 20th century. Parallels between the behaviour of a quantum system and its classical counterpart have been extensively studied in the context of few-body systems, such as the free propagating particle in a bounded domain (stadium billiard being a particular example), coupled harmonic oscillators, and a rotor subject to periodic driving 52 . In all these cases, the semiclassical limit can be taken by sending → 0, and the behaviour of the quantum system is often found to closely parallel its classical counterpart. For instance, quantum counterparts of classically chaotic systems can be typically described by random matrix theory 119 .
Intrinsically many-body systems, such as the PXP chain in Eq. (35), present challenges for the semiclassical description based on the conventional large-N limit or the mean-field approximation. The latter, for example, ignores any correlations between different spatial degrees of freedom, and thus violates the constraint in the PXP model in Eq. (35). The constraint introduces non-local correlations: if a given spin is in |• state, its neighbours must be in |• state. Such correlations build in quantum entanglement, which cannot be ignored unless one relaxes the constraint.
The MPS methods discussed in Sec. II allow to incorporate entanglement in a systematic way that goes beyond mean field 49 . For the PXP model, it has been first shown by Ho et al. 14 that this approach can be tailored to capture dynamics from the initial |Z 2 state. This state is a period-2 density wave, which can be described using a 2-atom unit cell. Hence, the variational MPS state is parametrised by two sets of angles, x i =(θ i , φ i ), i=1, 2, that describe the state of spins on odd and even sites in the lattice. The MPS matrices can be taken to be 97 where we notice that A • ∝σ + satisfies the condition A • (θ, φ)A • (θ , φ ) = 0, which effectively imposes the constraint that no two adjacent spins are in "•" state. The special points (θ 1 , θ 2 ) = (0, π/2) and (θ 1 , θ 2 ) = (π/2, 0) correspond to |Z 2 = |•••• . . . θ 1 = tan θ 2 sin θ 1 cos 2 θ 1 cos φ 1 + cos θ 2 cos φ 2 , (42a) which were derived for the PXP model in the presence of a chemical potential, µ z j n j , where n j ≡ 1 2 (1 + σ z j ). The equations for θ 2 , φ 2 can be obtained by substitution 1 ↔ 2.
Before discussing the solutions of the equations of motion (42), we note that the expectation value of the Hamiltonian is a conserved quantity, provided the the equations of motion are satisfied. This conservation law effectively reduces the dimensionality of the phase space from four down to three dimensions, as it restricts the dynamics to constant energy surfaces. However, when µ z =0 the system in Eqs. (42) actually has a class of solutions with two angles being stationary, φ 1,2 =0 14 . This class of solutions corresponds to a flow-invariant subspace, here arising due to particle-hole symmetry and time-reversal invariance of the PXP Hamiltonian in Eq. (35) 97 .
Restricting to the flow-invariant subspace φ i =0 (with µ z =0), we plot the flow diagram in θ 1 -θ 2 plane in Fig. 5(a). Shown in red colour is the periodic trajectory first identified in Ref. 14. This trajectory is intimately linked with quantum revivals for the initial Neél state |Z 2 in Fig. 4(a), since this state is represented by the point (θ 1 , θ 2 )=(π/2, 0) on the trajectory. We note that the variational description in Fig. 5(a) is not accurate at long times: the fidelity in the full quantum evolution in Fig. 4(a) shows a visible decay, which is in contrast with the variational description that implies perfect oscillations. As we mentioned in Sec. II B, this discrepancy between the exact quantum evolution and its projection onto the variational manifold can be quantified via quantum leakage in Eq. (14). The leakage measures the instantaneous disagreement between exact dynamics and its projection onto the variational manifold. Remarkably, the periodic trajectory is located entirely in the low leakage region 14 , thus providing a posteriori justification for the coherent quantum dynamics when the system is prepared in |Z 2 state. Two-dimensional phase space, however, is special as it precludes chaos. Dynamical systems with phase spaces higher than two-dimensional are known to display both classical chaos as well as regions of stable motion -the so-called Kolmogorov-Arnold-Moser (KAM) tori 120 . These features are prominent in Fig. 5(b) which shows the Poincaré section of the dynamical system in Eq. (42) when the chemical potential is present (µ z =0.325) and the variables φ i cannot be fixed to zero. To obtain such a graph, we follow the standard recipe: we choose some value for θ 1 and follow the successive intersections of the flow with the chosen section. This produces a discrete mapping known as the Poincaré map -this maps a given point (θ * 1 , φ 1 , θ 2 , φ 2 ) on the chosen hyperplane to a position (θ * 1 , φ 1 , θ 2 , φ 2 ) where the trajectory intersects this plane again. Periodic trajectories correspond to stationary points of the Poincaré map, while a chaotic system returns to the same plane at a location that is generally far away from the previous encounter. This procedure for evaluating the Poincaré sections can be conveniently implemented, e.g., in Mathematica using the built-in NDSolve routine. Fig. 5(b) illustrates a typical example of mixed phase space 52 : there exist (in general, multiple) KAM tori with stable periodic orbits in their centers, surrounded by the chaotic sea. Energy conservation results in complicated surfaces in the space (φ 1 , θ 2 , φ 2 ) that cannot be globally projected. Therefore, the figure shows a small part of the Poincaré section in the vicinity of the stationary point, (φ 1 , θ 2 ) ≈ (2.92, 2.98), corresponding to a stable periodic trajectory 97 . This periodic trajectory is surrounded by circle-shaped contours which are the intersections of KAM tori by our section plane. Unlike the pure PXP model (µ z =0) where the periodic orbit passes through a product state (|Z 2 ), in this case the states on the trajectory are entangled in general (though entanglement is necessarily small, since the orbit passes through MPS states of low bond dimension).

A. Discussion: benefits and pitfalls of TDVP
To summarise, the TDVP discussed above allows to obtain classical equations of motion and analyse parallels between the dynamics of a quantum system and its projection onto the variational manifold. The presence of (stable or unstable) classical periodic trajectories with small quantum leakage leads to slow relaxation in the full many-body dynamics.
There are several useful insights we glean from TDVP. (i) The KAM tori inform us about special features of quantum dynamics in the underlying model. For example, when dynamics is initialised precisely on the periodic orbit, there is slower relaxation of local observables and slower entanglement growth in quench dynamics 97 . Furthermore, initialising the dynamics in the centre of a big KAM torus results in slower entanglement growth. This highlights the role of quantum leakage in determining which of the trajectories give rise to slow thermalisation in exact quantum dynamics. (ii) Revivals in the PXP model are robust to various deformations of the Hamiltonian, provided that in the course of the deformation one follows the periodic trajectory that is influenced by the deformation 97 . (iii) The TDVP approach has also been fruitful in obtaining a more complete understanding of the PXP dynamics in manifolds with more degrees of freedom. In particular, this approach can describe the dynamics from the period-3 density wave state, |Z 3 , which also shows revivals of fidelity, as seen in Fig. 4(a). The corresponding variational ansatz for this state includes three angles θ i , i=1, 2, 3, which parametrize the state of three atoms within a unit cell (the corresponding φ i once again can be fixed to zero in the flow-invariant subspace). While the equations of motion are considerably more complicated than ones in Eq. (42), their numerical solution also generally gives rise to mixed phase space, similar to Fig. 5(b).
At the same time, one must beware some caveats to the TDVP approach. First, the variational manifold is not unique -one can always enlarge the MPS ansatz by including additional parameters that build in longer-range entanglement. This leads to potential ambiguities in the definition of the trajectory. While some TDVP ansätze may find unstable periodic trajectories, more generic embeddings may result in mixed phase space. Second, at present there is no understanding under which conditions the quantum leakage is small. Hence, even if one recovers a large number of periodic trajectories, they may all have strong leakage and not lead to any revivals in quantum dynamics, e.g. as in the chaotic Ising model in Ref. 97.
Finally, many questions remain about quantum-classical correspondence in TDVP as a function of bond dimension χ. Complementing the study of weak ergodicity breaking at small χ, TDVP has been used to study thermalisation at larger χ. In particular, Ref. 121 demonstrated that physical properties saturate in the case of transverse field Ising model when bond dimension χ>2. In addition, Ref. 122 studied the spectrum of Lyapunov exponents in the case of high bond dimension TDVP. The work in Ref. 97 shows that in the case of low bond dimensions one may encounter non-chaotic behaviour 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 123 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.

VI. QUANTUM MANY-BODY SCARS
In Sec. III we discussed different mechanisms that lead to the emergence of non-thermal eigenstates in highly-excited energy spectra of non-integrable systems. In the past few years, the name quantum many-body scars has been widely adopted as an umbrella term for such non-thermal eigenstates, regardless of the underlying mechanism of ergodicity breaking. In this section, we explain in more detail the reasons why this name was adopted by Turner et al. 13 . The inspiration was the PXP model and its many parallels with a single particle confined to a stadium-shaped billiard. In the latter case, the particle's eigenfuctions can concentrate around certain periodic orbits in the limit → 0, a phenomenon Heller 15 called quantum scarring in the early 1980s. In the remainder of this section, we discuss the semiclassical aspects of the PXP physics that make it the many-body analogue of quantum scarring.

A. Scars in few-body systems
What can we learn about the behaviour of a quantum system by looking at its classical counterpart? The Bohr-Sommerfeld quantisation demonstrates that classical integrable systems, such as the harmonic oscillator or hydrogen atom, have rather special quantum spectra and eigenstates. For classically chaotic systems, e.g., the Bunimovich stadium, traditional quantisation methods do not work. Typically, quantum counterparts of classically chaotic systems display level repulsion in their spectrum and their eigenstates look random. However, short unstable periodic orbits may leave a strong imprint on the system's quantum dynamics and eigenstate properties -this influence of classical periodic orbits is known as "quantum scarring" 15 : Definition (from Ref. 124): A quantum eigenstate of a classically chaotic system has a 'scar' of a periodic orbit if its density on the classical invariant manifolds near the periodic orbit is enhanced over the statistically expected density.
Utilising this definition, one can detect scars by visualising the wave function probability density, which can be seen to concentrate around certain periodic classical orbits, e.g., the diamond-shaped "bow tie" orbit in the Bunimovich stadium 124 . Moreover, scars leave an imprint on the dynamics: when a wave packet is launched in the vicinity of an unstable periodic orbit, it will tend to cluster around the orbit at later times, displaying a larger return probability than a wave packet launched elsewhere in the phase space. Furthermore, such a wave packet can be expanded over a small number of eigenstates that have approximately similar energy spacing, in contrast to an arbitrary wave packet. Much of this phenomenology was originally acquired through numerical solutions of quantum billliard problems, since the latter are non-integrable and do not admit analytical solution.
Nevertheless, rigorous proofs [125][126][127] have established that individual eigenstates of the billiard are "almost always" ergodic, since the phase space area affected by scarring vanishes in the semiclassical limit → 0, shrinking around the periodic orbit. These proofs of scarring have been based on constructions of approximate eigenstates or "quasimodes", i.e., special states designed to be strongly localised around a classical periodic orbit. For example, in the Bunimovich stadium, some of the quasimodes have the form ψ(x, y) = φ(x) sin(nπy), i.e., they represent a standing wave in one direction with a suitably chosen envelope function in the other direction. By carefully controlling the density of states, it was possible to show that there exist eigenstates with an anomalously high overlap with a small number of quasimodes, hence "inheriting" their scarring properties 124 .
Although the above suggests that scars are "rare" anomalies, they do have physical signifance. For example, scars provide an important counterexample to the intuitive expectation that every eigenstate of a classically chaotic system should locally look like a random superposition of plane waves 128 . Furthermore, also counterintuitively, a scarred quantum system appears more "regular" than its classical counterpart, since in the latter case there is no enhancement of density along the periodic orbit in the long-time limit. Finally, scars play a role in many experiments, including microwave cavities 129 , semiconductor quantum wells 130 , and the hydrogen atom in a magnetic field 131 .

B. Quasimodes in the PXP model
In few body systems, the understanding of quantum scars rests on two pillars: the existence of unstable classical periodic orbits and the quasimodes which leave a scar upon some of the eigenfunctions. The main difficulty in passing to the many-body case is how to give precise meaning to this terminology: the physics now takes place in the many-particle Hilbert space, precluding "easy" visualisation of the classical trajectory. In Secs. II B and V we introduced one way of defining the classical trajectory via TDVP, wherein many-body dynamics is projected onto a restricted space of MPS states. This approach systematically builds in entanglement, controlled by the MPS bond dimension, thus it forms a natural framework to define the semiclassical limit of the PXP model which does have an obvious → 0 limit.
The application of the TDVP method to the PXP model successfully captures the revivals following quenches from specific product states, hence it is natural to suggest that this approach defines an effective "semiclassical" description of the quantum dynamics, and holds the same relation with the exact PXP model as expected from the Bohr correspondence principle. However, in order to complete the analogy with single-particle scars, it is necessary to show that this limit also gives rise to quasimodes which, in turn, accurately approximate the quantum eigenstates of the PXP model.
In Ref.132 this correspondence was achieved by constructing a subspace K, which is fully symmetric over permutations within each of the two sublattices, comprising even and odd sites in the PXP chain -see Box 8. This forms a mean-field approximation for the PXP model wherein one describes the dynamics in terms of numbers of excitations on even or odd sites, while discarding any detailed local information (such as the precise positions of excitations within a sublattice). This approximation reduces the exponential complexity of the many-body problem down to polynomial in system size, and the quasimodes are obtained by simply diagonalising the Hamiltonian projected into the symmetric subspace, KH PXP K. As can be seen in Fig. 6(a), these quasimodes are excellent approximations to the exact eigenstates of the PXP model, in particular they capture the L+1 nonthermal eigenstates with the largest projection on |Z 2 state in Fig. 4(a). It turns out that the symmetric subspace K has an intimate relationship with the subspace M of MPS states discussed in Sec. V, which can be viewed as a Gutzwiller projection of spin coherent states with a unit cell of two sites 14 . The states in M are parametrised by two angles which represent probabilities for a site of each sublattice to hold an excitation, while K consists of states with definite occupations numbers for each sublattice. In this sense, there is a direct relationship between the two subspaces, reminiscent of the relationship between canonical and grand canonical ensembles in statistical mechanics. More formally, it can be proven that the linear span of M is equal to the symmetric subspace K, for every fixed system size L 132 . This implies that the dynamics within K necessarily includes all quantum fluctuations on top of the equations of motion in Eqs. (42).
The correspondence between the classical system M and the quantum system in its linear span (K) finallly allows us to view the quasimodes |ψ ∈ K as wavefunctions ψ(θ, φ): Here Ψ(θ, φ) ∈ M is the MPS state parametrised by Eq. (41), and S µ is the frame operator 132 . This operator is needed because we need to equip the classical state space with an integration measure µ in order to have a well-defined quantisation. In the more familiar manifolds of quantum states, such as spin coherent states 133 or unconstrained MPS 134 , the measure µ is simply given by the Haar measure for a transitive group action. However, for the subspace K with angles θ i , φ i having 2-site periodicity, the measure is non-trivial and one must use the frame operator.
Using the above correspondence Eq. (43), the real part of the quasimode wavefunctions is plotted in Fig. 6(b)-(c). For one of the top-band quasimodes near the middle of the spectrum [denoted by a cross in Fig. 6(a)], we observe enhanced concentration around the |Z 2 classical trajectory [shown by the black line in Fig. 6(b)]. This behaviour is strikingly reminiscent of wavefunction scarring in quantum billiards. In particular, following the classical trajectory, the phase of the wavefunction winds, and integrality of this winding number leads to the quasimodes' approximately equal spacing in energy. We recognise this behaviour from quantisation of regular trajectories 52 : in the old quantum theory, this underpins Sommerfeld-Wilson quantisation and leads to the de Broglie standing-wave condition for the Bohr model. These features are inherited by the exact scarred eigenstates given their strong overlap with the quasimodes. In contrast, Fig. 6(c) shows one of quasimodes deeper in the bulk of the spectrum. Such quasimodes are typically found around the orbits connecting corners of the manifold, where quantum leakage is large (see Sec. V). It is expected that these quasimodes would not survive the "injection" into the full Hilbert space, i.e., while the top-band quasimodes accurately approximate the true eigenstates of the full PXP model, this is not the case for the other quasimodes.

Box 8 | Mean field approximation for the PXP model
The essence of mean-field theory is the erasure of unnecessary local information. For example, to describe a free paramagnet, which is obtained by dropping the projectors in the PXP Hamiltonian in Eq. (35), we only need to know the total number of excitations relative to some reference state, not their precise positions in the lattice. This idea has been generalised to describe the reviving dynamics in the PXP model where the important information is the number of excitations on each of the two sublattices, comprising even and odd sites in the chain 132 . Such an approach is reminiscent of the symmetric subspaces in studies of fullyconnected models 135,136 . However, the key difference is that in the PXP model the permutation symmetry is broken to the sublattice level, and many of the permutation shuffles violate the constraint and therefore have to be excluded, which makes the analysis more involved.
To illustrate the approach, the symmetric subspace K is defined by forming a set of equivalence classes (n 1 , n 2 ), where integers n 1 , n 2 label the number of excitations on the two sublattices, encompassing the odd and even sites, respectively. Elements in these classes are equivalent under the action of the product of two symmetric groups S L/2 which "shuffle" the sites in each sublattice.  Fig. B8. The orthonormal basis for K is built from symmetric combinations of members of each class, |(n 1 , n 2 ) ∝ x∈(n1,n2) |x , where the sum runs over all product states |x that are members of the class (n 1 , n 2 ). Class sizes and matrix elements of the PXP Hamiltonian projected to K can be calculated using combinatorics 132 , allowing for a highly-efficient description of the PXP model and its dynamical properties. The results presented in this section and Sec. V provide justification for referring to the non-thermal eigenstates in the PXP model as a many-body analogue of single particle scars. One of the crucial ingredients that allows to make such a connection is the quantum-classical correspondence. At the same time, in Sec. III we introduced a number of other models that also feature non-thermal eigenstates, quite reminiscent of those in the PXP model. Here we briefly discuss the main similarities and differences of these models.
The PXP model, with its generalisations to higher spins 14 and clock models 118 , is broadly similar to models whose nonthermal eigenstates form an SGA. The main difference is that the SGA in the PXP model is only approximate, and the subspace has weak residual couplings to the thermal bulk. In addition, the PXP model also realises an approximate Krylov subspace built upon |Z 2 root state. The Krylov subspace becomes exact when generated by H + in Eq. (40). Finally, a single PXP eigenstate at zero energy in the middle of the spectrum can be viewed as projector-embedded AKLT ground state 86 . While similar connections between scarring mechanisms may be anticipated in other models mentioned in Sec. III, their demonstration remains an open problem.
In particular, the existence of classical periodic trajectories underlying quantum revivals in the PXP model calls for exploration of semiclassical dynamics in other scarred models discussed in Sec. III. For example, in the spin-1 XY model 53 in Eq. (20), the initial state with perfect revivals was identified to be This state is a superposition of the scar tower states |S XY n in Eq. (21) and it can be prepared as the ground state of the Hamiltonian H = Q † +Q ∝ J x . Since these states are equally spaced in energy with spacing 2h [see Eq. (20)], time evolving |ψ 0 gives perfect revivals with frequency 2h. For the AKLT tower of scars, such "rotations" in the space of scar states do not seem to produce a very simple initial state. However, as shown in Ref. 58 the AKLT scar tower can be compressed into a state with finite MPS bond dimension (e.g., bond dimension equal to 8 in a system with periodic boundary conditions).
We expect that similar constructions of reviving initial states can be performed for other models with an SGA using the properties of Q † . By contrast, for non-thermal eigenstates produced via the projector-based embedding scheme, we generally do not expect underlying periodic trajectories. Similarly, models with Krylov fragmentation (Sec. III B) are not guaranteed to have revivals, even if the root state |ψ 0 is experimentally preparable. This is because general tridiagonal matrices do not support revivals, unless their matrix elements are tuned to special values 137 . Thus, the study of classical periodic trajectories in a broader family of models could be helpful as a finer classification scheme for models displaying weak ETH violation in their eigenstate properties, in particular as a way of distinguishing many-body scarring from more generic embeddings of non-thermal eigenstates.

VII. WEAK ERGODICITY BREAKING IN EXPERIMENT
The experiments that opened the door to the investigation of weak ergodicity breaking discussed above were performed using quantum simulators based on Rydberg atom arrays 12,138 . Here we provide a brief overview of these experiments, focusing on the work by Bernien et al. 12 which reported the first experimental observation of quantum many-body scarring. Following this, we discuss more recent experiments realising the tilted Fermi-Hubbard in optical lattices, where signatures of non-ergodicity due to Hilbert space fragmentation have been detected 21 .

A. Rydberg atoms
An individual Rydberg atom may be viewed as an effective two-level system, where the two states |• , |• correspond, respectively, to an atom in the ground state and an atom in the so-called Rydberg state, i.e., with a highly excited electron in the outer shell, see Fig. 7(a). The state of an individual Rydberg atom may be manipulated via coupling to circularlypolarised radiation which causes Rabi oscillations with frequency Ω. The relative detuning ∆ of the Rydberg state off-resonance can also be controlled, leading to an effective single atom Rydberg Hamiltonian H = (Ω/2)σ x + ∆n, where the operator σ x = |• •| + |• •| is the Pauli matrix which describes Rabi oscillations, andn = |• •| measures the population of the Rydberg state.
Atoms in Rydberg states interact via dipole-dipole interactions, thus the many-body Hamiltonian describing this system is given by where V i,j = 1/r 6 ij and r ij is the distance between the atoms. By tuning r ij one can achieve the regime of the Rydberg blockade, as illustrated in Fig. 7(a). In this regime, the shift of an energy level with two excited Rydberg atoms is so strong that the state |•• is off-resonant and cannot be reached from the ground state. Rydberg blockades with varying radii were demonstrated experimentally 138,139 via the frequency renormalization of the Rabi oscillations. Coherent dynamics was observed in small arrays of Rydberg atoms in both 1D and 2D when the blockade radius exceeded the linear size of the chain. More recent experiments 12 have studied Rydberg arrays in a qualitatively different regime with a relatively short radius of the Rydberg blockade. In particular, long Rydberg chains, containing up to 51 atoms with the nearest-neighbour blockade, were quenched from a period-2 density wave initial state, |Z 2 ≡ |•••• . . . , corresponding to the state with the maximal possible number of Rydberg excitations allowed by the blockade. The dynamics of the domain wall density, where domain walls are defined as adjacent |•• or |•• configurations (the latter is excluded in the regime of perfect blockade), revealed long-time oscillations, similar to Fig. 4(a).
The density revivals observed in Ref. 12 were challenging to explain theoretically, for the following reasons. First, the density wave state |Z 2 forms an infinite-temperature ensemble for the atoms in the Rydberg blockade, thus it should display fast equilibration. Indeed, for other initial states at infinite temperature, such as |0 ≡ |••• . . . , the domain wall density was found to relax quickly to the steady-state value. Second, the persistence of crystalline order far from equilibrium was surprising because the Rydberg atom system was not known to have any conserved quantities other than the total energy. Finally, the frequency of domain wall density oscillations did not coincide with √ 2Ω, which would be naïvely expected for collective oscillations of two atoms in the Rydberg blockade regime, signalling that many-body effects played an important role. As we have already seen, these puzzles are now understood based on the studies of the PXP model summarised in Sec. IV. The PXP model arises as the effective model for the Hamiltonian in Eq. (45) when V i,i is infinite and interactions beyond nearest-neighbour are set to zero 90 . An important step towards the explanation of the observed oscillations was made by identifying the non-thermalising eigenstates in the spectrum of the PXP model in Fig. 4(b). The core phenomenology of the PXP model -the small number of ETH-violating eigenstates within the thermalising spectrum and the presence of many-body revivals and slow relaxation in quenches from specific initial states -led to the identification of the phenomenon in Ref. 12 with a many-body analogue of quantum scarring 15 .
We note that since the experiment in Ref. 12 there have been further advances in similar Rydberg platforms. In Ref. 140 it was shown that coherent revivals associated with quantum many-body scars can be enhanced via periodic driving, which generates a robust subharmonic response akin to discrete time-crystalline order 141 . The driving is realised by modulation of the detuning term, ∆(t) = ∆ 0 + ∆ m cos(ω m t), and the enhancement is found in the non-perturbative regime of ∆ m , ∆ 0 , ω m ∼ Ω. It was found that such time-dependent detuning can lead to a five-fold increase of scar lifetime compared to the fixed-detuning case 140 . Notably, parametric driving not only delays thermalisation, but also alters the actual trajectory being stabilised, which can be used to effectively "steer" complex dynamics in many-body systems 142 . This opens the door to robust creation and control of complex entangled states in the exponentially large Hilbert spaces of many-body systems, with intriguing potential applications in areas such as quantum metrology. Other fundamental phenomena recently realised in the same Rydberg platform include dynamical phase transitions and Kibble-Zurek scaling 143 and topological spin liquids 144 (see also Ref. 145 for a realisation in superconducting qubits.)

B. Tilted optical lattices
Another platform for probing weak ergodicity breaking are ultracold atoms in tilted optical lattices, which realise the 1D Fermi-Hubbard model 21 : whereĉ † j,σ denotes the usual electron creation operator on site j with spin projection σ,n j,σ ≡ĉ † j,σĉ j,σ , J and U are the hopping and on-site interaction terms, respectively. Compared to the standard Hubbard model, Eq. (46) includes the tilt potential ∆, which has the form of a dipole term, ∼ jn j .
The experimental setup in Ref. 21 consists of a degenerate Fermi gas of 40 K atoms that is prepared in an equal mixture of two spin components. The atoms are loaded into a 3D optical lattice from which 1D chains are isolated along x-direction with the length of about 290 lattice sites. The on-site interaction strength U is controlled by a Feshbach resonance and a magnetic field gradient is used to create the tilt ∆ (approximately independent of spin). We note, however, that linear potentials similar to the tilt one in Eq. (46) can also be realised in other platforms such as trapped ions 146 and superconducting qubits 147 , where signatures of strong ergodicity breaking have been recently observed and attributed to "Stark many-body localisation" 76,77 .
The model in Eq. (46) is a natural setting for exploring Hilbert space fragmentation. Ignoring the spin degree of freedom, in the limit of large tilt, ∆ U, J, the leading-order off-diagonal term 72 is precisely the pair-hopping Hamiltonian in Eq. (24). With spin included, the effective Hamiltonian in the limit of large tilt (at third order) comprises an off-diagonal term 21 whose strength is proportional to J (3) = J 2 U/∆ 2 (σ denotes the opposite spin of σ). At the same order, the effective Hamiltonian also contains another off-diagonal term, ∝ 2J (3)T XY whereT XY = i,σĉ † i,σĉ i+1,σĉ † i+1,σĉ i,σ , as well as two diagonal terms, U 1 − 4J 2 ∆ 2 in i,↑ni,↓ and 2J (3) i,σn i,σni+1,σ . The effective model in Eq. (47) conserves the dipole moment, i,σ in i,σ . The fact that the hopping rate J (3) is proportional to the interaction strength highlights that interactions are necessary to generate dipole-conserving processes 21 . For initial states that are products of singlons, the connected dynamical sector S only represents a vanishing fraction of the whole (effective) symmetry sector, thus severely restricting the dynamics of the system. The dipole-conserving processes in general involve the generation of doublons. This is, however, penalised by the Fermi-Hubbard on-site interaction and therefore, we expect a slowing down of the dipole-conserving dynamics. The additional spin-exchangeT XY increases the connectivity, but cannot fully connect the whole dipole symmetry sector and the system remains fragmented.
In experiment, ergodicity breaking was probed by quenching the system from an initial state chosen to be a charge-density wave of singlons on even sites (filling factor ν=1/2), which is prepared using a bichromatic optical superlattice. The initial state can be described as an incoherent mixture of site-localised particles with random spin configuration. The subsequent evolution is monitored by extracting the spin-resolved imbalance I σ = (N σ e − N σ o )/N σ , where N σ e(o) denotes the total number of spin-σ atoms on even (odd) sites and N σ = N σ e + N σ o . A non-zero steady-state imbalance signals a memory of the initial state, where I σ (t=0) = 1. In the regime ∆ J, |U |, the effective Hamiltonian is dipole-conserving up to third order in J/∆, which was argued 21 to be responsible for the non-ergodic behaviour of imbalance I, which did not decay to its thermal value at moderate times. At much later times, higher-order processes, beyond the ones in Eq. (47), are expected to reconnect the different Krylov sectors.
Tilted Fermi-Hubbard model is a promising platform to investigate a range of ergodicity breaking beyond fragmentation. Although experimentally challenging due to finite evolution times, it would be interesting to reconcile the phenomenon of Stark MBL and Hilbert-space fragmentation, by studying the impact of weak disorder or residual harmonic confinement on the long-time dynamics. Moreover, it would be interesting to explore the connection between lattice gauge theories and Hilbertspace fragmentation, which could be addressed experimentally in a similar system 148 . Finally, recent work 149 has proposed that quantum many-body scars could also be realised in the same model, albeit at a different filling factor compared to Ref. 21 -see Box 9.
Box 9 | Proposal for scars in the Fermi-Hubbard model In the regime ∆≈U J, the tilted Fermi-Hubbard model in Eq. (46) was argued to host many-body scars at filling factor ν=1 149 . In this case the sum of the dipole moment and the number of doublons is effectively conserved, giving rise to the effective Hamiltonian H eff = −J j,σĉ † j,σĉ j+1,σnj,σ (1 −n j+1,σ ) + h.c., (48) where hopping to the left (which decreases the total dipole moment by 1) is only allowed if it increases the number of doublons by the same amount (σ is opossite spin from σ).
The action of the Hamiltonian in Eq. (48) within the ν=1 sector fragments the Hilbert space beyond the simple conservation of U +∆. The largest connected component of the Hilbert space is the one containing the state with alternating ↑ and ↓ fermions. In such large sectors, after resolving all the symmetries, the model can be shown to be non-integrable 149 .
The proposal for many-body scars in the effective Hamiltonian in Eq. (48) is based on the existence of a regular subgraph which has the form of the hypergrid -a Cartesian product of line graphs (in our case, of length 3), i.e., the hypergrid is isomorphic to an adjacency graph of a free spin-1 paramagnet, as illustrated below. The mapping to the hypergrid follows from analysing pairs of sites, which can be in states " − " ≡ (↓↑), "2" ≡ ( 0) or " + " ≡ (↑↓), leading to a three level system. Note that the configuration (0 ) is omitted, as doublons can only be formed by hopping to the left. Inside the hypergrid, there are two states for which the cell alternates between − and +. These are the state |−+ ≡ |−+−+ . . . = |↓↑↑↓↓↑↑↓ . . . and its spin-inverted partner, |+− = |+−+− . . . ≡ |↑↓↓↑↑↓↓↑ . . . . The states |−+ and |+− are shown in black colour in figure below. These two states are the only corners of the hypergrid (state with only + and − cells) with no edges going out of it. When the system is initialised in either of these states, it was numerically shown to exhibit persistent oscillations in the quench dynamics and other scarring phenomenology similar to the PXP model discussed in Sec. IV.

VIII. CONCLUSIONS
The goal of this chapter was to present a pedagogical introduction to the new kind of dynamical behaviour in many-body systems -weak ergodicity breaking. While we refrained from giving a formal definition, weak ergodicity breaking was intuitively introduced as a strong dependence of relaxation dynamics on the system's initial configuration. We contrasted this behaviour against conventional ergodic systems, by showing that certain many-body states can be long-lived and they exhibit parametrically slow relaxation, unlike other initial configurations that quickly reach thermal equilibrium. Many recent studies discussed above, theoretical as well as experimental, have revealed that a number of familiar physical systems do not conform to the expectations of the ETH in its strong form, in that such systems can host non-thermal eigenstates and exhibit long-time coherent dynamics. The common pattern that arises in these diverse physical systems is an emergent non-thermalising subspace, which is dynamically decoupled from the rest of the many-body Hilbert space. Quantum entanglement, in particular, has been the vital tool for identifying such non-thermal subspaces in different models.
At this stage, a complete classification of non-thermal subspaces and their underlying mechanisms is still lacking. In particular, more work is needed to understand the connection between many-body scarring and two broad classes of weak ergodicity breaking phenomena: theories with confinement [150][151][152][153][154] and lattice gauge theories [155][156][157] . The latter, in particular the one-dimensional quantum link model -which has recently been realised in a Bose-Hubbard quantum simulator 148 -have intriguing connections with the PXP model 158 . It would be interesting to explore possible connections with PXP in higher dimensions 159 , as well as to understand the relation with kinetic constraints, which can generally lead to slow glassy-like dynamics 79,[160][161][162][163] .
On the theory side, many questions remain about the quantum-classical correspondence in the physical systems discussed above. The TDVP variational approach complements recent efforts in understanding parallels between classical chaos measures, on the one hand, and thermalising quantum dynamics and its underlying transport coefficients on the other hand 121,122 . While many-body scarred models are quantum chaotic, their properties deviate from other chaotic models such as the Sachdev-Ye-Kitaev (SYK) model 164,165 , which has attracted much attention as the fastest scrambler of quantum information 166 . Such deviations could be identified and studied more systematically using the variational approach and its resulting mixed phase portraits, which appear to be a generic feature of local Hamiltonians. One of the outstanding challenges is to bring together this approach to real-time dynamics with approaches that target the many-body eigenstates, possibly using quantum entanglement and other quantum information techniques as a way of linking quantum revivals and eigenstate properties 167 .
Finally, there is strong experimental and practical interest in weak ergodicity breaking. For example, many-body scarred revivals provide a mechanism for maintaining coherence, despite the presence of interactions which normally scramble local quantum information. In particular, scars in Rydberg chains have already been utilised for the preparation of specific entangled states 168 . This application made use of quantum control based on the variational TDVP approach and its identification of entangled periodic trajectories that simultaneously have small quantum leakage. Thus, scars may have a wider range of applications, for example in protected state transfer on quantum networks or in quantum sensing 169 . Such applications require deeper theoretical understanding of the effects that protect the coherence of scars, as well as the development of general experimental techniques for creating them on demand, e.g., using periodic driving in Rydberg arrays 140 , pumping protocols in dipolar Bose gases 170 or by tilting the Fermi-Hubbard model 21 .