Deep thermalization in constrained quantum systems

The concept of"deep thermalization"has recently been introduced to characterize moments of an ensemble of pure states, resulting from projective measurements on a subsystem, which lie beyond the purview of conventional Eigenstate Thermalization Hypothesis (ETH). In this work, we study deep thermalization in systems with kinetic constraints, such as the quantum East and the PXP models, which have been known to weakly break ETH by the slow dynamics and high sensitivity to the initial conditions. We demonstrate a sharp contrast in deep thermalization between the first and higher moments in these models by studying quench dynamics from initial product states in the computational basis: while the first moment shows good agreement with ETH, higher moments deviate from the uniform Haar ensemble at infinite temperature. We show that such behavior is caused by an interplay of time-reversal symmetry and an operator that anticommutes with the Hamiltonian. We formulate sufficient conditions for violating deep thermalization, even for systems that are otherwise"thermal"in the ETH sense. By appropriately breaking these properties, we illustrate how the PXP model fully deep-thermalizes for all initial product states in the thermodynamic limit. Our results highlight the sensitivity of deep thermalization as a probe of physics beyond ETH in kinetically-constrained systems.


I. INTRODUCTION
Thermal equilibrium is the eventual fate of an isolated chaotic quantum system in the absence of special mechanisms such as integrability [1] or localization [2,3].This behavior is encapsulated by the Eigenstate Thermalization Hypothesis (ETH) [4,5], which explains thermalization at the level of subsystems: under unitary dynamics, a subsystem entangles with its complement such that its reduced density matrix increasingly resembles a thermal Gibbs state.The complementary subsystem, whose states are traced over, assumes the role of thermal bath during the unitary dynamics.This scenario has been extensively tested numerically [6] and in experimental setups of quantum simulators with cold atoms, trapped ions and other engineered quantum systems [7][8][9].
The conventional picture of ETH, however, is blind to the microscopic details of the bath.Recent progress in the control and manipulation of individual degrees of freedom in quantum simulators [10][11][12][13] has brought in a refinement of this picture: states of the bath can be explicitly measured and the state of the subsystem can be studied conditional to a measurement outcome on the bath [14,15].Repeating such projective measurements in a fixed, local basis gives rise to an ensemble of pure states on the subsystem which, along with their corresponding Born probabilities, form the so-called "projected ensemble" [16,17].Through quench dynamics experiments on Rydberg atom arrays [14], numerical simulations of Hamiltonian models [15,18] and Floquet circuits [19,20], universal behavior has been found in the dynamics of the projected ensemble: under chaotic time evolution at infinite temperature, the statistical properties of the projected ensemble become indistinguishable from a uniform ensemble on the Hilbert space, i.e., the maximally entropic Haar ensemble [21].Remarkably, all moments of the projected ensemble were argued to be-come indistinguishable from those of the Haar ensemble.In the language of quantum information theory, the moments then furnish (approximate) quantum state designs [22,23].The picture summarized above, dubbed "deep thermalization" [14,15], illuminates a distinct role of measurements compared to a simple hindrance to thermalization.Accordingly, ETH is generalized to the statistical distribution of wave functions, instead of just expectation values of physical observables.This is part of a broader push to probe physics beyond ETH, complementing other approaches such as free probability [24,25], which make predictions about higher-order correlation functions.
The practical appeal of deep thermalization in chaotic systems lies in its universality: Haar-random ensembles can be generated under unitary dynamics starting from a simple initial state [14].An intriguing question therefore arises how this phenomenology changes in chaotic systems whose thermalization dynamics is strongly dependent on the initial state, such as systems displaying Hilbert space fragmentation, quantum many-body scars and other types of "weak" ergodicity breaking [26][27][28].
The paradigmatic examples include the experimental systems of Rydberg atoms [29,30] and ultracold bosons in a tilted optical lattice [31], which feature persistent quantum revivals when quenched from a special initial state, while they quickly thermalize for typical initial conditions.Such systems are described by an effective spin model called the "PXP model" [32,33], which imposes a kinetic constraint on simultaneous flips of neighboring spins.Similar types of constraints can give rise to slow glassy dynamics in the quantum East model [34,35].In this paper we explore the nature of deep thermalization and its sensitivity to initial conditions for such constrained systems.
In both the PXP and quantum East models studied below, we find deep thermalization to be absent, even at x-operators flipping a spin, subject to the state of its neighbor(s).In the quantum East model (c), the spin flip is allowed (green arrow) if the left neighbor is in the ↓ state, regardless of the state of right neighbor.In the PXP model (d), both neighbors must be in the ↓ state for a flip to occur.From these rules, it is clear that if the spin in B, adjacent to the bipartition, is in ↑ state, the state of A cannot be fully random.This must be taken into account with an appropriate postselection rule, as discussed in Sec.II B.
infinite temperature, while the first moment of the projected ensemble agrees well with ETH.We demonstrate this for a large class of time-evolved initial states as well as the eigenstates of the models.We elucidate the origin of this surprising behavior, finding that it does not stem from the constraints but rather from spectral properties and the existence of special operators that anticommute with the Hamiltonian.Once these are properly taken care of, we observe a restoration of deep thermalization in the thermodynamic limit.This is true even in the PXP model, where clear signs of ergodicity breaking are found in all accessible finite sizes.Our results illustrate that deep thermalization is not only a hallmark of "maximally chaotic" models but more broadly present in models that also display weak ergodicity breaking.The results furthermore highlight the sensitivity of deep thermalization framework compared to standard ETH, as the former re-quires more care beyond resolving the usual global symmetries of the model.This paper is organized as follows: In Sec.II, we review the construction of the projected ensemble and how we quantify its distance from the Haar ensemble.We also discuss the crucial role of symmetries, which are illustrated for the Sachdev-Ye-Kitaev (SYK) model in Sec.III.This maximally-chaotic model will be used as a benchmark for the rest of the paper.In Sec.IV, we show how time-reversal symmetry can prevent deep thermalization of eigenstates, and we illustrate this using the Ising model in mixed transverse and longitudinal fields.As our first kinetically constrained model, we consider the quantum East model in its thermalizing phase in Sec.V.In Sec.VI, we show how the interplay of timereversal invariance and "antisymmetries" can generally prevent deep thermalization for time-evolved states.Finally, in Sec.VII, we study deep thermalization in the constrained PXP model describing one-dimensional Rydberg atom arrays.Here we focus in particular on different types of initial states that have been known to give rise to anomalous dynamical behavior associated with quantum many-body scars [29,36].Our Conclusions are presented in Sec.VIII, while Appendixes contain further details of the analysis and effect of perturbations to the PXP model.

II. THE PROJECTED ENSEMBLE
Consider a quantum spin system in a pure state |Ψ⟩.The state can be prepared, e.g., via unitary time evolution from some product state, as sketched in Fig. 1(a).The system is bipartitioned into two contiguous regions: a subsystem A and its complement B, consisting of N A and N B spins, respectively.For simplicity, we assume the total Hilbert space is given by the tensor product H A ⊗ H B , although, as explained in Sec.II B, this assumption can be relaxed in certain cases.We consider performing projective measurements on the subsystem B in the local computational basis, |z⟩.Each such measurement outputs a classical string z B of length N B .According to the Born rule, such a measurement leaves the subsystem A in a pure state which is conditional to measuring the string z B in the complementary subsystem B, with a probability p z B given by: Note that each pure state on A is indexed by a measurement outcome z B on B. With the states ψ A z B and probabilities p z B , the projected ensemble is defined as the set: Since the measurement basis {|z B ⟩} is orthonormal, the probabilities {p z B } sum to unity.The projected ensemble then forms a probability distribution on H A , whose k th moment is given by: where the sum runs over all states in the projected ensemble E and ρ E acts as a density operator on the Hilbert space H ⊗k A , corresponding to k replicas of H A .The first moment of the distribution in Eq. ( 4) is exactly the reduced density matrix on H A : which is the central object of study in ETH.However, higher moments of the projected ensemble (k>1) are in general not equal to k−fold tensor products of the reduced density matrix.As such, the projected ensemble strictly encodes more information than the reduced density matrix on H A .

A. Haar ensemble
Consider the time evolution of an isolated quantum system initialized in a far-from-equilibrium state.ETH tells us that for ergodic systems in the absence of conservation laws, the reduced density matrix ρ A (t) of a subsystem A (much smaller than the rest of the system) relaxes to a thermal Gibbs state at late times, i.e., at times much larger than the inverse of the system's microscopic energy scales.This can be formally expressed as lim t→∞ ρ A (t) = exp −β Ĥ /Z, where β is the inverse temperature, set by the energy of the initial state, and Z is a normalizing factor.At finite times and finite system sizes, the equality is only exact up to fluctuations that typically subside exponentially with system size.During the dynamics, the microscopic details of the initial state are progressively lost and the system increasingly resembles a featureless thermal state.
Since the first moment of the projected ensemble is exactly the reduced density matrix, it is natural to ask whether this evolution to a featureless thermal state holds for all moments of the projected ensemble.We note that the evolution to a thermal Gibbs state is guided by the principle of maximization of entropy, i.e., the second law of thermodynamics.However, since we are dealing with a probability distribution, it is natural to consider maximizing the entropy of this distribution.This is accomplished by conjecturing that the projected ensemble dynamically evolves to the maximally entropic uniform distribution of pure states on a Hilbert space.For initial states with energies in the middle of the spectrum, i.e. states at infinite temperature, this maximally entropic ensemble is the Haar ensemble [21].One can then quantify the "depth" of thermalization by testing to what extent the moments of the projected ensemble agree with those of the Haar ensemble, see Fig. 1(b).
For pure states |ψ⟩ in a Hilbert space, k th -moments of the Haar ensemble can be constructed as: where the integral runs over the entire Hilbert space H.The k-th moment admits an analytical expression [15,21] ρ where Π k is the subspace of H ⊗k invariant under all permutations of the k copies and d is the dimension of H.
As a consistency check of the conjecture, we highlight that the first moment (mean) of the Haar ensemble is the Identity operator Î/d, which is exactly equal to the Gibbs ensemble at infinite temperature.At finite temperatures, the existence of a universal random ensemble has been hinted at, but its exact form remains unknown.[15].Here on, we restrict our focus to infinite temperature initial states.To quantify the degree to which a system deep thermalizes, the trace distance between the k th -moments has been used in the literature [14,15]: where ∥.∥ 1 denotes the trace norm.Note that ∆ (k) follows a monotonicity relation such that ∆ (k ′ ) < ∆ (k) ∀k ′ < k.This is closely linked to the concept of quantum k-designs, an important resource in quantum information theory.We say an ensemble E forms a kdesign if ∆ (k) = 0, which implies that ∆ (k ′ ) = 0, ∀k ′ < k .Then, a system deep thermalizes for a given k if it forms a k-design in the thermodynamic limit.The way this limit is taken is by fixing the subsystem A and then increasing the size of the "reservoir" subsystem B. Since we are limited to finite-sized systems in this study, we study the scaling of ∆ (k) as N B is increased.As shown in Ref. [15], for a chaotic system ∆ (k) is expected to decrease exponentially in N B for all k.The rate of the exponential decay, however, can be different for different k values and, typically, larger k are more strongly impacted by finite size effects.

B. Effect of symmetries and constraints
Until this point, we have considered the simple case where energy is the only conserved quantity.When additional conservation laws are present, generic states instead thermalize to a Generalized Gibbs Ensemble [37].This ensemble now features Lagrange multipliers enforcing the conservation of all charges.As a consequence, the projected ensemble no longer approaches the Haar ensemble, but a similar generalized Haar ensemble can be introduced using the same principle [18].This construction depends on the specifics of the model and its conserved charges, meaning that it generally has a much more complicated structure than the Haar ensemble.In order to simplify the study of deep thermalization with conservation laws, one can restrict to a single charge sector.However, this is usually not sufficient to recover thermalization to the Haar ensemble, as the conserved charge can introduce correlations between the measured string z B and the projected state ψ A z B .To recover the Haar ensemble, one then needs to use postselection on the measured strings z B [15].We present an example of this procedure in the following section.
The need for postselection also arises in a conceptually different case when dynamical constraints are present.The models we focus on in this work (Sec.V and VII) are defined in terms of Pauli x-matrices, which allow each spin to flip its state.However, the spin-flip operators are dressed with projectors, which encode dynamical constraints: a spin can only flip its state if its immediate neighbors to the left or right are in the ↓ state, see Figs. 1(c)-(d).While such a constraint cannot be expressed as a local operator that commutes with the Hamiltonian, it similarly places strong restrictions on the dynamics.For example, if the boundary spin in B subsystem (i.e., the one closest to the A subsystem) is in the ↑ state, this will impact the state of the spin in A closest to the bipartition.This will be true at all evolution times, hence the state of A will never reach the Haar ensemble.Nevertheless, since the constraints considered in this work are all local and affect only the nearest-neighbor spin pairs, it is easy to take care of them via postselection by keeping only the measurement outcomes in which the boundary spin in B is in the ↓ state.This will, a priori, not exclude any configurations in A, hence deep thermalization becomes possible, as we show in Secs.V-VII below.For a constrained system such as the PXP model in Fig. 1(d), the total Hilbert space does not have a tensor product structure, which can be seen from the fact that the total Hilbert space dimension grows as a Fibonacci number [38].Thus, with the help of postselection, we can lift the assumption that the total Hilbert space must be written in the form H A ⊗ H B .

III. MAXIMALLY CHAOTIC SYK MODEL
A natural toy model for illustrating the concept of deep thermalization is the Sachdev-Ye-Kitaev model [39][40][41].This model serves as a paradigm of quantum many-body chaos [42], thus we expect it to yield good agreement between the projected and Haar ensembles.The model is formulated in terms of 2N Majorana fermions with all-to-all interactions, ĤSYK = i<j<k<l J ijkl χi χj χk χl , where χ denote the Majorana operators and the summation indices take values between 1 and 2N .The couplings J ijkl are randomly drawn from a normal distribution with k) in the SYK model with NA=4 and various system sizes when starting from (a) state in Eq. ( 12) with random angles, and (b) a computational basis state.For each system size, the data is for a single random realization as there is no visible variation between them.
mean 0 and variance 6/(2N ) 3 .One of the many interesting features of the SYK model is that its random-matrix theory (RMT) class changes with system size [43].Indeed, for N =4n the model exhibits spectral statistics of the Gaussian Orthogonal Ensemble (GOE), N =4n + 1 and N =4n + 3 correspond instead to the Gaussian Unitary Ensemble (GUE), and for N =4n+2 to the Gaussian Symplectic Ensemble (GSE) [44].
For the numerical study, it is more convenient to work with spin-1/2 operators.We use the Jordan-Wigner transformation, to convert Eq. ( 9) to a model defined in terms of N spin-1/2 degrees of freedom.The resulting Hamiltonian has a rather complicated structure due to the Jordan-Wigner strings, hence we do not write it down explicitly.Nonetheless, we note that the various terms in the model can be either purely diagonal, flip two spins, or flip four spins.This depends on how many of the Majorana fermions correspond to the same spin site.As a consequence, the spin model conserves the parity of the number of ↓ spins: To probe deep thermalization in the SYK model, similar to Refs.[14,15], we study a global quench by preparing the system in different kinds of initial states.One natural choice is to use initial states belonging to the computational basis, which is also the measurement basis.As a second choice, we use product states where each N B,eff spin j points in a random direction on the Bloch sphere, where θ j and ϕ j are randomly drawn from the intervals [0, π] and [0, 2π], respectively.Below we show that these two classes of initial states lead to stark differences in the deep thermalization.Once we prepare the initial states, we evolve them using the SYK Hamiltonian and compute ∆ (k) at various time steps.Fig. 2 shows the averaged late-time value of ∆ (k) for both types of initial states and N A =4.To obtain a benchmark for ∆ (k) in a maximally chaotic state in the given Hilbert space, we also compute this value for "typical" states, which are simply complex random vectors where both the real and imaginary part of the wave function amplitude are drawn randomly from a normal distribution with mean 0 and variance unity.These are then pure states at infinite temperature.This benchmark value is also shown in Fig. 2, where we see it agrees closely with the results of quenches from states in Eq. ( 12) with angles sampled randomly.By contrast, for computational basis states the exponential decay with N B quickly reaches a plateau and completely stops, suggesting that the projected ensemble never reaches the Haar ensemble.
The origin of this difference between initial states is the symmetry ẐP .The computational basis states are eigenstates of this operator with eigenvalues α= ± 1, and so are the time-evolved states obtained from them.Furthermore, ẐP can be decomposed as As ẐP,B is diagonal in the measurement basis, the measured string z B has a definite value of α B = ± 1 under ẐP,B .Consequently, the corresponding state ψ A z B must be an eigenstate of ẐP,A with eigenvalue α A = α/α B .This additional constraint on the ψ A z B prevents them from uniformly exploring H A and reproducing the Haar ensemble.In contrast, states built from random angles are generically not eigenstates of ẐP and have expectation values of this operator close to 0. Thus, there is no individual constraint on the ψ A z B and we obtain good agreement with the Haar ensemble.
In order to recover deep thermalization for all initial states, a simple solution is to restrict to the symmetry sector with an even number of ↑ spins, or alternatively Z P = (−1) N .As discussed above, this introduces correlations between the measured string z B and the possible states in A. We thus post-select the z B strings to only keep those with Z P,B = (−1) N B .This implies that all ψ A z B will obey Z P,A = (−1) N A , effectively reducing the relevant Hilbert space of subsystem A. We now recover good agreement with the Haar ensemble, see Fig. 3.Note that this reduces the dimension of the subsystem B that acts as a reservoir for the states we are interested in.As such, we follow Ref. [15] and instead of N B we keep track of its effective counterpart N B,eff = log 2 (D B ), where D B is the number of postselected states in B (i.e., states with with Z P,B = (−1) N B in the symmetry sector Z P = (−1) N ).For the rest of this work, we will keep using the notation N B,eff for consistency even when no postselection is done.
Interestingly, for the time-evolved state we see no influence of the RMT class on deep thermalization, as different particle numbers in Fig. 3 all behave identically.The computed ∆ (k) for SYK eigenstates also shows no such dependence, Fig. 4. Furthermore, one can also see that, for all k investigated, the eigenstates at infinite temperature and time-evolved states show very close agreement, matching also the prediction for typical states.These results thus set our benchmark for deep thermalization in a fully chaotic system, which we will use for interpreting the results in the less generic models below.

IV. TIME REVERSAL
The SYK model in the previous section illustrated the effect of conventional symmetries on deep thermalization, as well as the apparent robustness of this phenomenon to the RMT class of the Hamiltonian studied.In this section, we show that for systems that have time-reversal symmetry (and so are in the GOE class), deep thermalization for eigenstates can be prevented and we formulate sufficient conditions for this to happen.
As before, we consider a quantum model defined on a system with N sites that can be divided in two subsystems with N A and N B sites, respectively.We will consider the case where the total Hilbert space is a tensor product H = H A ⊗ H B , as this is where we expect to see deep thermalization without any additional postselection procedure.Let us consider a Hamiltonian that is real in the measurement basis.Its eigenstates will also  3. In all cases, there is excellent agreement between typical states, time-evolved states and eigenstates near energy E=0.This is despite the level statistics being respectively GSE (N =14), GUE (N =15) and GOE (N =16).
be real in this basis and they can be expressed as with all ψ A z B real.This implies that even if the ψ A z B are maximally random, they will never be able to approximate the Haar ensemble, but will instead mimick the ensemble of states invariant under the application of orthogonal (instead of unitary) matrices.While this is trivial so far, we will show that the same effect can be obtained with Hamiltonians that are not real in the measurement basis.
Consider performing a basis change to the Hamiltonian while keeping the measurement basis the same.In particular, let us focus on the case where the change-of-basis matrix V obeys: 1. V can be decomposed into subsystems A and B as V = VA ⊗ VB , with VA and VB unitary; 2. VB is diagonal in the measurement basis.
The corresponding projected ensemble is obtained as where the term e iv(z B ) collects the phase change induced on each bitstring z B and adds an irrelevant overall phase term to each state in the projected ensemble.As VA is a unitary rotation on the subsystem A, it cannot change the agreement between the projected ensemble and the Haar ensemble.Indeed, when comparing the two state ensembles, applying the same unitary rotation to both of them will not modify the result.So we can keep the same projected ensemble and apply V † A to the Haar ensemble instead.But this ensemble is by definition invariant under any unitary rotation.Thus, in the end the basis change V will have no influence on ∆ (k) for the eigenstates.
To construct an example where ∆ (k) is not the same after a basis transformation, take a change-of-basis V that instead obeys: 3. All V j B are diagonal in the measurement basis.
In this case, we still have a departure from the Haar ensemble.Let us denote by K the antiunitary corresponding to complex conjugation.It is straightforward to see that K = KA ⊗ KB .As V |E⟩ is real, it follows that K V |E⟩ = V |E⟩.We can recast both sides of this equation using Eq. ( 14) to obtain where g j (z B ) are the eigenvalues of |z B ⟩ under V j B .As the rotation applied to { ψ A z B } can now depend on z B , this no longer corresponds to a global rotation of the basis in the subsystem A. As a consequence, ∆ (k) is generically not the same for |E⟩ and V |E⟩.Nonetheless, as all the |z B ⟩ are orthogonal, the only way to satisfy the equality between the second and last lines in Eq. ( 16) is if for all holds.This puts a constraint on the individual ψ A z B , with each possible value of z B restricting the corresponding states in A to a hyperplane in the Hilbert space.In specific cases, these hyperplanes can be identical or the dependence on z B can even vanish, leading to the same ∆ (k) as for real states.Appendix A illustrates several different possibilities, one of which will be used in Sec.VI below.

A. Ising model
The impact of basis changes can be showcased using the Ising model with both transverse and longitudinal fields: where we use the parameters h = (1 + √ 5)/4 and g = ( √ 5 + 5)/8 that were found in Ref. [45] to lead to strongly chaotic dynamics.In order to explore the different cases, we define various Hamiltonians obtained through the change of basis Ĥ′ = V ĤIsing V † .We consider three different change-of-basis matrices V : A few remarks are in order.While it is clear that both VS and VXY can be decomposed as VA ⊗ VB , VB is only diagonal for VS .VP cannot be decomposed in the same way, but it can be rewritten as VP = cos(π/8)1 − i sin(π/8) j σz j .This means it corresponds to the case V = j V j A ⊗ V j B with all V j B diagonal.Finally, we note that performing the transformation VXY ĤIsing V † XY leads to the same Ising Hamiltonian considered in Ref. [15], for which the eigenstates at infinite temperature display good agreement with the Haar ensemble.
The values of ∆ (k) with k = 1 to 3 for different bases are shown in Fig. 5.For k = 1, the eigenstates of the Hamiltonians generated by I, VS and VXY show the exact same values of ∆ (k) .However, this is no longer the case for k > 1, where I and VS are still identical but greatly differ from VXY due to the states in the projected ensemble being real.This illustrates the difference between the physical case of k = 1, which is directly linked to the expectation values of observables in subsystem A and thus to ETH, and k > 1.Meanwhile, for VP we see that ∆ (k) is different from the other cases, but the constraint on the states in the projected ensemble still prevents deep thermalization even at infinite temperature.In Appendix A, the dependence of deep thermalization on the angle used in VP is also explored.While these results on deep thermalization of eigenstates are not directly relevant for experiments since preparing high-energy eigenstates is a difficult task, we will show in Sec.VI that the same kind of mechanism can also prevent a time-evolved state from thermalizing.

V. QUANTUM EAST MODEL
Now that we have understood deep thermalization in the SYK model and the pitfalls associated with time reversal and basis transformations, we move on to our first kinetically-constrained model: the quantum East model [34,35].The model is inspired by classical structural glasses and features two different phases: one that has slow thermalization and localized eigenstates, and the other where local observables appear to thermalize from the point of view of ETH [46].We focus on the latter phase to see if deep thermalization can detect any hidden non-ergodicity in this regime.
The quantum East model is defined on a 1D lattice with spin-1/2 degrees of freedom and the Hamiltonian where Pj = (1 − σz j )/2 is a local projector on the spindown state at site j.The projector enforces the following constraint: a spin can only flip if its left neighbor is in the ↓ state, previously illustrated in Fig. 1(b).We will first consider open boundary conditions (OBCs) in order to get rid of lattice momentum as a conserved quantity.In this case, the first site does not have a neighbor to the left, thus the Hamiltonian term is only σx 1 , as indicated in Eq. (22).Note that ĤqEast commutes with the local operator σx N .However, as the latter is fully off-diagonal in the measurement basis, it should not introduce any correlation between subsystems A and B. As such we do not limit our study to a single sector of it.In addition, ĤqEast  is real and anticommutes with Ẑ = N j=1 σz j .Since this is not a conserved quantity, we take no additional step in the computation of ∆ (k) because of it.
However, when performing quenches from random basis states and states built from random angles (12), we see a stark difference between the two for k > 1, Fig. 6(a).For the initial states (12), here we set the ϕ angle of the last site to π/2 to always be equally split between the two symmetry sectors of σx N .While the random-angle initial states deep thermalize as expected, Fig. 6(a) shows that random basis states display the same value of ∆ (k) as real typical states.The source of this difference is not found in symmetries this time, but rather in the anticommuting operator Ẑ.Indeed, all the basis states considered are eigenstates of Ẑ with eigenvalue α = ±1.As they are also real, they are eigenstates of K Ẑ, where K is the complex conjugation operator.Finally, as the Hamiltonian is real, this allows to write implying that the time-evolved states are eigenstates of K Ẑ.The |z B ⟩ states are eigenstates of KB with eigenvalues 1 and eigenstates of Ẑ with eigenvalue α z B = ±1.This allows us to rewrite, using Eq. ( 23), (24) Since Ẑ can be rewritten as Ẑ = ẐA ⊗ ẐB , and the same  is trivially true for K, we have: As in the case of real eigenstates, this puts a constraint on the individual ψ A z B .To see the consequences of the constraint (26), let us consider the case of a single qubit in the projected ensemble with α=1 and N A =1, leading to ẐA =σ z 1 .If we parameterize the states in the projected ensemble as |θ, ϕ⟩ = cos(θ/2) |↓⟩ + e iϕ sin(θ/2) |↑⟩, we get We end up with θ free, but ϕ = ±π/2.If we plot the ψ A z B on the Bloch sphere, they will thus only lie in the Y Z plane where ϕ = ±π/2.This is confirmed in the numerical simulation in Fig. 6(b).Confining the projected states to a hyperplane is enough to reproduce the mean (first moment) of the Haar ensemble, but it is inadequate to reproduce the higher moments, which lies at the root of the disparity between k = 1 and k > 1. If, instead, the initial state is not an eigenstate of Ẑ, then Eq. ( 24) does not put an individual constraint on each state in the projected ensemble.As such, states in the projected ensemble are unconstrained (see Fig. 6) and we can expect deep thermalization for all k.To test this, we compare the late-time average of ∆ (k) of basis states with α=1, α= − 1, as well as symmetric superpositions (|α= + 1⟩ + |α=−1⟩)/ √ 2. This is shown in Fig. 7 (a), which reveals clear signatures of deep thermalization only for the latter states.In order to verify that the observed deviation from deep thermalization is not due to the conserved charge σx N , we can compare our results with the same model with periodic boundary conditions (PBCs).The Hamiltonian is the same as in Eq. ( 22), except for the first term that becomes PN σx 1 .Now, σx N no longer commutes with the Hamiltonian, however, we have a conservation of the lattice momentum.Note that with PBCs the state |↑↑ • • • ↑⟩ is also disconnected from the rest of the Hilbert space, and we explicitly discard it.The results in Fig. 7 (b) show that OBCs and PBCs give extremely similar results for ∆ (k) at late times.As they are also essentially identical to those of typical states, this shows that neither σx N nor lattice momentum affect deep thermalization.However, this is no longer true for eigenstates, as seen in Fig. 8. Momentum-resolved eigenstates are no longer real, except in the K = 0 and K = π sectors.Thus, we can see a clear difference between these sectors and and the rest.Interestingly, we also see a much narrower distribution of ∆ (k) , even in real sectors.
In summary, we find that the quantum East model shows clear signatures of deep thermalization for both eigenstates and time-evolved states, once the relevant symmetries and anticommuting operators are properly accounted for.We address the role of the anticommuting operators and the conditions for deep thermalization more generally in the subsequent section.

VI. ANTICOMMUTING OPERATORS AND TIME REVERSAL
While the Hamiltonian having time-reversal symmetry can change the ensemble to which the eigenstates deep-thermalize to, this is generally not the case for timeevolved state.Indeed, even if the initial state and Hamiltonian are real in the measurement basis, this will not be the case of the time-evolved state.However, for a Hamiltonian that is real in the canonical basis, we can recover the same constraints for the projected ensemble if there exists an observable Ẑ that: 1. anticommutes with the Hamiltonian, { Ẑ, Ĥ} = 0; 2. Can be decomposed into subsystems A and B as Ẑ = ẐA ⊗ ẐB with ẐA and ẐB hermitian; 3. ẐB is diagonal in the measurement basis; 4. Ẑ only has eigenvalues +1 and −1 and squares to the identity.
Note that these conditions imply that ẐA and ẐB both have eigenvalues ±1 and square to the identity.Let us assume that our initial state |ψ⟩ is real (up to a global phase) in the measurement basis and an eigenstate of Ẑ with eigenvalue α = ±1.Then one can perform a change of basis using The resulting time-evolved state will then be The expression in the last line is clearly real, and as such cannot lead to thermalization to the Haar ensemble.Now we can recognize that, if Ẑ = ẐA ⊗ ẐB , then V cannot be decomposed in the same way unless ẐB acts as the identity on B. However, due to the specific form of V , we can show that we still get convergence towards an ensemble similar to that of real states.Let us denote by R the real state such that |ψ(t)⟩ = V † |R⟩.We can then write where the ψA denote the states in the projected ensemble of |R⟩.These states all live in the real hyperplane of the Hilbert space of A. The action of exp −iα z B π 4 ẐA rotates this hyperplane along the axis specified by ẐA .As the angle of rotation is π/4 and the eigenvalues of ẐA are ±1, having α z B equal to plus or minus one ends up in the same hyperplane (see Appendix A for an example).This means that the ψ A z B also live in a hyperplane of the same dimensionality as the one of real states.As a consequence, ∆ (k) will have the same lower bound as for real states.Still, we emphasize here that as the direction of the rotation depends on α z B , the ∆ (k) for the ensemble { ψA z B } (corresponding to |R⟩) and { ψ A z B } (corresponding to |ψ(t)⟩) will not be identical but will only approach the same value as N B → ∞.
Some of the stated conditions can actually be relaxed.Indeed, the requirement that both the initial state and Hamiltonian are real can be changed by requiring that there exists a change of basis Û that satisfies Û = ÛA ⊗ ÛB with ÛB diagonal in the measurement basis and such that Û Ĥ Û † and Û |ψ⟩ are real (up to a global phase).One can then use Û V to turn |ψ(t)⟩ into a real vector.

VII. PXP MODEL
Finally, now that we have understood the different conditions that can prevent deep thermalization due to symmetries and anticommuting operators, we consider a model in which the constraint is stronger than in the East model and actually splits the Hilbert space into exponentially many, dynamically disconnected sectors.The model we study is the PXP model, which describes a 1D chain of Rydberg atoms [29,32,33]: Recall that Pj = (1 − σz j )/2 is a projector on the ↓ spin state on site j.In this model, the projectors physically originate from van der Waals interactions between Rydberg atoms, preventing the creation or destruction of neighboring ↑ spins, which is also known as the Rydberg blockade regime [47].In our study, we restrict to the largest connected sector, the one without any configurations that contain neighboring • • • ↑↑ • • • .In Eq. (31) we have assumed OBCs and separated out the boundary terms where the projectors falling outside the boundaries of the chain have been replaced by identity operators.
The PXP model is chaotic but it hosts nonthermal eigenstates known as quantum many-body scars (QMBSs) [29,48].Signatures of QMBSs are visible in global quenches from initial states that have high-overlap with these non-thermalizing QMBS eigenstates.Such quenches have been shown to lead to anomalous dynamics with long-lived coherent oscillations and slow thermalization, despite the model overall displaying chaotic level statistics [36].In particular, the states The origin of these revivals has been understood within a semiclassical approximation [49][50][51], which established a parallel between this many-body phenomenon and quantum scars of a single particle inside a stadium billiard [52].On the other hand, states  While for k=1 there is very little difference between the two, the contrast is strong for k>1 in large systems.31) is purely off-diagonal in the computational basis, all these initial states are effectively at infinite temperature.Therefore, according to strong ETH [9], they are expected to give rise to similar dynamics.This type of ETH violation therefore represents a form of "weak" ergodicity breaking [26].
We now probe deep thermalization in the PXP model.Due to the Rydberg blockade, as explained in Sec.II B, we perform postselection on the z B and only keep the strings where the first spin in B (the one next to subsystem A) is ↓.The same postselection was used in Ref. [14], which studied the full Rydberg model where the blockade is imperfect, i.e., it is enforced as a finite energetic penalty rather than as hard constraint like in our Eq.(31).Moreover, in light of Sec.VI, we must keep in mind that the PXP model in Eq. ( 31) anticommutes with Ẑ = N j=1 σz j , similar to the East model.This operator has also been referred to as "particle-hole transformation" in this context and it was shown to give rise to a large subspace of exact zero energy states [48,53,54].
From our discussion above, we expect a strong dependence of ∆ (k) on the initial state and not just on its energy.However, all the special states |Z n ⟩ are eigenstates of Ẑ with α= ± 1, hence we know they will never exhibit deep thermalization.We can remedy this by adding a small chemical potential in the form of where nj =1− Pj takes the value 1 if jth spin points ↑ (and 0 otherwise).As Ẑ anticommutes with ĤPXP but commutes with Ĥµ , it neither commutes nor anticommutes with Ĥ as long as µ̸ =0.We thus set µ=0.05, which causes all states to deep-thermalize but otherwise does not drastically affect the dynamics, as shown on Fig. 9.
In particular, the initial states considered remain close to infinite temperature and QMBS revivals are present when quenching from |Z 2 ⟩ and |Z 3 ⟩ states.Fig. 9 shows that even a small µ value is sufficient to make a generic state (even those with α= ± 1) deep-thermalize at late times.It also shows that while the agreement with typical states is not as good as for the models considered previously, we recover an exponential decay with N B .
We can now use the value of µ = 0.05 to study the dynamics from special initial states in Fig. 10.While for a given system-size the initial states display large variations in the late-time value of ∆ (k) , for all of them we recover a clear decay of that value with N B .This suggests that, in the thermodynamic limit, all states eventually deep-thermalize, regardless of scarring.Nonetheless, even for randomly sampled states, it is apparent that the time needed to deep-thermalize is remarkably long: the time is on the order of 10 3 , despite the microscopic energy scale of the Hamiltonian being O(1).We note, however, that this thermalization time scale is highly sensitive to the type of perturbation added to the PXP model to break the anticommutation with Ẑ.As shown in Appendix B, another perturbation with a similar strength leads to the a much shorter timescale for thermalization, in line with that in the East model.
The behavior of ∆ (k) for the eigenstates of the PXP model, shown in Fig. 11, also reveals some nonthermalizing features.We find a very broad distribution of values even at infinite temperature.This is similar to what is seen in the entanglement entropy of eigenstates in this model [36,55], due to the presence of many nonthermal eigenstates beyond the "obvious" QMBS outliers   For the latter two models, the initial state is a superposition of two states with α= ± 1, respectively.The difference in timescale as well as the variability between initial states is clearly visible between the three models.
with low entanglement entropy.In this sense, the signs of ergodicity breaking in higher moments (k>1) do not appear to provide additional information compared to the k=1 case, accessible in ETH.

VIII. CONCLUSIONS AND DISCUSSION
In this work, we have studied deep thermalization for several models ranging from the maximally chaotic SYK to the heavily constrained PXP.Our results shed light on how, beyond usual symmetries, the invariance under time reversal and the existence of operators that anticommute with the Hamiltonian can hinder deep thermalization in otherwise chaotic models.As these properties do not influence the first k=1 moment of the projected ensemble, they have no effect on conventional ETH and thermalization of expectation values of local observables.This highlights the sensitivity of deep thermalization to special properties of the model beyond global symmetries.However, once the symmetries and anticommuting operators are properly accounted for, we find an exponential decay of ∆ (k) at late times for all studied models and initial states in the numerically-accessible system sizes.Nonetheless, the rate of convergence to the Haar ensemble can vary greatly between different models or even between initial states in the same model, see Fig. 12 for a brief summary and comparison.We find that for the SYK model, the decay of ∆ (k) with time is well fitted by a power-law for all k, as found in the Ising model in Ref. [15].However, we find a decay as t −2.9 for k = 1 against t −1.2 in the Ising case, showcasing the faster thermalization in SYK.For the East and PXP models, while the decay also resembles a power-law, the variability between initial states implies that using a single exponent for all initial states is not meaningful.Overall, we find values similar to those of the Ising case for the East model, while for PXP the exponent is approximately halved.
The dependence of the exponent on the initial state is the most salient in the PXP model, where we also find that the time needed to reach the plateau is at least an order of magnitude larger than in other models considered.We note that this convergence rate has also been recently shown to sensitively depend on the boundary conditions [56].While the convergence to the Haar ensemble in the PXP model can be enhanced via weak perturbations, its slowness reveals the presence of anomalous dynamics which affects a large part of the Hilbert space (and not just a few specific initial states) and persists up to surprisingly long times.It would be interesting to see if these results could explain the anomalous energy transport in the PXP model at infinite temperature [57], where superdiffusion has been observed on time scales ∼ 10 2 accessible to large-scale tensor network simulations.The onset of deep thermalization at even later times ∼ 10 3 found here suggests that the observed superdiffusion may be transient before it gives way to diffusion.Nevertheless, this would still leave open the question of what physically sets such a long timescale for deep thermalization.The understanding of deep thermalization for initial conditions corresponding to a finite-temperature density matrix and the possible generalization of the projected ensemble to such cases [58] may shed light on this question.Furthermore, a more systematic investigation of models with variable constraints, proposed in Refs.[59,60], might yield more quantitative insights on the interplay between deep thermalization and dynamical connectivity of the Hilbert space, influenced by kinetic constraints.case the value of α z B still matters and thus this transformation is not equivalent to a global rotation.Indeed, let us consider a vector pointing in the +Y direction.In the case of ϕ = −π/4, for α z B = +1 it will end up pointing in the −Y direction, while for α z B = −1 it would be in the +Y direction.As a consequence, the ∆ (k) value will not be identical after this rotation.However, as the projected ensemble is still contained in a single hyperplane, ∆ (k) should fluctuate around the same lower limit as for ϕ = 0 for large enough N B .k) for time-evolved states in the PXP model with the PXPXP perturbation in Eq. (B2) and perturbation strength λ = 0.05.Even such a weak perturbation is sufficient to significantly lower the ∆ (k) plateau for all state while also shortening the time it takes to reach it.

Figure 1 .
Figure 1.(a)-(b): The projected ensemble.A typical setup consists of a spin system, prepared in a product state.The system then evolves under unitary dynamics generated by some Hamiltonian H, which keeps it in a pure state |ψ⟩.After time t, projective measurements are performed on each spin in the subsystem B. The output of this measurement, |zB⟩, indexes the resulting pure state on subsystem A, labeled by ψ A z B .For all the models studied in this work, we pick the subsystem A as the first NA spins, and subsystem B as the remaining N −NA spins.(b) The evolution of states in the projected ensemble for a single qubit.With time, the states uniformly occupy the Bloch sphere, becoming indistinguishable from the Haar ensemble.(c)-(d): The constrained models considered in this work are defined in terms of Paulix-operators flipping a spin, subject to the state of its neighbor(s).In the quantum East model (c), the spin flip is allowed (green arrow) if the left neighbor is in the ↓ state, regardless of the state of right neighbor.In the PXP model (d), both neighbors must be in the ↓ state for a flip to occur.From these rules, it is clear that if the spin in B, adjacent to the bipartition, is in ↑ state, the state of A cannot be fully random.This must be taken into account with an appropriate postselection rule, as discussed in Sec.II B.

Figure 3 .
Figure 3. Late-time average of ∆ (k) in the SYK model with NA=4, ẐP =(−1) N and postselection for (a) a computational basis state and (b) a product state (12) with random angles.

Figure 4 .
Figure 4. ∆ (k) for eigenstates of the SYK model for NA=4 and N =14, 15 and 16, with restriction to a single symmetry sector of ẐP and postselection.The red markers indicate the longtime average for the time-evolved states shown in Fig.3.In all cases, there is excellent agreement between typical states, time-evolved states and eigenstates near energy E=0.This is despite the level statistics being respectively GSE (N =14), GUE (N =15) and GOE (N =16).

Figure 5 .
Figure 5. ∆(k) for eigenstates of the quantum Ising model(18) for NA = 3 and N = 16 after various unitary transformations generated by matrices V .The dashed black (dash-dotted red) line indicates the average value for complex (real) typical states.The insets show zooms of the main panels.Only VXY leads to good agreement with typical complex states.

Figure 6 .
Figure 6.Deep thermalization in the quantum East model for N = 12.(a) ∆ (k) with NA = 3 for two different types of initial states at infinite temperature.The dashed black lines (dashdotted red lines) indicate the average ∆ (k) for (real) states.(b)-(c) For NA = 1, we plot 50 states from the projected ensemble of a time-evolved (b) basis state and (c) a state with random angles at t = 10000.For the former, all vectors are in the Y Z plane while for the latter they explore all sectors of the Bloch sphere.

Figure 7 .
Figure 7. Long-time average of ∆ (k) in the quantum East model with NA = 3 and various system sizes for (a) OBCs and (b) PBCs.Both types of boundary conditions show very similar results.The effect of the anticommuting operator Ẑ is clearly visible when comparing initial states.

Figure 8 .
Figure 8. ∆ (k) for the eigenstates in the quantum East model with NA = 3 and N = 14 with (a) OBCs and (b) PBCs.The effect of the eigenstates being real is clearly visible for all |E⟩ with OBCs, and for momenta K = 0 and K = π for PBCs.

Figure 9 .
Figure 9. Late-time average of ∆ (k) in the PXP model with NA=3 and various system sizes with (a) µ=0 and (b) µ=0.05.While for k=1 there is very little difference between the two, the contrast is strong for k>1 in large systems.

Figure 10 .
Figure 10.∆(k) for time-evolved states in the PXP model with NA = 3 and µ = 0.05.We contrast the states |Z2⟩ and |Z3⟩, which give to QMBS dynamics, with |Z0⟩ and |Z4⟩ that do not exhibit scarring revivals.We also consider the "010" state which contains a single ↑ spin in the middle of |Z0⟩.(a) Time series data for N = 22.(b) Late-time average for several system sizes.The strong influence on initial states is clearly visible, but all states still show a decrease of ∆ (k) with NB.However, even generic states and thermalizing states like |010⟩ show a very slow decay of ∆(k) , with the plateau only reached at times of order 10 3 .

Figure 11 .
Figure 11.∆ (k) for the eigenstates of the PXP model with NA=3 and N =22 with (a) µ=0 and (b) µ=0.05.The late-time averages of time-evolved states are also shown.The effect of eigenstates being real is clearly visible for all |E⟩, while µ has no strong influence for the eigenstates.Near E = 0, the range of values of ∆ (k) is much larger than for the other models considered in this work.

4 tFigure 12 .
Figure 12.A comparison of ∆(k) for (a) the SYK model with N = 16, (b) the East model with N = 16, and (c) the PXP model with N = 22.Panel (a) also shows a power-law fit for the SYK model of the form ∆ (k) ∼ t λ .For the latter two models, the initial state is a superposition of two states with α= ± 1, respectively.The difference in timescale as well as the variability between initial states is clearly visible between the three models.

Figure 15 .
Figure 15.∆(k) for time-evolved states in the PXP model with the PXPXP perturbation in Eq. (B2) and perturbation strength λ = 0.05.Even such a weak perturbation is sufficient to significantly lower the ∆ (k) plateau for all state while also shortening the time it takes to reach it.